Actual source code: borthog2.c
petsc-3.3-p7 2013-05-11
2: /*
3: Routines used for the orthogonalization of the Hessenberg matrix.
5: Note that for the complex numbers version, the VecDot() and
6: VecMDot() arguments within the code MUST remain in the order
7: given for correct computation of inner products.
8: */
9: #include <../src/ksp/ksp/impls/gmres/gmresimpl.h>
11: /*@C
12: KSPGMRESClassicalGramSchmidtOrthogonalization - This is the basic orthogonalization routine
13: using classical Gram-Schmidt with possible iterative refinement to improve the stability
15: Collective on KSP
17: Input Parameters:
18: + ksp - KSP object, must be associated with GMRES, FGMRES, or LGMRES Krylov method
19: - its - one less then the current GMRES restart iteration, i.e. the size of the Krylov space
21: Options Database Keys:
22: + -ksp_gmres_classicalgramschmidt - Activates KSPGMRESClassicalGramSchmidtOrthogonalization()
23: - -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - determine if iterative refinement is
24: used to increase the stability of the classical Gram-Schmidt orthogonalization.
26: Notes: Use KSPGMRESSetCGSRefinementType() to determine if iterative refinement is to be used
28: Level: intermediate
30: .seelaso: KSPGMRESSetOrthogonalization(), KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESSetCGSRefinementType(),
31: KSPGMRESGetCGSRefinementType(), KSPGMRESGetOrthogonalization()
33: @*/
36: PetscErrorCode KSPGMRESClassicalGramSchmidtOrthogonalization(KSP ksp,PetscInt it)
37: {
38: KSP_GMRES *gmres = (KSP_GMRES *)(ksp->data);
40: PetscInt j;
41: PetscScalar *hh,*hes,*lhh;
42: PetscReal hnrm, wnrm;
43: PetscBool refine = (PetscBool)(gmres->cgstype == KSP_GMRES_CGS_REFINE_ALWAYS);
46: PetscLogEventBegin(KSP_GMRESOrthogonalization,ksp,0,0,0);
47: if (!gmres->orthogwork) {
48: PetscMalloc((gmres->max_k + 2)*sizeof(PetscScalar),&gmres->orthogwork);
49: }
50: lhh = gmres->orthogwork;
51:
52: /* update Hessenberg matrix and do unmodified Gram-Schmidt */
53: hh = HH(0,it);
54: hes = HES(0,it);
56: /* Clear hh and hes since we will accumulate values into them */
57: for (j=0; j<=it; j++) {
58: hh[j] = 0.0;
59: hes[j] = 0.0;
60: }
62: /*
63: This is really a matrix-vector product, with the matrix stored
64: as pointer to rows
65: */
66: VecMDot(VEC_VV(it+1),it+1,&(VEC_VV(0)),lhh); /* <v,vnew> */
67: for (j=0; j<=it; j++) {
68: lhh[j] = - lhh[j];
69: }
71: /*
72: This is really a matrix vector product:
73: [h[0],h[1],...]*[ v[0]; v[1]; ...] subtracted from v[it+1].
74: */
75: VecMAXPY(VEC_VV(it+1),it+1,lhh,&VEC_VV(0));
76: /* note lhh[j] is -<v,vnew> , hence the subtraction */
77: for (j=0; j<=it; j++) {
78: hh[j] -= lhh[j]; /* hh += <v,vnew> */
79: hes[j] -= lhh[j]; /* hes += <v,vnew> */
80: }
82: /*
83: * the second step classical Gram-Schmidt is only necessary
84: * when a simple test criteria is not passed
85: */
86: if (gmres->cgstype == KSP_GMRES_CGS_REFINE_IFNEEDED) {
87: hnrm = 0.0;
88: for (j=0; j<=it; j++) {
89: hnrm += PetscRealPart(lhh[j] * PetscConj(lhh[j]));
90: }
91: hnrm = PetscSqrtReal(hnrm);
92: VecNorm(VEC_VV(it+1),NORM_2, &wnrm);
93: if (wnrm < 1.0286 * hnrm) {
94: refine = PETSC_TRUE;
95: PetscInfo2(ksp,"Performing iterative refinement wnorm %G hnorm %G\n",wnrm,hnrm);
96: }
97: }
99: if (refine) {
100: VecMDot(VEC_VV(it+1),it+1,&(VEC_VV(0)),lhh); /* <v,vnew> */
101: for (j=0; j<=it; j++) lhh[j] = - lhh[j];
102: VecMAXPY(VEC_VV(it+1),it+1,lhh,&VEC_VV(0));
103: /* note lhh[j] is -<v,vnew> , hence the subtraction */
104: for (j=0; j<=it; j++) {
105: hh[j] -= lhh[j]; /* hh += <v,vnew> */
106: hes[j] -= lhh[j]; /* hes += <v,vnew> */
107: }
108: }
109: PetscLogEventEnd(KSP_GMRESOrthogonalization,ksp,0,0,0);
110: return(0);
111: }