Actual source code: tr.c

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

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

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

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

 33: static PetscErrorCode SNESTR_KSPConverged_Destroy(void *cctx)
 34: {
 35:   SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx*)cctx;

 37:   (*ctx->convdestroy)(ctx->convctx);
 38:   PetscFree(ctx);
 39:   return 0;
 40: }

 42: /* ---------------------------------------------------------------- */
 43: /*
 44:    SNESTR_Converged_Private -test convergence JUST for
 45:    the trust region tolerance.

 47: */
 48: static PetscErrorCode SNESTR_Converged_Private(SNES snes,PetscInt it,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
 49: {
 50:   SNES_NEWTONTR  *neP = (SNES_NEWTONTR*)snes->data;

 52:   *reason = SNES_CONVERGED_ITERATING;
 53:   if (neP->delta < xnorm * snes->deltatol) {
 54:     PetscInfo(snes,"Converged due to trust region param %g<%g*%g\n",(double)neP->delta,(double)xnorm,(double)snes->deltatol);
 55:     *reason = SNES_DIVERGED_TR_DELTA;
 56:   } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
 57:     PetscInfo(snes,"Exceeded maximum number of function evaluations: %D\n",snes->max_funcs);
 58:     *reason = SNES_DIVERGED_FUNCTION_COUNT;
 59:   }
 60:   return 0;
 61: }

 63: /*@C
 64:    SNESNewtonTRSetPreCheck - Sets a user function that is called before the search step has been determined.
 65:        Allows the user a chance to change or override the decision of the line search routine.

 67:    Logically Collective on snes

 69:    Input Parameters:
 70: +  snes - the nonlinear solver object
 71: .  func - [optional] function evaluation routine, see SNESNewtonTRPreCheck()  for the calling sequence
 72: -  ctx  - [optional] user-defined context for private data for the function evaluation routine (may be NULL)

 74:    Level: intermediate

 76:    Note: This function is called BEFORE the function evaluation within the SNESNEWTONTR solver.

 78: .seealso: SNESNewtonTRPreCheck(), SNESNewtonTRGetPreCheck(), SNESNewtonTRSetPostCheck(), SNESNewtonTRGetPostCheck()
 79: @*/
 80: PetscErrorCode  SNESNewtonTRSetPreCheck(SNES snes, PetscErrorCode (*func)(SNES,Vec,Vec,PetscBool*,void*),void *ctx)
 81: {
 82:   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;

 85:   if (func) tr->precheck    = func;
 86:   if (ctx)  tr->precheckctx = ctx;
 87:   return 0;
 88: }

 90: /*@C
 91:    SNESNewtonTRGetPreCheck - Gets the pre-check function

 93:    Not collective

 95:    Input Parameter:
 96: .  snes - the nonlinear solver context

 98:    Output Parameters:
 99: +  func - [optional] function evaluation routine, see for the calling sequence SNESNewtonTRPreCheck()
100: -  ctx  - [optional] user-defined context for private data for the function evaluation routine (may be NULL)

102:    Level: intermediate

104: .seealso: SNESNewtonTRSetPreCheck(), SNESNewtonTRPreCheck()
105: @*/
106: PetscErrorCode  SNESNewtonTRGetPreCheck(SNES snes, PetscErrorCode (**func)(SNES,Vec,Vec,PetscBool*,void*),void **ctx)
107: {
108:   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;

111:   if (func) *func = tr->precheck;
112:   if (ctx)  *ctx  = tr->precheckctx;
113:   return 0;
114: }

116: /*@C
117:    SNESNewtonTRSetPostCheck - Sets a user function that is called after the search step has been determined but before the next
118:        function evaluation. Allows the user a chance to change or override the decision of the line search routine

120:    Logically Collective on snes

122:    Input Parameters:
123: +  snes - the nonlinear solver object
124: .  func - [optional] function evaluation routine, see SNESNewtonTRPostCheck()  for the calling sequence
125: -  ctx  - [optional] user-defined context for private data for the function evaluation routine (may be NULL)

127:    Level: intermediate

129:    Note: This function is called BEFORE the function evaluation within the SNESNEWTONTR solver while the function set in
130:    SNESLineSearchSetPostCheck() is called AFTER the function evaluation.

132: .seealso: SNESNewtonTRPostCheck(), SNESNewtonTRGetPostCheck()
133: @*/
134: PetscErrorCode  SNESNewtonTRSetPostCheck(SNES snes,PetscErrorCode (*func)(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*),void *ctx)
135: {
136:   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;

139:   if (func) tr->postcheck    = func;
140:   if (ctx)  tr->postcheckctx = ctx;
141:   return 0;
142: }

144: /*@C
145:    SNESNewtonTRGetPostCheck - Gets the post-check function

147:    Not collective

149:    Input Parameter:
150: .  snes - the nonlinear solver context

152:    Output Parameters:
153: +  func - [optional] function evaluation routine, see for the calling sequence SNESNewtonTRPostCheck()
154: -  ctx  - [optional] user-defined context for private data for the function evaluation routine (may be NULL)

156:    Level: intermediate

158: .seealso: SNESNewtonTRSetPostCheck(), SNESNewtonTRPostCheck()
159: @*/
160: PetscErrorCode  SNESNewtonTRGetPostCheck(SNES snes,PetscErrorCode (**func)(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*),void **ctx)
161: {
162:   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;

165:   if (func) *func = tr->postcheck;
166:   if (ctx)  *ctx  = tr->postcheckctx;
167:   return 0;
168: }

170: /*@C
171:    SNESNewtonTRPreCheck - Called before the step has been determined in SNESNEWTONTR

173:    Logically Collective on snes

175:    Input Parameters:
176: +  snes - the solver
177: .  X - The last solution
178: -  Y - The step direction

180:    Output Parameters:
181: .  changed_Y - Indicator that the step direction Y has been changed.

183:    Level: developer

185: .seealso: SNESNewtonTRSetPreCheck(), SNESNewtonTRGetPreCheck()
186: @*/
187: static PetscErrorCode SNESNewtonTRPreCheck(SNES snes,Vec X,Vec Y,PetscBool *changed_Y)
188: {
189:   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;

191:   *changed_Y = PETSC_FALSE;
192:   if (tr->precheck) {
193:     (*tr->precheck)(snes,X,Y,changed_Y,tr->precheckctx);
195:   }
196:   return 0;
197: }

199: /*@C
200:    SNESNewtonTRPostCheck - Called after the step has been determined in SNESNEWTONTR but before the function evaluation

202:    Logically Collective on snes

204:    Input Parameters:
205: +  snes - the solver
206: .  X - The last solution
207: .  Y - The full step direction
208: -  W - The updated solution, W = X - Y

210:    Output Parameters:
211: +  changed_Y - indicator if step has been changed
212: -  changed_W - Indicator if the new candidate solution W has been changed.

214:    Notes:
215:      If Y is changed then W is recomputed as X - Y

217:    Level: developer

219: .seealso: SNESNewtonTRSetPostCheck(), SNESNewtonTRGetPostCheck()
220: @*/
221: static PetscErrorCode SNESNewtonTRPostCheck(SNES snes,Vec X,Vec Y,Vec W,PetscBool *changed_Y,PetscBool *changed_W)
222: {
223:   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;

225:   *changed_Y = PETSC_FALSE;
226:   *changed_W = PETSC_FALSE;
227:   if (tr->postcheck) {
228:     (*tr->postcheck)(snes,X,Y,W,changed_Y,changed_W,tr->postcheckctx);
231:   }
232:   return 0;
233: }

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

239: */
240: static PetscErrorCode SNESSolve_NEWTONTR(SNES snes)
241: {
242:   SNES_NEWTONTR            *neP = (SNES_NEWTONTR*)snes->data;
243:   Vec                      X,F,Y,G,Ytmp,W;
244:   PetscInt                 maxits,i,lits;
245:   PetscReal                rho,fnorm,gnorm,gpnorm,xnorm=0,delta,nrm,ynorm,norm1;
246:   PetscScalar              cnorm;
247:   KSP                      ksp;
248:   SNESConvergedReason      reason = SNES_CONVERGED_ITERATING;
249:   PetscBool                breakout = PETSC_FALSE;
250:   SNES_TR_KSPConverged_Ctx *ctx;
251:   PetscErrorCode           (*convtest)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),(*convdestroy)(void*);
252:   void                     *convctx;


256:   maxits = snes->max_its;               /* maximum number of iterations */
257:   X      = snes->vec_sol;               /* solution vector */
258:   F      = snes->vec_func;              /* residual vector */
259:   Y      = snes->work[0];               /* work vectors */
260:   G      = snes->work[1];
261:   Ytmp   = snes->work[2];
262:   W      = snes->work[3];

264:   PetscObjectSAWsTakeAccess((PetscObject)snes);
265:   snes->iter = 0;
266:   PetscObjectSAWsGrantAccess((PetscObject)snes);

268:   /* Set the linear stopping criteria to use the More' trick. */
269:   SNESGetKSP(snes,&ksp);
270:   KSPGetConvergenceTest(ksp,&convtest,&convctx,&convdestroy);
271:   if (convtest != SNESTR_KSPConverged_Private) {
272:     PetscNew(&ctx);
273:     ctx->snes             = snes;
274:     KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy);
275:     KSPSetConvergenceTest(ksp,SNESTR_KSPConverged_Private,ctx,SNESTR_KSPConverged_Destroy);
276:     PetscInfo(snes,"Using Krylov convergence test SNESTR_KSPConverged_Private\n");
277:   }

279:   if (!snes->vec_func_init_set) {
280:     SNESComputeFunction(snes,X,F);          /* F(X) */
281:   } else snes->vec_func_init_set = PETSC_FALSE;

283:   VecNorm(F,NORM_2,&fnorm);             /* fnorm <- || F || */
284:   SNESCheckFunctionNorm(snes,fnorm);
285:   VecNorm(X,NORM_2,&xnorm);             /* xnorm <- || X || */
286:   PetscObjectSAWsTakeAccess((PetscObject)snes);
287:   snes->norm = fnorm;
288:   PetscObjectSAWsGrantAccess((PetscObject)snes);
289:   delta      = xnorm ? neP->delta0*xnorm : neP->delta0;
290:   neP->delta = delta;
291:   SNESLogConvergenceHistory(snes,fnorm,0);
292:   SNESMonitor(snes,0,fnorm);

294:   /* test convergence */
295:   (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
296:   if (snes->reason) return 0;

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

300:     /* Call general purpose update function */
301:     if (snes->ops->update) {
302:       (*snes->ops->update)(snes, snes->iter);
303:     }

305:     /* Solve J Y = F, where J is Jacobian matrix */
306:     SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
307:     SNESCheckJacobianDomainerror(snes);
308:     KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
309:     KSPSolve(snes->ksp,F,Ytmp);
310:     KSPGetIterationNumber(snes->ksp,&lits);
311:     snes->linear_its += lits;

313:     PetscInfo(snes,"iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n",snes->iter,lits);
314:     VecNorm(Ytmp,NORM_2,&nrm);
315:     norm1 = nrm;

317:     while (1) {
318:       PetscBool changed_y;
319:       PetscBool changed_w;
320:       VecCopy(Ytmp,Y);
321:       nrm  = norm1;

323:       /* Scale Y if need be and predict new value of F norm */
324:       if (nrm >= delta) {
325:         nrm    = delta/nrm;
326:         gpnorm = (1.0 - nrm)*fnorm;
327:         cnorm  = nrm;
328:         PetscInfo(snes,"Scaling direction by %g\n",(double)nrm);
329:         VecScale(Y,cnorm);
330:         nrm    = gpnorm;
331:         ynorm  = delta;
332:       } else {
333:         gpnorm = 0.0;
334:         PetscInfo(snes,"Direction is in Trust Region\n");
335:         ynorm  = nrm;
336:       }
337:       /* PreCheck() allows for updates to Y prior to W <- X - Y */
338:       SNESNewtonTRPreCheck(snes,X,Y,&changed_y);
339:       VecWAXPY(W,-1.0,Y,X);         /* W <- X - Y */
340:       SNESNewtonTRPostCheck(snes,X,Y,W,&changed_y,&changed_w);
341:       if (changed_y) VecWAXPY(W,-1.0,Y,X);
342:       VecCopy(Y,snes->vec_sol_update);
343:       SNESComputeFunction(snes,W,G); /*  F(X-Y) = G */
344:       VecNorm(G,NORM_2,&gnorm);      /* gnorm <- || g || */
345:       SNESCheckFunctionNorm(snes,gnorm);
346:       if (fnorm == gpnorm) rho = 0.0;
347:       else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm);

349:       /* Update size of trust region */
350:       if      (rho < neP->mu)  delta *= neP->delta1;
351:       else if (rho < neP->eta) delta *= neP->delta2;
352:       else                     delta *= neP->delta3;
353:       PetscInfo(snes,"fnorm=%g, gnorm=%g, ynorm=%g\n",(double)fnorm,(double)gnorm,(double)ynorm);
354:       PetscInfo(snes,"gpred=%g, rho=%g, delta=%g\n",(double)gpnorm,(double)rho,(double)delta);

356:       neP->delta = delta;
357:       if (rho > neP->sigma) break;
358:       PetscInfo(snes,"Trying again in smaller region\n");

360:       /* check to see if progress is hopeless */
361:       neP->itflag = PETSC_FALSE;
362:       SNESTR_Converged_Private(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
363:       if (!reason) (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
364:       if (reason == SNES_CONVERGED_SNORM_RELATIVE) reason = SNES_DIVERGED_INNER;
365:       if (reason) {
366:         /* We're not progressing, so return with the current iterate */
367:         SNESMonitor(snes,i+1,fnorm);
368:         breakout = PETSC_TRUE;
369:         break;
370:       }
371:       snes->numFailures++;
372:     }
373:     if (!breakout) {
374:       /* Update function and solution vectors */
375:       fnorm = gnorm;
376:       VecCopy(G,F);
377:       VecCopy(W,X);
378:       /* Monitor convergence */
379:       PetscObjectSAWsTakeAccess((PetscObject)snes);
380:       snes->iter = i+1;
381:       snes->norm = fnorm;
382:       snes->xnorm = xnorm;
383:       snes->ynorm = ynorm;
384:       PetscObjectSAWsGrantAccess((PetscObject)snes);
385:       SNESLogConvergenceHistory(snes,snes->norm,lits);
386:       SNESMonitor(snes,snes->iter,snes->norm);
387:       /* Test for convergence, xnorm = || X || */
388:       neP->itflag = PETSC_TRUE;
389:       if (snes->ops->converged != SNESConvergedSkip) VecNorm(X,NORM_2,&xnorm);
390:       (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
391:       if (reason) break;
392:     } else break;
393:   }

395:   if (i == maxits) {
396:     PetscInfo(snes,"Maximum number of iterations has been reached: %" PetscInt_FMT "\n",maxits);
397:     if (!reason) reason = SNES_DIVERGED_MAX_IT;
398:   }
399:   PetscObjectSAWsTakeAccess((PetscObject)snes);
400:   snes->reason = reason;
401:   PetscObjectSAWsGrantAccess((PetscObject)snes);
402:   if (convtest != SNESTR_KSPConverged_Private) {
403:     KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy);
404:     PetscFree(ctx);
405:     KSPSetConvergenceTest(ksp,convtest,convctx,convdestroy);
406:   }
407:   return 0;
408: }

410: /*------------------------------------------------------------*/
411: static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes)
412: {
413:   SNESSetWorkVecs(snes,4);
414:   SNESSetUpMatrices(snes);
415:   return 0;
416: }

418: PetscErrorCode SNESReset_NEWTONTR(SNES snes)
419: {
420:   return 0;
421: }

423: static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes)
424: {
425:   SNESReset_NEWTONTR(snes);
426:   PetscFree(snes->data);
427:   return 0;
428: }
429: /*------------------------------------------------------------*/

431: static PetscErrorCode SNESSetFromOptions_NEWTONTR(PetscOptionItems *PetscOptionsObject,SNES snes)
432: {
433:   SNES_NEWTONTR  *ctx = (SNES_NEWTONTR*)snes->data;

435:   PetscOptionsHead(PetscOptionsObject,"SNES trust region options for nonlinear equations");
436:   PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,NULL);
437:   PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,NULL);
438:   PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,NULL);
439:   PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,NULL);
440:   PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,NULL);
441:   PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,NULL);
442:   PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,NULL);
443:   PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,NULL);
444:   PetscOptionsTail();
445:   return 0;
446: }

448: static PetscErrorCode SNESView_NEWTONTR(SNES snes,PetscViewer viewer)
449: {
450:   SNES_NEWTONTR  *tr = (SNES_NEWTONTR*)snes->data;
451:   PetscBool      iascii;

453:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
454:   if (iascii) {
455:     PetscViewerASCIIPrintf(viewer,"  Trust region tolerance (-snes_trtol)\n",(double)snes->deltatol);
456:     PetscViewerASCIIPrintf(viewer,"  mu=%g, eta=%g, sigma=%g\n",(double)tr->mu,(double)tr->eta,(double)tr->sigma);
457:     PetscViewerASCIIPrintf(viewer,"  delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",(double)tr->delta0,(double)tr->delta1,(double)tr->delta2,(double)tr->delta3);
458:   }
459:   return 0;
460: }
461: /* ------------------------------------------------------------ */
462: /*MC
463:       SNESNEWTONTR - Newton based nonlinear solver that uses a trust region

465:    Options Database:
466: +    -snes_trtol <tol> - trust region tolerance
467: .    -snes_tr_mu <mu> - trust region parameter
468: .    -snes_tr_eta <eta> - trust region parameter
469: .    -snes_tr_sigma <sigma> - trust region parameter
470: .    -snes_tr_delta0 <delta0> -  initial size of the trust region is delta0*norm2(x)
471: .    -snes_tr_delta1 <delta1> - trust region parameter
472: .    -snes_tr_delta2 <delta2> - trust region parameter
473: -    -snes_tr_delta3 <delta3> - trust region parameter

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

479:    Level: intermediate

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

483: M*/
484: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes)
485: {
486:   SNES_NEWTONTR  *neP;

488:   snes->ops->setup          = SNESSetUp_NEWTONTR;
489:   snes->ops->solve          = SNESSolve_NEWTONTR;
490:   snes->ops->destroy        = SNESDestroy_NEWTONTR;
491:   snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR;
492:   snes->ops->view           = SNESView_NEWTONTR;
493:   snes->ops->reset          = SNESReset_NEWTONTR;

495:   snes->usesksp = PETSC_TRUE;
496:   snes->usesnpc = PETSC_FALSE;

498:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

500:   PetscNewLog(snes,&neP);
501:   snes->data  = (void*)neP;
502:   neP->mu     = 0.25;
503:   neP->eta    = 0.75;
504:   neP->delta  = 0.0;
505:   neP->delta0 = 0.2;
506:   neP->delta1 = 0.3;
507:   neP->delta2 = 0.75;
508:   neP->delta3 = 2.0;
509:   neP->sigma  = 0.0001;
510:   neP->itflag = PETSC_FALSE;
511:   neP->rnorm0 = 0.0;
512:   neP->ttol   = 0.0;
513:   return 0;
514: }