xref: /petsc/src/snes/impls/tr/tr.c (revision 4dfa11a44d5adf2389f1d3acbc8f3c1116dc6c3a)
1c6db04a5SJed Brown #include <../src/snes/impls/tr/trimpl.h> /*I   "petscsnes.h"   I*/
24800dd8cSBarry Smith 
3971273eeSBarry Smith typedef struct {
4971273eeSBarry Smith   SNES snes;
5df8705c3SBarry Smith   /*  Information on the regular SNES convergence test; which may have been user provided */
6df8705c3SBarry Smith   PetscErrorCode (*convtest)(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *);
7df8705c3SBarry Smith   PetscErrorCode (*convdestroy)(void *);
8df8705c3SBarry Smith   void *convctx;
9971273eeSBarry Smith } SNES_TR_KSPConverged_Ctx;
10971273eeSBarry Smith 
119371c9d4SSatish Balay static PetscErrorCode SNESTR_KSPConverged_Private(KSP ksp, PetscInt n, PetscReal rnorm, KSPConvergedReason *reason, void *cctx) {
12971273eeSBarry Smith   SNES_TR_KSPConverged_Ctx *ctx  = (SNES_TR_KSPConverged_Ctx *)cctx;
13971273eeSBarry Smith   SNES                      snes = ctx->snes;
1404d7464bSBarry Smith   SNES_NEWTONTR            *neP  = (SNES_NEWTONTR *)snes->data;
15df60cc22SBarry Smith   Vec                       x;
16064f8208SBarry Smith   PetscReal                 nrm;
17df60cc22SBarry Smith 
183a40ed3dSBarry Smith   PetscFunctionBegin;
199566063dSJacob Faibussowitsch   PetscCall((*ctx->convtest)(ksp, n, rnorm, reason, ctx->convctx));
2048a46eb9SPierre Jolivet   if (*reason) PetscCall(PetscInfo(snes, "Default or user provided convergence test KSP iterations=%" PetscInt_FMT ", rnorm=%g\n", n, (double)rnorm));
21a935fc98SLois Curfman McInnes   /* Determine norm of solution */
229566063dSJacob Faibussowitsch   PetscCall(KSPBuildSolution(ksp, NULL, &x));
239566063dSJacob Faibussowitsch   PetscCall(VecNorm(x, NORM_2, &nrm));
24064f8208SBarry Smith   if (nrm >= neP->delta) {
259566063dSJacob Faibussowitsch     PetscCall(PetscInfo(snes, "Ending linear iteration early, delta=%g, length=%g\n", (double)neP->delta, (double)nrm));
26329f5518SBarry Smith     *reason = KSP_CONVERGED_STEP_LENGTH;
27df60cc22SBarry Smith   }
283a40ed3dSBarry Smith   PetscFunctionReturn(0);
29df60cc22SBarry Smith }
3082bf6240SBarry Smith 
319371c9d4SSatish Balay static PetscErrorCode SNESTR_KSPConverged_Destroy(void *cctx) {
32971273eeSBarry Smith   SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx *)cctx;
33971273eeSBarry Smith 
34971273eeSBarry Smith   PetscFunctionBegin;
359566063dSJacob Faibussowitsch   PetscCall((*ctx->convdestroy)(ctx->convctx));
369566063dSJacob Faibussowitsch   PetscCall(PetscFree(ctx));
37971273eeSBarry Smith   PetscFunctionReturn(0);
38971273eeSBarry Smith }
39971273eeSBarry Smith 
4085385478SLisandro Dalcin /*
41470880abSPatrick Sanan    SNESTR_Converged_Private -test convergence JUST for
4285385478SLisandro Dalcin    the trust region tolerance.
4385385478SLisandro Dalcin */
449371c9d4SSatish Balay static PetscErrorCode SNESTR_Converged_Private(SNES snes, PetscInt it, PetscReal xnorm, PetscReal pnorm, PetscReal fnorm, SNESConvergedReason *reason, void *dummy) {
4504d7464bSBarry Smith   SNES_NEWTONTR *neP = (SNES_NEWTONTR *)snes->data;
4685385478SLisandro Dalcin 
4785385478SLisandro Dalcin   PetscFunctionBegin;
4885385478SLisandro Dalcin   *reason = SNES_CONVERGED_ITERATING;
4985385478SLisandro Dalcin   if (neP->delta < xnorm * snes->deltatol) {
509566063dSJacob Faibussowitsch     PetscCall(PetscInfo(snes, "Converged due to trust region param %g<%g*%g\n", (double)neP->delta, (double)xnorm, (double)snes->deltatol));
511c6b2ff8SBarry Smith     *reason = SNES_DIVERGED_TR_DELTA;
52e71169deSBarry Smith   } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
5363a3b9bcSJacob Faibussowitsch     PetscCall(PetscInfo(snes, "Exceeded maximum number of function evaluations: %" PetscInt_FMT "\n", snes->max_funcs));
5485385478SLisandro Dalcin     *reason = SNES_DIVERGED_FUNCTION_COUNT;
5585385478SLisandro Dalcin   }
5685385478SLisandro Dalcin   PetscFunctionReturn(0);
5785385478SLisandro Dalcin }
5885385478SLisandro Dalcin 
597cb011f5SBarry Smith /*@C
60c9368356SGlenn Hammond    SNESNewtonTRSetPreCheck - Sets a user function that is called before the search step has been determined.
61c9368356SGlenn Hammond        Allows the user a chance to change or override the decision of the line search routine.
62c9368356SGlenn Hammond 
63f6dfbefdSBarry Smith    Deprecated use `SNESNEWTONDCTRDC`
64f6dfbefdSBarry Smith 
65c9368356SGlenn Hammond    Logically Collective on snes
66c9368356SGlenn Hammond 
67c9368356SGlenn Hammond    Input Parameters:
68c9368356SGlenn Hammond +  snes - the nonlinear solver object
69f6dfbefdSBarry Smith .  func - [optional] function evaluation routine, see `SNESNewtonTRPreCheck()`  for the calling sequence
70c9368356SGlenn Hammond -  ctx  - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
71c9368356SGlenn Hammond 
72c9368356SGlenn Hammond    Level: intermediate
73c9368356SGlenn Hammond 
74f6dfbefdSBarry Smith    Note:
75f6dfbefdSBarry Smith    This function is called BEFORE the function evaluation within the `SNESNEWTONTR` solver.
76c9368356SGlenn Hammond 
77f6dfbefdSBarry Smith .seealso: `SNESNEWTONDCTRDC`, `SNESNEWTONDCTR`, `SNESNewtonTRPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`
78c9368356SGlenn Hammond @*/
799371c9d4SSatish Balay PetscErrorCode SNESNewtonTRSetPreCheck(SNES snes, PetscErrorCode (*func)(SNES, Vec, Vec, PetscBool *, void *), void *ctx) {
80c9368356SGlenn Hammond   SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
81c9368356SGlenn Hammond 
82c9368356SGlenn Hammond   PetscFunctionBegin;
83c9368356SGlenn Hammond   PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
84c9368356SGlenn Hammond   if (func) tr->precheck = func;
85c9368356SGlenn Hammond   if (ctx) tr->precheckctx = ctx;
86c9368356SGlenn Hammond   PetscFunctionReturn(0);
87c9368356SGlenn Hammond }
88c9368356SGlenn Hammond 
89c9368356SGlenn Hammond /*@C
90c9368356SGlenn Hammond    SNESNewtonTRGetPreCheck - Gets the pre-check function
91c9368356SGlenn Hammond 
92f6dfbefdSBarry Smith    Deprecated use `SNESNEWTONDCTRDC`
93f6dfbefdSBarry Smith 
94c9368356SGlenn Hammond    Not collective
95c9368356SGlenn Hammond 
96c9368356SGlenn Hammond    Input Parameter:
97c9368356SGlenn Hammond .  snes - the nonlinear solver context
98c9368356SGlenn Hammond 
99c9368356SGlenn Hammond    Output Parameters:
100f6dfbefdSBarry Smith +  func - [optional] function evaluation routine, see for the calling sequence `SNESNewtonTRPreCheck()`
101c9368356SGlenn Hammond -  ctx  - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
102c9368356SGlenn Hammond 
103c9368356SGlenn Hammond    Level: intermediate
104c9368356SGlenn Hammond 
105f6dfbefdSBarry Smith .seealso: `SNESNEWTONDCTRDC`, `SNESNEWTONDCTR`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRPreCheck()`
106c9368356SGlenn Hammond @*/
1079371c9d4SSatish Balay PetscErrorCode SNESNewtonTRGetPreCheck(SNES snes, PetscErrorCode (**func)(SNES, Vec, Vec, PetscBool *, void *), void **ctx) {
108c9368356SGlenn Hammond   SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
109c9368356SGlenn Hammond 
110c9368356SGlenn Hammond   PetscFunctionBegin;
111c9368356SGlenn Hammond   PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
112c9368356SGlenn Hammond   if (func) *func = tr->precheck;
113c9368356SGlenn Hammond   if (ctx) *ctx = tr->precheckctx;
114c9368356SGlenn Hammond   PetscFunctionReturn(0);
115c9368356SGlenn Hammond }
116c9368356SGlenn Hammond 
117c9368356SGlenn Hammond /*@C
1187cb011f5SBarry Smith    SNESNewtonTRSetPostCheck - Sets a user function that is called after the search step has been determined but before the next
1197cb011f5SBarry Smith        function evaluation. Allows the user a chance to change or override the decision of the line search routine
1207cb011f5SBarry Smith 
121f6dfbefdSBarry Smith    Deprecated use `SNESNEWTONDCTRDC`
122f6dfbefdSBarry Smith 
1237cb011f5SBarry Smith    Logically Collective on snes
1247cb011f5SBarry Smith 
1257cb011f5SBarry Smith    Input Parameters:
1267cb011f5SBarry Smith +  snes - the nonlinear solver object
127f6dfbefdSBarry Smith .  func - [optional] function evaluation routine, see `SNESNewtonTRPostCheck()`  for the calling sequence
1287cb011f5SBarry Smith -  ctx  - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
1297cb011f5SBarry Smith 
1307cb011f5SBarry Smith    Level: intermediate
1317cb011f5SBarry Smith 
132f6dfbefdSBarry Smith    Note:
133f6dfbefdSBarry Smith    This function is called BEFORE the function evaluation within the `SNESNEWTONTR` solver while the function set in
134f6dfbefdSBarry Smith    `SNESLineSearchSetPostCheck()` is called AFTER the function evaluation.
1357cb011f5SBarry Smith 
136f6dfbefdSBarry Smith .seealso: `SNESNEWTONDCTRDC`, `SNESNEWTONDCTR`, `SNESNewtonTRPostCheck()`, `SNESNewtonTRGetPostCheck()`
1377cb011f5SBarry Smith @*/
1389371c9d4SSatish Balay PetscErrorCode SNESNewtonTRSetPostCheck(SNES snes, PetscErrorCode (*func)(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *), void *ctx) {
1397cb011f5SBarry Smith   SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
1407cb011f5SBarry Smith 
1417cb011f5SBarry Smith   PetscFunctionBegin;
1427cb011f5SBarry Smith   PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
1437cb011f5SBarry Smith   if (func) tr->postcheck = func;
1447cb011f5SBarry Smith   if (ctx) tr->postcheckctx = ctx;
1457cb011f5SBarry Smith   PetscFunctionReturn(0);
1467cb011f5SBarry Smith }
1477cb011f5SBarry Smith 
1487cb011f5SBarry Smith /*@C
1497cb011f5SBarry Smith    SNESNewtonTRGetPostCheck - Gets the post-check function
1507cb011f5SBarry Smith 
151f6dfbefdSBarry Smith    Deprecated use `SNESNEWTONDCTRDC`
152f6dfbefdSBarry Smith 
1537cb011f5SBarry Smith    Not collective
1547cb011f5SBarry Smith 
1557cb011f5SBarry Smith    Input Parameter:
1567cb011f5SBarry Smith .  snes - the nonlinear solver context
1577cb011f5SBarry Smith 
1587cb011f5SBarry Smith    Output Parameters:
159f6dfbefdSBarry Smith +  func - [optional] function evaluation routine, see for the calling sequence `SNESNewtonTRPostCheck()`
1607cb011f5SBarry Smith -  ctx  - [optional] user-defined context for private data for the function evaluation routine (may be NULL)
1617cb011f5SBarry Smith 
1627cb011f5SBarry Smith    Level: intermediate
1637cb011f5SBarry Smith 
164f6dfbefdSBarry Smith .seealso: `SNESNEWTONDCTRDC`, `SNESNEWTONDCTR`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRPostCheck()`
1657cb011f5SBarry Smith @*/
1669371c9d4SSatish Balay PetscErrorCode SNESNewtonTRGetPostCheck(SNES snes, PetscErrorCode (**func)(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *), void **ctx) {
1677cb011f5SBarry Smith   SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
1687cb011f5SBarry Smith 
1697cb011f5SBarry Smith   PetscFunctionBegin;
1707cb011f5SBarry Smith   PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
1717cb011f5SBarry Smith   if (func) *func = tr->postcheck;
1727cb011f5SBarry Smith   if (ctx) *ctx = tr->postcheckctx;
1737cb011f5SBarry Smith   PetscFunctionReturn(0);
1747cb011f5SBarry Smith }
1757cb011f5SBarry Smith 
1767cb011f5SBarry Smith /*@C
177f6dfbefdSBarry Smith    SNESNewtonTRPreCheck - Called before the step has been determined in `SNESNEWTONTR`
178f6dfbefdSBarry Smith 
179f6dfbefdSBarry Smith    Deprecated use `SNESNEWTONDCTRDC`
180c9368356SGlenn Hammond 
181c9368356SGlenn Hammond    Logically Collective on snes
182c9368356SGlenn Hammond 
183c9368356SGlenn Hammond    Input Parameters:
184c9368356SGlenn Hammond +  snes - the solver
185c9368356SGlenn Hammond .  X - The last solution
186c9368356SGlenn Hammond -  Y - The step direction
187c9368356SGlenn Hammond 
188c9368356SGlenn Hammond    Output Parameters:
189c9368356SGlenn Hammond .  changed_Y - Indicator that the step direction Y has been changed.
190c9368356SGlenn Hammond 
191c9368356SGlenn Hammond    Level: developer
192c9368356SGlenn Hammond 
193f6dfbefdSBarry Smith .seealso: `SNESNEWTONDCTRDC`, `SNESNEWTONDCTR`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRGetPreCheck()`
194c9368356SGlenn Hammond @*/
1959371c9d4SSatish Balay static PetscErrorCode SNESNewtonTRPreCheck(SNES snes, Vec X, Vec Y, PetscBool *changed_Y) {
196c9368356SGlenn Hammond   SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
197c9368356SGlenn Hammond 
198c9368356SGlenn Hammond   PetscFunctionBegin;
199c9368356SGlenn Hammond   *changed_Y = PETSC_FALSE;
200c9368356SGlenn Hammond   if (tr->precheck) {
2019566063dSJacob Faibussowitsch     PetscCall((*tr->precheck)(snes, X, Y, changed_Y, tr->precheckctx));
202c9368356SGlenn Hammond     PetscValidLogicalCollectiveBool(snes, *changed_Y, 4);
203c9368356SGlenn Hammond   }
204c9368356SGlenn Hammond   PetscFunctionReturn(0);
205c9368356SGlenn Hammond }
206c9368356SGlenn Hammond 
207c9368356SGlenn Hammond /*@C
208f6dfbefdSBarry Smith    SNESNewtonTRPostCheck - Called after the step has been determined in `SNESNEWTONTR` but before the function evaluation
209f6dfbefdSBarry Smith 
210f6dfbefdSBarry Smith    Deprecated use `SNESNEWTONDCTRDC`
2117cb011f5SBarry Smith 
2127cb011f5SBarry Smith    Logically Collective on snes
2137cb011f5SBarry Smith 
2147cb011f5SBarry Smith    Input Parameters:
2156b867d5aSJose E. Roman +  snes - the solver
2166b867d5aSJose E. Roman .  X - The last solution
2177cb011f5SBarry Smith .  Y - The full step direction
2183312a946SBarry Smith -  W - The updated solution, W = X - Y
2197cb011f5SBarry Smith 
2207cb011f5SBarry Smith    Output Parameters:
2213312a946SBarry Smith +  changed_Y - indicator if step has been changed
2223312a946SBarry Smith -  changed_W - Indicator if the new candidate solution W has been changed.
2237cb011f5SBarry Smith 
224f6dfbefdSBarry Smith    Note:
2253312a946SBarry Smith      If Y is changed then W is recomputed as X - Y
2267cb011f5SBarry Smith 
2277cb011f5SBarry Smith    Level: developer
2287cb011f5SBarry Smith 
229f6dfbefdSBarry Smith .seealso: `SNESNEWTONDCTRDC`, `SNESNEWTONDCTR`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`
2307cb011f5SBarry Smith @*/
2319371c9d4SSatish Balay static PetscErrorCode SNESNewtonTRPostCheck(SNES snes, Vec X, Vec Y, Vec W, PetscBool *changed_Y, PetscBool *changed_W) {
2327cb011f5SBarry Smith   SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
2337cb011f5SBarry Smith 
2347cb011f5SBarry Smith   PetscFunctionBegin;
235c9368356SGlenn Hammond   *changed_Y = PETSC_FALSE;
2367cb011f5SBarry Smith   *changed_W = PETSC_FALSE;
2377cb011f5SBarry Smith   if (tr->postcheck) {
2389566063dSJacob Faibussowitsch     PetscCall((*tr->postcheck)(snes, X, Y, W, changed_Y, changed_W, tr->postcheckctx));
239c9368356SGlenn Hammond     PetscValidLogicalCollectiveBool(snes, *changed_Y, 5);
2407cb011f5SBarry Smith     PetscValidLogicalCollectiveBool(snes, *changed_W, 6);
2417cb011f5SBarry Smith   }
2427cb011f5SBarry Smith   PetscFunctionReturn(0);
2437cb011f5SBarry Smith }
24485385478SLisandro Dalcin 
245df60cc22SBarry Smith /*
24604d7464bSBarry Smith    SNESSolve_NEWTONTR - Implements Newton's Method with a very simple trust
247ddbbdb52SLois Curfman McInnes    region approach for solving systems of nonlinear equations.
2484800dd8cSBarry Smith 
2494800dd8cSBarry Smith */
2509371c9d4SSatish Balay static PetscErrorCode SNESSolve_NEWTONTR(SNES snes) {
25104d7464bSBarry Smith   SNES_NEWTONTR            *neP = (SNES_NEWTONTR *)snes->data;
252c9368356SGlenn Hammond   Vec                       X, F, Y, G, Ytmp, W;
253a7cc72afSBarry Smith   PetscInt                  maxits, i, lits;
25485385478SLisandro Dalcin   PetscReal                 rho, fnorm, gnorm, gpnorm, xnorm = 0, delta, nrm, ynorm, norm1;
25579f36460SBarry Smith   PetscScalar               cnorm;
256df60cc22SBarry Smith   KSP                       ksp;
2573f149594SLisandro Dalcin   SNESConvergedReason       reason   = SNES_CONVERGED_ITERATING;
258df8705c3SBarry Smith   PetscBool                 breakout = PETSC_FALSE;
259df8705c3SBarry Smith   SNES_TR_KSPConverged_Ctx *ctx;
2605e28bcb6Sprj-   PetscErrorCode (*convtest)(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *), (*convdestroy)(void *);
2615e28bcb6Sprj-   void *convctx;
2624800dd8cSBarry Smith 
2633a40ed3dSBarry Smith   PetscFunctionBegin;
2640b121fc5SBarry Smith   PetscCheck(!snes->xl && !snes->xu && !snes->ops->computevariablebounds, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
265c579b300SPatrick Farrell 
266fbe28522SBarry Smith   maxits = snes->max_its;  /* maximum number of iterations */
267fbe28522SBarry Smith   X      = snes->vec_sol;  /* solution vector */
26839e2f89bSBarry Smith   F      = snes->vec_func; /* residual vector */
269fbe28522SBarry Smith   Y      = snes->work[0];  /* work vectors */
270fbe28522SBarry Smith   G      = snes->work[1];
2716b5873e3SBarry Smith   Ytmp   = snes->work[2];
272c9368356SGlenn Hammond   W      = snes->work[3];
2734800dd8cSBarry Smith 
2749566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
2754c49b128SBarry Smith   snes->iter = 0;
2769566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
2774800dd8cSBarry Smith 
278df8705c3SBarry Smith   /* Set the linear stopping criteria to use the More' trick. */
2799566063dSJacob Faibussowitsch   PetscCall(SNESGetKSP(snes, &ksp));
2809566063dSJacob Faibussowitsch   PetscCall(KSPGetConvergenceTest(ksp, &convtest, &convctx, &convdestroy));
281df8705c3SBarry Smith   if (convtest != SNESTR_KSPConverged_Private) {
2829566063dSJacob Faibussowitsch     PetscCall(PetscNew(&ctx));
283df8705c3SBarry Smith     ctx->snes = snes;
2849566063dSJacob Faibussowitsch     PetscCall(KSPGetAndClearConvergenceTest(ksp, &ctx->convtest, &ctx->convctx, &ctx->convdestroy));
2859566063dSJacob Faibussowitsch     PetscCall(KSPSetConvergenceTest(ksp, SNESTR_KSPConverged_Private, ctx, SNESTR_KSPConverged_Destroy));
2869566063dSJacob Faibussowitsch     PetscCall(PetscInfo(snes, "Using Krylov convergence test SNESTR_KSPConverged_Private\n"));
287df8705c3SBarry Smith   }
288df8705c3SBarry Smith 
289e4ed7901SPeter Brune   if (!snes->vec_func_init_set) {
2909566063dSJacob Faibussowitsch     PetscCall(SNESComputeFunction(snes, X, F)); /* F(X) */
2911aa26658SKarl Rupp   } else snes->vec_func_init_set = PETSC_FALSE;
292e4ed7901SPeter Brune 
2939566063dSJacob Faibussowitsch   PetscCall(VecNorm(F, NORM_2, &fnorm)); /* fnorm <- || F || */
294422a814eSBarry Smith   SNESCheckFunctionNorm(snes, fnorm);
2959566063dSJacob Faibussowitsch   PetscCall(VecNorm(X, NORM_2, &xnorm)); /* xnorm <- || X || */
2969566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
297fbe28522SBarry Smith   snes->norm = fnorm;
2989566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
2998b84755bSBarry Smith   delta      = xnorm ? neP->delta0 * xnorm : neP->delta0;
3004800dd8cSBarry Smith   neP->delta = delta;
3019566063dSJacob Faibussowitsch   PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
3029566063dSJacob Faibussowitsch   PetscCall(SNESMonitor(snes, 0, fnorm));
303b37302c6SLois Curfman McInnes 
30485385478SLisandro Dalcin   /* test convergence */
305dbbe0bcdSBarry Smith   PetscUseTypeMethod(snes, converged, snes->iter, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
30685385478SLisandro Dalcin   if (snes->reason) PetscFunctionReturn(0);
3073f149594SLisandro Dalcin 
3084800dd8cSBarry Smith   for (i = 0; i < maxits; i++) {
30999a96b7cSMatthew Knepley     /* Call general purpose update function */
310dbbe0bcdSBarry Smith     PetscTryTypeMethod(snes, update, snes->iter);
31199a96b7cSMatthew Knepley 
31285385478SLisandro Dalcin     /* Solve J Y = F, where J is Jacobian matrix */
3139566063dSJacob Faibussowitsch     PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
31407b62357SFande Kong     SNESCheckJacobianDomainerror(snes);
3159566063dSJacob Faibussowitsch     PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
3169566063dSJacob Faibussowitsch     PetscCall(KSPSolve(snes->ksp, F, Ytmp));
3179566063dSJacob Faibussowitsch     PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
318329f5518SBarry Smith     snes->linear_its += lits;
3191aa26658SKarl Rupp 
3209566063dSJacob Faibussowitsch     PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
3219566063dSJacob Faibussowitsch     PetscCall(VecNorm(Ytmp, NORM_2, &nrm));
322064f8208SBarry Smith     norm1 = nrm;
323284fb49fSHeeho Park 
3246b5873e3SBarry Smith     while (1) {
325c9368356SGlenn Hammond       PetscBool changed_y;
3267cb011f5SBarry Smith       PetscBool changed_w;
3279566063dSJacob Faibussowitsch       PetscCall(VecCopy(Ytmp, Y));
328064f8208SBarry Smith       nrm = norm1;
3296b5873e3SBarry Smith 
3302deea200SLois Curfman McInnes       /* Scale Y if need be and predict new value of F norm */
331064f8208SBarry Smith       if (nrm >= delta) {
332064f8208SBarry Smith         nrm    = delta / nrm;
333064f8208SBarry Smith         gpnorm = (1.0 - nrm) * fnorm;
334064f8208SBarry Smith         cnorm  = nrm;
3359566063dSJacob Faibussowitsch         PetscCall(PetscInfo(snes, "Scaling direction by %g\n", (double)nrm));
3369566063dSJacob Faibussowitsch         PetscCall(VecScale(Y, cnorm));
337064f8208SBarry Smith         nrm   = gpnorm;
338fbe28522SBarry Smith         ynorm = delta;
339fbe28522SBarry Smith       } else {
340fbe28522SBarry Smith         gpnorm = 0.0;
3419566063dSJacob Faibussowitsch         PetscCall(PetscInfo(snes, "Direction is in Trust Region\n"));
342064f8208SBarry Smith         ynorm = nrm;
343fbe28522SBarry Smith       }
344c9368356SGlenn Hammond       /* PreCheck() allows for updates to Y prior to W <- X - Y */
3459566063dSJacob Faibussowitsch       PetscCall(SNESNewtonTRPreCheck(snes, X, Y, &changed_y));
3469566063dSJacob Faibussowitsch       PetscCall(VecWAXPY(W, -1.0, Y, X)); /* W <- X - Y */
3479566063dSJacob Faibussowitsch       PetscCall(SNESNewtonTRPostCheck(snes, X, Y, W, &changed_y, &changed_w));
3489566063dSJacob Faibussowitsch       if (changed_y) PetscCall(VecWAXPY(W, -1.0, Y, X));
3499566063dSJacob Faibussowitsch       PetscCall(VecCopy(Y, snes->vec_sol_update));
3509566063dSJacob Faibussowitsch       PetscCall(SNESComputeFunction(snes, W, G)); /*  F(X-Y) = G */
3519566063dSJacob Faibussowitsch       PetscCall(VecNorm(G, NORM_2, &gnorm));      /* gnorm <- || g || */
3527551ef16SBarry Smith       SNESCheckFunctionNorm(snes, gnorm);
3534800dd8cSBarry Smith       if (fnorm == gpnorm) rho = 0.0;
354fbe28522SBarry Smith       else rho = (fnorm * fnorm - gnorm * gnorm) / (fnorm * fnorm - gpnorm * gpnorm);
3554800dd8cSBarry Smith 
3564800dd8cSBarry Smith       /* Update size of trust region */
3574800dd8cSBarry Smith       if (rho < neP->mu) delta *= neP->delta1;
3584800dd8cSBarry Smith       else if (rho < neP->eta) delta *= neP->delta2;
3594800dd8cSBarry Smith       else delta *= neP->delta3;
3609566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes, "fnorm=%g, gnorm=%g, ynorm=%g\n", (double)fnorm, (double)gnorm, (double)ynorm));
3619566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes, "gpred=%g, rho=%g, delta=%g\n", (double)gpnorm, (double)rho, (double)delta));
3621aa26658SKarl Rupp 
3634800dd8cSBarry Smith       neP->delta = delta;
3644800dd8cSBarry Smith       if (rho > neP->sigma) break;
3659566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes, "Trying again in smaller region\n"));
366284fb49fSHeeho Park 
3676b5873e3SBarry Smith       /* check to see if progress is hopeless */
368ef87ac0dSKris Buschelman       neP->itflag = PETSC_FALSE;
3699566063dSJacob Faibussowitsch       PetscCall(SNESTR_Converged_Private(snes, snes->iter, xnorm, ynorm, fnorm, &reason, snes->cnvP));
370dbbe0bcdSBarry Smith       if (!reason) PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &reason, snes->cnvP);
3717cb011f5SBarry Smith       if (reason == SNES_CONVERGED_SNORM_RELATIVE) reason = SNES_DIVERGED_INNER;
372184914b5SBarry Smith       if (reason) {
37352392280SLois Curfman McInnes         /* We're not progressing, so return with the current iterate */
3749566063dSJacob Faibussowitsch         PetscCall(SNESMonitor(snes, i + 1, fnorm));
375a7cc72afSBarry Smith         breakout = PETSC_TRUE;
376454a90a3SBarry Smith         break;
37752392280SLois Curfman McInnes       }
37850ffb88aSMatthew Knepley       snes->numFailures++;
3794800dd8cSBarry Smith     }
3801acabf8cSLois Curfman McInnes     if (!breakout) {
38185385478SLisandro Dalcin       /* Update function and solution vectors */
3824800dd8cSBarry Smith       fnorm = gnorm;
3839566063dSJacob Faibussowitsch       PetscCall(VecCopy(G, F));
3849566063dSJacob Faibussowitsch       PetscCall(VecCopy(W, X));
38585385478SLisandro Dalcin       /* Monitor convergence */
3869566063dSJacob Faibussowitsch       PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
387c509a162SBarry Smith       snes->iter  = i + 1;
388fbe28522SBarry Smith       snes->norm  = fnorm;
389c1e67a49SFande Kong       snes->xnorm = xnorm;
390c1e67a49SFande Kong       snes->ynorm = ynorm;
3919566063dSJacob Faibussowitsch       PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
3929566063dSJacob Faibussowitsch       PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
3939566063dSJacob Faibussowitsch       PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
39485385478SLisandro Dalcin       /* Test for convergence, xnorm = || X || */
395ef87ac0dSKris Buschelman       neP->itflag = PETSC_TRUE;
3969566063dSJacob Faibussowitsch       if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X, NORM_2, &xnorm));
397dbbe0bcdSBarry Smith       PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &reason, snes->cnvP);
3983f149594SLisandro Dalcin       if (reason) break;
3991aa26658SKarl Rupp     } else break;
40038442cffSBarry Smith   }
401284fb49fSHeeho Park 
40252392280SLois Curfman McInnes   if (i == maxits) {
4039566063dSJacob Faibussowitsch     PetscCall(PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits));
40485385478SLisandro Dalcin     if (!reason) reason = SNES_DIVERGED_MAX_IT;
40552392280SLois Curfman McInnes   }
4069566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
407184914b5SBarry Smith   snes->reason = reason;
4089566063dSJacob Faibussowitsch   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
4095e28bcb6Sprj-   if (convtest != SNESTR_KSPConverged_Private) {
4109566063dSJacob Faibussowitsch     PetscCall(KSPGetAndClearConvergenceTest(ksp, &ctx->convtest, &ctx->convctx, &ctx->convdestroy));
4119566063dSJacob Faibussowitsch     PetscCall(PetscFree(ctx));
4129566063dSJacob Faibussowitsch     PetscCall(KSPSetConvergenceTest(ksp, convtest, convctx, convdestroy));
4135e28bcb6Sprj-   }
4143a40ed3dSBarry Smith   PetscFunctionReturn(0);
4154800dd8cSBarry Smith }
416284fb49fSHeeho Park 
4179371c9d4SSatish Balay static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes) {
4183a40ed3dSBarry Smith   PetscFunctionBegin;
4199566063dSJacob Faibussowitsch   PetscCall(SNESSetWorkVecs(snes, 4));
4209566063dSJacob Faibussowitsch   PetscCall(SNESSetUpMatrices(snes));
4213a40ed3dSBarry Smith   PetscFunctionReturn(0);
4224800dd8cSBarry Smith }
4236b8b9a38SLisandro Dalcin 
4249371c9d4SSatish Balay PetscErrorCode SNESReset_NEWTONTR(SNES snes) {
4256b8b9a38SLisandro Dalcin   PetscFunctionBegin;
4266b8b9a38SLisandro Dalcin   PetscFunctionReturn(0);
4276b8b9a38SLisandro Dalcin }
4286b8b9a38SLisandro Dalcin 
4299371c9d4SSatish Balay static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes) {
4303a40ed3dSBarry Smith   PetscFunctionBegin;
4319566063dSJacob Faibussowitsch   PetscCall(SNESReset_NEWTONTR(snes));
4329566063dSJacob Faibussowitsch   PetscCall(PetscFree(snes->data));
4333a40ed3dSBarry Smith   PetscFunctionReturn(0);
4344800dd8cSBarry Smith }
4354800dd8cSBarry Smith 
4369371c9d4SSatish Balay static PetscErrorCode SNESSetFromOptions_NEWTONTR(SNES snes, PetscOptionItems *PetscOptionsObject) {
43704d7464bSBarry Smith   SNES_NEWTONTR *ctx = (SNES_NEWTONTR *)snes->data;
4384800dd8cSBarry Smith 
4393a40ed3dSBarry Smith   PetscFunctionBegin;
440d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject, "SNES trust region options for nonlinear equations");
4419566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_trtol", "Trust region tolerance", "SNESSetTrustRegionTolerance", snes->deltatol, &snes->deltatol, NULL));
4429566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_tr_mu", "mu", "None", ctx->mu, &ctx->mu, NULL));
4439566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_tr_eta", "eta", "None", ctx->eta, &ctx->eta, NULL));
4449566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_tr_sigma", "sigma", "None", ctx->sigma, &ctx->sigma, NULL));
4459566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_tr_delta0", "delta0", "None", ctx->delta0, &ctx->delta0, NULL));
4469566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_tr_delta1", "delta1", "None", ctx->delta1, &ctx->delta1, NULL));
4479566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_tr_delta2", "delta2", "None", ctx->delta2, &ctx->delta2, NULL));
4489566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_tr_delta3", "delta3", "None", ctx->delta3, &ctx->delta3, NULL));
449d0609cedSBarry Smith   PetscOptionsHeadEnd();
4503a40ed3dSBarry Smith   PetscFunctionReturn(0);
4514800dd8cSBarry Smith }
4524800dd8cSBarry Smith 
4539371c9d4SSatish Balay static PetscErrorCode SNESView_NEWTONTR(SNES snes, PetscViewer viewer) {
45404d7464bSBarry Smith   SNES_NEWTONTR *tr = (SNES_NEWTONTR *)snes->data;
455ace3abfcSBarry Smith   PetscBool      iascii;
456a935fc98SLois Curfman McInnes 
4573a40ed3dSBarry Smith   PetscFunctionBegin;
4589566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
45932077d6dSBarry Smith   if (iascii) {
46063a3b9bcSJacob Faibussowitsch     PetscCall(PetscViewerASCIIPrintf(viewer, "  Trust region tolerance %g (-snes_trtol)\n", (double)snes->deltatol));
4619566063dSJacob Faibussowitsch     PetscCall(PetscViewerASCIIPrintf(viewer, "  mu=%g, eta=%g, sigma=%g\n", (double)tr->mu, (double)tr->eta, (double)tr->sigma));
4629566063dSJacob Faibussowitsch     PetscCall(PetscViewerASCIIPrintf(viewer, "  delta0=%g, delta1=%g, delta2=%g, delta3=%g\n", (double)tr->delta0, (double)tr->delta1, (double)tr->delta2, (double)tr->delta3));
46319bcc07fSBarry Smith   }
4643a40ed3dSBarry Smith   PetscFunctionReturn(0);
465a935fc98SLois Curfman McInnes }
466f6dfbefdSBarry Smith 
467ebe3b25bSBarry Smith /*MC
46804d7464bSBarry Smith       SNESNEWTONTR - Newton based nonlinear solver that uses a trust region
469ebe3b25bSBarry Smith 
470f6dfbefdSBarry Smith       Deprecated use `SNESNEWTONTRDC`
471f6dfbefdSBarry Smith 
472f6dfbefdSBarry Smith    Options Database Keys:
4738e434772SBarry Smith +    -snes_trtol <tol> - trust region tolerance
4748e434772SBarry Smith .    -snes_tr_mu <mu> - trust region parameter
4758e434772SBarry Smith .    -snes_tr_eta <eta> - trust region parameter
4768e434772SBarry Smith .    -snes_tr_sigma <sigma> - trust region parameter
4778b84755bSBarry Smith .    -snes_tr_delta0 <delta0> -  initial size of the trust region is delta0*norm2(x)
4788e434772SBarry Smith .    -snes_tr_delta1 <delta1> - trust region parameter
4798e434772SBarry Smith .    -snes_tr_delta2 <delta2> - trust region parameter
4808e434772SBarry Smith -    -snes_tr_delta3 <delta3> - trust region parameter
481ebe3b25bSBarry Smith 
482f6dfbefdSBarry Smith    Reference:
483f6dfbefdSBarry Smith .  - *  "The Minpack Project", by More', Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development
484ebe3b25bSBarry Smith    of Mathematical Software", Wayne Cowell, editor.
485ebe3b25bSBarry Smith 
486ee3001cbSBarry Smith    Level: intermediate
487ee3001cbSBarry Smith 
488f6dfbefdSBarry Smith .seealso: `SNESNEWTONTRDC`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESSetTrustRegionTolerance()`
489ebe3b25bSBarry Smith M*/
4909371c9d4SSatish Balay PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes) {
49104d7464bSBarry Smith   SNES_NEWTONTR *neP;
4924800dd8cSBarry Smith 
4933a40ed3dSBarry Smith   PetscFunctionBegin;
49404d7464bSBarry Smith   snes->ops->setup          = SNESSetUp_NEWTONTR;
49504d7464bSBarry Smith   snes->ops->solve          = SNESSolve_NEWTONTR;
49604d7464bSBarry Smith   snes->ops->destroy        = SNESDestroy_NEWTONTR;
49704d7464bSBarry Smith   snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR;
49804d7464bSBarry Smith   snes->ops->view           = SNESView_NEWTONTR;
49904d7464bSBarry Smith   snes->ops->reset          = SNESReset_NEWTONTR;
500fbe28522SBarry Smith 
50142f4f86dSBarry Smith   snes->usesksp = PETSC_TRUE;
502efd4aadfSBarry Smith   snes->usesnpc = PETSC_FALSE;
50342f4f86dSBarry Smith 
5044fc747eaSLawrence Mitchell   snes->alwayscomputesfinalresidual = PETSC_TRUE;
5054fc747eaSLawrence Mitchell 
506*4dfa11a4SJacob Faibussowitsch   PetscCall(PetscNew(&neP));
507fbe28522SBarry Smith   snes->data  = (void *)neP;
508fbe28522SBarry Smith   neP->mu     = 0.25;
509fbe28522SBarry Smith   neP->eta    = 0.75;
510fbe28522SBarry Smith   neP->delta  = 0.0;
511fbe28522SBarry Smith   neP->delta0 = 0.2;
512fbe28522SBarry Smith   neP->delta1 = 0.3;
513fbe28522SBarry Smith   neP->delta2 = 0.75;
514fbe28522SBarry Smith   neP->delta3 = 2.0;
515fbe28522SBarry Smith   neP->sigma  = 0.0001;
516ef87ac0dSKris Buschelman   neP->itflag = PETSC_FALSE;
51775567043SBarry Smith   neP->rnorm0 = 0.0;
51875567043SBarry Smith   neP->ttol   = 0.0;
5193a40ed3dSBarry Smith   PetscFunctionReturn(0);
5204800dd8cSBarry Smith }
521