Actual source code: gcr.c
petsc-3.3-p7 2013-05-11
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 */
14: void *modifypc_ctx; /* user defined data for the modifypc function */
15: } KSP_GCR;
19: PetscErrorCode KSPSolve_GCR_cycle( KSP ksp )
20: {
21: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
23: PetscScalar nrm,r_dot_v;
24: Mat A, B;
25: PC pc;
26: Vec s,v,r;
27: PetscReal norm_r;
28: PetscInt k, i, restart;
29: Vec x;
30: PetscReal res;
31:
33: restart = ctx->restart;
34: KSPGetPC( ksp, &pc );
35: KSPGetOperators( ksp, &A, &B, 0 );
36:
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: }
46:
47: PCApply( pc, r, s ); /* s = B^{-1} r */
48: MatMult( A, s, v ); /* v = A s */
49:
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 = PetscSqrtScalar(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 ) {
63: VecNorm( r, NORM_2, &norm_r );
64: }
65: /* update the local counter and the global counter */
66: ksp->its++;
67: res = norm_r;
68: ksp->rnorm = res;
69:
70: KSPLogResidualHistory(ksp,res);
71: KSPMonitor(ksp,ksp->its,res);
72:
73: if( ksp->its > ksp->chknorm ) {
74: (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);
75: if (ksp->reason) break;
76: }
77:
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: }
89: PetscErrorCode KSPSolve_GCR( KSP ksp )
90: {
91: KSP_GCR *ctx = (KSP_GCR*)ksp->data;
93: Mat A, B;
94: Vec r,b,x;
95: PetscReal norm_r;
96:
98: KSPGetOperators( ksp, &A, &B, PETSC_NULL );
99: x = ksp->vec_sol;
100: b = ksp->vec_rhs;
101: r = ctx->R;
102:
103: /* compute initial residual */
104: MatMult( A, x, r );
105: VecAYPX( r, -1.0, b ); /* r = b - A x */
106: VecNorm( r, NORM_2, &norm_r );
107:
108: ksp->its = 0;
109: ksp->rnorm0 = norm_r;
110:
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);
115:
116: do {
117: KSPSolve_GCR_cycle( ksp );
118: if (ksp->reason) break; /* catch case when convergence occurs inside the cycle */
119: } while( ksp->its < ksp->max_it );
120: if( ksp->its >= ksp->max_it) {
121: ksp->reason = KSP_DIVERGED_ITS;
122: }
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(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
157: KSPGetOperators( ksp, &A, 0, 0 );
158: MatGetVecs( A, &ctx->R, PETSC_NULL );
159: VecDuplicateVecs( ctx->R, ctx->restart, &ctx->VV );
160: VecDuplicateVecs( ctx->R, ctx->restart, &ctx->SS );
162: PetscMalloc( sizeof(PetscScalar)*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: KSPDefaultDestroy(ksp);
193: return(0);
194: }
198: PetscErrorCode KSPSetFromOptions_GCR(KSP ksp)
199: {
200: PetscErrorCode ierr;
201: KSP_GCR *ctx = (KSP_GCR *)ksp->data;
202: PetscInt restart;
203: PetscBool flg;
204:
206: PetscOptionsHead("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*);
217: EXTERN_C_BEGIN
220: PetscErrorCode KSPGCRSetModifyPC_GCR(KSP ksp,KSPGCRModifyPCFunction function,void *data,KSPGCRDestroyFunction destroy)
221: {
222: KSP_GCR *ctx = (KSP_GCR *)ksp->data;
223:
226: ctx->modifypc = function;
227: ctx->modifypc_destroy = destroy;
228: ctx->modifypc_ctx = data;
229: return(0);
230: }
231: EXTERN_C_END
235: /*@C
236: KSPGCRSetModifyPC - Sets the routine used by GCR to modify the preconditioner.
237:
238: Logically Collective on KSP
239:
240: Input Parameters:
241: + ksp - iterative context obtained from KSPCreate()
242: . function - user defined function to modify the preconditioner
243: . ctx - user provided contex for the modify preconditioner function
244: - destroy - the function to use to destroy the user provided application context.
245:
246: Calling Sequence of function:
247: PetscErrorCode function ( KSP ksp, PetscInt n, PetscReal rnorm, void *ctx )
248:
249: ksp - iterative context
250: n - the total number of GCR iterations that have occurred
251: rnorm - 2-norm residual value
252: ctx - the user provided application context
253:
254: Level: intermediate
255:
256: Notes:
257: The default modifypc routine is KSPGCRModifyPCNoChange()
258:
259: .seealso: KSPGCRModifyPCNoChange()
260:
261: @*/
262: PetscErrorCode KSPGCRSetModifyPC(KSP ksp,PetscErrorCode (*function)(KSP,PetscInt,PetscReal,void*),void *data,PetscErrorCode (*destroy)(void*))
263: {
265:
267: PetscUseMethod(ksp,"KSPGCRSetModifyPC_C",(KSP,PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*),void *data,PetscErrorCode (*)(void*)),(ksp,function,data,destroy));
268: return(0);
269: }
273: PetscErrorCode KSPGCRSetRestart_GCR(KSP ksp,PetscInt restart)
274: {
275: KSP_GCR *ctx;
276:
278: ctx = (KSP_GCR *)ksp->data;
279: ctx->restart = restart;
280: return(0);
281: }
285: PetscErrorCode KSPGCRSetRestart(KSP ksp, PetscInt restart)
286: {
290: PetscTryMethod(ksp,"KSPGCRSetRestart_C",(KSP,PetscInt),(ksp,restart));
291: return(0);
292: }
296: PetscErrorCode KSPBuildSolution_GCR(KSP ksp, Vec v, Vec *V)
297: {
299: Vec x;
300:
302: x = ksp->vec_sol;
303: if (v) {
304: VecCopy( x, v );
305: if (V) *V = v;
306: } else if (V) {
307: *V = ksp->vec_sol;
308: }
309: return(0);
310: }
314: PetscErrorCode KSPBuildResidual_GCR(KSP ksp, Vec t, Vec v, Vec *V)
315: {
317: KSP_GCR *ctx;
318:
320: ctx = (KSP_GCR *)ksp->data;
321: if (v) {
322: VecCopy( ctx->R, v );
323: if (V) *V = v;
324: } else if (V) {
325: *V = ctx->R;
326: }
327: return(0);
328: }
330: /*MC
331: KSPGCR - Implements the preconditioned Generalized Conjugate Residual method.
332:
333:
334: Options Database Keys:
335: + -ksp_gcr_restart <restart> - the number of stored vectors to orthogonalize against
336:
337: Level: beginner
338:
339: Notes: The GCR Krylov method supports non-symmetric matrices and permits the use of a preconditioner
340: which may vary from one iteration to the next. Users can can define a method to vary the
341: preconditioner between iterates via KSPGCRSetModifyPC().
342: Restarts are solves with x0 not equal to zero. When a restart occurs, the initial starting
343: solution is given by the current estimate for x which was obtained by the last restart
344: iterations of the GCR algorithm.
345: Unlike GMRES and FGMRES, when using GCR, the solution and residual vector can be directly accessed at any iterate,
346: 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().
350: The method implemented requires the storage of 2 x restart + 1 vectors, twice as much as GMRES.
351: Support only for right preconditioning.
353: Contributed by Dave May
354:
355: References:
356: S. C. Eisenstat, H. C. Elman, and H. C. Schultz. Variational iterative methods for
357: non-symmetric systems of linear equations. SIAM J. Numer. Anal., 20, 345-357, 1983
358:
359:
360: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
361: KSPGCRSetRestart(), KSPGCRSetModifyPC(), KSPGMRES, KSPFGMRES
362:
363: M*/
364: EXTERN_C_BEGIN
367: PetscErrorCode KSPCreate_GCR(KSP ksp)
368: {
370: KSP_GCR *ctx;
373: PetscNewLog(ksp,KSP_GCR,&ctx);
374: ctx->restart = 30;
375: ctx->n_restarts = 0;
376: ksp->data = (void*)ctx;
377: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);
379: ksp->ops->setup = KSPSetUp_GCR;
380: ksp->ops->solve = KSPSolve_GCR;
381: ksp->ops->reset = KSPReset_GCR;
382: ksp->ops->destroy = KSPDestroy_GCR;
383: ksp->ops->view = KSPView_GCR;
384: ksp->ops->setfromoptions = KSPSetFromOptions_GCR;
385: ksp->ops->buildsolution = KSPBuildSolution_GCR;
386: ksp->ops->buildresidual = KSPBuildResidual_GCR;
388: PetscObjectComposeFunctionDynamic( (PetscObject)ksp, "KSPGCRSetRestart_C",
389: "KSPGCRSetRestart_GCR",KSPGCRSetRestart_GCR );
390: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPGCRSetModifyPC_C",
391: "KSPGCRSetModifyPC_GCR",KSPGCRSetModifyPC_GCR);
392: return(0);
393: }
394: EXTERN_C_END