Actual source code: borthog2.c


  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:
 27:     Use KSPGMRESSetCGSRefinementType() to determine if iterative refinement is to be used.
 28:     This is much faster than KSPGMRESModifiedGramSchmidtOrthogonalization() but has the small possibility of stability issues
 29:     that can usually be handled by using a a single step of iterative refinement with KSPGMRESSetCGSRefinementType()

 31:    Level: intermediate

 33: .seelaso:  KSPGMRESSetOrthogonalization(), KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESSetCGSRefinementType(),
 34:            KSPGMRESGetCGSRefinementType(), KSPGMRESGetOrthogonalization(), KSPGMRESModifiedGramSchmidtOrthogonalization()

 36: @*/
 37: PetscErrorCode  KSPGMRESClassicalGramSchmidtOrthogonalization(KSP ksp,PetscInt it)
 38: {
 39:   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);

 45:   PetscLogEventBegin(KSP_GMRESOrthogonalization,ksp,0,0,0);
 46:   if (!gmres->orthogwork) {
 47:     PetscMalloc1(gmres->max_k + 2,&gmres->orthogwork);
 48:   }
 49:   lhh = gmres->orthogwork;

 51:   /* update Hessenberg matrix and do unmodified Gram-Schmidt */
 52:   hh  = HH(0,it);
 53:   hes = HES(0,it);

 55:   /* Clear hh and hes since we will accumulate values into them */
 56:   for (j=0; j<=it; j++) {
 57:     hh[j]  = 0.0;
 58:     hes[j] = 0.0;
 59:   }

 61:   /*
 62:      This is really a matrix-vector product, with the matrix stored
 63:      as pointer to rows
 64:   */
 65:   VecMDot(VEC_VV(it+1),it+1,&(VEC_VV(0)),lhh); /* <v,vnew> */
 66:   for (j=0; j<=it; j++) {
 67:     KSPCheckDot(ksp,lhh[j]);
 68:     if (ksp->reason) goto done;
 69:     lhh[j] = -lhh[j];
 70:   }

 72:   /*
 73:          This is really a matrix vector product:
 74:          [h[0],h[1],...]*[ v[0]; v[1]; ...] subtracted from v[it+1].
 75:   */
 76:   VecMAXPY(VEC_VV(it+1),it+1,lhh,&VEC_VV(0));
 77:   /* note lhh[j] is -<v,vnew> , hence the subtraction */
 78:   for (j=0; j<=it; j++) {
 79:     hh[j]  -= lhh[j];     /* hh += <v,vnew> */
 80:     hes[j] -= lhh[j];     /* hes += <v,vnew> */
 81:   }

 83:   /*
 84:      the second step classical Gram-Schmidt is only necessary
 85:      when a simple test criteria is not passed
 86:   */
 87:   if (gmres->cgstype == KSP_GMRES_CGS_REFINE_IFNEEDED) {
 88:     hnrm = 0.0;
 89:     for (j=0; j<=it; j++) hnrm +=  PetscRealPart(lhh[j] * PetscConj(lhh[j]));

 91:     hnrm = PetscSqrtReal(hnrm);
 92:     VecNorm(VEC_VV(it+1),NORM_2, &wnrm);
 93:     KSPCheckNorm(ksp,wnrm);
 94:     if (ksp->reason) goto done;
 95:     if (wnrm < hnrm) {
 96:       refine = PETSC_TRUE;
 97:       PetscInfo(ksp,"Performing iterative refinement wnorm %g hnorm %g\n",(double)wnrm,(double)hnrm);
 98:     }
 99:   }

101:   if (refine) {
102:     VecMDot(VEC_VV(it+1),it+1,&(VEC_VV(0)),lhh); /* <v,vnew> */
103:     for (j=0; j<=it; j++) {
104:        KSPCheckDot(ksp,lhh[j]);
105:        if (ksp->reason) goto done;
106:        lhh[j] = -lhh[j];
107:     }
108:     VecMAXPY(VEC_VV(it+1),it+1,lhh,&VEC_VV(0));
109:     /* note lhh[j] is -<v,vnew> , hence the subtraction */
110:     for (j=0; j<=it; j++) {
111:       hh[j]  -= lhh[j];     /* hh += <v,vnew> */
112:       hes[j] -= lhh[j];     /* hes += <v,vnew> */
113:     }
114:   }
115: done:
116:   PetscLogEventEnd(KSP_GMRESOrthogonalization,ksp,0,0,0);
117:   return 0;
118: }