1 /*
2 This implements Richardson Iteration.
3 */
4 #include <../src/ksp/ksp/impls/rich/richardsonimpl.h> /*I "petscksp.h" I*/
5
KSPSetUp_Richardson(KSP ksp)6 static PetscErrorCode KSPSetUp_Richardson(KSP ksp)
7 {
8 KSP_Richardson *richardsonP = (KSP_Richardson *)ksp->data;
9
10 PetscFunctionBegin;
11 if (richardsonP->selfscale) {
12 PetscCall(KSPSetWorkVecs(ksp, 4));
13 } else {
14 PetscCall(KSPSetWorkVecs(ksp, 2));
15 }
16 PetscFunctionReturn(PETSC_SUCCESS);
17 }
18
KSPSolve_Richardson(KSP ksp)19 static 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;
30
31 PetscFunctionBegin;
32 PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
33 PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
34
35 ksp->its = 0;
36
37 PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
38 x = ksp->vec_sol;
39 b = ksp->vec_rhs;
40 PetscCall(VecGetSize(x, &xs));
41 PetscCall(VecGetSize(ksp->work[0], &ws));
42 if (xs != ws) {
43 if (richardsonP->selfscale) {
44 PetscCall(KSPSetWorkVecs(ksp, 4));
45 } else {
46 PetscCall(KSPSetWorkVecs(ksp, 2));
47 }
48 }
49 r = ksp->work[0];
50 z = ksp->work[1];
51 if (richardsonP->selfscale) {
52 w = ksp->work[2];
53 y = ksp->work[3];
54 }
55 maxit = ksp->max_it;
56
57 /* if user has provided fast Richardson code use that */
58 PetscCall(PCApplyRichardsonExists(ksp->pc, &exists));
59 PetscCall(MatGetNullSpace(Pmat, &nullsp));
60 if (exists && maxit > 0 && richardsonP->scale == 1.0 && (ksp->converged == KSPConvergedDefault || ksp->converged == KSPConvergedSkip) && !ksp->numbermonitors && !ksp->transpose_solve && !nullsp) {
61 PCRichardsonConvergedReason reason;
62 PetscCall(PCApplyRichardson(ksp->pc, b, x, r, ksp->rtol, ksp->abstol, ksp->divtol, maxit, ksp->guess_zero, &ksp->its, &reason));
63 ksp->reason = (KSPConvergedReason)reason;
64 PetscFunctionReturn(PETSC_SUCCESS);
65 }
66
67 if (!ksp->guess_zero) { /* r <- b - A x */
68 PetscCall(KSP_MatMult(ksp, Amat, x, r));
69 PetscCall(VecAYPX(r, -1.0, b));
70 } else {
71 PetscCall(VecCopy(b, r));
72 }
73
74 ksp->its = 0;
75 if (richardsonP->selfscale) {
76 PetscCall(KSP_PCApply(ksp, r, z)); /* z <- B r */
77 for (i = 0; i < maxit; i++) {
78 if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
79 PetscCall(VecNorm(r, NORM_2, &rnorm)); /* rnorm <- r'*r */
80 } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
81 PetscCall(VecNorm(z, NORM_2, &rnorm)); /* rnorm <- z'*z */
82 } else rnorm = 0.0;
83
84 KSPCheckNorm(ksp, rnorm);
85 ksp->rnorm = rnorm;
86 PetscCall(KSPMonitor(ksp, i, rnorm));
87 PetscCall(KSPLogResidualHistory(ksp, rnorm));
88 PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
89 if (ksp->reason) break;
90 PetscCall(KSP_PCApplyBAorAB(ksp, z, y, w)); /* y = BAz = BABr */
91 PetscCall(VecDotNorm2(z, y, &rdot, &abr)); /* rdot = (Br)^T(BABR); abr = (BABr)^T (BABr) */
92 scale = rdot / abr;
93 PetscCall(PetscInfo(ksp, "Self-scale factor %g\n", (double)PetscRealPart(scale)));
94 PetscCall(VecAXPY(x, scale, z)); /* x <- x + scale z */
95 PetscCall(VecAXPY(r, -scale, w)); /* r <- r - scale*Az */
96 PetscCall(VecAXPY(z, -scale, y)); /* z <- z - scale*y */
97 ksp->its++;
98 }
99 } else {
100 for (i = 0; i < maxit; i++) {
101 if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
102 PetscCall(VecNorm(r, NORM_2, &rnorm)); /* rnorm <- r'*r */
103 } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
104 PetscCall(KSP_PCApply(ksp, r, z)); /* z <- B r */
105 PetscCall(VecNorm(z, NORM_2, &rnorm)); /* rnorm <- z'*z */
106 } else rnorm = 0.0;
107 ksp->rnorm = rnorm;
108 PetscCall(KSPMonitor(ksp, i, rnorm));
109 PetscCall(KSPLogResidualHistory(ksp, rnorm));
110 PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
111 if (ksp->reason) break;
112 if (ksp->normtype != KSP_NORM_PRECONDITIONED) PetscCall(KSP_PCApply(ksp, r, z)); /* z <- B r */
113
114 PetscCall(VecAXPY(x, richardsonP->scale, z)); /* x <- x + scale z */
115 ksp->its++;
116
117 if (i + 1 < maxit || ksp->normtype != KSP_NORM_NONE) {
118 PetscCall(KSP_MatMult(ksp, Amat, x, r)); /* r <- b - Ax */
119 PetscCall(VecAYPX(r, -1.0, b));
120 }
121 }
122 }
123 if (!ksp->reason) {
124 if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
125 PetscCall(VecNorm(r, NORM_2, &rnorm)); /* rnorm <- r'*r */
126 } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
127 PetscCall(KSP_PCApply(ksp, r, z)); /* z <- B r */
128 PetscCall(VecNorm(z, NORM_2, &rnorm)); /* rnorm <- z'*z */
129 } else rnorm = 0.0;
130
131 KSPCheckNorm(ksp, rnorm);
132 ksp->rnorm = rnorm;
133 PetscCall(KSPLogResidualHistory(ksp, rnorm));
134 PetscCall(KSPMonitor(ksp, i, rnorm));
135 if (ksp->its >= ksp->max_it) {
136 if (ksp->normtype != KSP_NORM_NONE) {
137 PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
138 if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
139 } else {
140 ksp->reason = KSP_CONVERGED_ITS;
141 }
142 }
143 }
144 PetscFunctionReturn(PETSC_SUCCESS);
145 }
146
KSPView_Richardson(KSP ksp,PetscViewer viewer)147 static PetscErrorCode KSPView_Richardson(KSP ksp, PetscViewer viewer)
148 {
149 KSP_Richardson *richardsonP = (KSP_Richardson *)ksp->data;
150 PetscBool isascii;
151
152 PetscFunctionBegin;
153 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
154 if (isascii) {
155 if (richardsonP->selfscale) {
156 PetscCall(PetscViewerASCIIPrintf(viewer, " using self-scale best computed damping factor\n"));
157 } else {
158 PetscCall(PetscViewerASCIIPrintf(viewer, " damping factor=%g\n", (double)richardsonP->scale));
159 }
160 }
161 PetscFunctionReturn(PETSC_SUCCESS);
162 }
163
KSPSetFromOptions_Richardson(KSP ksp,PetscOptionItems PetscOptionsObject)164 static PetscErrorCode KSPSetFromOptions_Richardson(KSP ksp, PetscOptionItems PetscOptionsObject)
165 {
166 KSP_Richardson *rich = (KSP_Richardson *)ksp->data;
167 PetscReal tmp;
168 PetscBool flg, flg2;
169
170 PetscFunctionBegin;
171 PetscOptionsHeadBegin(PetscOptionsObject, "KSP Richardson Options");
172 PetscCall(PetscOptionsReal("-ksp_richardson_scale", "damping factor", "KSPRichardsonSetScale", rich->scale, &tmp, &flg));
173 if (flg) PetscCall(KSPRichardsonSetScale(ksp, tmp));
174 PetscCall(PetscOptionsBool("-ksp_richardson_self_scale", "dynamically determine optimal damping factor", "KSPRichardsonSetSelfScale", rich->selfscale, &flg2, &flg));
175 if (flg) PetscCall(KSPRichardsonSetSelfScale(ksp, flg2));
176 PetscOptionsHeadEnd();
177 PetscFunctionReturn(PETSC_SUCCESS);
178 }
179
KSPDestroy_Richardson(KSP ksp)180 static PetscErrorCode KSPDestroy_Richardson(KSP ksp)
181 {
182 PetscFunctionBegin;
183 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetScale_C", NULL));
184 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetSelfScale_C", NULL));
185 PetscCall(KSPDestroyDefault(ksp));
186 PetscFunctionReturn(PETSC_SUCCESS);
187 }
188
KSPRichardsonSetScale_Richardson(KSP ksp,PetscReal scale)189 static PetscErrorCode KSPRichardsonSetScale_Richardson(KSP ksp, PetscReal scale)
190 {
191 KSP_Richardson *richardsonP;
192
193 PetscFunctionBegin;
194 richardsonP = (KSP_Richardson *)ksp->data;
195 richardsonP->scale = scale;
196 PetscFunctionReturn(PETSC_SUCCESS);
197 }
198
KSPRichardsonSetSelfScale_Richardson(KSP ksp,PetscBool selfscale)199 static PetscErrorCode KSPRichardsonSetSelfScale_Richardson(KSP ksp, PetscBool selfscale)
200 {
201 KSP_Richardson *richardsonP;
202
203 PetscFunctionBegin;
204 richardsonP = (KSP_Richardson *)ksp->data;
205 richardsonP->selfscale = selfscale;
206 PetscFunctionReturn(PETSC_SUCCESS);
207 }
208
KSPBuildResidual_Richardson(KSP ksp,Vec t,Vec v,Vec * V)209 static PetscErrorCode KSPBuildResidual_Richardson(KSP ksp, Vec t, Vec v, Vec *V)
210 {
211 PetscFunctionBegin;
212 if (ksp->normtype == KSP_NORM_NONE) {
213 PetscCall(KSPBuildResidualDefault(ksp, t, v, V));
214 } else {
215 PetscCall(VecCopy(ksp->work[0], v));
216 *V = v;
217 }
218 PetscFunctionReturn(PETSC_SUCCESS);
219 }
220
221 /*MC
222 KSPRICHARDSON - The preconditioned Richardson iterative method {cite}`richarson1911`
223
224 Options Database Key:
225 . -ksp_richardson_scale - damping factor on the correction (defaults to 1.0)
226
227 Level: beginner
228
229 Notes:
230 $ x^{n+1} = x^{n} + scale*B(b - A x^{n})$
231
232 Here B is the application of the preconditioner
233
234 This method often (usually) will not converge unless scale is very small.
235
236 For some preconditioners, currently `PCSOR`, the convergence test is skipped to improve speed,
237 thus it always iterates the maximum number of iterations you've selected. When -ksp_monitor
238 (or any other monitor) is turned on, the norm is computed at each iteration and so the convergence test is run unless
239 you specifically call `KSPSetNormType`(ksp,`KSP_NORM_NONE`);
240
241 For some preconditioners, currently `PCMG` and `PCHYPRE` with BoomerAMG if -ksp_monitor (and also
242 any other monitor) is not turned on then the convergence test is done by the preconditioner itself and
243 so the solver may run more or fewer iterations then if -ksp_monitor is selected.
244
245 Supports only left preconditioning
246
247 If using direct solvers such as `PCLU` and `PCCHOLESKY` one generally uses `KSPPREONLY` instead of this which uses exactly one iteration
248
249 `-ksp_type richardson -pc_type jacobi` gives one classical Jacobi preconditioning
250
251 .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
252 `KSPRichardsonSetScale()`, `KSPPREONLY`, `KSPRichardsonSetSelfScale()`
253 M*/
254
KSPCreate_Richardson(KSP ksp)255 PETSC_EXTERN PetscErrorCode KSPCreate_Richardson(KSP ksp)
256 {
257 KSP_Richardson *richardsonP;
258
259 PetscFunctionBegin;
260 PetscCall(PetscNew(&richardsonP));
261 ksp->data = (void *)richardsonP;
262
263 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
264 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
265 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
266
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;
274
275 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetScale_C", KSPRichardsonSetScale_Richardson));
276 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetSelfScale_C", KSPRichardsonSetSelfScale_Richardson));
277
278 richardsonP->scale = 1.0;
279 PetscFunctionReturn(PETSC_SUCCESS);
280 }
281