xref: /petsc/src/snes/linesearch/impls/bt/linesearchbt.c (revision 1d27aa22b2f6148b2c4e3f06a75e0638d6493e09)
1af0996ceSBarry Smith #include <petsc/private/linesearchimpl.h> /*I  "petscsnes.h"  I*/
2af0996ceSBarry Smith #include <petsc/private/snesimpl.h>
371b4ebd2SPeter Brune 
471b4ebd2SPeter Brune typedef struct {
571b4ebd2SPeter Brune   PetscReal alpha; /* sufficient decrease parameter */
6f1c6b773SPeter Brune } SNESLineSearch_BT;
771b4ebd2SPeter Brune 
82f4102e2SPeter Brune /*@
9420bcc1bSBarry Smith   SNESLineSearchBTSetAlpha - Sets the descent parameter, `alpha`, in the `SNESLINESEARCHBT` `SNESLineSearch` variant.
102f4102e2SPeter Brune 
112f4102e2SPeter Brune   Input Parameters:
122f4102e2SPeter Brune + linesearch - linesearch context
132f4102e2SPeter Brune - alpha      - The descent parameter
142f4102e2SPeter Brune 
152f4102e2SPeter Brune   Level: intermediate
162f4102e2SPeter Brune 
17420bcc1bSBarry Smith .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchSetLambda()`, `SNESLineSearchGetTolerances()`, `SNESLINESEARCHBT`, `SNESLineSearchBTGetAlpha()`
182f4102e2SPeter Brune @*/
19d71ae5a4SJacob Faibussowitsch PetscErrorCode SNESLineSearchBTSetAlpha(SNESLineSearch linesearch, PetscReal alpha)
20d71ae5a4SJacob Faibussowitsch {
21d7ca06c0SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
226e111a19SKarl Rupp 
232f4102e2SPeter Brune   PetscFunctionBegin;
242f4102e2SPeter Brune   PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 1);
252f4102e2SPeter Brune   bt->alpha = alpha;
263ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
272f4102e2SPeter Brune }
282f4102e2SPeter Brune 
292f4102e2SPeter Brune /*@
30420bcc1bSBarry Smith   SNESLineSearchBTGetAlpha - Gets the descent parameter, `alpha`, in the `SNESLINESEARCHBT` variant that was set with `SNESLineSearchBTSetAlpha()`
312f4102e2SPeter Brune 
32f6dfbefdSBarry Smith   Input Parameter:
332f4102e2SPeter Brune . linesearch - linesearch context
348e557f58SPeter Brune 
35f6dfbefdSBarry Smith   Output Parameter:
362f4102e2SPeter Brune . alpha - The descent parameter
372f4102e2SPeter Brune 
382f4102e2SPeter Brune   Level: intermediate
392f4102e2SPeter Brune 
40420bcc1bSBarry Smith .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchGetLambda()`, `SNESLineSearchGetTolerances()` `SNESLINESEARCHBT`, `SNESLineSearchBTSetAlpha()`
412f4102e2SPeter Brune @*/
42d71ae5a4SJacob Faibussowitsch PetscErrorCode SNESLineSearchBTGetAlpha(SNESLineSearch linesearch, PetscReal *alpha)
43d71ae5a4SJacob Faibussowitsch {
44d7ca06c0SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
456e111a19SKarl Rupp 
462f4102e2SPeter Brune   PetscFunctionBegin;
472f4102e2SPeter Brune   PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 1);
482f4102e2SPeter Brune   *alpha = bt->alpha;
493ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
502f4102e2SPeter Brune }
512f4102e2SPeter Brune 
52d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchApply_BT(SNESLineSearch linesearch)
53d71ae5a4SJacob Faibussowitsch {
5471b4ebd2SPeter Brune   PetscBool          changed_y, changed_w;
5571b4ebd2SPeter Brune   Vec                X, F, Y, W, G;
5671b4ebd2SPeter Brune   SNES               snes;
578a430afbSPeter Brune   PetscReal          fnorm, xnorm, ynorm, gnorm;
583b288129SPeter Brune   PetscReal          lambda, lambdatemp, lambdaprev, minlambda, maxstep, initslope, alpha, stol;
59dfcd3828SPeter Brune   PetscReal          t1, t2, a, b, d;
608a430afbSPeter Brune   PetscReal          f;
618a430afbSPeter Brune   PetscReal          g, gprev;
6271b4ebd2SPeter Brune   PetscViewer        monitor;
6371b4ebd2SPeter Brune   PetscInt           max_its, count;
64d7ca06c0SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
6571b4ebd2SPeter Brune   Mat                jac;
6632661dccSStefano Zampini   const char *const  ordStr[] = {"Linear", "Quadratic", "Cubic"};
67bd4a8a71SBarry Smith   PetscErrorCode (*objective)(SNES, Vec, PetscReal *, void *);
6871b4ebd2SPeter Brune 
6971b4ebd2SPeter Brune   PetscFunctionBegin;
709566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G));
71a68bbae5SBarry Smith   PetscCall(SNESLineSearchGetNorms(linesearch, NULL, &fnorm, NULL));
729566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetLambda(linesearch, &lambda));
739566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
749566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor));
759566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetTolerances(linesearch, &minlambda, &maxstep, NULL, NULL, NULL, &max_its));
769566063dSJacob Faibussowitsch   PetscCall(SNESGetTolerances(snes, NULL, NULL, &stol, NULL, NULL));
779566063dSJacob Faibussowitsch   PetscCall(SNESGetObjective(snes, &objective, NULL));
7871b4ebd2SPeter Brune   alpha = bt->alpha;
7971b4ebd2SPeter Brune 
809566063dSJacob Faibussowitsch   PetscCall(SNESGetJacobian(snes, &jac, NULL, NULL, NULL));
8108401ef6SPierre Jolivet   PetscCheck(jac || objective, PetscObjectComm((PetscObject)linesearch), PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix");
82c69d1a72SBarry Smith 
839566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y));
849566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
8571b4ebd2SPeter Brune 
869566063dSJacob Faibussowitsch   PetscCall(VecNormBegin(Y, NORM_2, &ynorm));
879566063dSJacob Faibussowitsch   PetscCall(VecNormBegin(X, NORM_2, &xnorm));
889566063dSJacob Faibussowitsch   PetscCall(VecNormEnd(Y, NORM_2, &ynorm));
899566063dSJacob Faibussowitsch   PetscCall(VecNormEnd(X, NORM_2, &xnorm));
903b288129SPeter Brune 
9171b4ebd2SPeter Brune   if (ynorm == 0.0) {
9271b4ebd2SPeter Brune     if (monitor) {
939566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
949566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Initial direction and size is 0\n"));
959566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
9671b4ebd2SPeter Brune     }
979566063dSJacob Faibussowitsch     PetscCall(VecCopy(X, W));
989566063dSJacob Faibussowitsch     PetscCall(VecCopy(F, G));
999566063dSJacob Faibussowitsch     PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm));
1009566063dSJacob Faibussowitsch     PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
1013ba16761SJacob Faibussowitsch     PetscFunctionReturn(PETSC_SUCCESS);
10271b4ebd2SPeter Brune   }
10371b4ebd2SPeter Brune   if (ynorm > maxstep) { /* Step too big, so scale back */
10471b4ebd2SPeter Brune     if (monitor) {
1059566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
1069566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Scaling step by %14.12e old ynorm %14.12e\n", (double)(maxstep / ynorm), (double)ynorm));
1079566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
10871b4ebd2SPeter Brune     }
1099566063dSJacob Faibussowitsch     PetscCall(VecScale(Y, maxstep / (ynorm)));
11071b4ebd2SPeter Brune     ynorm = maxstep;
11171b4ebd2SPeter Brune   }
1128a430afbSPeter Brune 
1138a430afbSPeter Brune   /* if the SNES has an objective set, use that instead of the function value */
114bd4a8a71SBarry Smith   if (objective) {
1159566063dSJacob Faibussowitsch     PetscCall(SNESComputeObjective(snes, X, &f));
1168a430afbSPeter Brune   } else {
117e3ce9940SStefano Zampini     f = 0.5 * PetscSqr(fnorm);
1188a430afbSPeter Brune   }
1198a430afbSPeter Brune 
1208a430afbSPeter Brune   /* compute the initial slope */
121bd4a8a71SBarry Smith   if (objective) {
12232661dccSStefano Zampini     /* slope comes from the function (assumed to be the gradient of the objective) */
1239566063dSJacob Faibussowitsch     PetscCall(VecDotRealPart(Y, F, &initslope));
1248a430afbSPeter Brune   } else {
1258a430afbSPeter Brune     /* slope comes from the normal equations */
1269566063dSJacob Faibussowitsch     PetscCall(MatMult(jac, Y, W));
1279566063dSJacob Faibussowitsch     PetscCall(VecDotRealPart(F, W, &initslope));
12871b4ebd2SPeter Brune     if (initslope > 0.0) initslope = -initslope;
12971b4ebd2SPeter Brune     if (initslope == 0.0) initslope = -1.0;
1308a430afbSPeter Brune   }
13171b4ebd2SPeter Brune 
132df019d78SBarry Smith   while (PETSC_TRUE) {
1339566063dSJacob Faibussowitsch     PetscCall(VecWAXPY(W, -lambda, Y, X));
1341baa6e33SBarry Smith     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
135e71169deSBarry Smith     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
1369566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while checking full step length!\n"));
13771b4ebd2SPeter Brune       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
1389566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
1393ba16761SJacob Faibussowitsch       PetscFunctionReturn(PETSC_SUCCESS);
14071b4ebd2SPeter Brune     }
1418a430afbSPeter Brune 
142bd4a8a71SBarry Smith     if (objective) {
1439566063dSJacob Faibussowitsch       PetscCall(SNESComputeObjective(snes, W, &g));
1448a430afbSPeter Brune     } else {
1459566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
1469bd66eb0SPeter Brune       if (linesearch->ops->vinorm) {
1479bd66eb0SPeter Brune         gnorm = fnorm;
1489566063dSJacob Faibussowitsch         PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
1499bd66eb0SPeter Brune       } else {
1509566063dSJacob Faibussowitsch         PetscCall(VecNorm(G, NORM_2, &gnorm));
1518a430afbSPeter Brune       }
152e3ce9940SStefano Zampini       g = 0.5 * PetscSqr(gnorm);
1539bd66eb0SPeter Brune     }
1549566063dSJacob Faibussowitsch     PetscCall(SNESLineSearchMonitor(linesearch));
1559bd66eb0SPeter Brune 
156df019d78SBarry Smith     if (!PetscIsInfOrNanReal(g)) break;
157df019d78SBarry Smith     if (monitor) {
1589566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
1599566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n", (double)lambda));
1609566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
161df019d78SBarry Smith     }
162ad540459SPierre Jolivet     if (lambda <= minlambda) SNESCheckFunctionNorm(snes, g);
16332661dccSStefano Zampini     lambda *= .5;
164df019d78SBarry Smith   }
165df019d78SBarry Smith 
16648a46eb9SPierre Jolivet   if (!objective) PetscCall(PetscInfo(snes, "Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm));
167e3ce9940SStefano Zampini   if (g <= f + lambda * alpha * initslope) { /* Sufficient reduction or step tolerance convergence */
16871b4ebd2SPeter Brune     if (monitor) {
1699566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
170bd4a8a71SBarry Smith       if (!objective) {
1719566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm));
1728a430afbSPeter Brune       } else {
173e3ce9940SStefano Zampini         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Using full step: old obj %14.12e new obj %14.12e\n", (double)f, (double)g));
1748a430afbSPeter Brune       }
1759566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
17671b4ebd2SPeter Brune     }
17771b4ebd2SPeter Brune   } else {
178c21ba15cSPeter Brune     if (stol * xnorm > ynorm) {
179a68bbae5SBarry Smith       /* Since the full step didn't give sufficient decrease and the step is tiny, exit */
1809566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm));
1819566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
182c61ba1e2SPeter Brune       if (monitor) {
1839566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
18463a3b9bcSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n", (double)ynorm, (double)(stol * xnorm)));
1859566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
186c61ba1e2SPeter Brune       }
1873ba16761SJacob Faibussowitsch       PetscFunctionReturn(PETSC_SUCCESS);
188c21ba15cSPeter Brune     }
18932661dccSStefano Zampini     /* Here to avoid -Wmaybe-uninitiliazed warnings */
19071b4ebd2SPeter Brune     lambdaprev = lambda;
1918a430afbSPeter Brune     gprev      = g;
19232661dccSStefano Zampini     if (linesearch->order != SNES_LINESEARCH_ORDER_LINEAR) {
19332661dccSStefano Zampini       /* Fit points with quadratic */
194e3ce9940SStefano Zampini       lambdatemp = -initslope * PetscSqr(lambda) / (2.0 * (g - f - lambda * initslope));
195e3ce9940SStefano Zampini       lambda     = PetscClipInterval(lambdatemp, .1 * lambda, .5 * lambda);
19671b4ebd2SPeter Brune 
1979566063dSJacob Faibussowitsch       PetscCall(VecWAXPY(W, -lambda, Y, X));
1981baa6e33SBarry Smith       if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
199e71169deSBarry Smith       if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
20063a3b9bcSJacob Faibussowitsch         PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while attempting quadratic backtracking! %" PetscInt_FMT " \n", snes->nfuncs));
20171b4ebd2SPeter Brune         snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
2029566063dSJacob Faibussowitsch         PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
2033ba16761SJacob Faibussowitsch         PetscFunctionReturn(PETSC_SUCCESS);
20471b4ebd2SPeter Brune       }
205bd4a8a71SBarry Smith       if (objective) {
2069566063dSJacob Faibussowitsch         PetscCall(SNESComputeObjective(snes, W, &g));
2078a430afbSPeter Brune       } else {
2089566063dSJacob Faibussowitsch         PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
2099bd66eb0SPeter Brune         if (linesearch->ops->vinorm) {
2109bd66eb0SPeter Brune           gnorm = fnorm;
2119566063dSJacob Faibussowitsch           PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
2129bd66eb0SPeter Brune         } else {
2139566063dSJacob Faibussowitsch           PetscCall(VecNorm(G, NORM_2, &gnorm));
2149bd66eb0SPeter Brune         }
215e3ce9940SStefano Zampini         g = 0.5 * PetscSqr(gnorm);
2168a430afbSPeter Brune       }
217c98378a5SBarry Smith       if (PetscIsInfOrNanReal(g)) {
2189566063dSJacob Faibussowitsch         PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
2199566063dSJacob Faibussowitsch         PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n"));
2203ba16761SJacob Faibussowitsch         PetscFunctionReturn(PETSC_SUCCESS);
221c98378a5SBarry Smith       }
22271b4ebd2SPeter Brune       if (monitor) {
2239566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
224bd4a8a71SBarry Smith         if (!objective) {
2259566063dSJacob Faibussowitsch           PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: gnorm after quadratic fit %14.12e\n", (double)gnorm));
2268a430afbSPeter Brune         } else {
2279566063dSJacob Faibussowitsch           PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: obj after quadratic fit %14.12e\n", (double)g));
2288a430afbSPeter Brune         }
2299566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
23071b4ebd2SPeter Brune       }
23132661dccSStefano Zampini     }
232e3ce9940SStefano Zampini     if (linesearch->order != SNES_LINESEARCH_ORDER_LINEAR && g <= f + lambda * alpha * initslope) { /* sufficient reduction */
23371b4ebd2SPeter Brune       if (monitor) {
2349566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
2359566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Quadratically determined step, lambda=%18.16e\n", (double)lambda));
2369566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
23771b4ebd2SPeter Brune       }
23871b4ebd2SPeter Brune     } else {
23971b4ebd2SPeter Brune       for (count = 0; count < max_its; count++) {
24071b4ebd2SPeter Brune         if (lambda <= minlambda) {
24171b4ebd2SPeter Brune           if (monitor) {
2429566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
24363a3b9bcSJacob Faibussowitsch             PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: unable to find good step length! After %" PetscInt_FMT " tries \n", count));
244bd4a8a71SBarry Smith             if (!objective) {
2459371c9d4SSatish Balay               PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope));
2468a430afbSPeter Brune             } else {
247e3ce9940SStefano Zampini               PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: obj(0)=%18.16e, obj=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope));
2488a430afbSPeter Brune             }
2499566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
25071b4ebd2SPeter Brune           }
2519566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
2523ba16761SJacob Faibussowitsch           PetscFunctionReturn(PETSC_SUCCESS);
25371b4ebd2SPeter Brune         }
254b000cd8dSPeter Brune         if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
255e3ce9940SStefano Zampini           /* Fit points with cubic */
256e3ce9940SStefano Zampini           t1 = g - f - lambda * initslope;
257e3ce9940SStefano Zampini           t2 = gprev - f - lambdaprev * initslope;
25871b4ebd2SPeter Brune           a  = (t1 / (lambda * lambda) - t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev);
25971b4ebd2SPeter Brune           b  = (-lambdaprev * t1 / (lambda * lambda) + lambda * t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev);
26071b4ebd2SPeter Brune           d  = b * b - 3 * a * initslope;
26171b4ebd2SPeter Brune           if (d < 0.0) d = 0.0;
262f5af7f23SKarl Rupp           if (a == 0.0) lambdatemp = -initslope / (2.0 * b);
263f5af7f23SKarl Rupp           else lambdatemp = (-b + PetscSqrtReal(d)) / (3.0 * a);
264b000cd8dSPeter Brune         } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
265e3ce9940SStefano Zampini           lambdatemp = -initslope * PetscSqr(lambda) / (2.0 * (g - f - lambda * initslope));
26632661dccSStefano Zampini         } else if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) { /* Just backtrack */
26732661dccSStefano Zampini           lambdatemp = .5 * lambda;
26832661dccSStefano Zampini         } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "Line search order %" PetscInt_FMT " for type bt", linesearch->order);
26971b4ebd2SPeter Brune         lambdaprev = lambda;
2708a430afbSPeter Brune         gprev      = g;
271e3ce9940SStefano Zampini 
272e3ce9940SStefano Zampini         lambda = PetscClipInterval(lambdatemp, .1 * lambda, .5 * lambda);
2739566063dSJacob Faibussowitsch         PetscCall(VecWAXPY(W, -lambda, Y, X));
2741baa6e33SBarry Smith         if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
275e71169deSBarry Smith         if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
27663a3b9bcSJacob Faibussowitsch           PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while looking for good step length! %" PetscInt_FMT " \n", count));
27748a46eb9SPierre Jolivet           if (!objective) PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", (double)fnorm, (double)gnorm, (double)ynorm, (double)lambda, (double)initslope));
2789566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
27971b4ebd2SPeter Brune           snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
2803ba16761SJacob Faibussowitsch           PetscFunctionReturn(PETSC_SUCCESS);
28171b4ebd2SPeter Brune         }
282bd4a8a71SBarry Smith         if (objective) {
2839566063dSJacob Faibussowitsch           PetscCall(SNESComputeObjective(snes, W, &g));
2848a430afbSPeter Brune         } else {
2859566063dSJacob Faibussowitsch           PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
2869bd66eb0SPeter Brune           if (linesearch->ops->vinorm) {
2879bd66eb0SPeter Brune             gnorm = fnorm;
2889566063dSJacob Faibussowitsch             PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
2899bd66eb0SPeter Brune           } else {
2909566063dSJacob Faibussowitsch             PetscCall(VecNorm(G, NORM_2, &gnorm));
2918a430afbSPeter Brune           }
292e3ce9940SStefano Zampini           g = 0.5 * PetscSqr(gnorm);
2939bd66eb0SPeter Brune         }
2944a2f8832SBarry Smith         if (PetscIsInfOrNanReal(g)) {
2959566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
2969566063dSJacob Faibussowitsch           PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n"));
2973ba16761SJacob Faibussowitsch           PetscFunctionReturn(PETSC_SUCCESS);
298c98378a5SBarry Smith         }
299e3ce9940SStefano Zampini         if (g <= f + lambda * alpha * initslope) { /* is reduction enough? */
30071b4ebd2SPeter Brune           if (monitor) {
3019566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
302bd4a8a71SBarry Smith             if (!objective) {
30332661dccSStefano Zampini               PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: %s step, current gnorm %14.12e lambda=%18.16e\n", ordStr[linesearch->order - 1], (double)gnorm, (double)lambda));
3049566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3058a430afbSPeter Brune             } else {
30632661dccSStefano Zampini               PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: %s step, obj %14.12e lambda=%18.16e\n", ordStr[linesearch->order - 1], (double)g, (double)lambda));
3079566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3088a430afbSPeter Brune             }
30971b4ebd2SPeter Brune           }
31071b4ebd2SPeter Brune           break;
311f5af7f23SKarl Rupp         } else if (monitor) {
3129566063dSJacob Faibussowitsch           PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
313bd4a8a71SBarry Smith           if (!objective) {
31432661dccSStefano Zampini             PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: %s step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n", ordStr[linesearch->order - 1], (double)gnorm, (double)lambda));
3159566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3168a430afbSPeter Brune           } else {
31732661dccSStefano Zampini             PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: %s step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n", ordStr[linesearch->order - 1], (double)g, (double)lambda));
3189566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3198a430afbSPeter Brune           }
32071b4ebd2SPeter Brune         }
32171b4ebd2SPeter Brune       }
32271b4ebd2SPeter Brune     }
32371b4ebd2SPeter Brune   }
32471b4ebd2SPeter Brune 
32571b4ebd2SPeter Brune   /* postcheck */
3260c1ac483SGlenn Hammond   /* update Y to lambda*Y so that W is consistent with  X - lambda*Y */
3279566063dSJacob Faibussowitsch   PetscCall(VecScale(Y, lambda));
3289566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
32971b4ebd2SPeter Brune   if (changed_y) {
3309566063dSJacob Faibussowitsch     PetscCall(VecWAXPY(W, -1.0, Y, X));
3311baa6e33SBarry Smith     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
33271b4ebd2SPeter Brune   }
333bd4a8a71SBarry Smith   if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
3349566063dSJacob Faibussowitsch     PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
3359bd66eb0SPeter Brune     if (linesearch->ops->vinorm) {
3369bd66eb0SPeter Brune       gnorm = fnorm;
3379566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
3389bd66eb0SPeter Brune     } else {
3399566063dSJacob Faibussowitsch       PetscCall(VecNorm(G, NORM_2, &gnorm));
3409bd66eb0SPeter Brune     }
3419566063dSJacob Faibussowitsch     PetscCall(VecNorm(Y, NORM_2, &ynorm));
342c98378a5SBarry Smith     if (PetscIsInfOrNanReal(gnorm)) {
3439566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
3449566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n"));
3453ba16761SJacob Faibussowitsch       PetscFunctionReturn(PETSC_SUCCESS);
346c98378a5SBarry Smith     }
347c427506bSPeter Brune   }
34871b4ebd2SPeter Brune 
34971b4ebd2SPeter Brune   /* copy the solution over */
3509566063dSJacob Faibussowitsch   PetscCall(VecCopy(W, X));
3519566063dSJacob Faibussowitsch   PetscCall(VecCopy(G, F));
3529566063dSJacob Faibussowitsch   PetscCall(VecNorm(X, NORM_2, &xnorm));
3539566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
3549566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm));
3553ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
35671b4ebd2SPeter Brune }
35771b4ebd2SPeter Brune 
35866976f2fSJacob Faibussowitsch static PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
359d71ae5a4SJacob Faibussowitsch {
3607f1410a3SPeter Brune   PetscBool          iascii;
361d7ca06c0SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
362075cc632SBarry Smith 
3637f1410a3SPeter Brune   PetscFunctionBegin;
3649566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
3657f1410a3SPeter Brune   if (iascii) {
366b000cd8dSPeter Brune     if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3679566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  interpolation: cubic\n"));
368b000cd8dSPeter Brune     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
3699566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  interpolation: quadratic\n"));
3707f1410a3SPeter Brune     }
3719566063dSJacob Faibussowitsch     PetscCall(PetscViewerASCIIPrintf(viewer, "  alpha=%e\n", (double)bt->alpha));
3727f1410a3SPeter Brune   }
3733ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
3747f1410a3SPeter Brune }
3757f1410a3SPeter Brune 
376d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
377d71ae5a4SJacob Faibussowitsch {
37871b4ebd2SPeter Brune   PetscFunctionBegin;
3799566063dSJacob Faibussowitsch   PetscCall(PetscFree(linesearch->data));
3803ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
38171b4ebd2SPeter Brune }
38271b4ebd2SPeter Brune 
383d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch, PetscOptionItems *PetscOptionsObject)
384d71ae5a4SJacob Faibussowitsch {
385eb23ba39SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
38671b4ebd2SPeter Brune 
3876e111a19SKarl Rupp   PetscFunctionBegin;
388d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject, "SNESLineSearch BT options");
3899566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL));
390d0609cedSBarry Smith   PetscOptionsHeadEnd();
3913ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
39271b4ebd2SPeter Brune }
39371b4ebd2SPeter Brune 
394954494b2SJed Brown /*MC
395*1d27aa22SBarry Smith    SNESLINESEARCHBT - Backtracking line search {cite}`dennis:83`.
396954494b2SJed Brown 
397db609ea7SPeter Brune    This line search finds the minimum of a polynomial fitting of the L2 norm of the
39832661dccSStefano Zampini    function or the objective function if it is provided with `SNESSetObjective()`.
39932661dccSStefano Zampini    If this fit does not satisfy the conditions for progress, the interval shrinks
400db609ea7SPeter Brune    and the fit is reattempted at most max_it times or until lambda is below minlambda.
401954494b2SJed Brown 
402cd7522eaSPeter Brune    Options Database Keys:
4033eec2f6aSPierre Jolivet +  -snes_linesearch_alpha <1e\-4>      - slope descent parameter
404db609ea7SPeter Brune .  -snes_linesearch_damping <1.0>      - initial step length
405eb23ba39SBarry Smith .  -snes_linesearch_maxstep <length>   - if the length the initial step is larger than this then the
406e1bc860dSBarry Smith                                          step is scaled back to be of this length at the beginning of the line search
407db609ea7SPeter Brune .  -snes_linesearch_max_it <40>        - maximum number of shrinking step
4083eec2f6aSPierre Jolivet .  -snes_linesearch_minlambda <1e\-12> - minimum step length allowed
40932661dccSStefano Zampini -  -snes_linesearch_order <1,2,3>      - order of the approximation. With order 1, it performs a simple backtracking without any curve fitting
410cd7522eaSPeter Brune 
411954494b2SJed Brown    Level: advanced
412954494b2SJed Brown 
413f6dfbefdSBarry Smith    Note:
414e55accfeSBarryFSmith    This line search will always produce a step that is less than or equal to, in length, the full step size.
415db609ea7SPeter Brune 
416420bcc1bSBarry Smith .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
417954494b2SJed Brown M*/
418d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
419d71ae5a4SJacob Faibussowitsch {
420f1c6b773SPeter Brune   SNESLineSearch_BT *bt;
42171b4ebd2SPeter Brune 
42271b4ebd2SPeter Brune   PetscFunctionBegin;
423f1c6b773SPeter Brune   linesearch->ops->apply          = SNESLineSearchApply_BT;
424f1c6b773SPeter Brune   linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
425f1c6b773SPeter Brune   linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
4260298fd71SBarry Smith   linesearch->ops->reset          = NULL;
4277f1410a3SPeter Brune   linesearch->ops->view           = SNESLineSearchView_BT;
4280298fd71SBarry Smith   linesearch->ops->setup          = NULL;
42971b4ebd2SPeter Brune 
4304dfa11a4SJacob Faibussowitsch   PetscCall(PetscNew(&bt));
431f5af7f23SKarl Rupp 
43271b4ebd2SPeter Brune   linesearch->data    = (void *)bt;
43371b4ebd2SPeter Brune   linesearch->max_its = 40;
434b000cd8dSPeter Brune   linesearch->order   = SNES_LINESEARCH_ORDER_CUBIC;
43571b4ebd2SPeter Brune   bt->alpha           = 1e-4;
4363ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
43771b4ebd2SPeter Brune }
438