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++) {
80: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
81: VecNorm(r,NORM_2,&rnorm); /* rnorm <- r'*r */
82: } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
83: VecNorm(z,NORM_2,&rnorm); /* rnorm <- z'*z */
84: } else rnorm = 0.0;
86: KSPCheckNorm(ksp,rnorm);
87: ksp->rnorm = rnorm;
88: KSPMonitor(ksp,i,rnorm);
89: KSPLogResidualHistory(ksp,rnorm);
90: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
91: if (ksp->reason) break;
92: KSP_PCApplyBAorAB(ksp,z,y,w); /* y = BAz = BABr */
93: VecDotNorm2(z,y,&rdot,&abr)); /* rdot = (Br)^T(BABR; abr = (BABr)^T (BABr) */
94: scale = rdot/abr;
95: PetscInfo(ksp,"Self-scale factor %g\n",(double)PetscRealPart(scale));
96: VecAXPY(x,scale,z); /* x <- x + scale z */
97: VecAXPY(r,-scale,w); /* r <- r - scale*Az */
98: VecAXPY(z,-scale,y); /* z <- z - scale*y */
99: ksp->its++;
100: }
101: } else {
102: for (i=0; i<maxit; i++) {
104: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
105: VecNorm(r,NORM_2,&rnorm); /* rnorm <- r'*r */
106: } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
107: KSP_PCApply(ksp,r,z); /* z <- B r */
108: VecNorm(z,NORM_2,&rnorm); /* rnorm <- z'*z */
109: } else rnorm = 0.0;
110: ksp->rnorm = rnorm;
111: KSPMonitor(ksp,i,rnorm);
112: KSPLogResidualHistory(ksp,rnorm);
113: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
114: if (ksp->reason) break;
115: if (ksp->normtype != KSP_NORM_PRECONDITIONED) {
116: KSP_PCApply(ksp,r,z); /* z <- B r */
117: }
119: VecAXPY(x,richardsonP->scale,z); /* x <- x + scale z */
120: ksp->its++;
122: if (i+1 < maxit || ksp->normtype != KSP_NORM_NONE) {
123: KSP_MatMult(ksp,Amat,x,r); /* r <- b - Ax */
124: VecAYPX(r,-1.0,b);
125: }
126: }
127: }
128: if (!ksp->reason) {
129: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
130: VecNorm(r,NORM_2,&rnorm); /* rnorm <- r'*r */
131: } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
132: KSP_PCApply(ksp,r,z); /* z <- B r */
133: VecNorm(z,NORM_2,&rnorm); /* rnorm <- z'*z */
134: } else rnorm = 0.0;
136: KSPCheckNorm(ksp,rnorm);
137: ksp->rnorm = rnorm;
138: KSPLogResidualHistory(ksp,rnorm);
139: KSPMonitor(ksp,i,rnorm);
140: if (ksp->its >= ksp->max_it) {
141: if (ksp->normtype != KSP_NORM_NONE) {
142: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
143: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
144: } else {
145: ksp->reason = KSP_CONVERGED_ITS;
146: }
147: }
148: }
149: return 0;
150: }
152: PetscErrorCode KSPView_Richardson(KSP ksp,PetscViewer viewer)
153: {
154: KSP_Richardson *richardsonP = (KSP_Richardson*)ksp->data;
155: PetscBool iascii;
157: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
158: if (iascii) {
159: if (richardsonP->selfscale) {
160: PetscViewerASCIIPrintf(viewer," using self-scale best computed damping factor\n");
161: } else {
162: PetscViewerASCIIPrintf(viewer," damping factor=%g\n",(double)richardsonP->scale);
163: }
164: }
165: return 0;
166: }
168: PetscErrorCode KSPSetFromOptions_Richardson(PetscOptionItems *PetscOptionsObject,KSP ksp)
169: {
170: KSP_Richardson *rich = (KSP_Richardson*)ksp->data;
171: PetscReal tmp;
172: PetscBool flg,flg2;
174: PetscOptionsHead(PetscOptionsObject,"KSP Richardson Options");
175: PetscOptionsReal("-ksp_richardson_scale","damping factor","KSPRichardsonSetScale",rich->scale,&tmp,&flg);
176: if (flg) KSPRichardsonSetScale(ksp,tmp);
177: PetscOptionsBool("-ksp_richardson_self_scale","dynamically determine optimal damping factor","KSPRichardsonSetSelfScale",rich->selfscale,&flg2,&flg);
178: if (flg) KSPRichardsonSetSelfScale(ksp,flg2);
179: PetscOptionsTail();
180: return 0;
181: }
183: PetscErrorCode KSPDestroy_Richardson(KSP ksp)
184: {
185: PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetScale_C",NULL);
186: KSPDestroyDefault(ksp);
187: return 0;
188: }
190: static PetscErrorCode KSPRichardsonSetScale_Richardson(KSP ksp,PetscReal scale)
191: {
192: KSP_Richardson *richardsonP;
194: richardsonP = (KSP_Richardson*)ksp->data;
195: richardsonP->scale = scale;
196: return 0;
197: }
199: static PetscErrorCode KSPRichardsonSetSelfScale_Richardson(KSP ksp,PetscBool selfscale)
200: {
201: KSP_Richardson *richardsonP;
203: richardsonP = (KSP_Richardson*)ksp->data;
204: richardsonP->selfscale = selfscale;
205: return 0;
206: }
208: static PetscErrorCode KSPBuildResidual_Richardson(KSP ksp,Vec t,Vec v,Vec *V)
209: {
210: if (ksp->normtype == KSP_NORM_NONE) {
211: KSPBuildResidualDefault(ksp,t,v,V);
212: } else {
213: VecCopy(ksp->work[0],v);
214: *V = v;
215: }
216: return 0;
217: }
219: /*MC
220: KSPRICHARDSON - The preconditioned Richardson iterative method
222: Options Database Keys:
223: . -ksp_richardson_scale - damping factor on the correction (defaults to 1.0)
225: Level: beginner
227: Notes:
228: x^{n+1} = x^{n} + scale*B(b - A x^{n})
230: Here B is the application of the preconditioner
232: This method often (usually) will not converge unless scale is very small.
234: Notes:
235: For some preconditioners, currently SOR, 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 which uses exactly one iteration
248: $ -ksp_type richardson -pc_type jacobi gives one classically Jacobi preconditioning
250: References:
251: . * - L. F. Richardson, "The Approximate Arithmetical Solution by Finite Differences of Physical Problems Involving
252: Differential Equations, with an Application to the Stresses in a Masonry Dam",
253: Philosophical Transactions of the Royal Society of London. Series A,
254: Containing Papers of a Mathematical or Physical Character, Vol. 210, 1911 (1911).
256: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
257: KSPRichardsonSetScale(), KSPPREONLY
259: M*/
261: PETSC_EXTERN PetscErrorCode KSPCreate_Richardson(KSP ksp)
262: {
263: KSP_Richardson *richardsonP;
265: PetscNewLog(ksp,&richardsonP);
266: ksp->data = (void*)richardsonP;
268: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
269: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
270: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
272: ksp->ops->setup = KSPSetUp_Richardson;
273: ksp->ops->solve = KSPSolve_Richardson;
274: ksp->ops->destroy = KSPDestroy_Richardson;
275: ksp->ops->buildsolution = KSPBuildSolutionDefault;
276: ksp->ops->buildresidual = KSPBuildResidual_Richardson;
277: ksp->ops->view = KSPView_Richardson;
278: ksp->ops->setfromoptions = KSPSetFromOptions_Richardson;
280: PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetScale_C",KSPRichardsonSetScale_Richardson);
281: PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetSelfScale_C",KSPRichardsonSetSelfScale_Richardson);
283: richardsonP->scale = 1.0;
284: return 0;
285: }