Actual source code: borthog2.c

      1: #define PETSCKSP_DLL

      3: /*
      4:     Routines used for the orthogonalization of the Hessenberg matrix.

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

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

     16:      Collective on KSP

     18:   Input Parameters:
     19: +   ksp - KSP object, must be associated with GMRES, FGMRES, or LGMRES Krylov method
     20: -   its - one less then the current GMRES restart iteration, i.e. the size of the Krylov space

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

     27:     Notes: Use KSPGMRESSetCGSRefinementType() to determine if iterative refinement is to be used

     29:    Level: intermediate

     31: .seelaso:  KSPGMRESSetOrthogonalization(), KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESSetCGSRefinementType()

     33: @*/
     36: PetscErrorCode PETSCKSP_DLLEXPORT KSPGMRESClassicalGramSchmidtOrthogonalization(KSP  ksp,PetscInt it)
     37: {
     38:   KSP_GMRES      *gmres = (KSP_GMRES *)(ksp->data);
     40:   PetscInt       j;
     41:   PetscScalar    *hh,*hes,shh[512],*lhh;
     42:   PetscReal      hnrm, wnrm;
     43:   PetscTruth     refine = (PetscTruth)(gmres->cgstype == KSP_GMRES_CGS_REFINE_ALWAYS);

     46:   PetscLogEventBegin(KSP_GMRESOrthogonalization,ksp,0,0,0);
     47:   /* Don't allocate small arrays */
     48:   if (it < 512) lhh = shh;
     49:   else {
     50:     PetscMalloc((it+1) * sizeof(PetscScalar),&lhh);
     51:   }
     52:
     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(it+1,VEC_VV(it+1),&(VEC_VV(0)),lhh); /* <v,vnew> */
     68:   for (j=0; j<=it; j++) {
     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++) {
     90:       hnrm  +=  PetscRealPart(lhh[j] * PetscConj(lhh[j]));
     91:     }
     92:     hnrm = sqrt(hnrm);
     93:     VecNorm(VEC_VV(it+1),NORM_2, &wnrm);
     94:     if (wnrm < 1.0286 * hnrm) {
     95:       refine = PETSC_TRUE;
     96:       PetscLogInfo((ksp,"KSPGMRESClassicalGramSchmidtOrthogonalization:Performing iterative refinement wnorm %g hnorm %g\n",wnrm,hnrm));
     97:     }
     98:   }

    100:   if (refine) {
    101:     VecMDot(it+1,VEC_VV(it+1),&(VEC_VV(0)),lhh); /* <v,vnew> */
    102:     for (j=0; j<=it; j++) lhh[j] = - lhh[j];
    103:     VecMAXPY(VEC_VV(it+1),it+1,lhh,&VEC_VV(0));
    104:     /* note lhh[j] is -<v,vnew> , hence the subtraction */
    105:     for (j=0; j<=it; j++) {
    106:       hh[j]  -= lhh[j];     /* hh += <v,vnew> */
    107:       hes[j] -= lhh[j];     /* hes += <v,vnew> */
    108:     }
    109:   }

    111:   if (it >= 512) {PetscFree(lhh);}
    112:   PetscLogEventEnd(KSP_GMRESOrthogonalization,ksp,0,0,0);
    113:   return(0);
    114: }