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