Actual source code: rich.c
petsc-3.5.4 2015-05-23
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;
39: PCGetDiagonalScale(ksp->pc,&diagonalscale);
40: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
42: ksp->its = 0;
44: PCGetOperators(ksp->pc,&Amat,&Pmat);
45: x = ksp->vec_sol;
46: b = ksp->vec_rhs;
47: VecGetSize(x,&xs);
48: VecGetSize(ksp->work[0],&ws);
49: if (xs != ws) {
50: if (richardsonP->selfscale) {
51: KSPSetWorkVecs(ksp,4);
52: } else {
53: KSPSetWorkVecs(ksp,2);
54: }
55: }
56: r = ksp->work[0];
57: z = ksp->work[1];
58: if (richardsonP->selfscale) {
59: w = ksp->work[2];
60: y = ksp->work[3];
61: }
62: maxit = ksp->max_it;
64: /* if user has provided fast Richardson code use that */
65: PCApplyRichardsonExists(ksp->pc,&exists);
66: if (exists && !ksp->numbermonitors && !ksp->transpose_solve & !ksp->nullsp) {
67: PCRichardsonConvergedReason reason;
68: PCApplyRichardson(ksp->pc,b,x,r,ksp->rtol,ksp->abstol,ksp->divtol,maxit,ksp->guess_zero,&ksp->its,&reason);
69: ksp->reason = (KSPConvergedReason)reason;
70: return(0);
71: }
73: scale = richardsonP->scale;
75: if (!ksp->guess_zero) { /* r <- b - A x */
76: KSP_MatMult(ksp,Amat,x,r);
77: VecAYPX(r,-1.0,b);
78: } else {
79: VecCopy(b,r);
80: }
82: ksp->its = 0;
83: if (richardsonP->selfscale) {
84: KSP_PCApply(ksp,r,z); /* z <- B r */
85: for (i=0; i<maxit; i++) {
87: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
88: VecNorm(r,NORM_2,&rnorm); /* rnorm <- r'*r */
89: KSPMonitor(ksp,i,rnorm);
90: ksp->rnorm = rnorm;
91: KSPLogResidualHistory(ksp,rnorm);
92: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
93: if (ksp->reason) break;
94: } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
95: VecNorm(z,NORM_2,&rnorm); /* rnorm <- z'*z */
96: KSPMonitor(ksp,i,rnorm);
97: ksp->rnorm = rnorm;
98: KSPLogResidualHistory(ksp,rnorm);
99: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
100: if (ksp->reason) break;
101: }
102: KSP_PCApplyBAorAB(ksp,z,y,w); /* y = BAz = BABr */
103: VecDotNorm2(z,y,&rdot,&abr); /* rdot = (Br)^T(BABR); abr = (BABr)^T (BABr) */
104: scale = rdot/abr;
105: PetscInfo1(ksp,"Self-scale factor %g\n",(double)PetscRealPart(scale));
106: VecAXPY(x,scale,z); /* x <- x + scale z */
107: VecAXPY(r,-scale,w); /* r <- r - scale*Az */
108: VecAXPY(z,-scale,y); /* z <- z - scale*y */
109: ksp->its++;
110: }
111: } else {
112: for (i=0; i<maxit; i++) {
114: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
115: VecNorm(r,NORM_2,&rnorm); /* rnorm <- r'*r */
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: KSPMonitor(ksp,i,rnorm);
128: ksp->rnorm = rnorm;
129: KSPLogResidualHistory(ksp,rnorm);
130: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
131: if (ksp->reason) break;
132: }
134: VecAXPY(x,scale,z); /* x <- x + scale z */
135: ksp->its++;
137: if (i+1 < maxit || ksp->normtype != KSP_NORM_NONE) {
138: KSP_MatMult(ksp,Amat,x,r); /* r <- b - Ax */
139: VecAYPX(r,-1.0,b);
140: }
141: }
142: }
143: if (!ksp->reason) {
144: if (ksp->normtype != KSP_NORM_NONE) {
145: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
146: VecNorm(r,NORM_2,&rnorm); /* rnorm <- r'*r */
147: } else {
148: KSP_PCApply(ksp,r,z); /* z <- B r */
149: VecNorm(z,NORM_2,&rnorm); /* rnorm <- z'*z */
150: }
151: ksp->rnorm = rnorm;
152: KSPLogResidualHistory(ksp,rnorm);
153: KSPMonitor(ksp,i,rnorm);
154: }
155: if (ksp->its >= ksp->max_it) {
156: if (ksp->normtype != KSP_NORM_NONE) {
157: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
158: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
159: } else {
160: ksp->reason = KSP_CONVERGED_ITS;
161: }
162: }
163: }
164: return(0);
165: }
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," Richardson: using self-scale best computed damping factor\n");
180: } else {
181: PetscViewerASCIIPrintf(viewer," Richardson: damping factor=%g\n",(double)richardsonP->scale);
182: }
183: }
184: return(0);
185: }
189: PetscErrorCode KSPSetFromOptions_Richardson(KSP ksp)
190: {
191: KSP_Richardson *rich = (KSP_Richardson*)ksp->data;
193: PetscReal tmp;
194: PetscBool flg,flg2;
197: PetscOptionsHead("KSP Richardson Options");
198: PetscOptionsReal("-ksp_richardson_scale","damping factor","KSPRichardsonSetScale",rich->scale,&tmp,&flg);
199: if (flg) { KSPRichardsonSetScale(ksp,tmp); }
200: PetscOptionsBool("-ksp_richardson_self_scale","dynamically determine optimal damping factor","KSPRichardsonSetSelfScale",rich->selfscale,&flg2,&flg);
201: if (flg) { KSPRichardsonSetSelfScale(ksp,flg2); }
202: PetscOptionsTail();
203: return(0);
204: }
208: PetscErrorCode KSPDestroy_Richardson(KSP ksp)
209: {
213: PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetScale_C",NULL);
214: KSPDestroyDefault(ksp);
215: return(0);
216: }
220: static PetscErrorCode KSPRichardsonSetScale_Richardson(KSP ksp,PetscReal scale)
221: {
222: KSP_Richardson *richardsonP;
225: richardsonP = (KSP_Richardson*)ksp->data;
226: richardsonP->scale = scale;
227: return(0);
228: }
232: static PetscErrorCode KSPRichardsonSetSelfScale_Richardson(KSP ksp,PetscBool selfscale)
233: {
234: KSP_Richardson *richardsonP;
237: richardsonP = (KSP_Richardson*)ksp->data;
238: richardsonP->selfscale = selfscale;
239: return(0);
240: }
242: /*MC
243: KSPRICHARDSON - The preconditioned Richardson iterative method
245: Options Database Keys:
246: . -ksp_richardson_scale - damping factor on the correction (defaults to 1.0)
248: Level: beginner
250: Notes: x^{n+1} = x^{n} + scale*B(b - A x^{n})
252: Here B is the application of the preconditioner
254: This method often (usually) will not converge unless scale is very small.
256: Notes: For some preconditioners, currently SOR, the convergence test is skipped to improve speed,
257: thus it always iterates the maximum number of iterations you've selected. When -ksp_monitor
258: (or any other monitor) is turned on, the norm is computed at each iteration and so the convergence test is run unless
259: you specifically call KSPSetNormType(ksp,KSP_NORM_NONE);
261: For some preconditioners, currently PCMG and PCHYPRE with BoomerAMG if -ksp_monitor (and also
262: any other monitor) is not turned on then the convergence test is done by the preconditioner itself and
263: so the solver may run more or fewer iterations then if -ksp_monitor is selected.
265: Supports only left preconditioning
267: References:
268: "The Approximate Arithmetical Solution by Finite Differences of Physical Problems Involving
269: Differential Equations, with an Application to the Stresses in a Masonry Dam",
270: L. F. Richardson, Philosophical Transactions of the Royal Society of London. Series A,
271: Containing Papers of a Mathematical or Physical Character, Vol. 210, 1911 (1911), pp. 307-357.
273: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
274: KSPRichardsonSetScale()
276: M*/
280: PETSC_EXTERN PetscErrorCode KSPCreate_Richardson(KSP ksp)
281: {
283: KSP_Richardson *richardsonP;
286: PetscNewLog(ksp,&richardsonP);
287: ksp->data = (void*)richardsonP;
289: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
290: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
292: ksp->ops->setup = KSPSetUp_Richardson;
293: ksp->ops->solve = KSPSolve_Richardson;
294: ksp->ops->destroy = KSPDestroy_Richardson;
295: ksp->ops->buildsolution = KSPBuildSolutionDefault;
296: ksp->ops->buildresidual = KSPBuildResidualDefault;
297: ksp->ops->view = KSPView_Richardson;
298: ksp->ops->setfromoptions = KSPSetFromOptions_Richardson;
300: PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetScale_C",KSPRichardsonSetScale_Richardson);
301: PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetSelfScale_C",KSPRichardsonSetSelfScale_Richardson);
303: richardsonP->scale = 1.0;
304: return(0);
305: }