Actual source code: borthog2.c

      1: /*
      2:     Routines used for the orthogonalization of the Hessenberg matrix.

      4:     Note that for the complex numbers version, the VecDot() and
      5:     VecMDot() arguments within the code MUST remain in the order
      6:     given for correct computation of inner products.
      7: */
      8: #include <../src/ksp/ksp/impls/gmres/gmresimpl.h>

     10: /*@C
     11:   KSPGMRESClassicalGramSchmidtOrthogonalization -  This is the basic orthogonalization routine
     12:   using classical Gram-Schmidt with possible iterative refinement to improve the stability

     14:   Collective, No Fortran Support

     16:   Input Parameters:
     17: + ksp - `KSP` object, must be associated with `KSPGMRES`, `KSPFGMRES`, or `KSPLGMRES` Krylov method
     18: - it  - one less than the current GMRES restart iteration, i.e. the size of the Krylov space

     20:   Options Database Keys:
     21: + -ksp_gmres_classicalgramschmidt                                             - Activates `KSPGMRESClassicalGramSchmidtOrthogonalization()`
     22: - -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - determine if iterative refinement is
     23:                                    used to increase the stability of the classical Gram-Schmidt  orthogonalization.

     25:   Level: intermediate

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

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

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

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

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

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

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

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

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

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