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: }