Actual source code: tr.c

petsc-3.11.3 2019-06-26
Report Typos and Errors

  2:  #include <../src/snes/impls/tr/trimpl.h>

  4: typedef struct {
  5:   SNES           snes;
  6:   /*  Information on the regular SNES convergence test; which may have been user provided */
  7:   PetscErrorCode (*convtest)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
  8:   PetscErrorCode (*convdestroy)(void*);
  9:   void           *convctx;
 10: } SNES_TR_KSPConverged_Ctx;

 12: static PetscErrorCode SNESTR_KSPConverged_Private(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *cctx)
 13: {
 14:   SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx*)cctx;
 15:   SNES                     snes = ctx->snes;
 16:   SNES_NEWTONTR            *neP = (SNES_NEWTONTR*)snes->data;
 17:   Vec                      x;
 18:   PetscReal                nrm;
 19:   PetscErrorCode           ierr;

 22:   (*ctx->convtest)(ksp,n,rnorm,reason,ctx->convctx);
 23:   if (*reason) {
 24:     PetscInfo2(snes,"Default or user provided convergence test KSP iterations=%D, rnorm=%g\n",n,(double)rnorm);
 25:   }
 26:   /* Determine norm of solution */
 27:   KSPBuildSolution(ksp,0,&x);
 28:   VecNorm(x,NORM_2,&nrm);
 29:   if (nrm >= neP->delta) {
 30:     PetscInfo2(snes,"Ending linear iteration early, delta=%g, length=%g\n",(double)neP->delta,(double)nrm);
 31:     *reason = KSP_CONVERGED_STEP_LENGTH;
 32:   }
 33:   return(0);
 34: }

 36: static PetscErrorCode SNESTR_KSPConverged_Destroy(void *cctx)
 37: {
 38:   SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx*)cctx;
 39:   PetscErrorCode           ierr;

 42:   (*ctx->convdestroy)(ctx->convctx);
 43:   PetscFree(ctx);
 44:   return(0);
 45: }

 47: /* ---------------------------------------------------------------- */
 48: /*
 49:    SNESTR_Converged_Private -test convergence JUST for
 50:    the trust region tolerance.

 52: */
 53: static PetscErrorCode SNESTR_Converged_Private(SNES snes,PetscInt it,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
 54: {
 55:   SNES_NEWTONTR  *neP = (SNES_NEWTONTR*)snes->data;

 59:   *reason = SNES_CONVERGED_ITERATING;
 60:   if (neP->delta < xnorm * snes->deltatol) {
 61:     PetscInfo3(snes,"Converged due to trust region param %g<%g*%g\n",(double)neP->delta,(double)xnorm,(double)snes->deltatol);
 62:     *reason = SNES_CONVERGED_TR_DELTA;
 63:   } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
 64:     PetscInfo1(snes,"Exceeded maximum number of function evaluations: %D\n",snes->max_funcs);
 65:     *reason = SNES_DIVERGED_FUNCTION_COUNT;
 66:   }
 67:   return(0);
 68: }


 71: /*
 72:    SNESSolve_NEWTONTR - Implements Newton's Method with a very simple trust
 73:    region approach for solving systems of nonlinear equations.


 76: */
 77: static PetscErrorCode SNESSolve_NEWTONTR(SNES snes)
 78: {
 79:   SNES_NEWTONTR            *neP = (SNES_NEWTONTR*)snes->data;
 80:   Vec                      X,F,Y,G,Ytmp;
 81:   PetscErrorCode           ierr;
 82:   PetscInt                 maxits,i,lits;
 83:   PetscReal                rho,fnorm,gnorm,gpnorm,xnorm=0,delta,nrm,ynorm,norm1;
 84:   PetscScalar              cnorm;
 85:   KSP                      ksp;
 86:   SNESConvergedReason      reason = SNES_CONVERGED_ITERATING;
 87:   PetscBool                breakout = PETSC_FALSE;
 88:   SNES_TR_KSPConverged_Ctx *ctx;
 89:   PetscErrorCode           (*convtest)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);

 92:   if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);

 94:   maxits = snes->max_its;               /* maximum number of iterations */
 95:   X      = snes->vec_sol;               /* solution vector */
 96:   F      = snes->vec_func;              /* residual vector */
 97:   Y      = snes->work[0];               /* work vectors */
 98:   G      = snes->work[1];
 99:   Ytmp   = snes->work[2];

101:   PetscObjectSAWsTakeAccess((PetscObject)snes);
102:   snes->iter = 0;
103:   PetscObjectSAWsGrantAccess((PetscObject)snes);

105:   /* Set the linear stopping criteria to use the More' trick. */
106:   SNESGetKSP(snes,&ksp);
107:   KSPGetConvergenceTest(ksp,&convtest,NULL,NULL);
108:   if (convtest != SNESTR_KSPConverged_Private) {
109:     PetscNew(&ctx);
110:     ctx->snes             = snes;
111:     KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy);
112:     KSPSetConvergenceTest(ksp,SNESTR_KSPConverged_Private,ctx,SNESTR_KSPConverged_Destroy);
113:     PetscInfo(snes,"Using Krylov convergence test SNESTR_KSPConverged_Private\n");
114:   }

116:   if (!snes->vec_func_init_set) {
117:     SNESComputeFunction(snes,X,F);          /* F(X) */
118:   } else snes->vec_func_init_set = PETSC_FALSE;

120:   VecNorm(F,NORM_2,&fnorm);             /* fnorm <- || F || */
121:   SNESCheckFunctionNorm(snes,fnorm);
122:   VecNorm(X,NORM_2,&xnorm);             /* fnorm <- || F || */
123:   PetscObjectSAWsTakeAccess((PetscObject)snes);
124:   snes->norm = fnorm;
125:   PetscObjectSAWsGrantAccess((PetscObject)snes);
126:   delta      = xnorm ? neP->delta0*xnorm : neP->delta0;
127:   neP->delta = delta;
128:   SNESLogConvergenceHistory(snes,fnorm,0);
129:   SNESMonitor(snes,0,fnorm);

131:   /* test convergence */
132:   (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
133:   if (snes->reason) return(0);


136:   for (i=0; i<maxits; i++) {

138:     /* Call general purpose update function */
139:     if (snes->ops->update) {
140:       (*snes->ops->update)(snes, snes->iter);
141:     }

143:     /* Solve J Y = F, where J is Jacobian matrix */
144:     SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
145:     SNESCheckJacobianDomainerror(snes);
146:     KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
147:     KSPSolve(snes->ksp,F,Ytmp);
148:     KSPGetIterationNumber(snes->ksp,&lits);

150:     snes->linear_its += lits;

152:     PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
153:     VecNorm(Ytmp,NORM_2,&nrm);
154:     norm1 = nrm;
155:     while (1) {
156:       VecCopy(Ytmp,Y);
157:       nrm  = norm1;

159:       /* Scale Y if need be and predict new value of F norm */
160:       if (nrm >= delta) {
161:         nrm    = delta/nrm;
162:         gpnorm = (1.0 - nrm)*fnorm;
163:         cnorm  = nrm;
164:         PetscInfo1(snes,"Scaling direction by %g\n",(double)nrm);
165:         VecScale(Y,cnorm);
166:         nrm    = gpnorm;
167:         ynorm  = delta;
168:       } else {
169:         gpnorm = 0.0;
170:         PetscInfo(snes,"Direction is in Trust Region\n");
171:         ynorm  = nrm;
172:       }
173:       VecCopy(Y,snes->vec_sol_update);
174:       VecAYPX(Y,-1.0,X);            /* Y <- X - Y */
175:       SNESComputeFunction(snes,Y,G); /*  F(X) */
176:       VecNorm(G,NORM_2,&gnorm);      /* gnorm <- || g || */
177:       SNESCheckFunctionNorm(snes,gnorm);
178:       if (fnorm == gpnorm) rho = 0.0;
179:       else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm);

181:       /* Update size of trust region */
182:       if      (rho < neP->mu)  delta *= neP->delta1;
183:       else if (rho < neP->eta) delta *= neP->delta2;
184:       else                     delta *= neP->delta3;
185:       PetscInfo3(snes,"fnorm=%g, gnorm=%g, ynorm=%g\n",(double)fnorm,(double)gnorm,(double)ynorm);
186:       PetscInfo3(snes,"gpred=%g, rho=%g, delta=%g\n",(double)gpnorm,(double)rho,(double)delta);

188:       neP->delta = delta;
189:       if (rho > neP->sigma) break;
190:       PetscInfo(snes,"Trying again in smaller region\n");
191:       /* check to see if progress is hopeless */
192:       neP->itflag = PETSC_FALSE;
193:       SNESTR_Converged_Private(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
194:       if (!reason) { (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP); }
195:       if (reason) {
196:         /* We're not progressing, so return with the current iterate */
197:         SNESMonitor(snes,i+1,fnorm);
198:         breakout = PETSC_TRUE;
199:         break;
200:       }
201:       snes->numFailures++;
202:     }
203:     if (!breakout) {
204:       /* Update function and solution vectors */
205:       fnorm = gnorm;
206:       VecCopy(G,F);
207:       VecCopy(Y,X);
208:       /* Monitor convergence */
209:       PetscObjectSAWsTakeAccess((PetscObject)snes);
210:       snes->iter = i+1;
211:       snes->norm = fnorm;
212:       snes->xnorm = xnorm;
213:       snes->ynorm = ynorm;
214:       PetscObjectSAWsGrantAccess((PetscObject)snes);
215:       SNESLogConvergenceHistory(snes,snes->norm,lits);
216:       SNESMonitor(snes,snes->iter,snes->norm);
217:       /* Test for convergence, xnorm = || X || */
218:       neP->itflag = PETSC_TRUE;
219:       if (snes->ops->converged != SNESConvergedSkip) { VecNorm(X,NORM_2,&xnorm); }
220:       (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
221:       if (reason) break;
222:     } else break;
223:   }
224:   if (i == maxits) {
225:     PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
226:     if (!reason) reason = SNES_DIVERGED_MAX_IT;
227:   }
228:   PetscObjectSAWsTakeAccess((PetscObject)snes);
229:   snes->reason = reason;
230:   PetscObjectSAWsGrantAccess((PetscObject)snes);
231:   return(0);
232: }
233: /*------------------------------------------------------------*/
234: static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes)
235: {

239:   SNESSetWorkVecs(snes,3);
240:   SNESSetUpMatrices(snes);
241:   return(0);
242: }

244: PetscErrorCode SNESReset_NEWTONTR(SNES snes)
245: {

248:   return(0);
249: }

251: static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes)
252: {

256:   SNESReset_NEWTONTR(snes);
257:   PetscFree(snes->data);
258:   return(0);
259: }
260: /*------------------------------------------------------------*/

262: static PetscErrorCode SNESSetFromOptions_NEWTONTR(PetscOptionItems *PetscOptionsObject,SNES snes)
263: {
264:   SNES_NEWTONTR  *ctx = (SNES_NEWTONTR*)snes->data;

268:   PetscOptionsHead(PetscOptionsObject,"SNES trust region options for nonlinear equations");
269:   PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,NULL);
270:   PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,NULL);
271:   PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,NULL);
272:   PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,NULL);
273:   PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,NULL);
274:   PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,NULL);
275:   PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,NULL);
276:   PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,NULL);
277:   PetscOptionsTail();
278:   return(0);
279: }

281: static PetscErrorCode SNESView_NEWTONTR(SNES snes,PetscViewer viewer)
282: {
283:   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;
285:   PetscBool      iascii;

288:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
289:   if (iascii) {
290:     PetscViewerASCIIPrintf(viewer,"  Trust region tolerance (-snes_trtol)\n",(double)snes->deltatol);
291:     PetscViewerASCIIPrintf(viewer,"  mu=%g, eta=%g, sigma=%g\n",(double)tr->mu,(double)tr->eta,(double)tr->sigma);
292:     PetscViewerASCIIPrintf(viewer,"  delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",(double)tr->delta0,(double)tr->delta1,(double)tr->delta2,(double)tr->delta3);
293:   }
294:   return(0);
295: }
296: /* ------------------------------------------------------------ */
297: /*MC
298:       SNESNEWTONTR - Newton based nonlinear solver that uses a trust region

300:    Options Database:
301: +    -snes_trtol <tol> - trust region tolerance
302: .    -snes_tr_mu <mu> - trust region parameter
303: .    -snes_tr_eta <eta> - trust region parameter
304: .    -snes_tr_sigma <sigma> - trust region parameter
305: .    -snes_tr_delta0 <delta0> -  initial size of the trust region is delta0*norm2(x)
306: .    -snes_tr_delta1 <delta1> - trust region parameter
307: .    -snes_tr_delta2 <delta2> - trust region parameter
308: -    -snes_tr_delta3 <delta3> - trust region parameter

310:    The basic algorithm is taken from "The Minpack Project", by More',
311:    Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development
312:    of Mathematical Software", Wayne Cowell, editor.

314:    Level: intermediate

316: .seealso:  SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESSetTrustRegionTolerance()

318: M*/
319: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes)
320: {
321:   SNES_NEWTONTR  *neP;

325:   snes->ops->setup          = SNESSetUp_NEWTONTR;
326:   snes->ops->solve          = SNESSolve_NEWTONTR;
327:   snes->ops->destroy        = SNESDestroy_NEWTONTR;
328:   snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR;
329:   snes->ops->view           = SNESView_NEWTONTR;
330:   snes->ops->reset          = SNESReset_NEWTONTR;

332:   snes->usesksp = PETSC_TRUE;
333:   snes->usesnpc = PETSC_FALSE;

335:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

337:   PetscNewLog(snes,&neP);
338:   snes->data  = (void*)neP;
339:   neP->mu     = 0.25;
340:   neP->eta    = 0.75;
341:   neP->delta  = 0.0;
342:   neP->delta0 = 0.2;
343:   neP->delta1 = 0.3;
344:   neP->delta2 = 0.75;
345:   neP->delta3 = 2.0;
346:   neP->sigma  = 0.0001;
347:   neP->itflag = PETSC_FALSE;
348:   neP->rnorm0 = 0.0;
349:   neP->ttol   = 0.0;
350:   return(0);
351: }