Actual source code: gcr.c

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

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

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

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

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

 30:   PetscFunctionBegin;
 31:   restart = ctx->restart;
 32:   PetscCall(KSPGetPC(ksp, &pc));
 33:   PetscCall(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) PetscCall((*ctx->modifypc)(ksp, ksp->its, ksp->rnorm, ctx->modifypc_ctx));

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

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

 51:     PetscCall(VecDotNorm2(r, v, &r_dot_v, &nrm));
 52:     nrm     = PetscSqrtReal(nrm);
 53:     r_dot_v = r_dot_v / nrm;
 54:     PetscCall(VecScale(v, 1.0 / nrm));
 55:     PetscCall(VecScale(s, 1.0 / nrm));
 56:     PetscCall(VecAXPY(x, r_dot_v, s));
 57:     PetscCall(VecAXPY(r, -r_dot_v, v));
 58:     if (ksp->its > ksp->chknorm && ksp->normtype != KSP_NORM_NONE) {
 59:       PetscCall(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:     PetscCall(KSPLogResidualHistory(ksp, norm_r));
 67:     PetscCall(KSPMonitor(ksp, ksp->its, norm_r));

 69:     if (ksp->its - 1 > ksp->chknorm) {
 70:       PetscCall((*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:   PetscFunctionReturn(PETSC_SUCCESS);
 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:   PetscFunctionBegin;
 91:   PetscCall(KSPGetOperators(ksp, &A, &B));
 92:   x = ksp->vec_sol;
 93:   b = ksp->vec_rhs;
 94:   r = ctx->R;

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

106:   PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm0));
107:   PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm0));
108:   PetscCall((*ksp->converged)(ksp, ksp->its, ksp->rnorm0, &ksp->reason, ksp->cnvP));
109:   if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

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

116:   if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
117:   PetscFunctionReturn(PETSC_SUCCESS);
118: }

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

125:   PetscFunctionBegin;
126:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
127:   if (iascii) {
128:     PetscCall(PetscViewerASCIIPrintf(viewer, "  restart = %" PetscInt_FMT " \n", ctx->restart));
129:     PetscCall(PetscViewerASCIIPrintf(viewer, "  restarts performed = %" PetscInt_FMT " \n", ctx->n_restarts));
130:   }
131:   PetscFunctionReturn(PETSC_SUCCESS);
132: }

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

140:   PetscFunctionBegin;
141:   PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
142:   PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);

144:   PetscCall(KSPGetOperators(ksp, &A, NULL));
145:   PetscCall(MatCreateVecs(A, &ctx->R, NULL));
146:   PetscCall(VecDuplicateVecs(ctx->R, ctx->restart, &ctx->VV));
147:   PetscCall(VecDuplicateVecs(ctx->R, ctx->restart, &ctx->SS));

149:   PetscCall(PetscMalloc1(ctx->restart, &ctx->val));
150:   PetscFunctionReturn(PETSC_SUCCESS);
151: }

153: static PetscErrorCode KSPReset_GCR(KSP ksp)
154: {
155:   KSP_GCR *ctx = (KSP_GCR *)ksp->data;

157:   PetscFunctionBegin;
158:   PetscCall(VecDestroy(&ctx->R));
159:   PetscCall(VecDestroyVecs(ctx->restart, &ctx->VV));
160:   PetscCall(VecDestroyVecs(ctx->restart, &ctx->SS));
161:   if (ctx->modifypc_destroy) PetscCall((*ctx->modifypc_destroy)(ctx->modifypc_ctx));
162:   PetscCall(PetscFree(ctx->val));
163:   PetscFunctionReturn(PETSC_SUCCESS);
164: }

166: static PetscErrorCode KSPDestroy_GCR(KSP ksp)
167: {
168:   PetscFunctionBegin;
169:   PetscCall(KSPReset_GCR(ksp));
170:   PetscCall(KSPDestroyDefault(ksp));
171:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRSetRestart_C", NULL));
172:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRGetRestart_C", NULL));
173:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRSetModifyPC_C", NULL));
174:   PetscFunctionReturn(PETSC_SUCCESS);
175: }

177: static PetscErrorCode KSPSetFromOptions_GCR(KSP ksp, PetscOptionItems *PetscOptionsObject)
178: {
179:   KSP_GCR  *ctx = (KSP_GCR *)ksp->data;
180:   PetscInt  restart;
181:   PetscBool flg;

183:   PetscFunctionBegin;
184:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP GCR options");
185:   PetscCall(PetscOptionsInt("-ksp_gcr_restart", "Number of Krylov search directions", "KSPGCRSetRestart", ctx->restart, &restart, &flg));
186:   if (flg) PetscCall(KSPGCRSetRestart(ksp, restart));
187:   PetscOptionsHeadEnd();
188:   PetscFunctionReturn(PETSC_SUCCESS);
189: }

192: typedef PetscErrorCode (*KSPGCRModifyPCFunction)(KSP, PetscInt, PetscReal, void *);
193: typedef PetscErrorCode (*KSPGCRDestroyFunction)(void *);

195: static PetscErrorCode KSPGCRSetModifyPC_GCR(KSP ksp, KSPGCRModifyPCFunction function, void *data, KSPGCRDestroyFunction destroy)
196: {
197:   KSP_GCR *ctx = (KSP_GCR *)ksp->data;

199:   PetscFunctionBegin;
201:   ctx->modifypc         = function;
202:   ctx->modifypc_destroy = destroy;
203:   ctx->modifypc_ctx     = data;
204:   PetscFunctionReturn(PETSC_SUCCESS);
205: }

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

210:   Logically Collective

212:   Input Parameters:
213: + ksp      - iterative context obtained from `KSPCreate()`
214: . function - user defined function to modify the preconditioner
215: . ctx      - user provided context for the modify preconditioner function
216: - destroy  - the function to use to destroy the user provided application context.

218:   Calling sequence of `function`:
219: + ksp   - iterative context
220: . n     - the total number of `PCGCR` iterations that have occurred
221: . rnorm - 2-norm residual value
222: - ctx   - the user provided application context

224:   Calling sequence of `destroy`:
225: . ctx - the user provided application context

227:   Level: intermediate

229:   Note:
230:   The default modifypc routine is `KSPGCRModifyPCNoChange()`

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

235: .seealso: [](ch_ksp), `KSP`, `KSPGCR`, `KSPGCRModifyPCNoChange()`, [](sec_flexibleksp)
236:  @*/
237: PetscErrorCode KSPGCRSetModifyPC(KSP ksp, PetscErrorCode (*function)(KSP ksp, PetscInt n, PetscReal rnorm, void *ctx), void *ctx, PetscErrorCode (*destroy)(void *ctx))
238: {
239:   PetscFunctionBegin;
240:   PetscUseMethod(ksp, "KSPGCRSetModifyPC_C", (KSP, PetscErrorCode(*)(KSP, PetscInt, PetscReal, void *), void *data, PetscErrorCode (*)(void *)), (ksp, function, ctx, destroy));
241:   PetscFunctionReturn(PETSC_SUCCESS);
242: }

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

248:   PetscFunctionBegin;
249:   ctx          = (KSP_GCR *)ksp->data;
250:   ctx->restart = restart;
251:   PetscFunctionReturn(PETSC_SUCCESS);
252: }

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

258:   PetscFunctionBegin;
259:   ctx      = (KSP_GCR *)ksp->data;
260:   *restart = ctx->restart;
261:   PetscFunctionReturn(PETSC_SUCCESS);
262: }

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

267:   Not Collective

269:   Input Parameters:
270: + ksp     - the Krylov space context
271: - restart - integer restart value

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

276:   Level: intermediate

278:   Note:
279:   The default value is 30.

281:   Developer Note:
282:   The API could be made uniform for all `KSP` methods that have have a restart.

284: .seealso: [](ch_ksp), `KSPGCR`, `KSPSetTolerances()`, `KSPGCRGetRestart()`, `KSPGMRESSetRestart()`
285: @*/
286: PetscErrorCode KSPGCRSetRestart(KSP ksp, PetscInt restart)
287: {
288:   PetscFunctionBegin;
289:   PetscTryMethod(ksp, "KSPGCRSetRestart_C", (KSP, PetscInt), (ksp, restart));
290:   PetscFunctionReturn(PETSC_SUCCESS);
291: }

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

296:   Not Collective

298:   Input Parameter:
299: . ksp - the Krylov space context

301:   Output Parameter:
302: . restart - integer restart value

304:   Level: intermediate

306: .seealso: [](ch_ksp), `KSPGCR`, `KSPSetTolerances()`, `KSPGCRSetRestart()`, `KSPGMRESGetRestart()`
307: @*/
308: PetscErrorCode KSPGCRGetRestart(KSP ksp, PetscInt *restart)
309: {
310:   PetscFunctionBegin;
311:   PetscTryMethod(ksp, "KSPGCRGetRestart_C", (KSP, PetscInt *), (ksp, restart));
312:   PetscFunctionReturn(PETSC_SUCCESS);
313: }

315: static PetscErrorCode KSPBuildSolution_GCR(KSP ksp, Vec v, Vec *V)
316: {
317:   Vec x;

319:   PetscFunctionBegin;
320:   x = ksp->vec_sol;
321:   if (v) {
322:     PetscCall(VecCopy(x, v));
323:     if (V) *V = v;
324:   } else if (V) {
325:     *V = ksp->vec_sol;
326:   }
327:   PetscFunctionReturn(PETSC_SUCCESS);
328: }

330: static PetscErrorCode KSPBuildResidual_GCR(KSP ksp, Vec t, Vec v, Vec *V)
331: {
332:   KSP_GCR *ctx;

334:   PetscFunctionBegin;
335:   ctx = (KSP_GCR *)ksp->data;
336:   if (v) {
337:     PetscCall(VecCopy(ctx->R, v));
338:     if (V) *V = v;
339:   } else if (V) {
340:     *V = ctx->R;
341:   }
342:   PetscFunctionReturn(PETSC_SUCCESS);
343: }

345: /*MC
346:      KSPGCR - Implements the preconditioned flexible Generalized Conjugate Residual method {cite}`eisenstat1983variational`. [](sec_flexibleksp),

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

351:    Level: beginner

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

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

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

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

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

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

375:     Contributed by:
376:     Dave May

378: .seealso: [](ch_ksp), [](sec_flexibleksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPGCRSetRestart()`, `KSPGCRGetRestart()`,
379:           `KSPGCRSetRestart()`, `KSPGCRSetModifyPC()`, `KSPGMRES`, `KSPFGMRES`
380: M*/
381: PETSC_EXTERN PetscErrorCode KSPCreate_GCR(KSP ksp)
382: {
383:   KSP_GCR *ctx;

385:   PetscFunctionBegin;
386:   PetscCall(PetscNew(&ctx));

388:   ctx->restart    = 30;
389:   ctx->n_restarts = 0;
390:   ksp->data       = (void *)ctx;

392:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
393:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 3));

395:   ksp->ops->setup          = KSPSetUp_GCR;
396:   ksp->ops->solve          = KSPSolve_GCR;
397:   ksp->ops->reset          = KSPReset_GCR;
398:   ksp->ops->destroy        = KSPDestroy_GCR;
399:   ksp->ops->view           = KSPView_GCR;
400:   ksp->ops->setfromoptions = KSPSetFromOptions_GCR;
401:   ksp->ops->buildsolution  = KSPBuildSolution_GCR;
402:   ksp->ops->buildresidual  = KSPBuildResidual_GCR;

404:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRSetRestart_C", KSPGCRSetRestart_GCR));
405:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRGetRestart_C", KSPGCRGetRestart_GCR));
406:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRSetModifyPC_C", KSPGCRSetModifyPC_GCR));
407:   PetscFunctionReturn(PETSC_SUCCESS);
408: }