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;
21: PetscScalar r_dot_v;
22: Mat A, B;
23: PC pc;
24: Vec s,v,r;
25: /*
26: The residual norm will not be computed when ksp->its > ksp->chknorm hence need to initialize norm_r with some dummy value
27: */
28: PetscReal norm_r = 0.0,nrm;
29: PetscInt k, i, restart;
30: Vec x;
33: restart = ctx->restart;
34: KSPGetPC(ksp, &pc);
35: KSPGetOperators(ksp, &A, &B);
37: x = ksp->vec_sol;
38: r = ctx->R;
40: for (k=0; k<restart; k++) {
41: v = ctx->VV[k];
42: s = ctx->SS[k];
43: if (ctx->modifypc) {
44: (*ctx->modifypc)(ksp,ksp->its,ksp->rnorm,ctx->modifypc_ctx);
45: }
47: KSP_PCApply(ksp, r, s); /* s = B^{-1} r */
48: KSP_MatMult(ksp,A, s, v); /* v = A s */
50: VecMDot(v,k, ctx->VV, ctx->val);
51: for (i=0; i<k; i++) ctx->val[i] = -ctx->val[i];
52: VecMAXPY(v,k,ctx->val,ctx->VV); /* v = v - sum_{i=0}^{k-1} alpha_i v_i */
53: VecMAXPY(s,k,ctx->val,ctx->SS); /* s = s - sum_{i=0}^{k-1} alpha_i s_i */
55: VecDotNorm2(r,v,&r_dot_v,&nrm);
56: nrm = PetscSqrtReal(nrm);
57: r_dot_v = r_dot_v/nrm;
58: VecScale(v, 1.0/nrm);
59: VecScale(s, 1.0/nrm);
60: VecAXPY(x, r_dot_v, s);
61: VecAXPY(r, -r_dot_v, v);
62: if (ksp->its > ksp->chknorm && ksp->normtype != KSP_NORM_NONE) {
63: VecNorm(r, NORM_2, &norm_r);
64: KSPCheckNorm(ksp,norm_r);
65: }
66: /* update the local counter and the global counter */
67: ksp->its++;
68: ksp->rnorm = norm_r;
70: KSPLogResidualHistory(ksp,norm_r);
71: KSPMonitor(ksp,ksp->its,norm_r);
73: if (ksp->its-1 > ksp->chknorm) {
74: (*ksp->converged)(ksp,ksp->its,norm_r,&ksp->reason,ksp->cnvP);
75: if (ksp->reason) break;
76: }
78: if (ksp->its >= ksp->max_it) {
79: ksp->reason = KSP_CONVERGED_ITS;
80: break;
81: }
82: }
83: ctx->n_restarts++;
84: return(0);
85: }
87: static PetscErrorCode KSPSolve_GCR(KSP ksp)
88: {
89: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
91: Mat A, B;
92: Vec r,b,x;
93: PetscReal norm_r = 0.0;
96: KSPGetOperators(ksp, &A, &B);
97: x = ksp->vec_sol;
98: b = ksp->vec_rhs;
99: r = ctx->R;
101: /* compute initial residual */
102: KSP_MatMult(ksp,A, x, r);
103: VecAYPX(r, -1.0, b); /* r = b - A x */
104: if (ksp->normtype != KSP_NORM_NONE) {
105: VecNorm(r, NORM_2, &norm_r);
106: KSPCheckNorm(ksp,norm_r);
107: }
108: ksp->its = 0;
109: ksp->rnorm0 = norm_r;
111: KSPLogResidualHistory(ksp,ksp->rnorm0);
112: KSPMonitor(ksp,ksp->its,ksp->rnorm0);
113: (*ksp->converged)(ksp,ksp->its,ksp->rnorm0,&ksp->reason,ksp->cnvP);
114: if (ksp->reason) return(0);
116: do {
117: KSPSolve_GCR_cycle(ksp);
118: if (ksp->reason) return(0); /* catch case when convergence occurs inside the cycle */
119: } while (ksp->its < ksp->max_it);
121: if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
122: return(0);
123: }
125: static PetscErrorCode KSPView_GCR(KSP ksp, PetscViewer viewer)
126: {
127: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
129: PetscBool iascii;
132: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
133: if (iascii) {
134: PetscViewerASCIIPrintf(viewer," restart = %D \n", ctx->restart);
135: PetscViewerASCIIPrintf(viewer," restarts performed = %D \n", ctx->n_restarts);
136: }
137: return(0);
138: }
141: static PetscErrorCode KSPSetUp_GCR(KSP ksp)
142: {
143: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
145: Mat A;
146: PetscBool diagonalscale;
149: PCGetDiagonalScale(ksp->pc,&diagonalscale);
150: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
152: KSPGetOperators(ksp, &A, NULL);
153: MatCreateVecs(A, &ctx->R, NULL);
154: VecDuplicateVecs(ctx->R, ctx->restart, &ctx->VV);
155: VecDuplicateVecs(ctx->R, ctx->restart, &ctx->SS);
157: PetscMalloc1(ctx->restart, &ctx->val);
158: return(0);
159: }
161: static PetscErrorCode KSPReset_GCR(KSP ksp)
162: {
164: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
167: VecDestroy(&ctx->R);
168: VecDestroyVecs(ctx->restart,&ctx->VV);
169: VecDestroyVecs(ctx->restart,&ctx->SS);
170: if (ctx->modifypc_destroy) {
171: (*ctx->modifypc_destroy)(ctx->modifypc_ctx);
172: }
173: PetscFree(ctx->val);
174: return(0);
175: }
177: static PetscErrorCode KSPDestroy_GCR(KSP ksp)
178: {
182: KSPReset_GCR(ksp);
183: KSPDestroyDefault(ksp);
184: PetscObjectComposeFunction((PetscObject)ksp,"KSPGCRSetRestart_C",NULL);
185: PetscObjectComposeFunction((PetscObject)ksp,"KSPGCRGetRestart_C",NULL);
186: PetscObjectComposeFunction((PetscObject)ksp,"KSPGCRSetModifyPC_C",NULL);
187: return(0);
188: }
190: static PetscErrorCode KSPSetFromOptions_GCR(PetscOptionItems *PetscOptionsObject,KSP ksp)
191: {
193: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
194: PetscInt restart;
195: PetscBool flg;
198: PetscOptionsHead(PetscOptionsObject,"KSP GCR options");
199: PetscOptionsInt("-ksp_gcr_restart","Number of Krylov search directions","KSPGCRSetRestart",ctx->restart,&restart,&flg);
200: if (flg) { KSPGCRSetRestart(ksp,restart); }
201: PetscOptionsTail();
202: return(0);
203: }
206: typedef PetscErrorCode (*KSPGCRModifyPCFunction)(KSP,PetscInt,PetscReal,void*);
207: typedef PetscErrorCode (*KSPGCRDestroyFunction)(void*);
209: static PetscErrorCode KSPGCRSetModifyPC_GCR(KSP ksp,KSPGCRModifyPCFunction function,void *data,KSPGCRDestroyFunction destroy)
210: {
211: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
215: ctx->modifypc = function;
216: ctx->modifypc_destroy = destroy;
217: ctx->modifypc_ctx = data;
218: return(0);
219: }
221: /*@C
222: KSPGCRSetModifyPC - Sets the routine used by GCR to modify the preconditioner.
224: Logically Collective on ksp
226: Input Parameters:
227: + ksp - iterative context obtained from KSPCreate()
228: . function - user defined function to modify the preconditioner
229: . ctx - user provided contex for the modify preconditioner function
230: - destroy - the function to use to destroy the user provided application context.
232: Calling Sequence of function:
233: PetscErrorCode function (KSP ksp, PetscInt n, PetscReal rnorm, void *ctx)
235: ksp - iterative context
236: n - the total number of GCR iterations that have occurred
237: rnorm - 2-norm residual value
238: ctx - the user provided application context
240: Level: intermediate
242: Notes:
243: The default modifypc routine is KSPGCRModifyPCNoChange()
245: .seealso: KSPGCRModifyPCNoChange()
247: @*/
248: PetscErrorCode KSPGCRSetModifyPC(KSP ksp,PetscErrorCode (*function)(KSP,PetscInt,PetscReal,void*),void *data,PetscErrorCode (*destroy)(void*))
249: {
253: PetscUseMethod(ksp,"KSPGCRSetModifyPC_C",(KSP,PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*),void *data,PetscErrorCode (*)(void*)),(ksp,function,data,destroy));
254: return(0);
255: }
257: static PetscErrorCode KSPGCRSetRestart_GCR(KSP ksp,PetscInt restart)
258: {
259: KSP_GCR *ctx;
262: ctx = (KSP_GCR*)ksp->data;
263: ctx->restart = restart;
264: return(0);
265: }
267: static PetscErrorCode KSPGCRGetRestart_GCR(KSP ksp,PetscInt *restart)
268: {
269: KSP_GCR *ctx;
272: ctx = (KSP_GCR*)ksp->data;
273: *restart = ctx->restart;
274: return(0);
275: }
277: /*@
278: KSPGCRSetRestart - Sets number of iterations at which GCR restarts.
280: Not Collective
282: Input Parameter:
283: + ksp - the Krylov space context
284: - restart - integer restart value
286: Note: The default value is 30.
288: Level: intermediate
290: .seealso: KSPSetTolerances(), KSPGCRGetRestart(), KSPGMRESSetRestart()
291: @*/
292: PetscErrorCode KSPGCRSetRestart(KSP ksp, PetscInt restart)
293: {
297: PetscTryMethod(ksp,"KSPGCRSetRestart_C",(KSP,PetscInt),(ksp,restart));
298: return(0);
299: }
301: /*@
302: KSPGCRGetRestart - Gets number of iterations at which GCR restarts.
304: Not Collective
306: Input Parameter:
307: . ksp - the Krylov space context
309: Output Parameter:
310: . restart - integer restart value
312: Note: The default value is 30.
314: Level: intermediate
316: .seealso: KSPSetTolerances(), KSPGCRSetRestart(), KSPGMRESGetRestart()
317: @*/
318: PetscErrorCode KSPGCRGetRestart(KSP ksp, PetscInt *restart)
319: {
323: PetscTryMethod(ksp,"KSPGCRGetRestart_C",(KSP,PetscInt*),(ksp,restart));
324: return(0);
325: }
327: static PetscErrorCode KSPBuildSolution_GCR(KSP ksp, Vec v, Vec *V)
328: {
330: Vec x;
333: x = ksp->vec_sol;
334: if (v) {
335: VecCopy(x, v);
336: if (V) *V = v;
337: } else if (V) {
338: *V = ksp->vec_sol;
339: }
340: return(0);
341: }
343: static PetscErrorCode KSPBuildResidual_GCR(KSP ksp, Vec t, Vec v, Vec *V)
344: {
346: KSP_GCR *ctx;
349: ctx = (KSP_GCR*)ksp->data;
350: if (v) {
351: VecCopy(ctx->R, v);
352: if (V) *V = v;
353: } else if (V) {
354: *V = ctx->R;
355: }
356: return(0);
357: }
359: /*MC
360: KSPGCR - Implements the preconditioned Generalized Conjugate Residual method.
362: Options Database Keys:
363: . -ksp_gcr_restart <restart> - the number of stored vectors to orthogonalize against
365: Level: beginner
367: Notes:
368: The GCR Krylov method supports non-symmetric matrices and permits the use of a preconditioner
369: which may vary from one iteration to the next. Users can can define a method to vary the
370: preconditioner between iterates via KSPGCRSetModifyPC().
372: Restarts are solves with x0 not equal to zero. When a restart occurs, the initial starting
373: solution is given by the current estimate for x which was obtained by the last restart
374: iterations of the GCR algorithm.
376: Unlike GMRES and FGMRES, when using GCR, the solution and residual vector can be directly accessed at any iterate,
377: with zero computational cost, via a call to KSPBuildSolution() and KSPBuildResidual() respectively.
379: This implementation of GCR will only apply the stopping condition test whenever ksp->its > ksp->chknorm,
380: where ksp->chknorm is specified via the command line argument -ksp_check_norm_iteration or via
381: the function KSPSetCheckNormIteration(). Hence the residual norm reported by the monitor and stored
382: in the residual history will be listed as 0.0 before this iteration. It is actually not 0.0; just not calculated.
384: The method implemented requires the storage of 2 x restart + 1 vectors, twice as much as GMRES.
385: Support only for right preconditioning.
387: Contributed by Dave May
389: References:
390: . 1. - S. C. Eisenstat, H. C. Elman, and H. C. Schultz. Variational iterative methods for
391: nonsymmetric systems of linear equations. SIAM J. Numer. Anal., 20, 1983
394: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
395: KSPGCRSetRestart(), KSPGCRSetModifyPC(), KSPGMRES, KSPFGMRES
397: M*/
398: PETSC_EXTERN PetscErrorCode KSPCreate_GCR(KSP ksp)
399: {
401: KSP_GCR *ctx;
404: PetscNewLog(ksp,&ctx);
406: ctx->restart = 30;
407: ctx->n_restarts = 0;
408: ksp->data = (void*)ctx;
410: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);
411: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,3);
413: ksp->ops->setup = KSPSetUp_GCR;
414: ksp->ops->solve = KSPSolve_GCR;
415: ksp->ops->reset = KSPReset_GCR;
416: ksp->ops->destroy = KSPDestroy_GCR;
417: ksp->ops->view = KSPView_GCR;
418: ksp->ops->setfromoptions = KSPSetFromOptions_GCR;
419: ksp->ops->buildsolution = KSPBuildSolution_GCR;
420: ksp->ops->buildresidual = KSPBuildResidual_GCR;
422: PetscObjectComposeFunction((PetscObject)ksp,"KSPGCRSetRestart_C",KSPGCRSetRestart_GCR);
423: PetscObjectComposeFunction((PetscObject)ksp,"KSPGCRGetRestart_C",KSPGCRGetRestart_GCR);
424: PetscObjectComposeFunction((PetscObject)ksp,"KSPGCRSetModifyPC_C",KSPGCRSetModifyPC_GCR);
425: return(0);
426: }