xref: /petsc/src/snes/linesearch/impls/bt/linesearchbt.c (revision 89dfbbd54135f7d16b8da6a6c7b6a7ee61018d1f)
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     }
1096b095a85SStefano Zampini     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 */
3266b095a85SStefano Zampini   PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
3279566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
32871b4ebd2SPeter Brune   if (changed_y) {
3296b095a85SStefano Zampini     if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X));
3301baa6e33SBarry Smith     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
33171b4ebd2SPeter Brune   }
332bd4a8a71SBarry Smith   if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
3339566063dSJacob Faibussowitsch     PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
3349bd66eb0SPeter Brune     if (linesearch->ops->vinorm) {
3359bd66eb0SPeter Brune       gnorm = fnorm;
3369566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
3379bd66eb0SPeter Brune     } else {
3389566063dSJacob Faibussowitsch       PetscCall(VecNorm(G, NORM_2, &gnorm));
3399bd66eb0SPeter Brune     }
3409566063dSJacob Faibussowitsch     PetscCall(VecNorm(Y, NORM_2, &ynorm));
341c98378a5SBarry Smith     if (PetscIsInfOrNanReal(gnorm)) {
3429566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
3439566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n"));
3443ba16761SJacob Faibussowitsch       PetscFunctionReturn(PETSC_SUCCESS);
345c98378a5SBarry Smith     }
346c427506bSPeter Brune   }
34771b4ebd2SPeter Brune 
34871b4ebd2SPeter Brune   /* copy the solution over */
3499566063dSJacob Faibussowitsch   PetscCall(VecCopy(W, X));
3509566063dSJacob Faibussowitsch   PetscCall(VecCopy(G, F));
3519566063dSJacob Faibussowitsch   PetscCall(VecNorm(X, NORM_2, &xnorm));
3529566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm));
3533ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
35471b4ebd2SPeter Brune }
35571b4ebd2SPeter Brune 
35666976f2fSJacob Faibussowitsch static PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
357d71ae5a4SJacob Faibussowitsch {
3587f1410a3SPeter Brune   PetscBool          iascii;
359d7ca06c0SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
360075cc632SBarry Smith 
3617f1410a3SPeter Brune   PetscFunctionBegin;
3629566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
3637f1410a3SPeter Brune   if (iascii) {
364b000cd8dSPeter Brune     if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3659566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  interpolation: cubic\n"));
366b000cd8dSPeter Brune     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
3679566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  interpolation: quadratic\n"));
3687f1410a3SPeter Brune     }
3699566063dSJacob Faibussowitsch     PetscCall(PetscViewerASCIIPrintf(viewer, "  alpha=%e\n", (double)bt->alpha));
3707f1410a3SPeter Brune   }
3713ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
3727f1410a3SPeter Brune }
3737f1410a3SPeter Brune 
374d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
375d71ae5a4SJacob Faibussowitsch {
37671b4ebd2SPeter Brune   PetscFunctionBegin;
3779566063dSJacob Faibussowitsch   PetscCall(PetscFree(linesearch->data));
3783ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
37971b4ebd2SPeter Brune }
38071b4ebd2SPeter Brune 
381d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch, PetscOptionItems *PetscOptionsObject)
382d71ae5a4SJacob Faibussowitsch {
383eb23ba39SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
38471b4ebd2SPeter Brune 
3856e111a19SKarl Rupp   PetscFunctionBegin;
386d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject, "SNESLineSearch BT options");
3879566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL));
388d0609cedSBarry Smith   PetscOptionsHeadEnd();
3893ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
39071b4ebd2SPeter Brune }
39171b4ebd2SPeter Brune 
392954494b2SJed Brown /*MC
3931d27aa22SBarry Smith    SNESLINESEARCHBT - Backtracking line search {cite}`dennis:83`.
394954494b2SJed Brown 
395db609ea7SPeter Brune    This line search finds the minimum of a polynomial fitting of the L2 norm of the
39632661dccSStefano Zampini    function or the objective function if it is provided with `SNESSetObjective()`.
39732661dccSStefano Zampini    If this fit does not satisfy the conditions for progress, the interval shrinks
398*89dfbbd5SBarry Smith    and the fit is reattempted at most `max_it` times or until $\lambda$ is below `minlambda`.
399954494b2SJed Brown 
400cd7522eaSPeter Brune    Options Database Keys:
4013eec2f6aSPierre Jolivet +  -snes_linesearch_alpha <1e\-4>      - slope descent parameter
402*89dfbbd5SBarry Smith .  -snes_linesearch_damping <1.0>      - scaling of initial step length on entry to the line search
403eb23ba39SBarry Smith .  -snes_linesearch_maxstep <length>   - if the length the initial step is larger than this then the
404e1bc860dSBarry Smith                                          step is scaled back to be of this length at the beginning of the line search
405*89dfbbd5SBarry Smith .  -snes_linesearch_max_it <40>        - maximum number of shrinking steps
4063eec2f6aSPierre Jolivet .  -snes_linesearch_minlambda <1e\-12> - minimum step length allowed
40732661dccSStefano Zampini -  -snes_linesearch_order <1,2,3>      - order of the approximation. With order 1, it performs a simple backtracking without any curve fitting
408cd7522eaSPeter Brune 
409954494b2SJed Brown    Level: advanced
410954494b2SJed Brown 
411f6dfbefdSBarry Smith    Note:
412e55accfeSBarryFSmith    This line search will always produce a step that is less than or equal to, in length, the full step size.
413db609ea7SPeter Brune 
414420bcc1bSBarry Smith .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
415954494b2SJed Brown M*/
416d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
417d71ae5a4SJacob Faibussowitsch {
418f1c6b773SPeter Brune   SNESLineSearch_BT *bt;
41971b4ebd2SPeter Brune 
42071b4ebd2SPeter Brune   PetscFunctionBegin;
421f1c6b773SPeter Brune   linesearch->ops->apply          = SNESLineSearchApply_BT;
422f1c6b773SPeter Brune   linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
423f1c6b773SPeter Brune   linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
4240298fd71SBarry Smith   linesearch->ops->reset          = NULL;
4257f1410a3SPeter Brune   linesearch->ops->view           = SNESLineSearchView_BT;
4260298fd71SBarry Smith   linesearch->ops->setup          = NULL;
42771b4ebd2SPeter Brune 
4284dfa11a4SJacob Faibussowitsch   PetscCall(PetscNew(&bt));
429f5af7f23SKarl Rupp 
43071b4ebd2SPeter Brune   linesearch->data    = (void *)bt;
43171b4ebd2SPeter Brune   linesearch->max_its = 40;
432b000cd8dSPeter Brune   linesearch->order   = SNES_LINESEARCH_ORDER_CUBIC;
43371b4ebd2SPeter Brune   bt->alpha           = 1e-4;
4343ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
43571b4ebd2SPeter Brune }
436