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) {
42: (*ctx->modifypc)(ksp,ksp->its,ksp->rnorm,ctx->modifypc_ctx);
43: }
45: KSP_PCApply(ksp, r, s); /* s = B^{-1} r */
46: KSP_MatMult(ksp,A, s, v); /* v = A s */
48: VecMDot(v,k, ctx->VV, ctx->val);
49: for (i=0; i<k; i++) ctx->val[i] = -ctx->val[i];
50: VecMAXPY(v,k,ctx->val,ctx->VV); /* v = v - sum_{i=0}^{k-1} alpha_i v_i */
51: VecMAXPY(s,k,ctx->val,ctx->SS); /* s = s - sum_{i=0}^{k-1} alpha_i s_i */
53: VecDotNorm2(r,v,&r_dot_v,&nrm);
54: nrm = PetscSqrtReal(nrm);
55: r_dot_v = r_dot_v/nrm;
56: VecScale(v, 1.0/nrm);
57: VecScale(s, 1.0/nrm);
58: VecAXPY(x, r_dot_v, s);
59: VecAXPY(r, -r_dot_v, v);
60: if (ksp->its > ksp->chknorm && ksp->normtype != KSP_NORM_NONE) {
61: VecNorm(r, NORM_2, &norm_r);
62: KSPCheckNorm(ksp,norm_r);
63: }
64: /* update the local counter and the global counter */
65: ksp->its++;
66: ksp->rnorm = norm_r;
68: KSPLogResidualHistory(ksp,norm_r);
69: KSPMonitor(ksp,ksp->its,norm_r);
71: if (ksp->its-1 > ksp->chknorm) {
72: (*ksp->converged)(ksp,ksp->its,norm_r,&ksp->reason,ksp->cnvP);
73: if (ksp->reason) break;
74: }
76: if (ksp->its >= ksp->max_it) {
77: ksp->reason = KSP_CONVERGED_ITS;
78: break;
79: }
80: }
81: ctx->n_restarts++;
82: return 0;
83: }
85: static PetscErrorCode KSPSolve_GCR(KSP ksp)
86: {
87: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
88: Mat A, B;
89: Vec r,b,x;
90: PetscReal norm_r = 0.0;
92: KSPGetOperators(ksp, &A, &B);
93: x = ksp->vec_sol;
94: b = ksp->vec_rhs;
95: r = ctx->R;
97: /* compute initial residual */
98: KSP_MatMult(ksp,A, x, r);
99: VecAYPX(r, -1.0, b); /* r = b - A x */
100: if (ksp->normtype != KSP_NORM_NONE) {
101: VecNorm(r, NORM_2, &norm_r);
102: KSPCheckNorm(ksp,norm_r);
103: }
104: ksp->its = 0;
105: ksp->rnorm0 = norm_r;
107: KSPLogResidualHistory(ksp,ksp->rnorm0);
108: KSPMonitor(ksp,ksp->its,ksp->rnorm0);
109: (*ksp->converged)(ksp,ksp->its,ksp->rnorm0,&ksp->reason,ksp->cnvP);
110: if (ksp->reason) return 0;
112: do {
113: KSPSolve_GCR_cycle(ksp);
114: if (ksp->reason) return 0; /* 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: return 0;
119: }
121: static PetscErrorCode KSPView_GCR(KSP ksp, PetscViewer viewer)
122: {
123: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
124: PetscBool iascii;
126: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
127: if (iascii) {
128: PetscViewerASCIIPrintf(viewer," restart = %D \n", ctx->restart);
129: PetscViewerASCIIPrintf(viewer," restarts performed = %D \n", ctx->n_restarts);
130: }
131: return 0;
132: }
134: static PetscErrorCode KSPSetUp_GCR(KSP ksp)
135: {
136: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
137: Mat A;
138: PetscBool diagonalscale;
140: PCGetDiagonalScale(ksp->pc,&diagonalscale);
143: KSPGetOperators(ksp, &A, NULL);
144: MatCreateVecs(A, &ctx->R, NULL);
145: VecDuplicateVecs(ctx->R, ctx->restart, &ctx->VV);
146: VecDuplicateVecs(ctx->R, ctx->restart, &ctx->SS);
148: PetscMalloc1(ctx->restart, &ctx->val);
149: return 0;
150: }
152: static PetscErrorCode KSPReset_GCR(KSP ksp)
153: {
154: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
156: VecDestroy(&ctx->R);
157: VecDestroyVecs(ctx->restart,&ctx->VV);
158: VecDestroyVecs(ctx->restart,&ctx->SS);
159: if (ctx->modifypc_destroy) {
160: (*ctx->modifypc_destroy)(ctx->modifypc_ctx);
161: }
162: PetscFree(ctx->val);
163: return 0;
164: }
166: static PetscErrorCode KSPDestroy_GCR(KSP ksp)
167: {
168: KSPReset_GCR(ksp);
169: KSPDestroyDefault(ksp);
170: PetscObjectComposeFunction((PetscObject)ksp,"KSPGCRSetRestart_C",NULL);
171: PetscObjectComposeFunction((PetscObject)ksp,"KSPGCRGetRestart_C",NULL);
172: PetscObjectComposeFunction((PetscObject)ksp,"KSPGCRSetModifyPC_C",NULL);
173: return 0;
174: }
176: static PetscErrorCode KSPSetFromOptions_GCR(PetscOptionItems *PetscOptionsObject,KSP ksp)
177: {
178: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
179: PetscInt restart;
180: PetscBool flg;
182: PetscOptionsHead(PetscOptionsObject,"KSP GCR options");
183: PetscOptionsInt("-ksp_gcr_restart","Number of Krylov search directions","KSPGCRSetRestart",ctx->restart,&restart,&flg);
184: if (flg) KSPGCRSetRestart(ksp,restart);
185: PetscOptionsTail();
186: return 0;
187: }
190: typedef PetscErrorCode (*KSPGCRModifyPCFunction)(KSP,PetscInt,PetscReal,void*);
191: typedef PetscErrorCode (*KSPGCRDestroyFunction)(void*);
193: static PetscErrorCode KSPGCRSetModifyPC_GCR(KSP ksp,KSPGCRModifyPCFunction function,void *data,KSPGCRDestroyFunction destroy)
194: {
195: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
198: ctx->modifypc = function;
199: ctx->modifypc_destroy = destroy;
200: ctx->modifypc_ctx = data;
201: return 0;
202: }
204: /*@C
205: KSPGCRSetModifyPC - Sets the routine used by GCR to modify the preconditioner.
207: Logically Collective on ksp
209: Input Parameters:
210: + ksp - iterative context obtained from KSPCreate()
211: . function - user defined function to modify the preconditioner
212: . ctx - user provided context for the modify preconditioner function
213: - destroy - the function to use to destroy the user provided application context.
215: Calling Sequence of function:
216: PetscErrorCode function (KSP ksp, PetscInt n, PetscReal rnorm, void *ctx)
218: ksp - iterative context
219: n - the total number of GCR iterations that have occurred
220: rnorm - 2-norm residual value
221: ctx - the user provided application context
223: Level: intermediate
225: Notes:
226: The default modifypc routine is KSPGCRModifyPCNoChange()
228: .seealso: KSPGCRModifyPCNoChange()
230: @*/
231: PetscErrorCode KSPGCRSetModifyPC(KSP ksp,PetscErrorCode (*function)(KSP,PetscInt,PetscReal,void*),void *data,PetscErrorCode (*destroy)(void*))
232: {
233: PetscUseMethod(ksp,"KSPGCRSetModifyPC_C",(KSP,PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*),void *data,PetscErrorCode (*)(void*)),(ksp,function,data,destroy));
234: return 0;
235: }
237: static PetscErrorCode KSPGCRSetRestart_GCR(KSP ksp,PetscInt restart)
238: {
239: KSP_GCR *ctx;
241: ctx = (KSP_GCR*)ksp->data;
242: ctx->restart = restart;
243: return 0;
244: }
246: static PetscErrorCode KSPGCRGetRestart_GCR(KSP ksp,PetscInt *restart)
247: {
248: KSP_GCR *ctx;
250: ctx = (KSP_GCR*)ksp->data;
251: *restart = ctx->restart;
252: return 0;
253: }
255: /*@
256: KSPGCRSetRestart - Sets number of iterations at which GCR restarts.
258: Not Collective
260: Input Parameters:
261: + ksp - the Krylov space context
262: - restart - integer restart value
264: Note: The default value is 30.
266: Level: intermediate
268: .seealso: KSPSetTolerances(), KSPGCRGetRestart(), KSPGMRESSetRestart()
269: @*/
270: PetscErrorCode KSPGCRSetRestart(KSP ksp, PetscInt restart)
271: {
272: PetscTryMethod(ksp,"KSPGCRSetRestart_C",(KSP,PetscInt),(ksp,restart));
273: return 0;
274: }
276: /*@
277: KSPGCRGetRestart - Gets number of iterations at which GCR restarts.
279: Not Collective
281: Input Parameter:
282: . ksp - the Krylov space context
284: Output Parameter:
285: . restart - integer restart value
287: Note: The default value is 30.
289: Level: intermediate
291: .seealso: KSPSetTolerances(), KSPGCRSetRestart(), KSPGMRESGetRestart()
292: @*/
293: PetscErrorCode KSPGCRGetRestart(KSP ksp, PetscInt *restart)
294: {
295: PetscTryMethod(ksp,"KSPGCRGetRestart_C",(KSP,PetscInt*),(ksp,restart));
296: return 0;
297: }
299: static PetscErrorCode KSPBuildSolution_GCR(KSP ksp, Vec v, Vec *V)
300: {
301: Vec x;
303: x = ksp->vec_sol;
304: if (v) {
305: VecCopy(x, v);
306: if (V) *V = v;
307: } else if (V) {
308: *V = ksp->vec_sol;
309: }
310: return 0;
311: }
313: static PetscErrorCode KSPBuildResidual_GCR(KSP ksp, Vec t, Vec v, Vec *V)
314: {
315: KSP_GCR *ctx;
317: ctx = (KSP_GCR*)ksp->data;
318: if (v) {
319: VecCopy(ctx->R, v);
320: if (V) *V = v;
321: } else if (V) {
322: *V = ctx->R;
323: }
324: return 0;
325: }
327: /*MC
328: KSPGCR - Implements the preconditioned Generalized Conjugate Residual method.
330: Options Database Keys:
331: . -ksp_gcr_restart <restart> - the number of stored vectors to orthogonalize against
333: Level: beginner
335: Notes:
336: The GCR Krylov method supports non-symmetric matrices and permits the use of a preconditioner
337: which may vary from one iteration to the next. Users can can define a method to vary the
338: preconditioner between iterates via KSPGCRSetModifyPC().
340: Restarts are solves with x0 not equal to zero. When a restart occurs, the initial starting
341: solution is given by the current estimate for x which was obtained by the last restart
342: iterations of the GCR algorithm.
344: Unlike GMRES and FGMRES, when using GCR, the solution and residual vector can be directly accessed at any iterate,
345: with zero computational cost, via a call to KSPBuildSolution() and KSPBuildResidual() respectively.
347: This implementation of GCR will only apply the stopping condition test whenever ksp->its > ksp->chknorm,
348: where ksp->chknorm is specified via the command line argument -ksp_check_norm_iteration or via
349: the function KSPSetCheckNormIteration(). Hence the residual norm reported by the monitor and stored
350: in the residual history will be listed as 0.0 before this iteration. It is actually not 0.0; just not calculated.
352: The method implemented requires the storage of 2 x restart + 1 vectors, twice as much as GMRES.
353: Support only for right preconditioning.
355: Contributed by Dave May
357: References:
358: . * - S. C. Eisenstat, H. C. Elman, and H. C. Schultz. Variational iterative methods for
359: nonsymmetric systems of linear equations. SIAM J. Numer. Anal., 20, 1983
361: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
362: KSPGCRSetRestart(), KSPGCRSetModifyPC(), KSPGMRES, KSPFGMRES
364: M*/
365: PETSC_EXTERN PetscErrorCode KSPCreate_GCR(KSP ksp)
366: {
367: KSP_GCR *ctx;
369: PetscNewLog(ksp,&ctx);
371: ctx->restart = 30;
372: ctx->n_restarts = 0;
373: ksp->data = (void*)ctx;
375: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);
376: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,3);
378: ksp->ops->setup = KSPSetUp_GCR;
379: ksp->ops->solve = KSPSolve_GCR;
380: ksp->ops->reset = KSPReset_GCR;
381: ksp->ops->destroy = KSPDestroy_GCR;
382: ksp->ops->view = KSPView_GCR;
383: ksp->ops->setfromoptions = KSPSetFromOptions_GCR;
384: ksp->ops->buildsolution = KSPBuildSolution_GCR;
385: ksp->ops->buildresidual = KSPBuildResidual_GCR;
387: PetscObjectComposeFunction((PetscObject)ksp,"KSPGCRSetRestart_C",KSPGCRSetRestart_GCR);
388: PetscObjectComposeFunction((PetscObject)ksp,"KSPGCRGetRestart_C",KSPGCRGetRestart_GCR);
389: PetscObjectComposeFunction((PetscObject)ksp,"KSPGCRSetModifyPC_C",KSPGCRSetModifyPC_GCR);
390: return 0;
391: }