Actual source code: rich.c

petsc-master 2020-08-06
Report Typos and Errors

  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->converged == KSPConvergedDefault || ksp->converged == KSPConvergedSkip) && !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 {
 69:     PetscInfo(ksp,"KSPSolve_Richardson: Warning, skipping optimized PCApplyRichardson()\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:       } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
 87:         VecNorm(z,NORM_2,&rnorm); /*   rnorm <- z'*z     */
 88:       } else rnorm = 0.0;

 90:       KSPCheckNorm(ksp,rnorm);
 91:       ksp->rnorm = rnorm;
 92:       KSPMonitor(ksp,i,rnorm);
 93:       KSPLogResidualHistory(ksp,rnorm);
 94:       (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
 95:       if (ksp->reason) break;
 96:       KSP_PCApplyBAorAB(ksp,z,y,w); /* y = BAz = BABr */
 97:       VecDotNorm2(z,y,&rdot,&abr);   /*   rdot = (Br)^T(BABR); abr = (BABr)^T (BABr) */
 98:       scale = rdot/abr;
 99:       PetscInfo1(ksp,"Self-scale factor %g\n",(double)PetscRealPart(scale));
100:       VecAXPY(x,scale,z);   /*   x  <- x + scale z */
101:       VecAXPY(r,-scale,w);  /*  r <- r - scale*Az */
102:       VecAXPY(z,-scale,y);  /*  z <- z - scale*y */
103:       ksp->its++;
104:     }
105:   } else {
106:     for (i=0; i<maxit; i++) {

108:       if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
109:         VecNorm(r,NORM_2,&rnorm); /*   rnorm <- r'*r     */
110:       } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
111:         KSP_PCApply(ksp,r,z);    /*   z <- B r          */
112:         VecNorm(z,NORM_2,&rnorm); /*   rnorm <- z'*z     */
113:       } else rnorm = 0.0;
114:       ksp->rnorm = rnorm;
115:       KSPMonitor(ksp,i,rnorm);
116:       KSPLogResidualHistory(ksp,rnorm);
117:       (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
118:       if (ksp->reason) break;
119:       if (ksp->normtype != KSP_NORM_PRECONDITIONED) {
120:         KSP_PCApply(ksp,r,z);    /*   z <- B r          */
121:       }

123:       VecAXPY(x,richardsonP->scale,z);    /*   x  <- x + scale z */
124:       ksp->its++;

126:       if (i+1 < maxit || ksp->normtype != KSP_NORM_NONE) {
127:         KSP_MatMult(ksp,Amat,x,r);      /*   r  <- b - Ax      */
128:         VecAYPX(r,-1.0,b);
129:       }
130:     }
131:   }
132:   if (!ksp->reason) {
133:     if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
134:       VecNorm(r,NORM_2,&rnorm);     /*   rnorm <- r'*r     */
135:     } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
136:       KSP_PCApply(ksp,r,z);   /*   z <- B r          */
137:       VecNorm(z,NORM_2,&rnorm);     /*   rnorm <- z'*z     */
138:     } else rnorm = 0.0;

140:     KSPCheckNorm(ksp,rnorm);
141:     ksp->rnorm = rnorm;
142:     KSPLogResidualHistory(ksp,rnorm);
143:     KSPMonitor(ksp,i,rnorm);
144:     if (ksp->its >= ksp->max_it) {
145:       if (ksp->normtype != KSP_NORM_NONE) {
146:         (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
147:         if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
148:       } else {
149:         ksp->reason = KSP_CONVERGED_ITS;
150:       }
151:     }
152:   }
153:   return(0);
154: }

156: PetscErrorCode KSPView_Richardson(KSP ksp,PetscViewer viewer)
157: {
158:   KSP_Richardson *richardsonP = (KSP_Richardson*)ksp->data;
160:   PetscBool      iascii;

163:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
164:   if (iascii) {
165:     if (richardsonP->selfscale) {
166:       PetscViewerASCIIPrintf(viewer,"  using self-scale best computed damping factor\n");
167:     } else {
168:       PetscViewerASCIIPrintf(viewer,"  damping factor=%g\n",(double)richardsonP->scale);
169:     }
170:   }
171:   return(0);
172: }

174: PetscErrorCode KSPSetFromOptions_Richardson(PetscOptionItems *PetscOptionsObject,KSP ksp)
175: {
176:   KSP_Richardson *rich = (KSP_Richardson*)ksp->data;
178:   PetscReal      tmp;
179:   PetscBool      flg,flg2;

182:   PetscOptionsHead(PetscOptionsObject,"KSP Richardson Options");
183:   PetscOptionsReal("-ksp_richardson_scale","damping factor","KSPRichardsonSetScale",rich->scale,&tmp,&flg);
184:   if (flg) { KSPRichardsonSetScale(ksp,tmp); }
185:   PetscOptionsBool("-ksp_richardson_self_scale","dynamically determine optimal damping factor","KSPRichardsonSetSelfScale",rich->selfscale,&flg2,&flg);
186:   if (flg) { KSPRichardsonSetSelfScale(ksp,flg2); }
187:   PetscOptionsTail();
188:   return(0);
189: }

191: PetscErrorCode KSPDestroy_Richardson(KSP ksp)
192: {

196:   PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetScale_C",NULL);
197:   KSPDestroyDefault(ksp);
198:   return(0);
199: }

201: static PetscErrorCode  KSPRichardsonSetScale_Richardson(KSP ksp,PetscReal scale)
202: {
203:   KSP_Richardson *richardsonP;

206:   richardsonP        = (KSP_Richardson*)ksp->data;
207:   richardsonP->scale = scale;
208:   return(0);
209: }

211: static PetscErrorCode  KSPRichardsonSetSelfScale_Richardson(KSP ksp,PetscBool selfscale)
212: {
213:   KSP_Richardson *richardsonP;

216:   richardsonP            = (KSP_Richardson*)ksp->data;
217:   richardsonP->selfscale = selfscale;
218:   return(0);
219: }

221: static PetscErrorCode KSPBuildResidual_Richardson(KSP ksp,Vec t,Vec v,Vec *V)
222: {

226:   if (ksp->normtype == KSP_NORM_NONE) {
227:     KSPBuildResidualDefault(ksp,t,v,V);
228:   } else {
229:     VecCopy(ksp->work[0],v);
230:     *V   = v;
231:   }
232:   return(0);
233: }

235: /*MC
236:      KSPRICHARDSON - The preconditioned Richardson iterative method

238:    Options Database Keys:
239: .   -ksp_richardson_scale - damping factor on the correction (defaults to 1.0)

241:    Level: beginner

243:    Notes:
244:     x^{n+1} = x^{n} + scale*B(b - A x^{n})

246:           Here B is the Section 1.5 Writing Application Codes with PETSc of the preconditioner

248:           This method often (usually) will not converge unless scale is very small.

250:    Notes:
251:     For some preconditioners, currently SOR, the convergence test is skipped to improve speed,
252:     thus it always iterates the maximum number of iterations you've selected. When -ksp_monitor
253:     (or any other monitor) is turned on, the norm is computed at each iteration and so the convergence test is run unless
254:     you specifically call KSPSetNormType(ksp,KSP_NORM_NONE);

256:          For some preconditioners, currently PCMG and PCHYPRE with BoomerAMG if -ksp_monitor (and also
257:     any other monitor) is not turned on then the convergence test is done by the preconditioner itself and
258:     so the solver may run more or fewer iterations then if -ksp_monitor is selected.

260:     Supports only left preconditioning

262:     If using direct solvers such as PCLU and PCCHOLESKY one generally uses KSPPREONLY which uses exactly one iteration

264: $    -ksp_type richardson -pc_type jacobi gives one classically Jacobi preconditioning

266:   References:
267: .  1. - L. F. Richardson, "The Approximate Arithmetical Solution by Finite Differences of Physical Problems Involving
268:    Differential Equations, with an Application to the Stresses in a Masonry Dam",
269:   Philosophical Transactions of the Royal Society of London. Series A,
270:   Containing Papers of a Mathematical or Physical Character, Vol. 210, 1911 (1911).

272: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
273:            KSPRichardsonSetScale(), KSPPREONLY

275: M*/

277: PETSC_EXTERN PetscErrorCode KSPCreate_Richardson(KSP ksp)
278: {
280:   KSP_Richardson *richardsonP;

283:   PetscNewLog(ksp,&richardsonP);
284:   ksp->data = (void*)richardsonP;

286:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
287:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
288:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);

290:   ksp->ops->setup          = KSPSetUp_Richardson;
291:   ksp->ops->solve          = KSPSolve_Richardson;
292:   ksp->ops->destroy        = KSPDestroy_Richardson;
293:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
294:   ksp->ops->buildresidual  = KSPBuildResidual_Richardson;
295:   ksp->ops->view           = KSPView_Richardson;
296:   ksp->ops->setfromoptions = KSPSetFromOptions_Richardson;

298:   PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetScale_C",KSPRichardsonSetScale_Richardson);
299:   PetscObjectComposeFunction((PetscObject)ksp,"KSPRichardsonSetSelfScale_C",KSPRichardsonSetSelfScale_Richardson);

301:   richardsonP->scale = 1.0;
302:   return(0);
303: }