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