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