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:   PetscFunctionBegin;
 12:   if (richardsonP->selfscale) {
 13:     PetscCall(KSPSetWorkVecs(ksp, 4));
 14:   } else {
 15:     PetscCall(KSPSetWorkVecs(ksp, 2));
 16:   }
 17:   PetscFunctionReturn(PETSC_SUCCESS);
 18: }

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

 32:   PetscFunctionBegin;
 33:   PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
 34:   PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);

 36:   ksp->its = 0;

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

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

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

 77:   ksp->its = 0;
 78:   if (richardsonP->selfscale) {
 79:     PetscCall(KSP_PCApply(ksp, r, z)); /*   z <- B r          */
 80:     for (i = 0; i < maxit; i++) {
 81:       if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
 82:         PetscCall(VecNorm(r, NORM_2, &rnorm)); /*   rnorm <- r'*r     */
 83:       } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
 84:         PetscCall(VecNorm(z, NORM_2, &rnorm)); /*   rnorm <- z'*z     */
 85:       } else rnorm = 0.0;

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

117:       PetscCall(VecAXPY(x, richardsonP->scale, z)); /*   x  <- x + scale z */
118:       ksp->its++;

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

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

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

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

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

173:   PetscFunctionBegin;
174:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP Richardson Options");
175:   PetscCall(PetscOptionsReal("-ksp_richardson_scale", "damping factor", "KSPRichardsonSetScale", rich->scale, &tmp, &flg));
176:   if (flg) PetscCall(KSPRichardsonSetScale(ksp, tmp));
177:   PetscCall(PetscOptionsBool("-ksp_richardson_self_scale", "dynamically determine optimal damping factor", "KSPRichardsonSetSelfScale", rich->selfscale, &flg2, &flg));
178:   if (flg) PetscCall(KSPRichardsonSetSelfScale(ksp, flg2));
179:   PetscOptionsHeadEnd();
180:   PetscFunctionReturn(PETSC_SUCCESS);
181: }

183: PetscErrorCode KSPDestroy_Richardson(KSP ksp)
184: {
185:   PetscFunctionBegin;
186:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetScale_C", NULL));
187:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetSelfScale_C", NULL));
188:   PetscCall(KSPDestroyDefault(ksp));
189:   PetscFunctionReturn(PETSC_SUCCESS);
190: }

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

196:   PetscFunctionBegin;
197:   richardsonP        = (KSP_Richardson *)ksp->data;
198:   richardsonP->scale = scale;
199:   PetscFunctionReturn(PETSC_SUCCESS);
200: }

202: static PetscErrorCode KSPRichardsonSetSelfScale_Richardson(KSP ksp, PetscBool selfscale)
203: {
204:   KSP_Richardson *richardsonP;

206:   PetscFunctionBegin;
207:   richardsonP            = (KSP_Richardson *)ksp->data;
208:   richardsonP->selfscale = selfscale;
209:   PetscFunctionReturn(PETSC_SUCCESS);
210: }

212: static PetscErrorCode KSPBuildResidual_Richardson(KSP ksp, Vec t, Vec v, Vec *V)
213: {
214:   PetscFunctionBegin;
215:   if (ksp->normtype == KSP_NORM_NONE) {
216:     PetscCall(KSPBuildResidualDefault(ksp, t, v, V));
217:   } else {
218:     PetscCall(VecCopy(ksp->work[0], v));
219:     *V = v;
220:   }
221:   PetscFunctionReturn(PETSC_SUCCESS);
222: }

224: /*MC
225:      KSPRICHARDSON - The preconditioned Richardson iterative method

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

230:    Level: beginner

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

235:           Here B is the application of the preconditioner

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

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

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

248:     Supports only left preconditioning

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

252:     `-ksp_type richardson -pc_type jacobi` gives one classical Jacobi preconditioning

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

260: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
261:           `KSPRichardsonSetScale()`, `KSPPREONLY`
262: M*/

264: PETSC_EXTERN PetscErrorCode KSPCreate_Richardson(KSP ksp)
265: {
266:   KSP_Richardson *richardsonP;

268:   PetscFunctionBegin;
269:   PetscCall(PetscNew(&richardsonP));
270:   ksp->data = (void *)richardsonP;

272:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
273:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
274:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));

276:   ksp->ops->setup          = KSPSetUp_Richardson;
277:   ksp->ops->solve          = KSPSolve_Richardson;
278:   ksp->ops->destroy        = KSPDestroy_Richardson;
279:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
280:   ksp->ops->buildresidual  = KSPBuildResidual_Richardson;
281:   ksp->ops->view           = KSPView_Richardson;
282:   ksp->ops->setfromoptions = KSPSetFromOptions_Richardson;

284:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetScale_C", KSPRichardsonSetScale_Richardson));
285:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetSelfScale_C", KSPRichardsonSetSelfScale_Richardson));

287:   richardsonP->scale = 1.0;
288:   PetscFunctionReturn(PETSC_SUCCESS);
289: }