Actual source code: rich.c

petsc-3.6.0 2015-06-09
Report Typos and Errors
  2: /*
  3:             This implements Richardson Iteration.
  4: */
  5: #include <petsc/private/kspimpl.h>              /*I "petscksp.h" I*/
  6: #include <../src/ksp/ksp/impls/rich/richardsonimpl.h>

 10: PetscErrorCode KSPSetUp_Richardson(KSP ksp)
 11: {
 13:   KSP_Richardson *richardsonP = (KSP_Richardson*)ksp->data;

 16:   if (richardsonP->selfscale) {
 17:     KSPSetWorkVecs(ksp,4);
 18:   } else {
 19:     KSPSetWorkVecs(ksp,2);
 20:   }
 21:   return(0);
 22: }

 26: PetscErrorCode  KSPSolve_Richardson(KSP ksp)
 27: {
 29:   PetscInt       i,maxit;
 30:   PetscReal      rnorm = 0.0,abr;
 31:   PetscScalar    scale,rdot;
 32:   Vec            x,b,r,z,w = NULL,y = NULL;
 33:   PetscInt       xs, ws;
 34:   Mat            Amat,Pmat;
 35:   KSP_Richardson *richardsonP = (KSP_Richardson*)ksp->data;
 36:   PetscBool      exists,diagonalscale;
 37:   MatNullSpace   nullsp;

 40:   PCGetDiagonalScale(ksp->pc,&diagonalscale);
 41:   if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);

 43:   ksp->its = 0;

 45:   PCGetOperators(ksp->pc,&Amat,&Pmat);
 46:   x    = ksp->vec_sol;
 47:   b    = ksp->vec_rhs;
 48:   VecGetSize(x,&xs);
 49:   VecGetSize(ksp->work[0],&ws);
 50:   if (xs != ws) {
 51:     if (richardsonP->selfscale) {
 52:       KSPSetWorkVecs(ksp,4);
 53:     } else {
 54:       KSPSetWorkVecs(ksp,2);
 55:     }
 56:   }
 57:   r = ksp->work[0];
 58:   z = ksp->work[1];
 59:   if (richardsonP->selfscale) {
 60:     w = ksp->work[2];
 61:     y = ksp->work[3];
 62:   }
 63:   maxit = ksp->max_it;

 65:   /* if user has provided fast Richardson code use that */
 66:   PCApplyRichardsonExists(ksp->pc,&exists);
 67:   MatGetNullSpace(Pmat,&nullsp);
 68:   if (exists && richardsonP->scale == 1.0 && !ksp->numbermonitors && !ksp->transpose_solve & !nullsp) {
 69:     PCRichardsonConvergedReason reason;
 70:     PCApplyRichardson(ksp->pc,b,x,r,ksp->rtol,ksp->abstol,ksp->divtol,maxit,ksp->guess_zero,&ksp->its,&reason);
 71:     ksp->reason = (KSPConvergedReason)reason;
 72:     return(0);
 73:   } else if (exists && !ksp->numbermonitors && !ksp->transpose_solve & !nullsp) {
 74:     PetscInfo(ksp,"KSPSolve_Richardson: Warning, skipping optimized PCApplyRichardson() because scale factor is not 1.0\n");
 75:   }

 77:   scale = richardsonP->scale;

 79:   if (!ksp->guess_zero) {                          /*   r <- b - A x     */
 80:     KSP_MatMult(ksp,Amat,x,r);
 81:     VecAYPX(r,-1.0,b);
 82:   } else {
 83:     VecCopy(b,r);
 84:   }

 86:   ksp->its = 0;
 87:   if (richardsonP->selfscale) {
 88:     KSP_PCApply(ksp,r,z);         /*   z <- B r          */
 89:     for (i=0; i<maxit; i++) {

 91:       if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
 92:         VecNorm(r,NORM_2,&rnorm); /*   rnorm <- r'*r     */
 93:         KSPMonitor(ksp,i,rnorm);
 94:         ksp->rnorm = rnorm;
 95:         KSPLogResidualHistory(ksp,rnorm);
 96:         (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
 97:         if (ksp->reason) break;
 98:       } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
 99:         VecNorm(z,NORM_2,&rnorm); /*   rnorm <- z'*z     */
100:         KSPMonitor(ksp,i,rnorm);
101:         ksp->rnorm = rnorm;
102:         KSPLogResidualHistory(ksp,rnorm);
103:         (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
104:         if (ksp->reason) break;
105:       }
106:       KSP_PCApplyBAorAB(ksp,z,y,w); /* y = BAz = BABr */
107:       VecDotNorm2(z,y,&rdot,&abr);   /*   rdot = (Br)^T(BABR); abr = (BABr)^T (BABr) */
108:       scale = rdot/abr;
109:       PetscInfo1(ksp,"Self-scale factor %g\n",(double)PetscRealPart(scale));
110:       VecAXPY(x,scale,z);   /*   x  <- x + scale z */
111:       VecAXPY(r,-scale,w);  /*  r <- r - scale*Az */
112:       VecAXPY(z,-scale,y);  /*  z <- z - scale*y */
113:       ksp->its++;
114:     }
115:   } else {
116:     for (i=0; i<maxit; i++) {

118:       if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
119:         VecNorm(r,NORM_2,&rnorm); /*   rnorm <- r'*r     */
120:         KSPMonitor(ksp,i,rnorm);
121:         ksp->rnorm = rnorm;
122:         KSPLogResidualHistory(ksp,rnorm);
123:         (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
124:         if (ksp->reason) break;
125:       }

127:       KSP_PCApply(ksp,r,z);    /*   z <- B r          */

129:       if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
130:         VecNorm(z,NORM_2,&rnorm); /*   rnorm <- z'*z     */
131:         KSPMonitor(ksp,i,rnorm);
132:         ksp->rnorm = rnorm;
133:         KSPLogResidualHistory(ksp,rnorm);
134:         (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
135:         if (ksp->reason) break;
136:       }

138:       VecAXPY(x,scale,z);    /*   x  <- x + scale z */
139:       ksp->its++;

141:       if (i+1 < maxit || ksp->normtype != KSP_NORM_NONE) {
142:         KSP_MatMult(ksp,Amat,x,r);      /*   r  <- b - Ax      */
143:         VecAYPX(r,-1.0,b);
144:       }
145:     }
146:   }
147:   if (!ksp->reason) {
148:     if (ksp->normtype != KSP_NORM_NONE) {
149:       if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
150:         VecNorm(r,NORM_2,&rnorm);     /*   rnorm <- r'*r     */
151:       } else {
152:         KSP_PCApply(ksp,r,z);   /*   z <- B r          */
153:         VecNorm(z,NORM_2,&rnorm);     /*   rnorm <- z'*z     */
154:       }
155:       ksp->rnorm = rnorm;
156:       KSPLogResidualHistory(ksp,rnorm);
157:       KSPMonitor(ksp,i,rnorm);
158:     }
159:     if (ksp->its >= ksp->max_it) {
160:       if (ksp->normtype != KSP_NORM_NONE) {
161:         (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
162:         if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
163:       } else {
164:         ksp->reason = KSP_CONVERGED_ITS;
165:       }
166:     }
167:   }
168:   return(0);
169: }

173: PetscErrorCode KSPView_Richardson(KSP ksp,PetscViewer viewer)
174: {
175:   KSP_Richardson *richardsonP = (KSP_Richardson*)ksp->data;
177:   PetscBool      iascii;

180:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
181:   if (iascii) {
182:     if (richardsonP->selfscale) {
183:       PetscViewerASCIIPrintf(viewer,"  Richardson: using self-scale best computed damping factor\n");
184:     } else {
185:       PetscViewerASCIIPrintf(viewer,"  Richardson: damping factor=%g\n",(double)richardsonP->scale);
186:     }
187:   }
188:   return(0);
189: }

193: PetscErrorCode KSPSetFromOptions_Richardson(PetscOptions *PetscOptionsObject,KSP ksp)
194: {
195:   KSP_Richardson *rich = (KSP_Richardson*)ksp->data;
197:   PetscReal      tmp;
198:   PetscBool      flg,flg2;

201:   PetscOptionsHead(PetscOptionsObject,"KSP Richardson Options");
202:   PetscOptionsReal("-ksp_richardson_scale","damping factor","KSPRichardsonSetScale",rich->scale,&tmp,&flg);
203:   if (flg) { KSPRichardsonSetScale(ksp,tmp); }
204:   PetscOptionsBool("-ksp_richardson_self_scale","dynamically determine optimal damping factor","KSPRichardsonSetSelfScale",rich->selfscale,&flg2,&flg);
205:   if (flg) { KSPRichardsonSetSelfScale(ksp,flg2); }
206:   PetscOptionsTail();
207:   return(0);
208: }

212: PetscErrorCode KSPDestroy_Richardson(KSP ksp)
213: {

217:   PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetScale_C",NULL);
218:   KSPDestroyDefault(ksp);
219:   return(0);
220: }

224: static PetscErrorCode  KSPRichardsonSetScale_Richardson(KSP ksp,PetscReal scale)
225: {
226:   KSP_Richardson *richardsonP;

229:   richardsonP        = (KSP_Richardson*)ksp->data;
230:   richardsonP->scale = scale;
231:   return(0);
232: }

236: static PetscErrorCode  KSPRichardsonSetSelfScale_Richardson(KSP ksp,PetscBool selfscale)
237: {
238:   KSP_Richardson *richardsonP;

241:   richardsonP            = (KSP_Richardson*)ksp->data;
242:   richardsonP->selfscale = selfscale;
243:   return(0);
244: }

246: /*MC
247:      KSPRICHARDSON - The preconditioned Richardson iterative method

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

252:    Level: beginner

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

256:           Here B is the application of the preconditioner

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

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

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

269:     Supports only left preconditioning

271:   References:
272:   "The Approximate Arithmetical Solution by Finite Differences of Physical Problems Involving
273:    Differential Equations, with an Application to the Stresses in a Masonry Dam",
274:   L. F. Richardson, Philosophical Transactions of the Royal Society of London. Series A,
275:   Containing Papers of a Mathematical or Physical Character, Vol. 210, 1911 (1911), pp. 307-357.

277: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
278:            KSPRichardsonSetScale()

280: M*/

284: PETSC_EXTERN PetscErrorCode KSPCreate_Richardson(KSP ksp)
285: {
287:   KSP_Richardson *richardsonP;

290:   PetscNewLog(ksp,&richardsonP);
291:   ksp->data = (void*)richardsonP;

293:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
294:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);

296:   ksp->ops->setup          = KSPSetUp_Richardson;
297:   ksp->ops->solve          = KSPSolve_Richardson;
298:   ksp->ops->destroy        = KSPDestroy_Richardson;
299:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
300:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
301:   ksp->ops->view           = KSPView_Richardson;
302:   ksp->ops->setfromoptions = KSPSetFromOptions_Richardson;

304:   PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetScale_C",KSPRichardsonSetScale_Richardson);
305:   PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetSelfScale_C",KSPRichardsonSetSelfScale_Richardson);

307:   richardsonP->scale = 1.0;
308:   return(0);
309: }