Actual source code: rich.c
1: /*
2: This implements Richardson Iteration.
3: */
4: #include <../src/ksp/ksp/impls/rich/richardsonimpl.h>
6: static PetscErrorCode KSPSetUp_Richardson(KSP ksp)
7: {
8: KSP_Richardson *richardsonP = (KSP_Richardson *)ksp->data;
10: PetscFunctionBegin;
11: if (richardsonP->selfscale) {
12: PetscCall(KSPSetWorkVecs(ksp, 4));
13: } else {
14: PetscCall(KSPSetWorkVecs(ksp, 2));
15: }
16: PetscFunctionReturn(PETSC_SUCCESS);
17: }
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;
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);
35: ksp->its = 0;
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;
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: } else {
66: PetscCall(PetscInfo(ksp, "KSPSolve_Richardson: Warning, skipping optimized PCApplyRichardson()\n"));
67: }
69: if (!ksp->guess_zero) { /* r <- b - A x */
70: PetscCall(KSP_MatMult(ksp, Amat, x, r));
71: PetscCall(VecAYPX(r, -1.0, b));
72: } else {
73: PetscCall(VecCopy(b, r));
74: }
76: ksp->its = 0;
77: if (richardsonP->selfscale) {
78: PetscCall(KSP_PCApply(ksp, r, z)); /* z <- B r */
79: for (i = 0; i < maxit; i++) {
80: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
81: PetscCall(VecNorm(r, NORM_2, &rnorm)); /* rnorm <- r'*r */
82: } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
83: PetscCall(VecNorm(z, NORM_2, &rnorm)); /* rnorm <- z'*z */
84: } else rnorm = 0.0;
86: KSPCheckNorm(ksp, rnorm);
87: ksp->rnorm = rnorm;
88: PetscCall(KSPMonitor(ksp, i, rnorm));
89: PetscCall(KSPLogResidualHistory(ksp, rnorm));
90: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
91: if (ksp->reason) break;
92: PetscCall(KSP_PCApplyBAorAB(ksp, z, y, w)); /* y = BAz = BABr */
93: PetscCall(VecDotNorm2(z, y, &rdot, &abr)); /* rdot = (Br)^T(BABR); abr = (BABr)^T (BABr) */
94: scale = rdot / abr;
95: PetscCall(PetscInfo(ksp, "Self-scale factor %g\n", (double)PetscRealPart(scale)));
96: PetscCall(VecAXPY(x, scale, z)); /* x <- x + scale z */
97: PetscCall(VecAXPY(r, -scale, w)); /* r <- r - scale*Az */
98: PetscCall(VecAXPY(z, -scale, y)); /* z <- z - scale*y */
99: ksp->its++;
100: }
101: } else {
102: for (i = 0; i < maxit; i++) {
103: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
104: PetscCall(VecNorm(r, NORM_2, &rnorm)); /* rnorm <- r'*r */
105: } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
106: PetscCall(KSP_PCApply(ksp, r, z)); /* z <- B r */
107: PetscCall(VecNorm(z, NORM_2, &rnorm)); /* rnorm <- z'*z */
108: } else rnorm = 0.0;
109: ksp->rnorm = rnorm;
110: PetscCall(KSPMonitor(ksp, i, rnorm));
111: PetscCall(KSPLogResidualHistory(ksp, rnorm));
112: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
113: if (ksp->reason) break;
114: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, z)); /* z <- B r */ }
116: PetscCall(VecAXPY(x, richardsonP->scale, z)); /* x <- x + scale z */
117: ksp->its++;
119: if (i + 1 < maxit || ksp->normtype != KSP_NORM_NONE) {
120: PetscCall(KSP_MatMult(ksp, Amat, x, r)); /* r <- b - Ax */
121: PetscCall(VecAYPX(r, -1.0, b));
122: }
123: }
124: }
125: if (!ksp->reason) {
126: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
127: PetscCall(VecNorm(r, NORM_2, &rnorm)); /* rnorm <- r'*r */
128: } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
129: PetscCall(KSP_PCApply(ksp, r, z)); /* z <- B r */
130: PetscCall(VecNorm(z, NORM_2, &rnorm)); /* rnorm <- z'*z */
131: } else rnorm = 0.0;
133: KSPCheckNorm(ksp, rnorm);
134: ksp->rnorm = rnorm;
135: PetscCall(KSPLogResidualHistory(ksp, rnorm));
136: PetscCall(KSPMonitor(ksp, i, rnorm));
137: if (ksp->its >= ksp->max_it) {
138: if (ksp->normtype != KSP_NORM_NONE) {
139: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
140: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
141: } else {
142: ksp->reason = KSP_CONVERGED_ITS;
143: }
144: }
145: }
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: static PetscErrorCode KSPView_Richardson(KSP ksp, PetscViewer viewer)
150: {
151: KSP_Richardson *richardsonP = (KSP_Richardson *)ksp->data;
152: PetscBool iascii;
154: PetscFunctionBegin;
155: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
156: if (iascii) {
157: if (richardsonP->selfscale) {
158: PetscCall(PetscViewerASCIIPrintf(viewer, " using self-scale best computed damping factor\n"));
159: } else {
160: PetscCall(PetscViewerASCIIPrintf(viewer, " damping factor=%g\n", (double)richardsonP->scale));
161: }
162: }
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: static PetscErrorCode KSPSetFromOptions_Richardson(KSP ksp, PetscOptionItems *PetscOptionsObject)
167: {
168: KSP_Richardson *rich = (KSP_Richardson *)ksp->data;
169: PetscReal tmp;
170: PetscBool flg, flg2;
172: PetscFunctionBegin;
173: PetscOptionsHeadBegin(PetscOptionsObject, "KSP Richardson Options");
174: PetscCall(PetscOptionsReal("-ksp_richardson_scale", "damping factor", "KSPRichardsonSetScale", rich->scale, &tmp, &flg));
175: if (flg) PetscCall(KSPRichardsonSetScale(ksp, tmp));
176: PetscCall(PetscOptionsBool("-ksp_richardson_self_scale", "dynamically determine optimal damping factor", "KSPRichardsonSetSelfScale", rich->selfscale, &flg2, &flg));
177: if (flg) PetscCall(KSPRichardsonSetSelfScale(ksp, flg2));
178: PetscOptionsHeadEnd();
179: PetscFunctionReturn(PETSC_SUCCESS);
180: }
182: static PetscErrorCode KSPDestroy_Richardson(KSP ksp)
183: {
184: PetscFunctionBegin;
185: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetScale_C", NULL));
186: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetSelfScale_C", NULL));
187: PetscCall(KSPDestroyDefault(ksp));
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: static PetscErrorCode KSPRichardsonSetScale_Richardson(KSP ksp, PetscReal scale)
192: {
193: KSP_Richardson *richardsonP;
195: PetscFunctionBegin;
196: richardsonP = (KSP_Richardson *)ksp->data;
197: richardsonP->scale = scale;
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }
201: static PetscErrorCode KSPRichardsonSetSelfScale_Richardson(KSP ksp, PetscBool selfscale)
202: {
203: KSP_Richardson *richardsonP;
205: PetscFunctionBegin;
206: richardsonP = (KSP_Richardson *)ksp->data;
207: richardsonP->selfscale = selfscale;
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: static PetscErrorCode KSPBuildResidual_Richardson(KSP ksp, Vec t, Vec v, Vec *V)
212: {
213: PetscFunctionBegin;
214: if (ksp->normtype == KSP_NORM_NONE) {
215: PetscCall(KSPBuildResidualDefault(ksp, t, v, V));
216: } else {
217: PetscCall(VecCopy(ksp->work[0], v));
218: *V = v;
219: }
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: /*MC
224: KSPRICHARDSON - The preconditioned Richardson iterative method
226: Options Database Key:
227: . -ksp_richardson_scale - damping factor on the correction (defaults to 1.0)
229: Level: beginner
231: Notes:
232: x^{n+1} = x^{n} + scale*B(b - A x^{n})
234: Here B is the application of the preconditioner
236: This method often (usually) will not converge unless scale is very small.
238: For some preconditioners, currently `PCSOR`, the convergence test is skipped to improve speed,
239: thus it always iterates the maximum number of iterations you've selected. When -ksp_monitor
240: (or any other monitor) is turned on, the norm is computed at each iteration and so the convergence test is run unless
241: you specifically call `KSPSetNormType`(ksp,`KSP_NORM_NONE`);
243: For some preconditioners, currently `PCMG` and `PCHYPRE` with BoomerAMG if -ksp_monitor (and also
244: any other monitor) is not turned on then the convergence test is done by the preconditioner itself and
245: so the solver may run more or fewer iterations then if -ksp_monitor is selected.
247: Supports only left preconditioning
249: If using direct solvers such as `PCLU` and `PCCHOLESKY` one generally uses `KSPPREONLY` which uses exactly one iteration
251: `-ksp_type richardson -pc_type jacobi` gives one classical Jacobi preconditioning
253: Reference:
254: . * - L. F. Richardson, "The Approximate Arithmetical Solution by Finite Differences of Physical Problems Involving
255: Differential Equations, with an Application to the Stresses in a Masonry Dam",
256: Philosophical Transactions of the Royal Society of London. Series A,
257: Containing Papers of a Mathematical or Physical Character, Vol. 210, 1911 (1911).
259: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
260: `KSPRichardsonSetScale()`, `KSPPREONLY`
261: M*/
263: PETSC_EXTERN PetscErrorCode KSPCreate_Richardson(KSP ksp)
264: {
265: KSP_Richardson *richardsonP;
267: PetscFunctionBegin;
268: PetscCall(PetscNew(&richardsonP));
269: ksp->data = (void *)richardsonP;
271: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
272: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
273: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
275: ksp->ops->setup = KSPSetUp_Richardson;
276: ksp->ops->solve = KSPSolve_Richardson;
277: ksp->ops->destroy = KSPDestroy_Richardson;
278: ksp->ops->buildsolution = KSPBuildSolutionDefault;
279: ksp->ops->buildresidual = KSPBuildResidual_Richardson;
280: ksp->ops->view = KSPView_Richardson;
281: ksp->ops->setfromoptions = KSPSetFromOptions_Richardson;
283: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetScale_C", KSPRichardsonSetScale_Richardson));
284: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetSelfScale_C", KSPRichardsonSetSelfScale_Richardson));
286: richardsonP->scale = 1.0;
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }