Actual source code: gcr.c

petsc-master 2020-08-06
Report Typos and Errors

  2:  #include <petsc/private/kspimpl.h>

  4: typedef struct {
  5:   PetscInt    restart;
  6:   PetscInt    n_restarts;
  7:   PetscScalar *val;
  8:   Vec         *VV, *SS;
  9:   Vec         R;

 11:   PetscErrorCode (*modifypc)(KSP,PetscInt,PetscReal,void*);  /* function to modify the preconditioner*/
 12:   PetscErrorCode (*modifypc_destroy)(void*);                 /* function to destroy the user context for the modifypc function */

 14:   void *modifypc_ctx;                                        /* user defined data for the modifypc function */
 15: } KSP_GCR;

 17: static PetscErrorCode KSPSolve_GCR_cycle(KSP ksp)
 18: {
 19:   KSP_GCR        *ctx = (KSP_GCR*)ksp->data;
 21:   PetscScalar    r_dot_v;
 22:   Mat            A, B;
 23:   PC             pc;
 24:   Vec            s,v,r;
 25:   /*
 26:      The residual norm will not be computed when ksp->its > ksp->chknorm hence need to initialize norm_r with some dummy value
 27:   */
 28:   PetscReal      norm_r = 0.0,nrm;
 29:   PetscInt       k, i, restart;
 30:   Vec            x;

 33:   restart = ctx->restart;
 34:   KSPGetPC(ksp, &pc);
 35:   KSPGetOperators(ksp, &A, &B);

 37:   x = ksp->vec_sol;
 38:   r = ctx->R;

 40:   for (k=0; k<restart; k++) {
 41:     v = ctx->VV[k];
 42:     s = ctx->SS[k];
 43:     if (ctx->modifypc) {
 44:       (*ctx->modifypc)(ksp,ksp->its,ksp->rnorm,ctx->modifypc_ctx);
 45:     }

 47:     KSP_PCApply(ksp, r, s); /* s = B^{-1} r */
 48:     KSP_MatMult(ksp,A, s, v);  /* v = A s */

 50:     VecMDot(v,k, ctx->VV, ctx->val);
 51:     for (i=0; i<k; i++) ctx->val[i] = -ctx->val[i];
 52:     VecMAXPY(v,k,ctx->val,ctx->VV); /* v = v - sum_{i=0}^{k-1} alpha_i v_i */
 53:     VecMAXPY(s,k,ctx->val,ctx->SS); /* s = s - sum_{i=0}^{k-1} alpha_i s_i */

 55:     VecDotNorm2(r,v,&r_dot_v,&nrm);
 56:     nrm     = PetscSqrtReal(nrm);
 57:     r_dot_v = r_dot_v/nrm;
 58:     VecScale(v, 1.0/nrm);
 59:     VecScale(s, 1.0/nrm);
 60:     VecAXPY(x,  r_dot_v, s);
 61:     VecAXPY(r, -r_dot_v, v);
 62:     if (ksp->its > ksp->chknorm && ksp->normtype != KSP_NORM_NONE) {
 63:       VecNorm(r, NORM_2, &norm_r);
 64:       KSPCheckNorm(ksp,norm_r);
 65:     }
 66:     /* update the local counter and the global counter */
 67:     ksp->its++;
 68:     ksp->rnorm = norm_r;

 70:     KSPLogResidualHistory(ksp,norm_r);
 71:     KSPMonitor(ksp,ksp->its,norm_r);

 73:     if (ksp->its-1 > ksp->chknorm) {
 74:       (*ksp->converged)(ksp,ksp->its,norm_r,&ksp->reason,ksp->cnvP);
 75:       if (ksp->reason) break;
 76:     }

 78:     if (ksp->its >= ksp->max_it) {
 79:       ksp->reason = KSP_CONVERGED_ITS;
 80:       break;
 81:     }
 82:   }
 83:   ctx->n_restarts++;
 84:   return(0);
 85: }

 87: static PetscErrorCode KSPSolve_GCR(KSP ksp)
 88: {
 89:   KSP_GCR        *ctx = (KSP_GCR*)ksp->data;
 91:   Mat            A, B;
 92:   Vec            r,b,x;
 93:   PetscReal      norm_r = 0.0;

 96:   KSPGetOperators(ksp, &A, &B);
 97:   x    = ksp->vec_sol;
 98:   b    = ksp->vec_rhs;
 99:   r    = ctx->R;

101:   /* compute initial residual */
102:   KSP_MatMult(ksp,A, x, r);
103:   VecAYPX(r, -1.0, b); /* r = b - A x  */
104:   if (ksp->normtype != KSP_NORM_NONE) {
105:     VecNorm(r, NORM_2, &norm_r);
106:     KSPCheckNorm(ksp,norm_r);
107:   }
108:   ksp->its    = 0;
109:   ksp->rnorm0 = norm_r;

111:   KSPLogResidualHistory(ksp,ksp->rnorm0);
112:   KSPMonitor(ksp,ksp->its,ksp->rnorm0);
113:   (*ksp->converged)(ksp,ksp->its,ksp->rnorm0,&ksp->reason,ksp->cnvP);
114:   if (ksp->reason) return(0);

116:   do {
117:     KSPSolve_GCR_cycle(ksp);
118:     if (ksp->reason) return(0); /* catch case when convergence occurs inside the cycle */
119:   } while (ksp->its < ksp->max_it);

121:   if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
122:   return(0);
123: }

125: static PetscErrorCode KSPView_GCR(KSP ksp, PetscViewer viewer)
126: {
127:   KSP_GCR        *ctx = (KSP_GCR*)ksp->data;
129:   PetscBool      iascii;

132:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
133:   if (iascii) {
134:     PetscViewerASCIIPrintf(viewer,"  restart = %D \n", ctx->restart);
135:     PetscViewerASCIIPrintf(viewer,"  restarts performed = %D \n", ctx->n_restarts);
136:   }
137:   return(0);
138: }


141: static PetscErrorCode KSPSetUp_GCR(KSP ksp)
142: {
143:   KSP_GCR        *ctx = (KSP_GCR*)ksp->data;
145:   Mat            A;
146:   PetscBool      diagonalscale;

149:   PCGetDiagonalScale(ksp->pc,&diagonalscale);
150:   if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);

152:   KSPGetOperators(ksp, &A, NULL);
153:   MatCreateVecs(A, &ctx->R, NULL);
154:   VecDuplicateVecs(ctx->R, ctx->restart, &ctx->VV);
155:   VecDuplicateVecs(ctx->R, ctx->restart, &ctx->SS);

157:   PetscMalloc1(ctx->restart, &ctx->val);
158:   return(0);
159: }

161: static PetscErrorCode KSPReset_GCR(KSP ksp)
162: {
164:   KSP_GCR        *ctx = (KSP_GCR*)ksp->data;

167:   VecDestroy(&ctx->R);
168:   VecDestroyVecs(ctx->restart,&ctx->VV);
169:   VecDestroyVecs(ctx->restart,&ctx->SS);
170:   if (ctx->modifypc_destroy) {
171:     (*ctx->modifypc_destroy)(ctx->modifypc_ctx);
172:   }
173:   PetscFree(ctx->val);
174:   return(0);
175: }

177: static PetscErrorCode KSPDestroy_GCR(KSP ksp)
178: {

182:   KSPReset_GCR(ksp);
183:   KSPDestroyDefault(ksp);
184:   return(0);
185: }

187: static PetscErrorCode KSPSetFromOptions_GCR(PetscOptionItems *PetscOptionsObject,KSP ksp)
188: {
190:   KSP_GCR        *ctx = (KSP_GCR*)ksp->data;
191:   PetscInt       restart;
192:   PetscBool      flg;

195:   PetscOptionsHead(PetscOptionsObject,"KSP GCR options");
196:   PetscOptionsInt("-ksp_gcr_restart","Number of Krylov search directions","KSPGCRSetRestart",ctx->restart,&restart,&flg);
197:   if (flg) { KSPGCRSetRestart(ksp,restart); }
198:   PetscOptionsTail();
199:   return(0);
200: }

203: typedef PetscErrorCode (*KSPGCRModifyPCFunction)(KSP,PetscInt,PetscReal,void*);
204: typedef PetscErrorCode (*KSPGCRDestroyFunction)(void*);

206: static PetscErrorCode  KSPGCRSetModifyPC_GCR(KSP ksp,KSPGCRModifyPCFunction function,void *data,KSPGCRDestroyFunction destroy)
207: {
208:   KSP_GCR *ctx = (KSP_GCR*)ksp->data;

212:   ctx->modifypc         = function;
213:   ctx->modifypc_destroy = destroy;
214:   ctx->modifypc_ctx     = data;
215:   return(0);
216: }

218: /*@C
219:  KSPGCRSetModifyPC - Sets the routine used by GCR to modify the preconditioner.

221:  Logically Collective on ksp

223:  Input Parameters:
224:  +  ksp      - iterative context obtained from KSPCreate()
225:  .  function - user defined function to modify the preconditioner
226:  .  ctx      - user provided contex for the modify preconditioner function
227:  -  destroy  - the function to use to destroy the user provided Section 1.5 Writing Application Codes with PETSc context.

229:  Calling Sequence of function:
230:   PetscErrorCode function (KSP ksp, PetscInt n, PetscReal rnorm, void *ctx)

232:  ksp   - iterative context
233:  n     - the total number of GCR iterations that have occurred
234:  rnorm - 2-norm residual value
235:  ctx   - the user provided Section 1.5 Writing Application Codes with PETSc context

237:  Level: intermediate

239:  Notes:
240:  The default modifypc routine is KSPGCRModifyPCNoChange()

242:  .seealso: KSPGCRModifyPCNoChange()

244:  @*/
245: PetscErrorCode  KSPGCRSetModifyPC(KSP ksp,PetscErrorCode (*function)(KSP,PetscInt,PetscReal,void*),void *data,PetscErrorCode (*destroy)(void*))
246: {

250:   PetscUseMethod(ksp,"KSPGCRSetModifyPC_C",(KSP,PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*),void *data,PetscErrorCode (*)(void*)),(ksp,function,data,destroy));
251:   return(0);
252: }

254: static PetscErrorCode KSPGCRSetRestart_GCR(KSP ksp,PetscInt restart)
255: {
256:   KSP_GCR *ctx;

259:   ctx          = (KSP_GCR*)ksp->data;
260:   ctx->restart = restart;
261:   return(0);
262: }

264: PetscErrorCode  KSPGCRSetRestart(KSP ksp, PetscInt restart)
265: {

269:   PetscTryMethod(ksp,"KSPGCRSetRestart_C",(KSP,PetscInt),(ksp,restart));
270:   return(0);
271: }

273: static PetscErrorCode  KSPBuildSolution_GCR(KSP ksp, Vec v, Vec *V)
274: {
276:   Vec            x;

279:   x = ksp->vec_sol;
280:   if (v) {
281:     VecCopy(x, v);
282:     if (V) *V = v;
283:   } else if (V) {
284:     *V = ksp->vec_sol;
285:   }
286:   return(0);
287: }

289: static PetscErrorCode  KSPBuildResidual_GCR(KSP ksp, Vec t, Vec v, Vec *V)
290: {
292:   KSP_GCR        *ctx;

295:   ctx = (KSP_GCR*)ksp->data;
296:   if (v) {
297:     VecCopy(ctx->R, v);
298:     if (V) *V = v;
299:   } else if (V) {
300:     *V = ctx->R;
301:   }
302:   return(0);
303: }

305: /*MC
306:      KSPGCR - Implements the preconditioned Generalized Conjugate Residual method.

308:    Options Database Keys:
309: .   -ksp_gcr_restart <restart> - the number of stored vectors to orthogonalize against

311:    Level: beginner

313:     Notes:
314:     The GCR Krylov method supports non-symmetric matrices and permits the use of a preconditioner
315:            which may vary from one iteration to the next. Users can can define a method to vary the
316:            preconditioner between iterates via KSPGCRSetModifyPC().

318:            Restarts are solves with x0 not equal to zero. When a restart occurs, the initial starting
319:            solution is given by the current estimate for x which was obtained by the last restart
320:            iterations of the GCR algorithm.

322:            Unlike GMRES and FGMRES, when using GCR, the solution and residual vector can be directly accessed at any iterate,
323:            with zero computational cost, via a call to KSPBuildSolution() and KSPBuildResidual() respectively.

325:            This implementation of GCR will only apply the stopping condition test whenever ksp->its > ksp->chknorm,
326:            where ksp->chknorm is specified via the command line argument -ksp_check_norm_iteration or via
327:            the function KSPSetCheckNormIteration(). Hence the residual norm reported by the monitor and stored
328:            in the residual history will be listed as 0.0 before this iteration. It is actually not 0.0; just not calculated.

330:            The method implemented requires the storage of 2 x restart + 1 vectors, twice as much as GMRES.
331:            Support only for right preconditioning.

333:     Contributed by Dave May

335:     References:
336: .          1. - S. C. Eisenstat, H. C. Elman, and H. C. Schultz. Variational iterative methods for
337:            nonsymmetric systems of linear equations. SIAM J. Numer. Anal., 20, 1983


340: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
341:            KSPGCRSetRestart(), KSPGCRSetModifyPC(), KSPGMRES, KSPFGMRES

343: M*/
344: PETSC_EXTERN PetscErrorCode KSPCreate_GCR(KSP ksp)
345: {
347:   KSP_GCR        *ctx;

350:   PetscNewLog(ksp,&ctx);

352:   ctx->restart    = 30;
353:   ctx->n_restarts = 0;
354:   ksp->data       = (void*)ctx;

356:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);
357:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,3);

359:   ksp->ops->setup          = KSPSetUp_GCR;
360:   ksp->ops->solve          = KSPSolve_GCR;
361:   ksp->ops->reset          = KSPReset_GCR;
362:   ksp->ops->destroy        = KSPDestroy_GCR;
363:   ksp->ops->view           = KSPView_GCR;
364:   ksp->ops->setfromoptions = KSPSetFromOptions_GCR;
365:   ksp->ops->buildsolution  = KSPBuildSolution_GCR;
366:   ksp->ops->buildresidual  = KSPBuildResidual_GCR;

368:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGCRSetRestart_C",KSPGCRSetRestart_GCR);
369:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGCRSetModifyPC_C",KSPGCRSetModifyPC_GCR);
370:   return(0);
371: }