Actual source code: rich.c
petsc-3.13.6 2020-09-29
2: /*
3: This implements Richardson Iteration.
4: */
5: #include <../src/ksp/ksp/impls/rich/richardsonimpl.h>
7: PetscErrorCode KSPSetUp_Richardson(KSP ksp)
8: {
10: KSP_Richardson *richardsonP = (KSP_Richardson*)ksp->data;
13: if (richardsonP->selfscale) {
14: KSPSetWorkVecs(ksp,4);
15: } else {
16: KSPSetWorkVecs(ksp,2);
17: }
18: return(0);
19: }
21: PetscErrorCode KSPSolve_Richardson(KSP ksp)
22: {
24: PetscInt i,maxit;
25: PetscReal rnorm = 0.0,abr;
26: PetscScalar scale,rdot;
27: Vec x,b,r,z,w = NULL,y = NULL;
28: PetscInt xs, ws;
29: Mat Amat,Pmat;
30: KSP_Richardson *richardsonP = (KSP_Richardson*)ksp->data;
31: PetscBool exists,diagonalscale;
32: MatNullSpace nullsp;
35: PCGetDiagonalScale(ksp->pc,&diagonalscale);
36: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
38: ksp->its = 0;
40: PCGetOperators(ksp->pc,&Amat,&Pmat);
41: x = ksp->vec_sol;
42: b = ksp->vec_rhs;
43: VecGetSize(x,&xs);
44: VecGetSize(ksp->work[0],&ws);
45: if (xs != ws) {
46: if (richardsonP->selfscale) {
47: KSPSetWorkVecs(ksp,4);
48: } else {
49: KSPSetWorkVecs(ksp,2);
50: }
51: }
52: r = ksp->work[0];
53: z = ksp->work[1];
54: if (richardsonP->selfscale) {
55: w = ksp->work[2];
56: y = ksp->work[3];
57: }
58: maxit = ksp->max_it;
60: /* if user has provided fast Richardson code use that */
61: PCApplyRichardsonExists(ksp->pc,&exists);
62: MatGetNullSpace(Pmat,&nullsp);
63: if (exists && maxit > 0 && richardsonP->scale == 1.0 && !ksp->numbermonitors && !ksp->transpose_solve && !nullsp) {
64: PCRichardsonConvergedReason reason;
65: PCApplyRichardson(ksp->pc,b,x,r,ksp->rtol,ksp->abstol,ksp->divtol,maxit,ksp->guess_zero,&ksp->its,&reason);
66: ksp->reason = (KSPConvergedReason)reason;
67: return(0);
68: } else if (exists && maxit > 0 && !ksp->numbermonitors && !ksp->transpose_solve && !nullsp) {
69: PetscInfo(ksp,"KSPSolve_Richardson: Warning, skipping optimized PCApplyRichardson() because scale factor is not 1.0\n");
70: }
72: if (!ksp->guess_zero) { /* r <- b - A x */
73: KSP_MatMult(ksp,Amat,x,r);
74: VecAYPX(r,-1.0,b);
75: } else {
76: VecCopy(b,r);
77: }
79: ksp->its = 0;
80: if (richardsonP->selfscale) {
81: KSP_PCApply(ksp,r,z); /* z <- B r */
82: for (i=0; i<maxit; i++) {
84: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
85: VecNorm(r,NORM_2,&rnorm); /* rnorm <- r'*r */
86: KSPCheckNorm(ksp,rnorm);
87: KSPMonitor(ksp,i,rnorm);
88: ksp->rnorm = rnorm;
89: KSPLogResidualHistory(ksp,rnorm);
90: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
91: if (ksp->reason) break;
92: } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
93: VecNorm(z,NORM_2,&rnorm); /* rnorm <- z'*z */
94: KSPCheckNorm(ksp,rnorm);
95: KSPMonitor(ksp,i,rnorm);
96: ksp->rnorm = rnorm;
97: KSPLogResidualHistory(ksp,rnorm);
98: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
99: if (ksp->reason) break;
100: }
101: KSP_PCApplyBAorAB(ksp,z,y,w); /* y = BAz = BABr */
102: VecDotNorm2(z,y,&rdot,&abr); /* rdot = (Br)^T(BABR); abr = (BABr)^T (BABr) */
103: scale = rdot/abr;
104: PetscInfo1(ksp,"Self-scale factor %g\n",(double)PetscRealPart(scale));
105: VecAXPY(x,scale,z); /* x <- x + scale z */
106: VecAXPY(r,-scale,w); /* r <- r - scale*Az */
107: VecAXPY(z,-scale,y); /* z <- z - scale*y */
108: ksp->its++;
109: }
110: } else {
111: for (i=0; i<maxit; i++) {
113: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
114: VecNorm(r,NORM_2,&rnorm); /* rnorm <- r'*r */
115: KSPCheckNorm(ksp,rnorm);
116: KSPMonitor(ksp,i,rnorm);
117: ksp->rnorm = rnorm;
118: KSPLogResidualHistory(ksp,rnorm);
119: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
120: if (ksp->reason) break;
121: }
123: KSP_PCApply(ksp,r,z); /* z <- B r */
125: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
126: VecNorm(z,NORM_2,&rnorm); /* rnorm <- z'*z */
127: KSPCheckNorm(ksp,rnorm);
128: KSPMonitor(ksp,i,rnorm);
129: ksp->rnorm = rnorm;
130: KSPLogResidualHistory(ksp,rnorm);
131: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
132: if (ksp->reason) break;
133: }
135: VecAXPY(x,richardsonP->scale,z); /* x <- x + scale z */
136: ksp->its++;
138: if (i+1 < maxit || ksp->normtype != KSP_NORM_NONE) {
139: KSP_MatMult(ksp,Amat,x,r); /* r <- b - Ax */
140: VecAYPX(r,-1.0,b);
141: }
142: }
143: }
144: if (!ksp->reason) {
145: if (ksp->normtype != KSP_NORM_NONE) {
146: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
147: VecNorm(r,NORM_2,&rnorm); /* rnorm <- r'*r */
148: } else {
149: KSP_PCApply(ksp,r,z); /* z <- B r */
150: VecNorm(z,NORM_2,&rnorm); /* rnorm <- z'*z */
151: }
152: KSPCheckNorm(ksp,rnorm);
153: ksp->rnorm = rnorm;
154: KSPLogResidualHistory(ksp,rnorm);
155: KSPMonitor(ksp,i,rnorm);
156: }
157: if (ksp->its >= ksp->max_it) {
158: if (ksp->normtype != KSP_NORM_NONE) {
159: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
160: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
161: } else {
162: ksp->reason = KSP_CONVERGED_ITS;
163: }
164: }
165: }
166: return(0);
167: }
169: PetscErrorCode KSPView_Richardson(KSP ksp,PetscViewer viewer)
170: {
171: KSP_Richardson *richardsonP = (KSP_Richardson*)ksp->data;
173: PetscBool iascii;
176: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
177: if (iascii) {
178: if (richardsonP->selfscale) {
179: PetscViewerASCIIPrintf(viewer," using self-scale best computed damping factor\n");
180: } else {
181: PetscViewerASCIIPrintf(viewer," damping factor=%g\n",(double)richardsonP->scale);
182: }
183: }
184: return(0);
185: }
187: PetscErrorCode KSPSetFromOptions_Richardson(PetscOptionItems *PetscOptionsObject,KSP ksp)
188: {
189: KSP_Richardson *rich = (KSP_Richardson*)ksp->data;
191: PetscReal tmp;
192: PetscBool flg,flg2;
195: PetscOptionsHead(PetscOptionsObject,"KSP Richardson Options");
196: PetscOptionsReal("-ksp_richardson_scale","damping factor","KSPRichardsonSetScale",rich->scale,&tmp,&flg);
197: if (flg) { KSPRichardsonSetScale(ksp,tmp); }
198: PetscOptionsBool("-ksp_richardson_self_scale","dynamically determine optimal damping factor","KSPRichardsonSetSelfScale",rich->selfscale,&flg2,&flg);
199: if (flg) { KSPRichardsonSetSelfScale(ksp,flg2); }
200: PetscOptionsTail();
201: return(0);
202: }
204: PetscErrorCode KSPDestroy_Richardson(KSP ksp)
205: {
209: PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetScale_C",NULL);
210: KSPDestroyDefault(ksp);
211: return(0);
212: }
214: static PetscErrorCode KSPRichardsonSetScale_Richardson(KSP ksp,PetscReal scale)
215: {
216: KSP_Richardson *richardsonP;
219: richardsonP = (KSP_Richardson*)ksp->data;
220: richardsonP->scale = scale;
221: return(0);
222: }
224: static PetscErrorCode KSPRichardsonSetSelfScale_Richardson(KSP ksp,PetscBool selfscale)
225: {
226: KSP_Richardson *richardsonP;
229: richardsonP = (KSP_Richardson*)ksp->data;
230: richardsonP->selfscale = selfscale;
231: return(0);
232: }
234: /*MC
235: KSPRICHARDSON - The preconditioned Richardson iterative method
237: Options Database Keys:
238: . -ksp_richardson_scale - damping factor on the correction (defaults to 1.0)
240: Level: beginner
242: Notes:
243: x^{n+1} = x^{n} + scale*B(b - A x^{n})
245: Here B is the Section 1.5 Writing Application Codes with PETSc of the preconditioner
247: This method often (usually) will not converge unless scale is very small.
249: Notes:
250: For some preconditioners, currently SOR, the convergence test is skipped to improve speed,
251: thus it always iterates the maximum number of iterations you've selected. When -ksp_monitor
252: (or any other monitor) is turned on, the norm is computed at each iteration and so the convergence test is run unless
253: you specifically call KSPSetNormType(ksp,KSP_NORM_NONE);
255: For some preconditioners, currently PCMG and PCHYPRE with BoomerAMG if -ksp_monitor (and also
256: any other monitor) is not turned on then the convergence test is done by the preconditioner itself and
257: so the solver may run more or fewer iterations then if -ksp_monitor is selected.
259: Supports only left preconditioning
261: If using direct solvers such as PCLU and PCCHOLESKY one generally uses KSPPREONLY which uses exactly one iteration
263: $ -ksp_type richardson -pc_type jacobi gives one classically Jacobi preconditioning
265: References:
266: . 1. - L. F. Richardson, "The Approximate Arithmetical Solution by Finite Differences of Physical Problems Involving
267: Differential Equations, with an Application to the Stresses in a Masonry Dam",
268: Philosophical Transactions of the Royal Society of London. Series A,
269: Containing Papers of a Mathematical or Physical Character, Vol. 210, 1911 (1911).
271: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
272: KSPRichardsonSetScale(), KSPPREONLY
274: M*/
276: PETSC_EXTERN PetscErrorCode KSPCreate_Richardson(KSP ksp)
277: {
279: KSP_Richardson *richardsonP;
282: PetscNewLog(ksp,&richardsonP);
283: ksp->data = (void*)richardsonP;
285: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
286: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
287: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
289: ksp->ops->setup = KSPSetUp_Richardson;
290: ksp->ops->solve = KSPSolve_Richardson;
291: ksp->ops->destroy = KSPDestroy_Richardson;
292: ksp->ops->buildsolution = KSPBuildSolutionDefault;
293: ksp->ops->buildresidual = KSPBuildResidualDefault;
294: ksp->ops->view = KSPView_Richardson;
295: ksp->ops->setfromoptions = KSPSetFromOptions_Richardson;
297: PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetScale_C",KSPRichardsonSetScale_Richardson);
298: PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetSelfScale_C",KSPRichardsonSetSelfScale_Richardson);
300: richardsonP->scale = 1.0;
301: return(0);
302: }