Actual source code: gcr.c
petsc-3.7.5 2017-01-01
2: #include <petscksp.h>
3: #include <petsc/private/kspimpl.h>
5: typedef struct {
6: PetscInt restart;
7: PetscInt n_restarts;
8: PetscScalar *val;
9: Vec *VV, *SS;
10: Vec R;
12: PetscErrorCode (*modifypc)(KSP,PetscInt,PetscReal,void*); /* function to modify the preconditioner*/
13: PetscErrorCode (*modifypc_destroy)(void*); /* function to destroy the user context for the modifypc function */
15: void *modifypc_ctx; /* user defined data for the modifypc function */
16: } KSP_GCR;
20: PetscErrorCode KSPSolve_GCR_cycle(KSP ksp)
21: {
22: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
24: PetscScalar r_dot_v;
25: Mat A, B;
26: PC pc;
27: Vec s,v,r;
28: PetscReal norm_r,nrm;
29: PetscInt k, i, restart;
30: Vec x;
31: PetscReal res;
34: restart = ctx->restart;
35: KSPGetPC(ksp, &pc);
36: KSPGetOperators(ksp, &A, &B);
38: x = ksp->vec_sol;
39: r = ctx->R;
41: for (k=0; k<restart; k++) {
42: v = ctx->VV[k];
43: s = ctx->SS[k];
44: if (ctx->modifypc) {
45: (*ctx->modifypc)(ksp,ksp->its,ksp->rnorm,ctx->modifypc_ctx);
46: }
48: KSP_PCApply(ksp, r, s); /* s = B^{-1} r */
49: KSP_MatMult(ksp,A, s, v); /* v = A s */
51: VecMDot(v,k, ctx->VV, ctx->val);
52: for (i=0; i<k; i++) ctx->val[i] = -ctx->val[i];
53: VecMAXPY(v,k,ctx->val,ctx->VV); /* v = v - sum_{i=0}^{k-1} alpha_i v_i */
54: VecMAXPY(s,k,ctx->val,ctx->SS); /* s = s - sum_{i=0}^{k-1} alpha_i s_i */
56: VecDotNorm2(r,v,&r_dot_v,&nrm);
57: nrm = PetscSqrtReal(nrm);
58: r_dot_v = r_dot_v/nrm;
59: VecScale(v, 1.0/nrm);
60: VecScale(s, 1.0/nrm);
61: VecAXPY(x, r_dot_v, s);
62: VecAXPY(r, -r_dot_v, v);
63: if (ksp->its > ksp->chknorm) {
64: VecNorm(r, NORM_2, &norm_r);
65: }
66: /* update the local counter and the global counter */
67: ksp->its++;
68: res = norm_r;
69: ksp->rnorm = res;
71: KSPLogResidualHistory(ksp,res);
72: KSPMonitor(ksp,ksp->its,res);
74: if (ksp->its-1 > ksp->chknorm) {
75: (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);
76: if (ksp->reason) break;
77: }
79: if (ksp->its >= ksp->max_it) {
80: ksp->reason = KSP_CONVERGED_ITS;
81: break;
82: }
83: }
84: ctx->n_restarts++;
85: return(0);
86: }
90: PetscErrorCode KSPSolve_GCR(KSP ksp)
91: {
92: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
94: Mat A, B;
95: Vec r,b,x;
96: PetscReal norm_r;
99: KSPGetOperators(ksp, &A, &B);
100: x = ksp->vec_sol;
101: b = ksp->vec_rhs;
102: r = ctx->R;
104: /* compute initial residual */
105: KSP_MatMult(ksp,A, x, r);
106: VecAYPX(r, -1.0, b); /* r = b - A x */
107: VecNorm(r, NORM_2, &norm_r);
109: ksp->its = 0;
110: ksp->rnorm0 = norm_r;
112: KSPLogResidualHistory(ksp,ksp->rnorm0);
113: KSPMonitor(ksp,ksp->its,ksp->rnorm0);
114: (*ksp->converged)(ksp,ksp->its,ksp->rnorm0,&ksp->reason,ksp->cnvP);
115: if (ksp->reason) return(0);
117: do {
118: KSPSolve_GCR_cycle(ksp);
119: if (ksp->reason) break; /* catch case when convergence occurs inside the cycle */
120: } while (ksp->its < ksp->max_it);
122: if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
123: return(0);
124: }
128: PetscErrorCode KSPView_GCR(KSP ksp, PetscViewer viewer)
129: {
130: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
132: PetscBool iascii;
135: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
136: if (iascii) {
137: PetscViewerASCIIPrintf(viewer," GCR: restart = %D \n", ctx->restart);
138: PetscViewerASCIIPrintf(viewer," GCR: restarts performed = %D \n", ctx->n_restarts);
139: }
140: return(0);
141: }
146: PetscErrorCode KSPSetUp_GCR(KSP ksp)
147: {
148: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
150: Mat A;
151: PetscBool diagonalscale;
154: PCGetDiagonalScale(ksp->pc,&diagonalscale);
155: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
157: KSPGetOperators(ksp, &A, NULL);
158: MatCreateVecs(A, &ctx->R, NULL);
159: VecDuplicateVecs(ctx->R, ctx->restart, &ctx->VV);
160: VecDuplicateVecs(ctx->R, ctx->restart, &ctx->SS);
162: PetscMalloc1(ctx->restart, &ctx->val);
163: return(0);
164: }
168: PetscErrorCode KSPReset_GCR(KSP ksp)
169: {
171: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
174: VecDestroy(&ctx->R);
175: VecDestroyVecs(ctx->restart,&ctx->VV);
176: VecDestroyVecs(ctx->restart,&ctx->SS);
177: if (ctx->modifypc_destroy) {
178: (*ctx->modifypc_destroy)(ctx->modifypc_ctx);
179: }
180: PetscFree(ctx->val);
181: return(0);
182: }
186: PetscErrorCode KSPDestroy_GCR(KSP ksp)
187: {
191: KSPReset_GCR(ksp);
192: KSPDestroyDefault(ksp);
193: return(0);
194: }
198: PetscErrorCode KSPSetFromOptions_GCR(PetscOptionItems *PetscOptionsObject,KSP ksp)
199: {
201: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
202: PetscInt restart;
203: PetscBool flg;
206: PetscOptionsHead(PetscOptionsObject,"KSP GCR options");
207: PetscOptionsInt("-ksp_gcr_restart","Number of Krylov search directions","KSPGCRSetRestart",ctx->restart,&restart,&flg);
208: if (flg) { KSPGCRSetRestart(ksp,restart); }
209: PetscOptionsTail();
210: return(0);
211: }
214: typedef PetscErrorCode (*KSPGCRModifyPCFunction)(KSP,PetscInt,PetscReal,void*);
215: typedef PetscErrorCode (*KSPGCRDestroyFunction)(void*);
219: static PetscErrorCode KSPGCRSetModifyPC_GCR(KSP ksp,KSPGCRModifyPCFunction function,void *data,KSPGCRDestroyFunction destroy)
220: {
221: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
225: ctx->modifypc = function;
226: ctx->modifypc_destroy = destroy;
227: ctx->modifypc_ctx = data;
228: return(0);
229: }
233: /*@C
234: KSPGCRSetModifyPC - Sets the routine used by GCR to modify the preconditioner.
236: Logically Collective on KSP
238: Input Parameters:
239: + ksp - iterative context obtained from KSPCreate()
240: . function - user defined function to modify the preconditioner
241: . ctx - user provided contex for the modify preconditioner function
242: - destroy - the function to use to destroy the user provided application context.
244: Calling Sequence of function:
245: PetscErrorCode function (KSP ksp, PetscInt n, PetscReal rnorm, void *ctx)
247: ksp - iterative context
248: n - the total number of GCR iterations that have occurred
249: rnorm - 2-norm residual value
250: ctx - the user provided application context
252: Level: intermediate
254: Notes:
255: The default modifypc routine is KSPGCRModifyPCNoChange()
257: .seealso: KSPGCRModifyPCNoChange()
259: @*/
260: PetscErrorCode KSPGCRSetModifyPC(KSP ksp,PetscErrorCode (*function)(KSP,PetscInt,PetscReal,void*),void *data,PetscErrorCode (*destroy)(void*))
261: {
265: PetscUseMethod(ksp,"KSPGCRSetModifyPC_C",(KSP,PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*),void *data,PetscErrorCode (*)(void*)),(ksp,function,data,destroy));
266: return(0);
267: }
271: static PetscErrorCode KSPGCRSetRestart_GCR(KSP ksp,PetscInt restart)
272: {
273: KSP_GCR *ctx;
276: ctx = (KSP_GCR*)ksp->data;
277: ctx->restart = restart;
278: return(0);
279: }
283: PetscErrorCode KSPGCRSetRestart(KSP ksp, PetscInt restart)
284: {
288: PetscTryMethod(ksp,"KSPGCRSetRestart_C",(KSP,PetscInt),(ksp,restart));
289: return(0);
290: }
294: PetscErrorCode KSPBuildSolution_GCR(KSP ksp, Vec v, Vec *V)
295: {
297: Vec x;
300: x = ksp->vec_sol;
301: if (v) {
302: VecCopy(x, v);
303: if (V) *V = v;
304: } else if (V) {
305: *V = ksp->vec_sol;
306: }
307: return(0);
308: }
312: PetscErrorCode KSPBuildResidual_GCR(KSP ksp, Vec t, Vec v, Vec *V)
313: {
315: KSP_GCR *ctx;
318: ctx = (KSP_GCR*)ksp->data;
319: if (v) {
320: VecCopy(ctx->R, v);
321: if (V) *V = v;
322: } else if (V) {
323: *V = ctx->R;
324: }
325: return(0);
326: }
328: /*MC
329: KSPGCR - Implements the preconditioned Generalized Conjugate Residual method.
331: Options Database Keys:
332: . -ksp_gcr_restart <restart> - the number of stored vectors to orthogonalize against
334: Level: beginner
336: Notes: 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().
339: Restarts are solves with x0 not equal to zero. When a restart occurs, the initial starting
340: solution is given by the current estimate for x which was obtained by the last restart
341: iterations of the GCR algorithm.
342: Unlike GMRES and FGMRES, when using GCR, the solution and residual vector can be directly accessed at any iterate,
343: with zero computational cost, via a call to KSPBuildSolution() and KSPBuildResidual() respectively.
344: This implementation of GCR will only apply the stopping condition test whenever ksp->its > ksp->chknorm,
345: where ksp->chknorm is specified via the command line argument -ksp_check_norm_iteration or via
346: the function KSPSetCheckNormIteration().
347: The method implemented requires the storage of 2 x restart + 1 vectors, twice as much as GMRES.
348: Support only for right preconditioning.
350: Contributed by Dave May
352: References:
353: . 1. - S. C. Eisenstat, H. C. Elman, and H. C. Schultz. Variational iterative methods for
354: nonsymmetric systems of linear equations. SIAM J. Numer. Anal., 20, 1983
357: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
358: KSPGCRSetRestart(), KSPGCRSetModifyPC(), KSPGMRES, KSPFGMRES
360: M*/
363: PETSC_EXTERN PetscErrorCode KSPCreate_GCR(KSP ksp)
364: {
366: 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_UNPRECONDITIONED,PC_RIGHT,3);
377: ksp->ops->setup = KSPSetUp_GCR;
378: ksp->ops->solve = KSPSolve_GCR;
379: ksp->ops->reset = KSPReset_GCR;
380: ksp->ops->destroy = KSPDestroy_GCR;
381: ksp->ops->view = KSPView_GCR;
382: ksp->ops->setfromoptions = KSPSetFromOptions_GCR;
383: ksp->ops->buildsolution = KSPBuildSolution_GCR;
384: ksp->ops->buildresidual = KSPBuildResidual_GCR;
386: PetscObjectComposeFunction((PetscObject)ksp, "KSPGCRSetRestart_C",KSPGCRSetRestart_GCR);
387: PetscObjectComposeFunction((PetscObject)ksp,"KSPGCRSetModifyPC_C",KSPGCRSetModifyPC_GCR);
388: return(0);
389: }