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