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);
     41:   PetscInt       j;
     42:   PetscScalar    *hh,*hes,*lhh;
     43:   PetscReal      hnrm, wnrm;
     44:   PetscBool      refine = (PetscBool)(gmres->cgstype == KSP_GMRES_CGS_REFINE_ALWAYS);

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

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

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

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

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

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

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

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