xref: /petsc/src/snes/linesearch/impls/bt/linesearchbt.c (revision 48a46eb9bd028bec07ec0f396b1a3abb43f14558)
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 /*@
92f4102e2SPeter Brune    SNESLineSearchBTSetAlpha - Sets the descent parameter, alpha, in the BT linesearch 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 
17db781477SPatrick Sanan .seealso: `SNESLineSearchSetLambda()`, `SNESLineSearchGetTolerances()` `SNESLINESEARCHBT`
182f4102e2SPeter Brune @*/
199371c9d4SSatish Balay PetscErrorCode SNESLineSearchBTSetAlpha(SNESLineSearch linesearch, PetscReal alpha) {
20d7ca06c0SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
216e111a19SKarl Rupp 
222f4102e2SPeter Brune   PetscFunctionBegin;
232f4102e2SPeter Brune   PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 1);
242f4102e2SPeter Brune   bt->alpha = alpha;
252f4102e2SPeter Brune   PetscFunctionReturn(0);
262f4102e2SPeter Brune }
272f4102e2SPeter Brune 
282f4102e2SPeter Brune /*@
292f4102e2SPeter Brune    SNESLineSearchBTGetAlpha - Gets the descent parameter, alpha, in the BT linesearch variant.
302f4102e2SPeter Brune 
312f4102e2SPeter Brune    Input Parameters:
322f4102e2SPeter Brune .  linesearch - linesearch context
338e557f58SPeter Brune 
348e557f58SPeter Brune    Output Parameters:
352f4102e2SPeter Brune .  alpha - The descent parameter
362f4102e2SPeter Brune 
372f4102e2SPeter Brune    Level: intermediate
382f4102e2SPeter Brune 
39db781477SPatrick Sanan .seealso: `SNESLineSearchGetLambda()`, `SNESLineSearchGetTolerances()` `SNESLINESEARCHBT`
402f4102e2SPeter Brune @*/
419371c9d4SSatish Balay PetscErrorCode SNESLineSearchBTGetAlpha(SNESLineSearch linesearch, PetscReal *alpha) {
42d7ca06c0SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
436e111a19SKarl Rupp 
442f4102e2SPeter Brune   PetscFunctionBegin;
452f4102e2SPeter Brune   PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 1);
462f4102e2SPeter Brune   *alpha = bt->alpha;
472f4102e2SPeter Brune   PetscFunctionReturn(0);
482f4102e2SPeter Brune }
492f4102e2SPeter Brune 
509371c9d4SSatish Balay static PetscErrorCode SNESLineSearchApply_BT(SNESLineSearch linesearch) {
5171b4ebd2SPeter Brune   PetscBool          changed_y, changed_w;
5271b4ebd2SPeter Brune   Vec                X, F, Y, W, G;
5371b4ebd2SPeter Brune   SNES               snes;
548a430afbSPeter Brune   PetscReal          fnorm, xnorm, ynorm, gnorm;
553b288129SPeter Brune   PetscReal          lambda, lambdatemp, lambdaprev, minlambda, maxstep, initslope, alpha, stol;
56dfcd3828SPeter Brune   PetscReal          t1, t2, a, b, d;
578a430afbSPeter Brune   PetscReal          f;
588a430afbSPeter Brune   PetscReal          g, gprev;
5971b4ebd2SPeter Brune   PetscViewer        monitor;
6071b4ebd2SPeter Brune   PetscInt           max_its, count;
61d7ca06c0SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
6271b4ebd2SPeter Brune   Mat                jac;
63bd4a8a71SBarry Smith   PetscErrorCode (*objective)(SNES, Vec, PetscReal *, void *);
6471b4ebd2SPeter Brune 
6571b4ebd2SPeter Brune   PetscFunctionBegin;
669566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G));
679566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm));
689566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetLambda(linesearch, &lambda));
699566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
709566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor));
719566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetTolerances(linesearch, &minlambda, &maxstep, NULL, NULL, NULL, &max_its));
729566063dSJacob Faibussowitsch   PetscCall(SNESGetTolerances(snes, NULL, NULL, &stol, NULL, NULL));
739566063dSJacob Faibussowitsch   PetscCall(SNESGetObjective(snes, &objective, NULL));
7471b4ebd2SPeter Brune   alpha = bt->alpha;
7571b4ebd2SPeter Brune 
769566063dSJacob Faibussowitsch   PetscCall(SNESGetJacobian(snes, &jac, NULL, NULL, NULL));
7708401ef6SPierre Jolivet   PetscCheck(jac || objective, PetscObjectComm((PetscObject)linesearch), PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix");
78c69d1a72SBarry Smith 
799566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y));
809566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
8171b4ebd2SPeter Brune 
829566063dSJacob Faibussowitsch   PetscCall(VecNormBegin(Y, NORM_2, &ynorm));
839566063dSJacob Faibussowitsch   PetscCall(VecNormBegin(X, NORM_2, &xnorm));
849566063dSJacob Faibussowitsch   PetscCall(VecNormEnd(Y, NORM_2, &ynorm));
859566063dSJacob Faibussowitsch   PetscCall(VecNormEnd(X, NORM_2, &xnorm));
863b288129SPeter Brune 
8771b4ebd2SPeter Brune   if (ynorm == 0.0) {
8871b4ebd2SPeter Brune     if (monitor) {
899566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
909566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Initial direction and size is 0\n"));
919566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
9271b4ebd2SPeter Brune     }
939566063dSJacob Faibussowitsch     PetscCall(VecCopy(X, W));
949566063dSJacob Faibussowitsch     PetscCall(VecCopy(F, G));
959566063dSJacob Faibussowitsch     PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm));
969566063dSJacob Faibussowitsch     PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
9771b4ebd2SPeter Brune     PetscFunctionReturn(0);
9871b4ebd2SPeter Brune   }
9971b4ebd2SPeter Brune   if (ynorm > maxstep) { /* Step too big, so scale back */
10071b4ebd2SPeter Brune     if (monitor) {
1019566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
1029566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Scaling step by %14.12e old ynorm %14.12e\n", (double)(maxstep / ynorm), (double)ynorm));
1039566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
10471b4ebd2SPeter Brune     }
1059566063dSJacob Faibussowitsch     PetscCall(VecScale(Y, maxstep / (ynorm)));
10671b4ebd2SPeter Brune     ynorm = maxstep;
10771b4ebd2SPeter Brune   }
1088a430afbSPeter Brune 
1098a430afbSPeter Brune   /* if the SNES has an objective set, use that instead of the function value */
110bd4a8a71SBarry Smith   if (objective) {
1119566063dSJacob Faibussowitsch     PetscCall(SNESComputeObjective(snes, X, &f));
1128a430afbSPeter Brune   } else {
1138a430afbSPeter Brune     f = fnorm * fnorm;
1148a430afbSPeter Brune   }
1158a430afbSPeter Brune 
1168a430afbSPeter Brune   /* compute the initial slope */
117bd4a8a71SBarry Smith   if (objective) {
1188a430afbSPeter Brune     /* slope comes from the function (assumed to be the gradient of the objective */
1199566063dSJacob Faibussowitsch     PetscCall(VecDotRealPart(Y, F, &initslope));
1208a430afbSPeter Brune   } else {
1218a430afbSPeter Brune     /* slope comes from the normal equations */
1229566063dSJacob Faibussowitsch     PetscCall(MatMult(jac, Y, W));
1239566063dSJacob Faibussowitsch     PetscCall(VecDotRealPart(F, W, &initslope));
12471b4ebd2SPeter Brune     if (initslope > 0.0) initslope = -initslope;
12571b4ebd2SPeter Brune     if (initslope == 0.0) initslope = -1.0;
1268a430afbSPeter Brune   }
12771b4ebd2SPeter Brune 
128df019d78SBarry Smith   while (PETSC_TRUE) {
1299566063dSJacob Faibussowitsch     PetscCall(VecWAXPY(W, -lambda, Y, X));
1301baa6e33SBarry Smith     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
131e71169deSBarry Smith     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
1329566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while checking full step length!\n"));
13371b4ebd2SPeter Brune       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
1349566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
13571b4ebd2SPeter Brune       PetscFunctionReturn(0);
13671b4ebd2SPeter Brune     }
1378a430afbSPeter Brune 
138bd4a8a71SBarry Smith     if (objective) {
1399566063dSJacob Faibussowitsch       PetscCall(SNESComputeObjective(snes, W, &g));
1408a430afbSPeter Brune     } else {
1419566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
1429bd66eb0SPeter Brune       if (linesearch->ops->vinorm) {
1439bd66eb0SPeter Brune         gnorm = fnorm;
1449566063dSJacob Faibussowitsch         PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
1459bd66eb0SPeter Brune       } else {
1469566063dSJacob Faibussowitsch         PetscCall(VecNorm(G, NORM_2, &gnorm));
1478a430afbSPeter Brune       }
14879e7e81bSJed Brown       g = PetscSqr(gnorm);
1499bd66eb0SPeter Brune     }
1509566063dSJacob Faibussowitsch     PetscCall(SNESLineSearchMonitor(linesearch));
1519bd66eb0SPeter Brune 
152df019d78SBarry Smith     if (!PetscIsInfOrNanReal(g)) break;
153df019d78SBarry Smith     if (monitor) {
1549566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
1559566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n", (double)lambda));
1569566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
157df019d78SBarry Smith     }
1589371c9d4SSatish Balay     if (lambda <= minlambda) { SNESCheckFunctionNorm(snes, g); }
159df019d78SBarry Smith     lambda = .5 * lambda;
160df019d78SBarry Smith   }
161df019d78SBarry Smith 
162*48a46eb9SPierre Jolivet   if (!objective) PetscCall(PetscInfo(snes, "Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm));
1638a430afbSPeter Brune   if (.5 * g <= .5 * f + lambda * alpha * initslope) { /* Sufficient reduction or step tolerance convergence */
16471b4ebd2SPeter Brune     if (monitor) {
1659566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
166bd4a8a71SBarry Smith       if (!objective) {
1679566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm));
1688a430afbSPeter Brune       } else {
1699566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Using full step: obj %14.12e obj %14.12e\n", (double)f, (double)g));
1708a430afbSPeter Brune       }
1719566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
17271b4ebd2SPeter Brune     }
17371b4ebd2SPeter Brune   } else {
174c21ba15cSPeter Brune     /* Since the full step didn't work and the step is tiny, quit */
175c21ba15cSPeter Brune     if (stol * xnorm > ynorm) {
1769566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm));
1779566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
178c61ba1e2SPeter Brune       if (monitor) {
1799566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
18063a3b9bcSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n", (double)ynorm, (double)(stol * xnorm)));
1819566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
182c61ba1e2SPeter Brune       }
183c21ba15cSPeter Brune       PetscFunctionReturn(0);
184c21ba15cSPeter Brune     }
18571b4ebd2SPeter Brune     /* Fit points with quadratic */
1868a430afbSPeter Brune     lambdatemp = -initslope / (g - f - 2.0 * lambda * initslope);
18771b4ebd2SPeter Brune     lambdaprev = lambda;
1888a430afbSPeter Brune     gprev      = g;
18971b4ebd2SPeter Brune     if (lambdatemp > .5 * lambda) lambdatemp = .5 * lambda;
19071b4ebd2SPeter Brune     if (lambdatemp <= .1 * lambda) lambda = .1 * lambda;
19171b4ebd2SPeter Brune     else lambda = lambdatemp;
19271b4ebd2SPeter Brune 
1939566063dSJacob Faibussowitsch     PetscCall(VecWAXPY(W, -lambda, Y, X));
1941baa6e33SBarry Smith     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
195e71169deSBarry Smith     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
19663a3b9bcSJacob Faibussowitsch       PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while attempting quadratic backtracking! %" PetscInt_FMT " \n", snes->nfuncs));
19771b4ebd2SPeter Brune       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
1989566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
19971b4ebd2SPeter Brune       PetscFunctionReturn(0);
20071b4ebd2SPeter Brune     }
201bd4a8a71SBarry Smith     if (objective) {
2029566063dSJacob Faibussowitsch       PetscCall(SNESComputeObjective(snes, W, &g));
2038a430afbSPeter Brune     } else {
2049566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
2059bd66eb0SPeter Brune       if (linesearch->ops->vinorm) {
2069bd66eb0SPeter Brune         gnorm = fnorm;
2079566063dSJacob Faibussowitsch         PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
2089bd66eb0SPeter Brune       } else {
2099566063dSJacob Faibussowitsch         PetscCall(VecNorm(G, NORM_2, &gnorm));
2109bd66eb0SPeter Brune       }
211f1ab69e3SPeter Brune       g = PetscSqr(gnorm);
2128a430afbSPeter Brune     }
213c98378a5SBarry Smith     if (PetscIsInfOrNanReal(g)) {
2149566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
2159566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n"));
216c98378a5SBarry Smith       PetscFunctionReturn(0);
217c98378a5SBarry Smith     }
21871b4ebd2SPeter Brune     if (monitor) {
2199566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
220bd4a8a71SBarry Smith       if (!objective) {
2219566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: gnorm after quadratic fit %14.12e\n", (double)gnorm));
2228a430afbSPeter Brune       } else {
2239566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: obj after quadratic fit %14.12e\n", (double)g));
2248a430afbSPeter Brune       }
2259566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
22671b4ebd2SPeter Brune     }
2278a430afbSPeter Brune     if (.5 * g < .5 * f + lambda * alpha * initslope) { /* sufficient reduction */
22871b4ebd2SPeter Brune       if (monitor) {
2299566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
2309566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Quadratically determined step, lambda=%18.16e\n", (double)lambda));
2319566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
23271b4ebd2SPeter Brune       }
23371b4ebd2SPeter Brune     } else {
23471b4ebd2SPeter Brune       /* Fit points with cubic */
23571b4ebd2SPeter Brune       for (count = 0; count < max_its; count++) {
23671b4ebd2SPeter Brune         if (lambda <= minlambda) {
23771b4ebd2SPeter Brune           if (monitor) {
2389566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
23963a3b9bcSJacob Faibussowitsch             PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: unable to find good step length! After %" PetscInt_FMT " tries \n", count));
240bd4a8a71SBarry Smith             if (!objective) {
2419371c9d4SSatish 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));
2428a430afbSPeter Brune             } else {
2439371c9d4SSatish Balay               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));
2448a430afbSPeter Brune             }
2459566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
24671b4ebd2SPeter Brune           }
2479566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
24871b4ebd2SPeter Brune           PetscFunctionReturn(0);
24971b4ebd2SPeter Brune         }
250b000cd8dSPeter Brune         if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
2518a430afbSPeter Brune           t1 = .5 * (g - f) - lambda * initslope;
2528a430afbSPeter Brune           t2 = .5 * (gprev - f) - lambdaprev * initslope;
25371b4ebd2SPeter Brune           a  = (t1 / (lambda * lambda) - t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev);
25471b4ebd2SPeter Brune           b  = (-lambdaprev * t1 / (lambda * lambda) + lambda * t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev);
25571b4ebd2SPeter Brune           d  = b * b - 3 * a * initslope;
25671b4ebd2SPeter Brune           if (d < 0.0) d = 0.0;
257f5af7f23SKarl Rupp           if (a == 0.0) lambdatemp = -initslope / (2.0 * b);
258f5af7f23SKarl Rupp           else lambdatemp = (-b + PetscSqrtReal(d)) / (3.0 * a);
259f5af7f23SKarl Rupp 
260b000cd8dSPeter Brune         } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
2618a430afbSPeter Brune           lambdatemp = -initslope / (g - f - 2.0 * initslope);
262ce94432eSBarry Smith         } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt");
26371b4ebd2SPeter Brune         lambdaprev = lambda;
2648a430afbSPeter Brune         gprev      = g;
26571b4ebd2SPeter Brune         if (lambdatemp > .5 * lambda) lambdatemp = .5 * lambda;
26671b4ebd2SPeter Brune         if (lambdatemp <= .1 * lambda) lambda = .1 * lambda;
26771b4ebd2SPeter Brune         else lambda = lambdatemp;
2689566063dSJacob Faibussowitsch         PetscCall(VecWAXPY(W, -lambda, Y, X));
2691baa6e33SBarry Smith         if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
270e71169deSBarry Smith         if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
27163a3b9bcSJacob Faibussowitsch           PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while looking for good step length! %" PetscInt_FMT " \n", count));
272*48a46eb9SPierre 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));
2739566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
27471b4ebd2SPeter Brune           snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
27571b4ebd2SPeter Brune           PetscFunctionReturn(0);
27671b4ebd2SPeter Brune         }
277bd4a8a71SBarry Smith         if (objective) {
2789566063dSJacob Faibussowitsch           PetscCall(SNESComputeObjective(snes, W, &g));
2798a430afbSPeter Brune         } else {
2809566063dSJacob Faibussowitsch           PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
2819bd66eb0SPeter Brune           if (linesearch->ops->vinorm) {
2829bd66eb0SPeter Brune             gnorm = fnorm;
2839566063dSJacob Faibussowitsch             PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
2849bd66eb0SPeter Brune           } else {
2859566063dSJacob Faibussowitsch             PetscCall(VecNorm(G, NORM_2, &gnorm));
2868a430afbSPeter Brune           }
287f1ab69e3SPeter Brune           g = PetscSqr(gnorm);
2889bd66eb0SPeter Brune         }
2894a2f8832SBarry Smith         if (PetscIsInfOrNanReal(g)) {
2909566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
2919566063dSJacob Faibussowitsch           PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n"));
292c98378a5SBarry Smith           PetscFunctionReturn(0);
293c98378a5SBarry Smith         }
2948a430afbSPeter Brune         if (.5 * g < .5 * f + lambda * alpha * initslope) { /* is reduction enough? */
29571b4ebd2SPeter Brune           if (monitor) {
2969566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
297bd4a8a71SBarry Smith             if (!objective) {
298b000cd8dSPeter Brune               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
2999566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n", (double)gnorm, (double)lambda));
30071b4ebd2SPeter Brune               } else {
3019566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n", (double)gnorm, (double)lambda));
30271b4ebd2SPeter Brune               }
3039566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3048a430afbSPeter Brune             } else {
3058a430afbSPeter Brune               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3069566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n", (double)g, (double)lambda));
3078a430afbSPeter Brune               } else {
3089566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n", (double)g, (double)lambda));
3098a430afbSPeter Brune               }
3109566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3118a430afbSPeter Brune             }
31271b4ebd2SPeter Brune           }
31371b4ebd2SPeter Brune           break;
314f5af7f23SKarl Rupp         } else if (monitor) {
3159566063dSJacob Faibussowitsch           PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
316bd4a8a71SBarry Smith           if (!objective) {
317b000cd8dSPeter Brune             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3189566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n", (double)gnorm, (double)lambda));
31971b4ebd2SPeter Brune             } else {
3209566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n", (double)gnorm, (double)lambda));
32171b4ebd2SPeter Brune             }
3229566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3238a430afbSPeter Brune           } else {
3248a430afbSPeter Brune             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3259566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n", (double)g, (double)lambda));
3268a430afbSPeter Brune             } else {
3279566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n", (double)g, (double)lambda));
3288a430afbSPeter Brune             }
3299566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3308a430afbSPeter Brune           }
33171b4ebd2SPeter Brune         }
33271b4ebd2SPeter Brune       }
33371b4ebd2SPeter Brune     }
33471b4ebd2SPeter Brune   }
33571b4ebd2SPeter Brune 
33671b4ebd2SPeter Brune   /* postcheck */
3370c1ac483SGlenn Hammond   /* update Y to lambda*Y so that W is consistent with  X - lambda*Y */
3389566063dSJacob Faibussowitsch   PetscCall(VecScale(Y, lambda));
3399566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
34071b4ebd2SPeter Brune   if (changed_y) {
3419566063dSJacob Faibussowitsch     PetscCall(VecWAXPY(W, -1.0, Y, X));
3421baa6e33SBarry Smith     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
34371b4ebd2SPeter Brune   }
344bd4a8a71SBarry Smith   if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
3459566063dSJacob Faibussowitsch     PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
3469bd66eb0SPeter Brune     if (linesearch->ops->vinorm) {
3479bd66eb0SPeter Brune       gnorm = fnorm;
3489566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
3499bd66eb0SPeter Brune     } else {
3509566063dSJacob Faibussowitsch       PetscCall(VecNorm(G, NORM_2, &gnorm));
3519bd66eb0SPeter Brune     }
3529566063dSJacob Faibussowitsch     PetscCall(VecNorm(Y, NORM_2, &ynorm));
353c98378a5SBarry Smith     if (PetscIsInfOrNanReal(gnorm)) {
3549566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
3559566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n"));
356c98378a5SBarry Smith       PetscFunctionReturn(0);
357c98378a5SBarry Smith     }
358c427506bSPeter Brune   }
35971b4ebd2SPeter Brune 
36071b4ebd2SPeter Brune   /* copy the solution over */
3619566063dSJacob Faibussowitsch   PetscCall(VecCopy(W, X));
3629566063dSJacob Faibussowitsch   PetscCall(VecCopy(G, F));
3639566063dSJacob Faibussowitsch   PetscCall(VecNorm(X, NORM_2, &xnorm));
3649566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
3659566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm));
36671b4ebd2SPeter Brune   PetscFunctionReturn(0);
36771b4ebd2SPeter Brune }
36871b4ebd2SPeter Brune 
3699371c9d4SSatish Balay PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer) {
3707f1410a3SPeter Brune   PetscBool          iascii;
371d7ca06c0SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
372075cc632SBarry Smith 
3737f1410a3SPeter Brune   PetscFunctionBegin;
3749566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
3757f1410a3SPeter Brune   if (iascii) {
376b000cd8dSPeter Brune     if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3779566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  interpolation: cubic\n"));
378b000cd8dSPeter Brune     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
3799566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  interpolation: quadratic\n"));
3807f1410a3SPeter Brune     }
3819566063dSJacob Faibussowitsch     PetscCall(PetscViewerASCIIPrintf(viewer, "  alpha=%e\n", (double)bt->alpha));
3827f1410a3SPeter Brune   }
3837f1410a3SPeter Brune   PetscFunctionReturn(0);
3847f1410a3SPeter Brune }
3857f1410a3SPeter Brune 
3869371c9d4SSatish Balay static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch) {
38771b4ebd2SPeter Brune   PetscFunctionBegin;
3889566063dSJacob Faibussowitsch   PetscCall(PetscFree(linesearch->data));
38971b4ebd2SPeter Brune   PetscFunctionReturn(0);
39071b4ebd2SPeter Brune }
39171b4ebd2SPeter Brune 
3929371c9d4SSatish Balay static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch, PetscOptionItems *PetscOptionsObject) {
393eb23ba39SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
39471b4ebd2SPeter Brune 
3956e111a19SKarl Rupp   PetscFunctionBegin;
396d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject, "SNESLineSearch BT options");
3979566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL));
398d0609cedSBarry Smith   PetscOptionsHeadEnd();
39971b4ebd2SPeter Brune   PetscFunctionReturn(0);
40071b4ebd2SPeter Brune }
40171b4ebd2SPeter Brune 
402954494b2SJed Brown /*MC
403db609ea7SPeter Brune    SNESLINESEARCHBT - Backtracking line search.
404954494b2SJed Brown 
405db609ea7SPeter Brune    This line search finds the minimum of a polynomial fitting of the L2 norm of the
406eb23ba39SBarry Smith    function or the objective function if it is provided with SNESSetObjective(). If this fit does not satisfy the conditions for progress, the interval shrinks
407db609ea7SPeter Brune    and the fit is reattempted at most max_it times or until lambda is below minlambda.
408954494b2SJed Brown 
409cd7522eaSPeter Brune    Options Database Keys:
4103eec2f6aSPierre Jolivet +  -snes_linesearch_alpha <1e\-4> - slope descent parameter
411db609ea7SPeter Brune .  -snes_linesearch_damping <1.0> - initial step length
412eb23ba39SBarry Smith .  -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the
413e1bc860dSBarry Smith                                        step is scaled back to be of this length at the beginning of the line search
414db609ea7SPeter Brune .  -snes_linesearch_max_it <40> - maximum number of shrinking step
4153eec2f6aSPierre Jolivet .  -snes_linesearch_minlambda <1e\-12> - minimum step length allowed
416cd7522eaSPeter Brune -  -snes_linesearch_order <cubic,quadratic> - order of the approximation
417cd7522eaSPeter Brune 
418954494b2SJed Brown    Level: advanced
419954494b2SJed Brown 
420db609ea7SPeter Brune    Notes:
421db609ea7SPeter Brune    This line search is taken from "Numerical Methods for Unconstrained
422db609ea7SPeter Brune    Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
423e55accfeSBarryFSmith 
424e55accfeSBarryFSmith    This line search will always produce a step that is less than or equal to, in length, the full step size.
425db609ea7SPeter Brune 
426db781477SPatrick Sanan .seealso: `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
427954494b2SJed Brown M*/
4289371c9d4SSatish Balay PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch) {
429f1c6b773SPeter Brune   SNESLineSearch_BT *bt;
43071b4ebd2SPeter Brune 
43171b4ebd2SPeter Brune   PetscFunctionBegin;
432f1c6b773SPeter Brune   linesearch->ops->apply          = SNESLineSearchApply_BT;
433f1c6b773SPeter Brune   linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
434f1c6b773SPeter Brune   linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
4350298fd71SBarry Smith   linesearch->ops->reset          = NULL;
4367f1410a3SPeter Brune   linesearch->ops->view           = SNESLineSearchView_BT;
4370298fd71SBarry Smith   linesearch->ops->setup          = NULL;
43871b4ebd2SPeter Brune 
4399566063dSJacob Faibussowitsch   PetscCall(PetscNewLog(linesearch, &bt));
440f5af7f23SKarl Rupp 
44171b4ebd2SPeter Brune   linesearch->data    = (void *)bt;
44271b4ebd2SPeter Brune   linesearch->max_its = 40;
443b000cd8dSPeter Brune   linesearch->order   = SNES_LINESEARCH_ORDER_CUBIC;
44471b4ebd2SPeter Brune   bt->alpha           = 1e-4;
44571b4ebd2SPeter Brune   PetscFunctionReturn(0);
44671b4ebd2SPeter Brune }
447