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