Actual source code: tr.c
petsc-3.14.6 2021-03-30
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,NULL,&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_DIVERGED_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: }
70: /*@C
71: SNESNewtonTRSetPreCheck - Sets a user function that is called before the search step has been determined.
72: Allows the user a chance to change or override the decision of the line search routine.
74: Logically Collective on snes
76: Input Parameters:
77: + snes - the nonlinear solver object
78: . func - [optional] function evaluation routine, see SNESNewtonTRPreCheck() for the calling sequence
79: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
81: Level: intermediate
83: Note: This function is called BEFORE the function evaluation within the SNESNEWTONTR solver.
85: .seealso: SNESNewtonTRPreCheck(), SNESNewtonTRGetPreCheck(), SNESNewtonTRSetPostCheck(), SNESNewtonTRGetPostCheck()
86: @*/
87: PetscErrorCode SNESNewtonTRSetPreCheck(SNES snes, PetscErrorCode (*func)(SNES,Vec,Vec,PetscBool*,void*),void *ctx)
88: {
89: SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data;
93: if (func) tr->precheck = func;
94: if (ctx) tr->precheckctx = ctx;
95: return(0);
96: }
98: /*@C
99: SNESNewtonTRGetPreCheck - Gets the pre-check function
101: Not collective
103: Input Parameter:
104: . snes - the nonlinear solver context
106: Output Parameters:
107: + func - [optional] function evaluation routine, see for the calling sequence SNESNewtonTRPreCheck()
108: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
110: Level: intermediate
112: .seealso: SNESNewtonTRSetPreCheck(), SNESNewtonTRPreCheck()
113: @*/
114: PetscErrorCode SNESNewtonTRGetPreCheck(SNES snes, PetscErrorCode (**func)(SNES,Vec,Vec,PetscBool*,void*),void **ctx)
115: {
116: SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data;
120: if (func) *func = tr->precheck;
121: if (ctx) *ctx = tr->precheckctx;
122: return(0);
123: }
125: /*@C
126: SNESNewtonTRSetPostCheck - Sets a user function that is called after the search step has been determined but before the next
127: function evaluation. Allows the user a chance to change or override the decision of the line search routine
129: Logically Collective on snes
131: Input Parameters:
132: + snes - the nonlinear solver object
133: . func - [optional] function evaluation routine, see SNESNewtonTRPostCheck() for the calling sequence
134: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
136: Level: intermediate
138: Note: This function is called BEFORE the function evaluation within the SNESNEWTONTR solver while the function set in
139: SNESLineSearchSetPostCheck() is called AFTER the function evaluation.
141: .seealso: SNESNewtonTRPostCheck(), SNESNewtonTRGetPostCheck()
142: @*/
143: PetscErrorCode SNESNewtonTRSetPostCheck(SNES snes, PetscErrorCode (*func)(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*),void *ctx)
144: {
145: SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data;
149: if (func) tr->postcheck = func;
150: if (ctx) tr->postcheckctx = ctx;
151: return(0);
152: }
154: /*@C
155: SNESNewtonTRGetPostCheck - Gets the post-check function
157: Not collective
159: Input Parameter:
160: . snes - the nonlinear solver context
162: Output Parameters:
163: + func - [optional] function evaluation routine, see for the calling sequence SNESNewtonTRPostCheck()
164: - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
166: Level: intermediate
168: .seealso: SNESNewtonTRSetPostCheck(), SNESNewtonTRPostCheck()
169: @*/
170: PetscErrorCode SNESNewtonTRGetPostCheck(SNES snes, PetscErrorCode (**func)(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*),void **ctx)
171: {
172: SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data;
176: if (func) *func = tr->postcheck;
177: if (ctx) *ctx = tr->postcheckctx;
178: return(0);
179: }
181: /*@C
182: SNESNewtonTRPreCheck - Called before the step has been determined in SNESNEWTONTR
184: Logically Collective on snes
186: Input Parameters:
187: + snes - the solver
188: . X - The last solution
189: - Y - The step direction
191: Output Parameters:
192: . changed_Y - Indicator that the step direction Y has been changed.
194: Level: developer
196: .seealso: SNESNewtonTRSetPreCheck(), SNESNewtonTRGetPreCheck()
197: @*/
198: static PetscErrorCode SNESNewtonTRPreCheck(SNES snes,Vec X,Vec Y,PetscBool *changed_Y)
199: {
200: SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data;
204: *changed_Y = PETSC_FALSE;
205: if (tr->precheck) {
206: (*tr->precheck)(snes,X,Y,changed_Y,tr->precheckctx);
208: }
209: return(0);
210: }
212: /*@C
213: SNESNewtonTRPostCheck - Called after the step has been determined in SNESNEWTONTR but before the function evaluation
215: Logically Collective on snes
217: Input Parameters:
218: + snes - the solver. X - The last solution
219: . Y - The full step direction
220: - W - The updated solution, W = X - Y
222: Output Parameters:
223: + changed_Y - indicator if step has been changed
224: - changed_W - Indicator if the new candidate solution W has been changed.
226: Notes:
227: If Y is changed then W is recomputed as X - Y
229: Level: developer
231: .seealso: SNESNewtonTRSetPostCheck(), SNESNewtonTRGetPostCheck()
232: @*/
233: static PetscErrorCode SNESNewtonTRPostCheck(SNES snes,Vec X,Vec Y,Vec W,PetscBool *changed_Y,PetscBool *changed_W)
234: {
235: SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data;
239: *changed_Y = PETSC_FALSE;
240: *changed_W = PETSC_FALSE;
241: if (tr->postcheck) {
242: (*tr->postcheck)(snes,X,Y,W,changed_Y,changed_W,tr->postcheckctx);
245: }
246: return(0);
247: }
249: /*
250: SNESSolve_NEWTONTR - Implements Newton's Method with a very simple trust
251: region approach for solving systems of nonlinear equations.
254: */
255: static PetscErrorCode SNESSolve_NEWTONTR(SNES snes)
256: {
257: SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data;
258: Vec X,F,Y,G,Ytmp,W;
259: PetscErrorCode ierr;
260: PetscInt maxits,i,lits;
261: PetscReal rho,fnorm,gnorm,gpnorm,xnorm=0,delta,nrm,ynorm,norm1;
262: PetscScalar cnorm;
263: KSP ksp;
264: SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
265: PetscBool breakout = PETSC_FALSE;
266: SNES_TR_KSPConverged_Ctx *ctx;
267: PetscErrorCode (*convtest)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),(*convdestroy)(void*);
268: void *convctx;
271: 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);
273: maxits = snes->max_its; /* maximum number of iterations */
274: X = snes->vec_sol; /* solution vector */
275: F = snes->vec_func; /* residual vector */
276: Y = snes->work[0]; /* work vectors */
277: G = snes->work[1];
278: Ytmp = snes->work[2];
279: W = snes->work[3];
281: PetscObjectSAWsTakeAccess((PetscObject)snes);
282: snes->iter = 0;
283: PetscObjectSAWsGrantAccess((PetscObject)snes);
285: /* Set the linear stopping criteria to use the More' trick. */
286: SNESGetKSP(snes,&ksp);
287: KSPGetConvergenceTest(ksp,&convtest,&convctx,&convdestroy);
288: if (convtest != SNESTR_KSPConverged_Private) {
289: PetscNew(&ctx);
290: ctx->snes = snes;
291: KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy);
292: KSPSetConvergenceTest(ksp,SNESTR_KSPConverged_Private,ctx,SNESTR_KSPConverged_Destroy);
293: PetscInfo(snes,"Using Krylov convergence test SNESTR_KSPConverged_Private\n");
294: }
296: if (!snes->vec_func_init_set) {
297: SNESComputeFunction(snes,X,F); /* F(X) */
298: } else snes->vec_func_init_set = PETSC_FALSE;
300: VecNorm(F,NORM_2,&fnorm); /* fnorm <- || F || */
301: SNESCheckFunctionNorm(snes,fnorm);
302: VecNorm(X,NORM_2,&xnorm); /* fnorm <- || F || */
303: PetscObjectSAWsTakeAccess((PetscObject)snes);
304: snes->norm = fnorm;
305: PetscObjectSAWsGrantAccess((PetscObject)snes);
306: delta = xnorm ? neP->delta0*xnorm : neP->delta0;
307: neP->delta = delta;
308: SNESLogConvergenceHistory(snes,fnorm,0);
309: SNESMonitor(snes,0,fnorm);
311: /* test convergence */
312: (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
313: if (snes->reason) return(0);
316: for (i=0; i<maxits; i++) {
318: /* Call general purpose update function */
319: if (snes->ops->update) {
320: (*snes->ops->update)(snes, snes->iter);
321: }
323: /* Solve J Y = F, where J is Jacobian matrix */
324: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
325: SNESCheckJacobianDomainerror(snes);
326: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
327: KSPSolve(snes->ksp,F,Ytmp);
328: KSPGetIterationNumber(snes->ksp,&lits);
330: snes->linear_its += lits;
332: PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
333: VecNorm(Ytmp,NORM_2,&nrm);
334: norm1 = nrm;
335: while (1) {
336: PetscBool changed_y;
337: PetscBool changed_w;
338: VecCopy(Ytmp,Y);
339: nrm = norm1;
341: /* Scale Y if need be and predict new value of F norm */
342: if (nrm >= delta) {
343: nrm = delta/nrm;
344: gpnorm = (1.0 - nrm)*fnorm;
345: cnorm = nrm;
346: PetscInfo1(snes,"Scaling direction by %g\n",(double)nrm);
347: VecScale(Y,cnorm);
348: nrm = gpnorm;
349: ynorm = delta;
350: } else {
351: gpnorm = 0.0;
352: PetscInfo(snes,"Direction is in Trust Region\n");
353: ynorm = nrm;
354: }
355: /* PreCheck() allows for updates to Y prior to W <- X - Y */
356: SNESNewtonTRPreCheck(snes,X,Y,&changed_y);
357: VecWAXPY(W,-1.0,Y,X); /* W <- X - Y */
358: SNESNewtonTRPostCheck(snes,X,Y,W,&changed_y,&changed_w);
359: if (changed_y) VecWAXPY(W,-1.0,Y,X);
360: VecCopy(Y,snes->vec_sol_update);
361: SNESComputeFunction(snes,W,G); /* F(X) */
362: VecNorm(G,NORM_2,&gnorm); /* gnorm <- || g || */
363: SNESCheckFunctionNorm(snes,gnorm);
364: if (fnorm == gpnorm) rho = 0.0;
365: else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm);
367: /* Update size of trust region */
368: if (rho < neP->mu) delta *= neP->delta1;
369: else if (rho < neP->eta) delta *= neP->delta2;
370: else delta *= neP->delta3;
371: PetscInfo3(snes,"fnorm=%g, gnorm=%g, ynorm=%g\n",(double)fnorm,(double)gnorm,(double)ynorm);
372: PetscInfo3(snes,"gpred=%g, rho=%g, delta=%g\n",(double)gpnorm,(double)rho,(double)delta);
374: neP->delta = delta;
375: if (rho > neP->sigma) break;
376: PetscInfo(snes,"Trying again in smaller region\n");
377: /* check to see if progress is hopeless */
378: neP->itflag = PETSC_FALSE;
379: SNESTR_Converged_Private(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
380: if (!reason) {(*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);}
381: if (reason == SNES_CONVERGED_SNORM_RELATIVE) reason = SNES_DIVERGED_INNER;
382: if (reason) {
383: /* We're not progressing, so return with the current iterate */
384: SNESMonitor(snes,i+1,fnorm);
385: breakout = PETSC_TRUE;
386: break;
387: }
388: snes->numFailures++;
389: }
390: if (!breakout) {
391: /* Update function and solution vectors */
392: fnorm = gnorm;
393: VecCopy(G,F);
394: VecCopy(W,X);
395: /* Monitor convergence */
396: PetscObjectSAWsTakeAccess((PetscObject)snes);
397: snes->iter = i+1;
398: snes->norm = fnorm;
399: snes->xnorm = xnorm;
400: snes->ynorm = ynorm;
401: PetscObjectSAWsGrantAccess((PetscObject)snes);
402: SNESLogConvergenceHistory(snes,snes->norm,lits);
403: SNESMonitor(snes,snes->iter,snes->norm);
404: /* Test for convergence, xnorm = || X || */
405: neP->itflag = PETSC_TRUE;
406: if (snes->ops->converged != SNESConvergedSkip) { VecNorm(X,NORM_2,&xnorm); }
407: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
408: if (reason) break;
409: } else break;
410: }
411: if (i == maxits) {
412: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
413: if (!reason) reason = SNES_DIVERGED_MAX_IT;
414: }
415: PetscObjectSAWsTakeAccess((PetscObject)snes);
416: snes->reason = reason;
417: PetscObjectSAWsGrantAccess((PetscObject)snes);
418: if (convtest != SNESTR_KSPConverged_Private) {
419: KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy);
420: PetscFree(ctx);
421: KSPSetConvergenceTest(ksp,convtest,convctx,convdestroy);
422: }
423: return(0);
424: }
425: /*------------------------------------------------------------*/
426: static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes)
427: {
431: SNESSetWorkVecs(snes,4);
432: SNESSetUpMatrices(snes);
433: return(0);
434: }
436: PetscErrorCode SNESReset_NEWTONTR(SNES snes)
437: {
440: return(0);
441: }
443: static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes)
444: {
448: SNESReset_NEWTONTR(snes);
449: PetscFree(snes->data);
450: return(0);
451: }
452: /*------------------------------------------------------------*/
454: static PetscErrorCode SNESSetFromOptions_NEWTONTR(PetscOptionItems *PetscOptionsObject,SNES snes)
455: {
456: SNES_NEWTONTR *ctx = (SNES_NEWTONTR*)snes->data;
460: PetscOptionsHead(PetscOptionsObject,"SNES trust region options for nonlinear equations");
461: PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,NULL);
462: PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,NULL);
463: PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,NULL);
464: PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,NULL);
465: PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,NULL);
466: PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,NULL);
467: PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,NULL);
468: PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,NULL);
469: PetscOptionsTail();
470: return(0);
471: }
473: static PetscErrorCode SNESView_NEWTONTR(SNES snes,PetscViewer viewer)
474: {
475: SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data;
477: PetscBool iascii;
480: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
481: if (iascii) {
482: PetscViewerASCIIPrintf(viewer," Trust region tolerance (-snes_trtol)\n",(double)snes->deltatol);
483: PetscViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",(double)tr->mu,(double)tr->eta,(double)tr->sigma);
484: PetscViewerASCIIPrintf(viewer," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",(double)tr->delta0,(double)tr->delta1,(double)tr->delta2,(double)tr->delta3);
485: }
486: return(0);
487: }
488: /* ------------------------------------------------------------ */
489: /*MC
490: SNESNEWTONTR - Newton based nonlinear solver that uses a trust region
492: Options Database:
493: + -snes_trtol <tol> - trust region tolerance
494: . -snes_tr_mu <mu> - trust region parameter
495: . -snes_tr_eta <eta> - trust region parameter
496: . -snes_tr_sigma <sigma> - trust region parameter
497: . -snes_tr_delta0 <delta0> - initial size of the trust region is delta0*norm2(x)
498: . -snes_tr_delta1 <delta1> - trust region parameter
499: . -snes_tr_delta2 <delta2> - trust region parameter
500: - -snes_tr_delta3 <delta3> - trust region parameter
502: The basic algorithm is taken from "The Minpack Project", by More',
503: Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development
504: of Mathematical Software", Wayne Cowell, editor.
506: Level: intermediate
508: .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESSetTrustRegionTolerance()
510: M*/
511: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes)
512: {
513: SNES_NEWTONTR *neP;
517: snes->ops->setup = SNESSetUp_NEWTONTR;
518: snes->ops->solve = SNESSolve_NEWTONTR;
519: snes->ops->destroy = SNESDestroy_NEWTONTR;
520: snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR;
521: snes->ops->view = SNESView_NEWTONTR;
522: snes->ops->reset = SNESReset_NEWTONTR;
524: snes->usesksp = PETSC_TRUE;
525: snes->usesnpc = PETSC_FALSE;
527: snes->alwayscomputesfinalresidual = PETSC_TRUE;
529: PetscNewLog(snes,&neP);
530: snes->data = (void*)neP;
531: neP->mu = 0.25;
532: neP->eta = 0.75;
533: neP->delta = 0.0;
534: neP->delta0 = 0.2;
535: neP->delta1 = 0.3;
536: neP->delta2 = 0.75;
537: neP->delta3 = 2.0;
538: neP->sigma = 0.0001;
539: neP->itflag = PETSC_FALSE;
540: neP->rnorm0 = 0.0;
541: neP->ttol = 0.0;
542: return(0);
543: }