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++) {
 79:       if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
 80:         VecNorm(r, NORM_2, &rnorm); /*   rnorm <- r'*r     */
 81:       } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
 82:         VecNorm(z, NORM_2, &rnorm); /*   rnorm <- z'*z     */
 83:       } else rnorm = 0.0;

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

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

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

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

148: PetscErrorCode KSPView_Richardson(KSP ksp, PetscViewer viewer)
149: {
150:   KSP_Richardson *richardsonP = (KSP_Richardson *)ksp->data;
151:   PetscBool       iascii;

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

164: PetscErrorCode KSPSetFromOptions_Richardson(KSP ksp, PetscOptionItems *PetscOptionsObject)
165: {
166:   KSP_Richardson *rich = (KSP_Richardson *)ksp->data;
167:   PetscReal       tmp;
168:   PetscBool       flg, flg2;

170:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP Richardson Options");
171:   PetscOptionsReal("-ksp_richardson_scale", "damping factor", "KSPRichardsonSetScale", rich->scale, &tmp, &flg);
172:   if (flg) KSPRichardsonSetScale(ksp, tmp);
173:   PetscOptionsBool("-ksp_richardson_self_scale", "dynamically determine optimal damping factor", "KSPRichardsonSetSelfScale", rich->selfscale, &flg2, &flg);
174:   if (flg) KSPRichardsonSetSelfScale(ksp, flg2);
175:   PetscOptionsHeadEnd();
176:   return 0;
177: }

179: PetscErrorCode KSPDestroy_Richardson(KSP ksp)
180: {
181:   PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetScale_C", NULL);
182:   PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetSelfScale_C", NULL);
183:   KSPDestroyDefault(ksp);
184:   return 0;
185: }

187: static PetscErrorCode KSPRichardsonSetScale_Richardson(KSP ksp, PetscReal scale)
188: {
189:   KSP_Richardson *richardsonP;

191:   richardsonP        = (KSP_Richardson *)ksp->data;
192:   richardsonP->scale = scale;
193:   return 0;
194: }

196: static PetscErrorCode KSPRichardsonSetSelfScale_Richardson(KSP ksp, PetscBool selfscale)
197: {
198:   KSP_Richardson *richardsonP;

200:   richardsonP            = (KSP_Richardson *)ksp->data;
201:   richardsonP->selfscale = selfscale;
202:   return 0;
203: }

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

216: /*MC
217:      KSPRICHARDSON - The preconditioned Richardson iterative method

219:    Options Database Key:
220: .   -ksp_richardson_scale - damping factor on the correction (defaults to 1.0)

222:    Level: beginner

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

227:           Here B is the application of the preconditioner

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

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

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

240:     Supports only left preconditioning

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

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

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

252: .seealso: [](chapter_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
253:           `KSPRichardsonSetScale()`, `KSPPREONLY`
254: M*/

256: PETSC_EXTERN PetscErrorCode KSPCreate_Richardson(KSP ksp)
257: {
258:   KSP_Richardson *richardsonP;

260:   PetscNew(&richardsonP);
261:   ksp->data = (void *)richardsonP;

263:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3);
264:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2);
265:   KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1);

267:   ksp->ops->setup          = KSPSetUp_Richardson;
268:   ksp->ops->solve          = KSPSolve_Richardson;
269:   ksp->ops->destroy        = KSPDestroy_Richardson;
270:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
271:   ksp->ops->buildresidual  = KSPBuildResidual_Richardson;
272:   ksp->ops->view           = KSPView_Richardson;
273:   ksp->ops->setfromoptions = KSPSetFromOptions_Richardson;

275:   PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetScale_C", KSPRichardsonSetScale_Richardson);
276:   PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetSelfScale_C", KSPRichardsonSetSelfScale_Richardson);

278:   richardsonP->scale = 1.0;
279:   return 0;
280: }