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:   PetscFunctionBegin;
 32:   restart = ctx->restart;
 33:   PetscCall(KSPGetPC(ksp, &pc));
 34:   PetscCall(KSPGetOperators(ksp, &A, &B));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

211:  Logically Collective

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

219:  Calling Sequence of `function`:
220: $  PetscErrorCode function(KSP ksp, PetscInt n, PetscReal rnorm, void *ctx)
221: +  ksp   - iterative context
222: .  n     - the total number of GCR iterations that have occurred
223: .  rnorm - 2-norm residual value
224: -  ctx   - the user provided application context

226:  Calling Sequence of `destroy`:
227: $ PetscErrorCode destroy(void *ctx)

229:  Level: intermediate

231:  Note:
232:  The default modifypc routine is `KSPGCRModifyPCNoChange()`

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

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

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

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

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

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

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

269:    Not Collective

271:    Input Parameters:
272: +  ksp - the Krylov space context
273: -  restart - integer restart value

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

278:    Level: intermediate

280:    Note:
281:    The default value is 30.

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

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

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

298:    Not Collective

300:    Input Parameter:
301: .  ksp - the Krylov space context

303:    Output Parameter:
304: .   restart - integer restart value

306:    Level: intermediate

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

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

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

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

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

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

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

353:    Level: beginner

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

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

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

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

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

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

377:     Contributed by:
378:     Dave May

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

384: .seealso: [](ch_ksp), [](sec_flexibleksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPGCRSetRestart()`, `KSPGCRGetRestart()`,
385:           `KSPGCRSetRestart()`, `KSPGCRSetModifyPC()`, `KSPGMRES`, `KSPFGMRES`
386: M*/
387: PETSC_EXTERN PetscErrorCode KSPCreate_GCR(KSP ksp)
388: {
389:   KSP_GCR *ctx;

391:   PetscFunctionBegin;
392:   PetscCall(PetscNew(&ctx));

394:   ctx->restart    = 30;
395:   ctx->n_restarts = 0;
396:   ksp->data       = (void *)ctx;

398:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
399:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 3));

401:   ksp->ops->setup          = KSPSetUp_GCR;
402:   ksp->ops->solve          = KSPSolve_GCR;
403:   ksp->ops->reset          = KSPReset_GCR;
404:   ksp->ops->destroy        = KSPDestroy_GCR;
405:   ksp->ops->view           = KSPView_GCR;
406:   ksp->ops->setfromoptions = KSPSetFromOptions_GCR;
407:   ksp->ops->buildsolution  = KSPBuildSolution_GCR;
408:   ksp->ops->buildresidual  = KSPBuildResidual_GCR;

410:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRSetRestart_C", KSPGCRSetRestart_GCR));
411:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRGetRestart_C", KSPGCRGetRestart_GCR));
412:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRSetModifyPC_C", KSPGCRSetModifyPC_GCR));
413:   PetscFunctionReturn(PETSC_SUCCESS);
414: }