Actual source code: rich.c


  2: /*
  3:             This implements Richardson Iteration.
  4: */
  5: #include <../src/ksp/ksp/impls/rich/richardsonimpl.h>

  7: PetscErrorCode KSPSetUp_Richardson(KSP ksp)
  8: {
  9:   KSP_Richardson *richardsonP = (KSP_Richardson*)ksp->data;

 11:   if (richardsonP->selfscale) {
 12:     KSPSetWorkVecs(ksp,4);
 13:   } else {
 14:     KSPSetWorkVecs(ksp,2);
 15:   }
 16:   return 0;
 17: }

 19: PetscErrorCode  KSPSolve_Richardson(KSP ksp)
 20: {
 21:   PetscInt       i,maxit;
 22:   PetscReal      rnorm = 0.0,abr;
 23:   PetscScalar    scale,rdot;
 24:   Vec            x,b,r,z,w = NULL,y = NULL;
 25:   PetscInt       xs, ws;
 26:   Mat            Amat,Pmat;
 27:   KSP_Richardson *richardsonP = (KSP_Richardson*)ksp->data;
 28:   PetscBool      exists,diagonalscale;
 29:   MatNullSpace   nullsp;

 31:   PCGetDiagonalScale(ksp->pc,&diagonalscale);

 34:   ksp->its = 0;

 36:   PCGetOperators(ksp->pc,&Amat,&Pmat);
 37:   x    = ksp->vec_sol;
 38:   b    = ksp->vec_rhs;
 39:   VecGetSize(x,&xs);
 40:   VecGetSize(ksp->work[0],&ws);
 41:   if (xs != ws) {
 42:     if (richardsonP->selfscale) {
 43:       KSPSetWorkVecs(ksp,4);
 44:     } else {
 45:       KSPSetWorkVecs(ksp,2);
 46:     }
 47:   }
 48:   r = ksp->work[0];
 49:   z = ksp->work[1];
 50:   if (richardsonP->selfscale) {
 51:     w = ksp->work[2];
 52:     y = ksp->work[3];
 53:   }
 54:   maxit = ksp->max_it;

 56:   /* if user has provided fast Richardson code use that */
 57:   PCApplyRichardsonExists(ksp->pc,&exists);
 58:   MatGetNullSpace(Pmat,&nullsp);
 59:   if (exists && maxit > 0 && richardsonP->scale == 1.0 && (ksp->converged == KSPConvergedDefault || ksp->converged == KSPConvergedSkip) && !ksp->numbermonitors && !ksp->transpose_solve && !nullsp) {
 60:     PCRichardsonConvergedReason reason;
 61:     PCApplyRichardson(ksp->pc,b,x,r,ksp->rtol,ksp->abstol,ksp->divtol,maxit,ksp->guess_zero,&ksp->its,&reason);
 62:     ksp->reason = (KSPConvergedReason)reason;
 63:     return 0;
 64:   } else {
 65:     PetscInfo(ksp,"KSPSolve_Richardson: Warning, skipping optimized PCApplyRichardson()\n");
 66:   }

 68:   if (!ksp->guess_zero) {                          /*   r <- b - A x     */
 69:     KSP_MatMult(ksp,Amat,x,r);
 70:     VecAYPX(r,-1.0,b);
 71:   } else {
 72:     VecCopy(b,r);
 73:   }

 75:   ksp->its = 0;
 76:   if (richardsonP->selfscale) {
 77:     KSP_PCApply(ksp,r,z);         /*   z <- B r          */
 78:     for (i=0; i<maxit; i++) {

 80:       if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
 81:         VecNorm(r,NORM_2,&rnorm); /*   rnorm <- r'*r     */
 82:       } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
 83:         VecNorm(z,NORM_2,&rnorm); /*   rnorm <- z'*z     */
 84:       } else rnorm = 0.0;

 86:       KSPCheckNorm(ksp,rnorm);
 87:       ksp->rnorm = rnorm;
 88:       KSPMonitor(ksp,i,rnorm);
 89:       KSPLogResidualHistory(ksp,rnorm);
 90:       (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
 91:       if (ksp->reason) break;
 92:       KSP_PCApplyBAorAB(ksp,z,y,w); /* y = BAz = BABr */
 93:       VecDotNorm2(z,y,&rdot,&abr));   /*   rdot = (Br)^T(BABR; abr = (BABr)^T (BABr) */
 94:       scale = rdot/abr;
 95:       PetscInfo(ksp,"Self-scale factor %g\n",(double)PetscRealPart(scale));
 96:       VecAXPY(x,scale,z);   /*   x  <- x + scale z */
 97:       VecAXPY(r,-scale,w);  /*  r <- r - scale*Az */
 98:       VecAXPY(z,-scale,y);  /*  z <- z - scale*y */
 99:       ksp->its++;
100:     }
101:   } else {
102:     for (i=0; i<maxit; i++) {

104:       if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
105:         VecNorm(r,NORM_2,&rnorm); /*   rnorm <- r'*r     */
106:       } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
107:         KSP_PCApply(ksp,r,z);    /*   z <- B r          */
108:         VecNorm(z,NORM_2,&rnorm); /*   rnorm <- z'*z     */
109:       } else rnorm = 0.0;
110:       ksp->rnorm = rnorm;
111:       KSPMonitor(ksp,i,rnorm);
112:       KSPLogResidualHistory(ksp,rnorm);
113:       (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
114:       if (ksp->reason) break;
115:       if (ksp->normtype != KSP_NORM_PRECONDITIONED) {
116:         KSP_PCApply(ksp,r,z);    /*   z <- B r          */
117:       }

119:       VecAXPY(x,richardsonP->scale,z);    /*   x  <- x + scale z */
120:       ksp->its++;

122:       if (i+1 < maxit || ksp->normtype != KSP_NORM_NONE) {
123:         KSP_MatMult(ksp,Amat,x,r);      /*   r  <- b - Ax      */
124:         VecAYPX(r,-1.0,b);
125:       }
126:     }
127:   }
128:   if (!ksp->reason) {
129:     if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
130:       VecNorm(r,NORM_2,&rnorm);     /*   rnorm <- r'*r     */
131:     } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
132:       KSP_PCApply(ksp,r,z);   /*   z <- B r          */
133:       VecNorm(z,NORM_2,&rnorm);     /*   rnorm <- z'*z     */
134:     } else rnorm = 0.0;

136:     KSPCheckNorm(ksp,rnorm);
137:     ksp->rnorm = rnorm;
138:     KSPLogResidualHistory(ksp,rnorm);
139:     KSPMonitor(ksp,i,rnorm);
140:     if (ksp->its >= ksp->max_it) {
141:       if (ksp->normtype != KSP_NORM_NONE) {
142:         (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
143:         if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
144:       } else {
145:         ksp->reason = KSP_CONVERGED_ITS;
146:       }
147:     }
148:   }
149:   return 0;
150: }

152: PetscErrorCode KSPView_Richardson(KSP ksp,PetscViewer viewer)
153: {
154:   KSP_Richardson *richardsonP = (KSP_Richardson*)ksp->data;
155:   PetscBool      iascii;

157:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
158:   if (iascii) {
159:     if (richardsonP->selfscale) {
160:       PetscViewerASCIIPrintf(viewer,"  using self-scale best computed damping factor\n");
161:     } else {
162:       PetscViewerASCIIPrintf(viewer,"  damping factor=%g\n",(double)richardsonP->scale);
163:     }
164:   }
165:   return 0;
166: }

168: PetscErrorCode KSPSetFromOptions_Richardson(PetscOptionItems *PetscOptionsObject,KSP ksp)
169: {
170:   KSP_Richardson *rich = (KSP_Richardson*)ksp->data;
171:   PetscReal      tmp;
172:   PetscBool      flg,flg2;

174:   PetscOptionsHead(PetscOptionsObject,"KSP Richardson Options");
175:   PetscOptionsReal("-ksp_richardson_scale","damping factor","KSPRichardsonSetScale",rich->scale,&tmp,&flg);
176:   if (flg) KSPRichardsonSetScale(ksp,tmp);
177:   PetscOptionsBool("-ksp_richardson_self_scale","dynamically determine optimal damping factor","KSPRichardsonSetSelfScale",rich->selfscale,&flg2,&flg);
178:   if (flg) KSPRichardsonSetSelfScale(ksp,flg2);
179:   PetscOptionsTail();
180:   return 0;
181: }

183: PetscErrorCode KSPDestroy_Richardson(KSP ksp)
184: {
185:   PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetScale_C",NULL);
186:   KSPDestroyDefault(ksp);
187:   return 0;
188: }

190: static PetscErrorCode  KSPRichardsonSetScale_Richardson(KSP ksp,PetscReal scale)
191: {
192:   KSP_Richardson *richardsonP;

194:   richardsonP        = (KSP_Richardson*)ksp->data;
195:   richardsonP->scale = scale;
196:   return 0;
197: }

199: static PetscErrorCode  KSPRichardsonSetSelfScale_Richardson(KSP ksp,PetscBool selfscale)
200: {
201:   KSP_Richardson *richardsonP;

203:   richardsonP            = (KSP_Richardson*)ksp->data;
204:   richardsonP->selfscale = selfscale;
205:   return 0;
206: }

208: static PetscErrorCode KSPBuildResidual_Richardson(KSP ksp,Vec t,Vec v,Vec *V)
209: {
210:   if (ksp->normtype == KSP_NORM_NONE) {
211:     KSPBuildResidualDefault(ksp,t,v,V);
212:   } else {
213:     VecCopy(ksp->work[0],v);
214:     *V   = v;
215:   }
216:   return 0;
217: }

219: /*MC
220:      KSPRICHARDSON - The preconditioned Richardson iterative method

222:    Options Database Keys:
223: .   -ksp_richardson_scale - damping factor on the correction (defaults to 1.0)

225:    Level: beginner

227:    Notes:
228:     x^{n+1} = x^{n} + scale*B(b - A x^{n})

230:           Here B is the application of the preconditioner

232:           This method often (usually) will not converge unless scale is very small.

234:    Notes:
235:     For some preconditioners, currently SOR, the convergence test is skipped to improve speed,
236:     thus it always iterates the maximum number of iterations you've selected. When -ksp_monitor
237:     (or any other monitor) is turned on, the norm is computed at each iteration and so the convergence test is run unless
238:     you specifically call KSPSetNormType(ksp,KSP_NORM_NONE);

240:          For some preconditioners, currently PCMG and PCHYPRE with BoomerAMG if -ksp_monitor (and also
241:     any other monitor) is not turned on then the convergence test is done by the preconditioner itself and
242:     so the solver may run more or fewer iterations then if -ksp_monitor is selected.

244:     Supports only left preconditioning

246:     If using direct solvers such as PCLU and PCCHOLESKY one generally uses KSPPREONLY which uses exactly one iteration

248: $    -ksp_type richardson -pc_type jacobi gives one classically Jacobi preconditioning

250:   References:
251: . * - L. F. Richardson, "The Approximate Arithmetical Solution by Finite Differences of Physical Problems Involving
252:    Differential Equations, with an Application to the Stresses in a Masonry Dam",
253:   Philosophical Transactions of the Royal Society of London. Series A,
254:   Containing Papers of a Mathematical or Physical Character, Vol. 210, 1911 (1911).

256: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
257:            KSPRichardsonSetScale(), KSPPREONLY

259: M*/

261: PETSC_EXTERN PetscErrorCode KSPCreate_Richardson(KSP ksp)
262: {
263:   KSP_Richardson *richardsonP;

265:   PetscNewLog(ksp,&richardsonP);
266:   ksp->data = (void*)richardsonP;

268:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
269:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
270:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);

272:   ksp->ops->setup          = KSPSetUp_Richardson;
273:   ksp->ops->solve          = KSPSolve_Richardson;
274:   ksp->ops->destroy        = KSPDestroy_Richardson;
275:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
276:   ksp->ops->buildresidual  = KSPBuildResidual_Richardson;
277:   ksp->ops->view           = KSPView_Richardson;
278:   ksp->ops->setfromoptions = KSPSetFromOptions_Richardson;

280:   PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetScale_C",KSPRichardsonSetScale_Richardson);
281:   PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetSelfScale_C",KSPRichardsonSetSelfScale_Richardson);

283:   richardsonP->scale = 1.0;
284:   return 0;
285: }