Actual source code: rich.c
petsc-3.7.5 2017-01-01
2: /*
3: This implements Richardson Iteration.
4: */
5: #include <petsc/private/kspimpl.h> /*I "petscksp.h" I*/
6: #include <../src/ksp/ksp/impls/rich/richardsonimpl.h>
10: PetscErrorCode KSPSetUp_Richardson(KSP ksp)
11: {
13: KSP_Richardson *richardsonP = (KSP_Richardson*)ksp->data;
16: if (richardsonP->selfscale) {
17: KSPSetWorkVecs(ksp,4);
18: } else {
19: KSPSetWorkVecs(ksp,2);
20: }
21: return(0);
22: }
26: PetscErrorCode KSPSolve_Richardson(KSP ksp)
27: {
29: PetscInt i,maxit;
30: PetscReal rnorm = 0.0,abr;
31: PetscScalar scale,rdot;
32: Vec x,b,r,z,w = NULL,y = NULL;
33: PetscInt xs, ws;
34: Mat Amat,Pmat;
35: KSP_Richardson *richardsonP = (KSP_Richardson*)ksp->data;
36: PetscBool exists,diagonalscale;
37: MatNullSpace nullsp;
40: PCGetDiagonalScale(ksp->pc,&diagonalscale);
41: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
43: ksp->its = 0;
45: PCGetOperators(ksp->pc,&Amat,&Pmat);
46: x = ksp->vec_sol;
47: b = ksp->vec_rhs;
48: VecGetSize(x,&xs);
49: VecGetSize(ksp->work[0],&ws);
50: if (xs != ws) {
51: if (richardsonP->selfscale) {
52: KSPSetWorkVecs(ksp,4);
53: } else {
54: KSPSetWorkVecs(ksp,2);
55: }
56: }
57: r = ksp->work[0];
58: z = ksp->work[1];
59: if (richardsonP->selfscale) {
60: w = ksp->work[2];
61: y = ksp->work[3];
62: }
63: maxit = ksp->max_it;
65: /* if user has provided fast Richardson code use that */
66: PCApplyRichardsonExists(ksp->pc,&exists);
67: MatGetNullSpace(Pmat,&nullsp);
68: if (exists && richardsonP->scale == 1.0 && !ksp->numbermonitors && !ksp->transpose_solve & !nullsp) {
69: PCRichardsonConvergedReason reason;
70: PCApplyRichardson(ksp->pc,b,x,r,ksp->rtol,ksp->abstol,ksp->divtol,maxit,ksp->guess_zero,&ksp->its,&reason);
71: ksp->reason = (KSPConvergedReason)reason;
72: return(0);
73: } else if (exists && !ksp->numbermonitors && !ksp->transpose_solve & !nullsp) {
74: PetscInfo(ksp,"KSPSolve_Richardson: Warning, skipping optimized PCApplyRichardson() because scale factor is not 1.0\n");
75: }
77: scale = richardsonP->scale;
79: if (!ksp->guess_zero) { /* r <- b - A x */
80: KSP_MatMult(ksp,Amat,x,r);
81: VecAYPX(r,-1.0,b);
82: } else {
83: VecCopy(b,r);
84: }
86: ksp->its = 0;
87: if (richardsonP->selfscale) {
88: KSP_PCApply(ksp,r,z); /* z <- B r */
89: for (i=0; i<maxit; i++) {
91: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
92: VecNorm(r,NORM_2,&rnorm); /* rnorm <- r'*r */
93: KSPMonitor(ksp,i,rnorm);
94: ksp->rnorm = rnorm;
95: KSPLogResidualHistory(ksp,rnorm);
96: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
97: if (ksp->reason) break;
98: } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
99: VecNorm(z,NORM_2,&rnorm); /* rnorm <- z'*z */
100: KSPMonitor(ksp,i,rnorm);
101: ksp->rnorm = rnorm;
102: KSPLogResidualHistory(ksp,rnorm);
103: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
104: if (ksp->reason) break;
105: }
106: KSP_PCApplyBAorAB(ksp,z,y,w); /* y = BAz = BABr */
107: VecDotNorm2(z,y,&rdot,&abr); /* rdot = (Br)^T(BABR); abr = (BABr)^T (BABr) */
108: scale = rdot/abr;
109: PetscInfo1(ksp,"Self-scale factor %g\n",(double)PetscRealPart(scale));
110: VecAXPY(x,scale,z); /* x <- x + scale z */
111: VecAXPY(r,-scale,w); /* r <- r - scale*Az */
112: VecAXPY(z,-scale,y); /* z <- z - scale*y */
113: ksp->its++;
114: }
115: } else {
116: for (i=0; i<maxit; i++) {
118: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
119: VecNorm(r,NORM_2,&rnorm); /* rnorm <- r'*r */
120: KSPMonitor(ksp,i,rnorm);
121: ksp->rnorm = rnorm;
122: KSPLogResidualHistory(ksp,rnorm);
123: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
124: if (ksp->reason) break;
125: }
127: KSP_PCApply(ksp,r,z); /* z <- B r */
129: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
130: VecNorm(z,NORM_2,&rnorm); /* rnorm <- z'*z */
131: KSPMonitor(ksp,i,rnorm);
132: ksp->rnorm = rnorm;
133: KSPLogResidualHistory(ksp,rnorm);
134: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
135: if (ksp->reason) break;
136: }
138: VecAXPY(x,richardsonP->scale,z); /* x <- x + scale z */
139: ksp->its++;
141: if (i+1 < maxit || ksp->normtype != KSP_NORM_NONE) {
142: KSP_MatMult(ksp,Amat,x,r); /* r <- b - Ax */
143: VecAYPX(r,-1.0,b);
144: }
145: }
146: }
147: if (!ksp->reason) {
148: if (ksp->normtype != KSP_NORM_NONE) {
149: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
150: VecNorm(r,NORM_2,&rnorm); /* rnorm <- r'*r */
151: } else {
152: KSP_PCApply(ksp,r,z); /* z <- B r */
153: VecNorm(z,NORM_2,&rnorm); /* rnorm <- z'*z */
154: }
155: ksp->rnorm = rnorm;
156: KSPLogResidualHistory(ksp,rnorm);
157: KSPMonitor(ksp,i,rnorm);
158: }
159: if (ksp->its >= ksp->max_it) {
160: if (ksp->normtype != KSP_NORM_NONE) {
161: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
162: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
163: } else {
164: ksp->reason = KSP_CONVERGED_ITS;
165: }
166: }
167: }
168: return(0);
169: }
173: PetscErrorCode KSPView_Richardson(KSP ksp,PetscViewer viewer)
174: {
175: KSP_Richardson *richardsonP = (KSP_Richardson*)ksp->data;
177: PetscBool iascii;
180: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
181: if (iascii) {
182: if (richardsonP->selfscale) {
183: PetscViewerASCIIPrintf(viewer," Richardson: using self-scale best computed damping factor\n");
184: } else {
185: PetscViewerASCIIPrintf(viewer," Richardson: damping factor=%g\n",(double)richardsonP->scale);
186: }
187: }
188: return(0);
189: }
193: PetscErrorCode KSPSetFromOptions_Richardson(PetscOptionItems *PetscOptionsObject,KSP ksp)
194: {
195: KSP_Richardson *rich = (KSP_Richardson*)ksp->data;
197: PetscReal tmp;
198: PetscBool flg,flg2;
201: PetscOptionsHead(PetscOptionsObject,"KSP Richardson Options");
202: PetscOptionsReal("-ksp_richardson_scale","damping factor","KSPRichardsonSetScale",rich->scale,&tmp,&flg);
203: if (flg) { KSPRichardsonSetScale(ksp,tmp); }
204: PetscOptionsBool("-ksp_richardson_self_scale","dynamically determine optimal damping factor","KSPRichardsonSetSelfScale",rich->selfscale,&flg2,&flg);
205: if (flg) { KSPRichardsonSetSelfScale(ksp,flg2); }
206: PetscOptionsTail();
207: return(0);
208: }
212: PetscErrorCode KSPDestroy_Richardson(KSP ksp)
213: {
217: PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetScale_C",NULL);
218: KSPDestroyDefault(ksp);
219: return(0);
220: }
224: static PetscErrorCode KSPRichardsonSetScale_Richardson(KSP ksp,PetscReal scale)
225: {
226: KSP_Richardson *richardsonP;
229: richardsonP = (KSP_Richardson*)ksp->data;
230: richardsonP->scale = scale;
231: return(0);
232: }
236: static PetscErrorCode KSPRichardsonSetSelfScale_Richardson(KSP ksp,PetscBool selfscale)
237: {
238: KSP_Richardson *richardsonP;
241: richardsonP = (KSP_Richardson*)ksp->data;
242: richardsonP->selfscale = selfscale;
243: return(0);
244: }
246: /*MC
247: KSPRICHARDSON - The preconditioned Richardson iterative method
249: Options Database Keys:
250: . -ksp_richardson_scale - damping factor on the correction (defaults to 1.0)
252: Level: beginner
254: Notes: x^{n+1} = x^{n} + scale*B(b - A x^{n})
256: Here B is the application of the preconditioner
258: This method often (usually) will not converge unless scale is very small.
260: Notes: For some preconditioners, currently SOR, the convergence test is skipped to improve speed,
261: thus it always iterates the maximum number of iterations you've selected. When -ksp_monitor
262: (or any other monitor) is turned on, the norm is computed at each iteration and so the convergence test is run unless
263: you specifically call KSPSetNormType(ksp,KSP_NORM_NONE);
265: For some preconditioners, currently PCMG and PCHYPRE with BoomerAMG if -ksp_monitor (and also
266: any other monitor) is not turned on then the convergence test is done by the preconditioner itself and
267: so the solver may run more or fewer iterations then if -ksp_monitor is selected.
269: Supports only left preconditioning
271: References:
272: . 1. - L. F. Richardson, "The Approximate Arithmetical Solution by Finite Differences of Physical Problems Involving
273: Differential Equations, with an Application to the Stresses in a Masonry Dam",
274: Philosophical Transactions of the Royal Society of London. Series A,
275: Containing Papers of a Mathematical or Physical Character, Vol. 210, 1911 (1911).
277: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
278: KSPRichardsonSetScale()
280: M*/
284: PETSC_EXTERN PetscErrorCode KSPCreate_Richardson(KSP ksp)
285: {
287: KSP_Richardson *richardsonP;
290: PetscNewLog(ksp,&richardsonP);
291: ksp->data = (void*)richardsonP;
293: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
294: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
296: ksp->ops->setup = KSPSetUp_Richardson;
297: ksp->ops->solve = KSPSolve_Richardson;
298: ksp->ops->destroy = KSPDestroy_Richardson;
299: ksp->ops->buildsolution = KSPBuildSolutionDefault;
300: ksp->ops->buildresidual = KSPBuildResidualDefault;
301: ksp->ops->view = KSPView_Richardson;
302: ksp->ops->setfromoptions = KSPSetFromOptions_Richardson;
304: PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetScale_C",KSPRichardsonSetScale_Richardson);
305: PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetSelfScale_C",KSPRichardsonSetSelfScale_Richardson);
307: richardsonP->scale = 1.0;
308: return(0);
309: }