Actual source code: gcr.c


  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;
 20:   PetscScalar r_dot_v;
 21:   Mat         A, B;
 22:   PC          pc;
 23:   Vec         s, v, r;
 24:   /*
 25:      The residual norm will not be computed when ksp->its > ksp->chknorm hence need to initialize norm_r with some dummy value
 26:   */
 27:   PetscReal norm_r = 0.0, nrm;
 28:   PetscInt  k, i, restart;
 29:   Vec       x;

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

 35:   x = ksp->vec_sol;
 36:   r = ctx->R;

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

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

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

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

 66:     KSPLogResidualHistory(ksp, norm_r);
 67:     KSPMonitor(ksp, ksp->its, norm_r);

 69:     if (ksp->its - 1 > ksp->chknorm) {
 70:       (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
 71:       if (ksp->reason) break;
 72:     }

 74:     if (ksp->its >= ksp->max_it) {
 75:       ksp->reason = KSP_CONVERGED_ITS;
 76:       break;
 77:     }
 78:   }
 79:   ctx->n_restarts++;
 80:   return 0;
 81: }

 83: static PetscErrorCode KSPSolve_GCR(KSP ksp)
 84: {
 85:   KSP_GCR  *ctx = (KSP_GCR *)ksp->data;
 86:   Mat       A, B;
 87:   Vec       r, b, x;
 88:   PetscReal norm_r = 0.0;

 90:   KSPGetOperators(ksp, &A, &B);
 91:   x = ksp->vec_sol;
 92:   b = ksp->vec_rhs;
 93:   r = ctx->R;

 95:   /* compute initial residual */
 96:   KSP_MatMult(ksp, A, x, r);
 97:   VecAYPX(r, -1.0, b); /* r = b - A x  */
 98:   if (ksp->normtype != KSP_NORM_NONE) {
 99:     VecNorm(r, NORM_2, &norm_r);
100:     KSPCheckNorm(ksp, norm_r);
101:   }
102:   ksp->its    = 0;
103:   ksp->rnorm0 = norm_r;

105:   KSPLogResidualHistory(ksp, ksp->rnorm0);
106:   KSPMonitor(ksp, ksp->its, ksp->rnorm0);
107:   (*ksp->converged)(ksp, ksp->its, ksp->rnorm0, &ksp->reason, ksp->cnvP);
108:   if (ksp->reason) return 0;

110:   do {
111:     KSPSolve_GCR_cycle(ksp);
112:     if (ksp->reason) return 0; /* catch case when convergence occurs inside the cycle */
113:   } while (ksp->its < ksp->max_it);

115:   if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
116:   return 0;
117: }

119: static PetscErrorCode KSPView_GCR(KSP ksp, PetscViewer viewer)
120: {
121:   KSP_GCR  *ctx = (KSP_GCR *)ksp->data;
122:   PetscBool iascii;

124:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
125:   if (iascii) {
126:     PetscViewerASCIIPrintf(viewer, "  restart = %" PetscInt_FMT " \n", ctx->restart);
127:     PetscViewerASCIIPrintf(viewer, "  restarts performed = %" PetscInt_FMT " \n", ctx->n_restarts);
128:   }
129:   return 0;
130: }

132: static PetscErrorCode KSPSetUp_GCR(KSP ksp)
133: {
134:   KSP_GCR  *ctx = (KSP_GCR *)ksp->data;
135:   Mat       A;
136:   PetscBool diagonalscale;

138:   PCGetDiagonalScale(ksp->pc, &diagonalscale);

141:   KSPGetOperators(ksp, &A, NULL);
142:   MatCreateVecs(A, &ctx->R, NULL);
143:   VecDuplicateVecs(ctx->R, ctx->restart, &ctx->VV);
144:   VecDuplicateVecs(ctx->R, ctx->restart, &ctx->SS);

146:   PetscMalloc1(ctx->restart, &ctx->val);
147:   return 0;
148: }

150: static PetscErrorCode KSPReset_GCR(KSP ksp)
151: {
152:   KSP_GCR *ctx = (KSP_GCR *)ksp->data;

154:   VecDestroy(&ctx->R);
155:   VecDestroyVecs(ctx->restart, &ctx->VV);
156:   VecDestroyVecs(ctx->restart, &ctx->SS);
157:   if (ctx->modifypc_destroy) (*ctx->modifypc_destroy)(ctx->modifypc_ctx);
158:   PetscFree(ctx->val);
159:   return 0;
160: }

162: static PetscErrorCode KSPDestroy_GCR(KSP ksp)
163: {
164:   KSPReset_GCR(ksp);
165:   KSPDestroyDefault(ksp);
166:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRSetRestart_C", NULL);
167:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRGetRestart_C", NULL);
168:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRSetModifyPC_C", NULL);
169:   return 0;
170: }

172: static PetscErrorCode KSPSetFromOptions_GCR(KSP ksp, PetscOptionItems *PetscOptionsObject)
173: {
174:   KSP_GCR  *ctx = (KSP_GCR *)ksp->data;
175:   PetscInt  restart;
176:   PetscBool flg;

178:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP GCR options");
179:   PetscOptionsInt("-ksp_gcr_restart", "Number of Krylov search directions", "KSPGCRSetRestart", ctx->restart, &restart, &flg);
180:   if (flg) KSPGCRSetRestart(ksp, restart);
181:   PetscOptionsHeadEnd();
182:   return 0;
183: }

186: typedef PetscErrorCode (*KSPGCRModifyPCFunction)(KSP, PetscInt, PetscReal, void *);
187: typedef PetscErrorCode (*KSPGCRDestroyFunction)(void *);

189: static PetscErrorCode KSPGCRSetModifyPC_GCR(KSP ksp, KSPGCRModifyPCFunction function, void *data, KSPGCRDestroyFunction destroy)
190: {
191:   KSP_GCR *ctx = (KSP_GCR *)ksp->data;

194:   ctx->modifypc         = function;
195:   ctx->modifypc_destroy = destroy;
196:   ctx->modifypc_ctx     = data;
197:   return 0;
198: }

200: /*@C
201:  KSPGCRSetModifyPC - Sets the routine used by `KSPGCR` to modify the preconditioner for each iteration

203:  Logically Collective

205:  Input Parameters:
206:  +  ksp      - iterative context obtained from KSPCreate()
207:  .  function - user defined function to modify the preconditioner
208:  .  ctx      - user provided context for the modify preconditioner function
209:  -  destroy  - the function to use to destroy the user provided application context.

211:  Calling Sequence of function:
212:   PetscErrorCode function (KSP ksp, PetscInt n, PetscReal rnorm, void *ctx)
213: +  ksp   - iterative context
214: .  n     - the total number of GCR iterations that have occurred
215: .  rnorm - 2-norm residual value
216: -  ctx   - the user provided application context

218:  Level: intermediate

220:  Note:
221:  The default modifypc routine is `KSPGCRModifyPCNoChange()`

223:  Developer Note:
224:  The API should make uniform for all flexible types, [](sec_flexibleksp), and not have separate function calls for each type.

226: .seealso: [](chapter_ksp), `KSP`, `KSPGCR`, `KSPGCRModifyPCNoChange()`, [](sec_flexibleksp)
227:  @*/
228: PetscErrorCode KSPGCRSetModifyPC(KSP ksp, PetscErrorCode (*function)(KSP, PetscInt, PetscReal, void *), void *data, PetscErrorCode (*destroy)(void *))
229: {
230:   PetscUseMethod(ksp, "KSPGCRSetModifyPC_C", (KSP, PetscErrorCode(*)(KSP, PetscInt, PetscReal, void *), void *data, PetscErrorCode (*)(void *)), (ksp, function, data, destroy));
231:   return 0;
232: }

234: static PetscErrorCode KSPGCRSetRestart_GCR(KSP ksp, PetscInt restart)
235: {
236:   KSP_GCR *ctx;

238:   ctx          = (KSP_GCR *)ksp->data;
239:   ctx->restart = restart;
240:   return 0;
241: }

243: static PetscErrorCode KSPGCRGetRestart_GCR(KSP ksp, PetscInt *restart)
244: {
245:   KSP_GCR *ctx;

247:   ctx      = (KSP_GCR *)ksp->data;
248:   *restart = ctx->restart;
249:   return 0;
250: }

252: /*@
253:    KSPGCRSetRestart - Sets number of iterations at which `KSPGCR` restarts.

255:    Not Collective

257:    Input Parameters:
258: +  ksp - the Krylov space context
259: -  restart - integer restart value

261:    Options Database Key:
262: .   -ksp_gcr_restart <restart> - the number of stored vectors to orthogonalize against

264:    Level: intermediate

266:    Note:
267:    The default value is 30.

269:    Developer Note:
270:    The API could be made uniform for all `KSP` methods have have a restart.

272: .seealso: [](chapter_ksp), `KSPGCR`, `KSPSetTolerances()`, `KSPGCRGetRestart()`, `KSPGMRESSetRestart()`
273: @*/
274: PetscErrorCode KSPGCRSetRestart(KSP ksp, PetscInt restart)
275: {
276:   PetscTryMethod(ksp, "KSPGCRSetRestart_C", (KSP, PetscInt), (ksp, restart));
277:   return 0;
278: }

280: /*@
281:    KSPGCRGetRestart - Gets number of iterations at which `KSPGCR` restarts.

283:    Not Collective

285:    Input Parameter:
286: .  ksp - the Krylov space context

288:    Output Parameter:
289: .   restart - integer restart value

291:    Level: intermediate

293: .seealso: [](chapter_ksp), `KSPGCR`, `KSPSetTolerances()`, `KSPGCRSetRestart()`, `KSPGMRESGetRestart()`
294: @*/
295: PetscErrorCode KSPGCRGetRestart(KSP ksp, PetscInt *restart)
296: {
297:   PetscTryMethod(ksp, "KSPGCRGetRestart_C", (KSP, PetscInt *), (ksp, restart));
298:   return 0;
299: }

301: static PetscErrorCode KSPBuildSolution_GCR(KSP ksp, Vec v, Vec *V)
302: {
303:   Vec x;

305:   x = ksp->vec_sol;
306:   if (v) {
307:     VecCopy(x, v);
308:     if (V) *V = v;
309:   } else if (V) {
310:     *V = ksp->vec_sol;
311:   }
312:   return 0;
313: }

315: static PetscErrorCode KSPBuildResidual_GCR(KSP ksp, Vec t, Vec v, Vec *V)
316: {
317:   KSP_GCR *ctx;

319:   ctx = (KSP_GCR *)ksp->data;
320:   if (v) {
321:     VecCopy(ctx->R, v);
322:     if (V) *V = v;
323:   } else if (V) {
324:     *V = ctx->R;
325:   }
326:   return 0;
327: }

329: /*MC
330:      KSPGCR - Implements the preconditioned flexible Generalized Conjugate Residual method. [](sec_flexibleksp),

332:    Options Database Key:
333: .   -ksp_gcr_restart <restart> - the number of stored vectors to orthogonalize against

335:    Level: beginner

337:     Notes:
338:     The GCR Krylov method supports non-symmetric matrices and permits the use of a preconditioner
339:     which may vary from one iteration to the next.

341:     Users can can define a method to vary the
342:            preconditioner between iterates via `KSPGCRSetModifyPC()`.

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

348:            Unlike `KSPGMRES` and `KSPFGMRES`, when using GCR, the solution and residual vector can be directly accessed at any iterate,
349:            with zero computational cost, via a call to `KSPBuildSolution()` and `KSPBuildResidual()` respectively.

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

356:            The method implemented requires the storage of 2 x restart + 1 vectors, twice as much as `KSPGMRES`.
357:            Support only for right preconditioning.

359:     Contributed by:
360:     Dave May

362:     References:
363: .   * - S. C. Eisenstat, H. C. Elman, and H. C. Schultz. Variational iterative methods for
364:            nonsymmetric systems of linear equations. SIAM J. Numer. Anal., 20, 1983

366: .seealso: [](chapter_ksp), [](sec_flexibleksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPGCRSetRestart()`, `KSPGCRGetRestart()`,
367:           `KSPGCRSetRestart()`, `KSPGCRSetModifyPC()`, `KSPGMRES`, `KSPFGMRES`
368: M*/
369: PETSC_EXTERN PetscErrorCode KSPCreate_GCR(KSP ksp)
370: {
371:   KSP_GCR *ctx;

373:   PetscNew(&ctx);

375:   ctx->restart    = 30;
376:   ctx->n_restarts = 0;
377:   ksp->data       = (void *)ctx;

379:   KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1);
380:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 3);

382:   ksp->ops->setup          = KSPSetUp_GCR;
383:   ksp->ops->solve          = KSPSolve_GCR;
384:   ksp->ops->reset          = KSPReset_GCR;
385:   ksp->ops->destroy        = KSPDestroy_GCR;
386:   ksp->ops->view           = KSPView_GCR;
387:   ksp->ops->setfromoptions = KSPSetFromOptions_GCR;
388:   ksp->ops->buildsolution  = KSPBuildSolution_GCR;
389:   ksp->ops->buildresidual  = KSPBuildResidual_GCR;

391:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRSetRestart_C", KSPGCRSetRestart_GCR);
392:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRGetRestart_C", KSPGCRGetRestart_GCR);
393:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRSetModifyPC_C", KSPGCRSetModifyPC_GCR);
394:   return 0;
395: }