Actual source code: gcr.c

petsc-3.5.1 2014-08-06
Report Typos and Errors
  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:     PCApply(pc, r, s); /* s = B^{-1} r */
 49:     MatMult(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 > 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:   MatMult(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:   MatGetVecs(A, &ctx->R, 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:   KSPDestroyDefault(ksp);
193:   return(0);
194: }

198: PetscErrorCode KSPSetFromOptions_GCR(KSP ksp)
199: {
201:   KSP_GCR        *ctx = (KSP_GCR*)ksp->data;
202:   PetscInt       restart;
203:   PetscBool      flg;

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*);

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.


332:    Options Database Keys:
333: +   -ksp_gcr_restart <restart> - the number of stored vectors to orthogonalize against

335:    Level: beginner

337:     Notes: The GCR Krylov method supports non-symmetric matrices and permits the use of a preconditioner
338:            which may vary from one iteration to the next. Users can can define a method to vary the
339:            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.
343:            Unlike GMRES and FGMRES, when using GCR, the solution and residual vector can be directly accessed at any iterate,
344:            with zero computational cost, via a call to KSPBuildSolution() and KSPBuildResidual() respectively.
345:            This implementation of GCR will only apply the stopping condition test whenever ksp->its > ksp->chknorm,
346:            where ksp->chknorm is specified via the command line argument -ksp_check_norm_iteration or via
347:            the function KSPSetCheckNormIteration().
348:            The method implemented requires the storage of 2 x restart + 1 vectors, twice as much as GMRES.
349:            Support only for right preconditioning.

351:     Contributed by Dave May

353:     References:
354:            S. C. Eisenstat, H. C. Elman, and H. C. Schultz. Variational iterative methods for
355:            non-symmetric systems of linear equations. SIAM J. Numer. Anal., 20, 345-357, 1983


358: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
359:            KSPGCRSetRestart(), KSPGCRSetModifyPC(), KSPGMRES, KSPFGMRES

361: M*/
364: PETSC_EXTERN PetscErrorCode KSPCreate_GCR(KSP ksp)
365: {
367:   KSP_GCR        *ctx;

370:   PetscNewLog(ksp,&ctx);

372:   ctx->restart    = 30;
373:   ctx->n_restarts = 0;
374:   ksp->data       = (void*)ctx;

376:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);

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,"KSPGCRSetModifyPC_C",KSPGCRSetModifyPC_GCR);
389:   return(0);
390: }