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

 17:   Input Parameters:
 18: +   ksp - KSP object, must be associated with `KSPGMRES`, `KSPFGMRES`, or `KSPLGMRES` 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:     Level: intermediate

 28:     Notes:
 29:     Use `KSPGMRESSetCGSRefinementType()` to determine if iterative refinement is to be used.
 30:     This is much faster than `KSPGMRESModifiedGramSchmidtOrthogonalization()` but has the small possibility of stability issues
 31:     that can usually be handled by using a a single step of iterative refinement with `KSPGMRESSetCGSRefinementType()`

 33: .seealso: [](ch_ksp), `KSPGMRESCGSRefinementType`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESSetCGSRefinementType()`,
 34:            `KSPGMRESGetCGSRefinementType()`, `KSPGMRESGetOrthogonalization()`, `KSPGMRESModifiedGramSchmidtOrthogonalization()`
 35: @*/
 36: PetscErrorCode KSPGMRESClassicalGramSchmidtOrthogonalization(KSP ksp, PetscInt it)
 37: {
 38:   KSP_GMRES   *gmres = (KSP_GMRES *)(ksp->data);
 39:   PetscInt     j;
 40:   PetscScalar *hh, *hes, *lhh;
 41:   PetscReal    hnrm, wnrm;
 42:   PetscBool    refine = (PetscBool)(gmres->cgstype == KSP_GMRES_CGS_REFINE_ALWAYS);

 44:   PetscFunctionBegin;
 45:   PetscCall(PetscLogEventBegin(KSP_GMRESOrthogonalization, ksp, 0, 0, 0));
 46:   if (!gmres->orthogwork) PetscCall(PetscMalloc1(gmres->max_k + 2, &gmres->orthogwork));
 47:   lhh = gmres->orthogwork;

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

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

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

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

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

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

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