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: }