xref: /petsc/src/snes/linesearch/impls/bt/linesearchbt.c (revision a68bbae58a07f2fb515cab24a67de1159d72e8a2)
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 /*@
9f6dfbefdSBarry Smith    SNESLineSearchBTSetAlpha - Sets the descent parameter, alpha, in the `SNESLINESEARCHBT` 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 
17f6dfbefdSBarry Smith .seealso: `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 /*@
30f6dfbefdSBarry Smith    SNESLineSearchBTGetAlpha - Gets the descent parameter, alpha, in the `SNESLINESEARCHBT` variant.
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 
40f6dfbefdSBarry Smith .seealso: `SNESLineSearchGetLambda()`, `SNESLineSearchGetTolerances()` `SNESLINESEARCHBT`, `SNESLineSearchBTGetAlpha()`
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));
71*a68bbae5SBarry 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 {
11732661dccSStefano Zampini     f = 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));
12432661dccSStefano Zampini     initslope *= 2.0; /* all the fitting assumes f = ||grad||^2, not f = 0.5 ||grad||^2 */
1258a430afbSPeter Brune   } else {
1268a430afbSPeter Brune     /* slope comes from the normal equations */
1279566063dSJacob Faibussowitsch     PetscCall(MatMult(jac, Y, W));
1289566063dSJacob Faibussowitsch     PetscCall(VecDotRealPart(F, W, &initslope));
12971b4ebd2SPeter Brune     if (initslope > 0.0) initslope = -initslope;
13071b4ebd2SPeter Brune     if (initslope == 0.0) initslope = -1.0;
1318a430afbSPeter Brune   }
13271b4ebd2SPeter Brune 
133df019d78SBarry Smith   while (PETSC_TRUE) {
1349566063dSJacob Faibussowitsch     PetscCall(VecWAXPY(W, -lambda, Y, X));
1351baa6e33SBarry Smith     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
136e71169deSBarry Smith     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
1379566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while checking full step length!\n"));
13871b4ebd2SPeter Brune       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
1399566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
1403ba16761SJacob Faibussowitsch       PetscFunctionReturn(PETSC_SUCCESS);
14171b4ebd2SPeter Brune     }
1428a430afbSPeter Brune 
143bd4a8a71SBarry Smith     if (objective) {
1449566063dSJacob Faibussowitsch       PetscCall(SNESComputeObjective(snes, W, &g));
1458a430afbSPeter Brune     } else {
1469566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
1479bd66eb0SPeter Brune       if (linesearch->ops->vinorm) {
1489bd66eb0SPeter Brune         gnorm = fnorm;
1499566063dSJacob Faibussowitsch         PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
1509bd66eb0SPeter Brune       } else {
1519566063dSJacob Faibussowitsch         PetscCall(VecNorm(G, NORM_2, &gnorm));
1528a430afbSPeter Brune       }
15379e7e81bSJed Brown       g = PetscSqr(gnorm);
1549bd66eb0SPeter Brune     }
1559566063dSJacob Faibussowitsch     PetscCall(SNESLineSearchMonitor(linesearch));
1569bd66eb0SPeter Brune 
157df019d78SBarry Smith     if (!PetscIsInfOrNanReal(g)) break;
158df019d78SBarry Smith     if (monitor) {
1599566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
1609566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n", (double)lambda));
1619566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
162df019d78SBarry Smith     }
163ad540459SPierre Jolivet     if (lambda <= minlambda) SNESCheckFunctionNorm(snes, g);
16432661dccSStefano Zampini     lambda *= .5;
165df019d78SBarry Smith   }
166df019d78SBarry Smith 
16748a46eb9SPierre Jolivet   if (!objective) PetscCall(PetscInfo(snes, "Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm));
1688a430afbSPeter Brune   if (.5 * g <= .5 * f + lambda * alpha * initslope) { /* Sufficient reduction or step tolerance convergence */
16971b4ebd2SPeter Brune     if (monitor) {
1709566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
171bd4a8a71SBarry Smith       if (!objective) {
1729566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm));
1738a430afbSPeter Brune       } else {
1749566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Using full step: obj %14.12e obj %14.12e\n", (double)f, (double)g));
1758a430afbSPeter Brune       }
1769566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
17771b4ebd2SPeter Brune     }
17871b4ebd2SPeter Brune   } else {
179c21ba15cSPeter Brune     if (stol * xnorm > ynorm) {
180*a68bbae5SBarry Smith       /* Since the full step didn't give sufficient decrease and the step is tiny, exit */
1819566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm));
1829566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
183c61ba1e2SPeter Brune       if (monitor) {
1849566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
18563a3b9bcSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n", (double)ynorm, (double)(stol * xnorm)));
1869566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
187c61ba1e2SPeter Brune       }
1883ba16761SJacob Faibussowitsch       PetscFunctionReturn(PETSC_SUCCESS);
189c21ba15cSPeter Brune     }
19032661dccSStefano Zampini     /* Here to avoid -Wmaybe-uninitiliazed warnings */
19171b4ebd2SPeter Brune     lambdaprev = lambda;
1928a430afbSPeter Brune     gprev      = g;
19332661dccSStefano Zampini     if (linesearch->order != SNES_LINESEARCH_ORDER_LINEAR) {
19432661dccSStefano Zampini       /* Fit points with quadratic */
19532661dccSStefano Zampini       lambdatemp = -initslope / (g - f - 2.0 * lambda * initslope);
19671b4ebd2SPeter Brune       if (lambdatemp > .5 * lambda) lambdatemp = .5 * lambda;
19771b4ebd2SPeter Brune       if (lambdatemp <= .1 * lambda) lambda = .1 * lambda;
19871b4ebd2SPeter Brune       else lambda = lambdatemp;
19971b4ebd2SPeter Brune 
2009566063dSJacob Faibussowitsch       PetscCall(VecWAXPY(W, -lambda, Y, X));
2011baa6e33SBarry Smith       if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
202e71169deSBarry Smith       if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
20363a3b9bcSJacob Faibussowitsch         PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while attempting quadratic backtracking! %" PetscInt_FMT " \n", snes->nfuncs));
20471b4ebd2SPeter Brune         snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
2059566063dSJacob Faibussowitsch         PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
2063ba16761SJacob Faibussowitsch         PetscFunctionReturn(PETSC_SUCCESS);
20771b4ebd2SPeter Brune       }
208bd4a8a71SBarry Smith       if (objective) {
2099566063dSJacob Faibussowitsch         PetscCall(SNESComputeObjective(snes, W, &g));
2108a430afbSPeter Brune       } else {
2119566063dSJacob Faibussowitsch         PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
2129bd66eb0SPeter Brune         if (linesearch->ops->vinorm) {
2139bd66eb0SPeter Brune           gnorm = fnorm;
2149566063dSJacob Faibussowitsch           PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
2159bd66eb0SPeter Brune         } else {
2169566063dSJacob Faibussowitsch           PetscCall(VecNorm(G, NORM_2, &gnorm));
2179bd66eb0SPeter Brune         }
218f1ab69e3SPeter Brune         g = PetscSqr(gnorm);
2198a430afbSPeter Brune       }
220c98378a5SBarry Smith       if (PetscIsInfOrNanReal(g)) {
2219566063dSJacob Faibussowitsch         PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
2229566063dSJacob Faibussowitsch         PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n"));
2233ba16761SJacob Faibussowitsch         PetscFunctionReturn(PETSC_SUCCESS);
224c98378a5SBarry Smith       }
22571b4ebd2SPeter Brune       if (monitor) {
2269566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
227bd4a8a71SBarry Smith         if (!objective) {
2289566063dSJacob Faibussowitsch           PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: gnorm after quadratic fit %14.12e\n", (double)gnorm));
2298a430afbSPeter Brune         } else {
2309566063dSJacob Faibussowitsch           PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: obj after quadratic fit %14.12e\n", (double)g));
2318a430afbSPeter Brune         }
2329566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
23371b4ebd2SPeter Brune       }
23432661dccSStefano Zampini     }
23532661dccSStefano Zampini     if (linesearch->order != SNES_LINESEARCH_ORDER_LINEAR && .5 * g < .5 * f + lambda * alpha * initslope) { /* sufficient reduction */
23671b4ebd2SPeter Brune       if (monitor) {
2379566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
2389566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Quadratically determined step, lambda=%18.16e\n", (double)lambda));
2399566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
24071b4ebd2SPeter Brune       }
24171b4ebd2SPeter Brune     } else {
24271b4ebd2SPeter Brune       /* Fit points with cubic */
24371b4ebd2SPeter Brune       for (count = 0; count < max_its; count++) {
24471b4ebd2SPeter Brune         if (lambda <= minlambda) {
24571b4ebd2SPeter Brune           if (monitor) {
2469566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
24763a3b9bcSJacob Faibussowitsch             PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: unable to find good step length! After %" PetscInt_FMT " tries \n", count));
248bd4a8a71SBarry Smith             if (!objective) {
2499371c9d4SSatish 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));
2508a430afbSPeter Brune             } else {
25132661dccSStefano 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 / 2.0));
2528a430afbSPeter Brune             }
2539566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
25471b4ebd2SPeter Brune           }
2559566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
2563ba16761SJacob Faibussowitsch           PetscFunctionReturn(PETSC_SUCCESS);
25771b4ebd2SPeter Brune         }
258b000cd8dSPeter Brune         if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
2598a430afbSPeter Brune           t1 = .5 * (g - f) - lambda * initslope;
2608a430afbSPeter Brune           t2 = .5 * (gprev - f) - lambdaprev * initslope;
26171b4ebd2SPeter Brune           a  = (t1 / (lambda * lambda) - t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev);
26271b4ebd2SPeter Brune           b  = (-lambdaprev * t1 / (lambda * lambda) + lambda * t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev);
26371b4ebd2SPeter Brune           d  = b * b - 3 * a * initslope;
26471b4ebd2SPeter Brune           if (d < 0.0) d = 0.0;
265f5af7f23SKarl Rupp           if (a == 0.0) lambdatemp = -initslope / (2.0 * b);
266f5af7f23SKarl Rupp           else lambdatemp = (-b + PetscSqrtReal(d)) / (3.0 * a);
267b000cd8dSPeter Brune         } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
2688a430afbSPeter Brune           lambdatemp = -initslope / (g - f - 2.0 * initslope);
26932661dccSStefano Zampini         } else if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) { /* Just backtrack */
27032661dccSStefano Zampini           lambdatemp = .5 * lambda;
27132661dccSStefano Zampini         } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "Line search order %" PetscInt_FMT " for type bt", linesearch->order);
27271b4ebd2SPeter Brune         lambdaprev = lambda;
2738a430afbSPeter Brune         gprev      = g;
27471b4ebd2SPeter Brune         if (lambdatemp > .5 * lambda) lambdatemp = .5 * lambda;
27571b4ebd2SPeter Brune         if (lambdatemp <= .1 * lambda) lambda = .1 * lambda;
27671b4ebd2SPeter Brune         else lambda = lambdatemp;
2779566063dSJacob Faibussowitsch         PetscCall(VecWAXPY(W, -lambda, Y, X));
2781baa6e33SBarry Smith         if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
279e71169deSBarry Smith         if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
28063a3b9bcSJacob Faibussowitsch           PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while looking for good step length! %" PetscInt_FMT " \n", count));
28148a46eb9SPierre 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));
2829566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
28371b4ebd2SPeter Brune           snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
2843ba16761SJacob Faibussowitsch           PetscFunctionReturn(PETSC_SUCCESS);
28571b4ebd2SPeter Brune         }
286bd4a8a71SBarry Smith         if (objective) {
2879566063dSJacob Faibussowitsch           PetscCall(SNESComputeObjective(snes, W, &g));
2888a430afbSPeter Brune         } else {
2899566063dSJacob Faibussowitsch           PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
2909bd66eb0SPeter Brune           if (linesearch->ops->vinorm) {
2919bd66eb0SPeter Brune             gnorm = fnorm;
2929566063dSJacob Faibussowitsch             PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
2939bd66eb0SPeter Brune           } else {
2949566063dSJacob Faibussowitsch             PetscCall(VecNorm(G, NORM_2, &gnorm));
2958a430afbSPeter Brune           }
296f1ab69e3SPeter Brune           g = PetscSqr(gnorm);
2979bd66eb0SPeter Brune         }
2984a2f8832SBarry Smith         if (PetscIsInfOrNanReal(g)) {
2999566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
3009566063dSJacob Faibussowitsch           PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n"));
3013ba16761SJacob Faibussowitsch           PetscFunctionReturn(PETSC_SUCCESS);
302c98378a5SBarry Smith         }
3038a430afbSPeter Brune         if (.5 * g < .5 * f + lambda * alpha * initslope) { /* is reduction enough? */
30471b4ebd2SPeter Brune           if (monitor) {
3059566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
306bd4a8a71SBarry Smith             if (!objective) {
30732661dccSStefano Zampini               PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: %s step, current gnorm %14.12e lambda=%18.16e\n", ordStr[linesearch->order - 1], (double)gnorm, (double)lambda));
3089566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3098a430afbSPeter Brune             } else {
31032661dccSStefano Zampini               PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: %s step, obj %14.12e lambda=%18.16e\n", ordStr[linesearch->order - 1], (double)g, (double)lambda));
3119566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3128a430afbSPeter Brune             }
31371b4ebd2SPeter Brune           }
31471b4ebd2SPeter Brune           break;
315f5af7f23SKarl Rupp         } else if (monitor) {
3169566063dSJacob Faibussowitsch           PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
317bd4a8a71SBarry Smith           if (!objective) {
31832661dccSStefano 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));
3199566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3208a430afbSPeter Brune           } else {
32132661dccSStefano 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));
3229566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3238a430afbSPeter Brune           }
32471b4ebd2SPeter Brune         }
32571b4ebd2SPeter Brune       }
32671b4ebd2SPeter Brune     }
32771b4ebd2SPeter Brune   }
32871b4ebd2SPeter Brune 
32971b4ebd2SPeter Brune   /* postcheck */
3300c1ac483SGlenn Hammond   /* update Y to lambda*Y so that W is consistent with  X - lambda*Y */
3319566063dSJacob Faibussowitsch   PetscCall(VecScale(Y, lambda));
3329566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
33371b4ebd2SPeter Brune   if (changed_y) {
3349566063dSJacob Faibussowitsch     PetscCall(VecWAXPY(W, -1.0, Y, X));
3351baa6e33SBarry Smith     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
33671b4ebd2SPeter Brune   }
337bd4a8a71SBarry Smith   if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
3389566063dSJacob Faibussowitsch     PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
3399bd66eb0SPeter Brune     if (linesearch->ops->vinorm) {
3409bd66eb0SPeter Brune       gnorm = fnorm;
3419566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
3429bd66eb0SPeter Brune     } else {
3439566063dSJacob Faibussowitsch       PetscCall(VecNorm(G, NORM_2, &gnorm));
3449bd66eb0SPeter Brune     }
3459566063dSJacob Faibussowitsch     PetscCall(VecNorm(Y, NORM_2, &ynorm));
346c98378a5SBarry Smith     if (PetscIsInfOrNanReal(gnorm)) {
3479566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
3489566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n"));
3493ba16761SJacob Faibussowitsch       PetscFunctionReturn(PETSC_SUCCESS);
350c98378a5SBarry Smith     }
351c427506bSPeter Brune   }
35271b4ebd2SPeter Brune 
35371b4ebd2SPeter Brune   /* copy the solution over */
3549566063dSJacob Faibussowitsch   PetscCall(VecCopy(W, X));
3559566063dSJacob Faibussowitsch   PetscCall(VecCopy(G, F));
3569566063dSJacob Faibussowitsch   PetscCall(VecNorm(X, NORM_2, &xnorm));
3579566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
3589566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm));
3593ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
36071b4ebd2SPeter Brune }
36171b4ebd2SPeter Brune 
362d71ae5a4SJacob Faibussowitsch PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
363d71ae5a4SJacob Faibussowitsch {
3647f1410a3SPeter Brune   PetscBool          iascii;
365d7ca06c0SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
366075cc632SBarry Smith 
3677f1410a3SPeter Brune   PetscFunctionBegin;
3689566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
3697f1410a3SPeter Brune   if (iascii) {
370b000cd8dSPeter Brune     if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3719566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  interpolation: cubic\n"));
372b000cd8dSPeter Brune     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
3739566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  interpolation: quadratic\n"));
3747f1410a3SPeter Brune     }
3759566063dSJacob Faibussowitsch     PetscCall(PetscViewerASCIIPrintf(viewer, "  alpha=%e\n", (double)bt->alpha));
3767f1410a3SPeter Brune   }
3773ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
3787f1410a3SPeter Brune }
3797f1410a3SPeter Brune 
380d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
381d71ae5a4SJacob Faibussowitsch {
38271b4ebd2SPeter Brune   PetscFunctionBegin;
3839566063dSJacob Faibussowitsch   PetscCall(PetscFree(linesearch->data));
3843ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
38571b4ebd2SPeter Brune }
38671b4ebd2SPeter Brune 
387d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch, PetscOptionItems *PetscOptionsObject)
388d71ae5a4SJacob Faibussowitsch {
389eb23ba39SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
39071b4ebd2SPeter Brune 
3916e111a19SKarl Rupp   PetscFunctionBegin;
392d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject, "SNESLineSearch BT options");
3939566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL));
394d0609cedSBarry Smith   PetscOptionsHeadEnd();
3953ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
39671b4ebd2SPeter Brune }
39771b4ebd2SPeter Brune 
398954494b2SJed Brown /*MC
399db609ea7SPeter Brune    SNESLINESEARCHBT - Backtracking line search.
400954494b2SJed Brown 
401db609ea7SPeter Brune    This line search finds the minimum of a polynomial fitting of the L2 norm of the
40232661dccSStefano Zampini    function or the objective function if it is provided with `SNESSetObjective()`.
40332661dccSStefano Zampini    If this fit does not satisfy the conditions for progress, the interval shrinks
404db609ea7SPeter Brune    and the fit is reattempted at most max_it times or until lambda is below minlambda.
405954494b2SJed Brown 
406cd7522eaSPeter Brune    Options Database Keys:
4073eec2f6aSPierre Jolivet +  -snes_linesearch_alpha <1e\-4> - slope descent parameter
408db609ea7SPeter Brune .  -snes_linesearch_damping <1.0> - initial step length
409eb23ba39SBarry Smith .  -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the
410e1bc860dSBarry Smith                                        step is scaled back to be of this length at the beginning of the line search
411db609ea7SPeter Brune .  -snes_linesearch_max_it <40> - maximum number of shrinking step
4123eec2f6aSPierre Jolivet .  -snes_linesearch_minlambda <1e\-12> - minimum step length allowed
41332661dccSStefano Zampini -  -snes_linesearch_order <1,2,3> - order of the approximation. With order 1, it performs a simple backtracking without any curve fitting
414cd7522eaSPeter Brune 
415954494b2SJed Brown    Level: advanced
416954494b2SJed Brown 
417f6dfbefdSBarry Smith    Note:
418e55accfeSBarryFSmith    This line search will always produce a step that is less than or equal to, in length, the full step size.
419db609ea7SPeter Brune 
420f6dfbefdSBarry Smith    Reference:
421f6dfbefdSBarry Smith .  - * - "Numerical Methods for Unconstrained Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
422f6dfbefdSBarry Smith 
423f6dfbefdSBarry Smith .seealso: `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
424954494b2SJed Brown M*/
425d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
426d71ae5a4SJacob Faibussowitsch {
427f1c6b773SPeter Brune   SNESLineSearch_BT *bt;
42871b4ebd2SPeter Brune 
42971b4ebd2SPeter Brune   PetscFunctionBegin;
430f1c6b773SPeter Brune   linesearch->ops->apply          = SNESLineSearchApply_BT;
431f1c6b773SPeter Brune   linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
432f1c6b773SPeter Brune   linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
4330298fd71SBarry Smith   linesearch->ops->reset          = NULL;
4347f1410a3SPeter Brune   linesearch->ops->view           = SNESLineSearchView_BT;
4350298fd71SBarry Smith   linesearch->ops->setup          = NULL;
43671b4ebd2SPeter Brune 
4374dfa11a4SJacob Faibussowitsch   PetscCall(PetscNew(&bt));
438f5af7f23SKarl Rupp 
43971b4ebd2SPeter Brune   linesearch->data    = (void *)bt;
44071b4ebd2SPeter Brune   linesearch->max_its = 40;
441b000cd8dSPeter Brune   linesearch->order   = SNES_LINESEARCH_ORDER_CUBIC;
44271b4ebd2SPeter Brune   bt->alpha           = 1e-4;
4433ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
44471b4ebd2SPeter Brune }
445