xref: /petsc/src/snes/linesearch/impls/bt/linesearchbt.c (revision 3ba1676111f5c958fe6c2729b46ca4d523958bb3)
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;
26*3ba16761SJacob 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;
49*3ba16761SJacob 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;
66bd4a8a71SBarry Smith   PetscErrorCode (*objective)(SNES, Vec, PetscReal *, void *);
6771b4ebd2SPeter Brune 
6871b4ebd2SPeter Brune   PetscFunctionBegin;
699566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G));
709566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm));
719566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetLambda(linesearch, &lambda));
729566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
739566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor));
749566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetTolerances(linesearch, &minlambda, &maxstep, NULL, NULL, NULL, &max_its));
759566063dSJacob Faibussowitsch   PetscCall(SNESGetTolerances(snes, NULL, NULL, &stol, NULL, NULL));
769566063dSJacob Faibussowitsch   PetscCall(SNESGetObjective(snes, &objective, NULL));
7771b4ebd2SPeter Brune   alpha = bt->alpha;
7871b4ebd2SPeter Brune 
799566063dSJacob Faibussowitsch   PetscCall(SNESGetJacobian(snes, &jac, NULL, NULL, NULL));
8008401ef6SPierre Jolivet   PetscCheck(jac || objective, PetscObjectComm((PetscObject)linesearch), PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix");
81c69d1a72SBarry Smith 
829566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y));
839566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
8471b4ebd2SPeter Brune 
859566063dSJacob Faibussowitsch   PetscCall(VecNormBegin(Y, NORM_2, &ynorm));
869566063dSJacob Faibussowitsch   PetscCall(VecNormBegin(X, NORM_2, &xnorm));
879566063dSJacob Faibussowitsch   PetscCall(VecNormEnd(Y, NORM_2, &ynorm));
889566063dSJacob Faibussowitsch   PetscCall(VecNormEnd(X, NORM_2, &xnorm));
893b288129SPeter Brune 
9071b4ebd2SPeter Brune   if (ynorm == 0.0) {
9171b4ebd2SPeter Brune     if (monitor) {
929566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
939566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Initial direction and size is 0\n"));
949566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
9571b4ebd2SPeter Brune     }
969566063dSJacob Faibussowitsch     PetscCall(VecCopy(X, W));
979566063dSJacob Faibussowitsch     PetscCall(VecCopy(F, G));
989566063dSJacob Faibussowitsch     PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm));
999566063dSJacob Faibussowitsch     PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
100*3ba16761SJacob Faibussowitsch     PetscFunctionReturn(PETSC_SUCCESS);
10171b4ebd2SPeter Brune   }
10271b4ebd2SPeter Brune   if (ynorm > maxstep) { /* Step too big, so scale back */
10371b4ebd2SPeter Brune     if (monitor) {
1049566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
1059566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Scaling step by %14.12e old ynorm %14.12e\n", (double)(maxstep / ynorm), (double)ynorm));
1069566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
10771b4ebd2SPeter Brune     }
1089566063dSJacob Faibussowitsch     PetscCall(VecScale(Y, maxstep / (ynorm)));
10971b4ebd2SPeter Brune     ynorm = maxstep;
11071b4ebd2SPeter Brune   }
1118a430afbSPeter Brune 
1128a430afbSPeter Brune   /* if the SNES has an objective set, use that instead of the function value */
113bd4a8a71SBarry Smith   if (objective) {
1149566063dSJacob Faibussowitsch     PetscCall(SNESComputeObjective(snes, X, &f));
1158a430afbSPeter Brune   } else {
1168a430afbSPeter Brune     f = fnorm * fnorm;
1178a430afbSPeter Brune   }
1188a430afbSPeter Brune 
1198a430afbSPeter Brune   /* compute the initial slope */
120bd4a8a71SBarry Smith   if (objective) {
1218a430afbSPeter Brune     /* slope comes from the function (assumed to be the gradient of the objective */
1229566063dSJacob Faibussowitsch     PetscCall(VecDotRealPart(Y, F, &initslope));
1238a430afbSPeter Brune   } else {
1248a430afbSPeter Brune     /* slope comes from the normal equations */
1259566063dSJacob Faibussowitsch     PetscCall(MatMult(jac, Y, W));
1269566063dSJacob Faibussowitsch     PetscCall(VecDotRealPart(F, W, &initslope));
12771b4ebd2SPeter Brune     if (initslope > 0.0) initslope = -initslope;
12871b4ebd2SPeter Brune     if (initslope == 0.0) initslope = -1.0;
1298a430afbSPeter Brune   }
13071b4ebd2SPeter Brune 
131df019d78SBarry Smith   while (PETSC_TRUE) {
1329566063dSJacob Faibussowitsch     PetscCall(VecWAXPY(W, -lambda, Y, X));
1331baa6e33SBarry Smith     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
134e71169deSBarry Smith     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
1359566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while checking full step length!\n"));
13671b4ebd2SPeter Brune       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
1379566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
138*3ba16761SJacob Faibussowitsch       PetscFunctionReturn(PETSC_SUCCESS);
13971b4ebd2SPeter Brune     }
1408a430afbSPeter Brune 
141bd4a8a71SBarry Smith     if (objective) {
1429566063dSJacob Faibussowitsch       PetscCall(SNESComputeObjective(snes, W, &g));
1438a430afbSPeter Brune     } else {
1449566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
1459bd66eb0SPeter Brune       if (linesearch->ops->vinorm) {
1469bd66eb0SPeter Brune         gnorm = fnorm;
1479566063dSJacob Faibussowitsch         PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
1489bd66eb0SPeter Brune       } else {
1499566063dSJacob Faibussowitsch         PetscCall(VecNorm(G, NORM_2, &gnorm));
1508a430afbSPeter Brune       }
15179e7e81bSJed Brown       g = PetscSqr(gnorm);
1529bd66eb0SPeter Brune     }
1539566063dSJacob Faibussowitsch     PetscCall(SNESLineSearchMonitor(linesearch));
1549bd66eb0SPeter Brune 
155df019d78SBarry Smith     if (!PetscIsInfOrNanReal(g)) break;
156df019d78SBarry Smith     if (monitor) {
1579566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
1589566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n", (double)lambda));
1599566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
160df019d78SBarry Smith     }
161ad540459SPierre Jolivet     if (lambda <= minlambda) SNESCheckFunctionNorm(snes, g);
162df019d78SBarry Smith     lambda = .5 * lambda;
163df019d78SBarry Smith   }
164df019d78SBarry Smith 
16548a46eb9SPierre Jolivet   if (!objective) PetscCall(PetscInfo(snes, "Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm));
1668a430afbSPeter Brune   if (.5 * g <= .5 * f + lambda * alpha * initslope) { /* Sufficient reduction or step tolerance convergence */
16771b4ebd2SPeter Brune     if (monitor) {
1689566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
169bd4a8a71SBarry Smith       if (!objective) {
1709566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm));
1718a430afbSPeter Brune       } else {
1729566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Using full step: obj %14.12e obj %14.12e\n", (double)f, (double)g));
1738a430afbSPeter Brune       }
1749566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
17571b4ebd2SPeter Brune     }
17671b4ebd2SPeter Brune   } else {
177c21ba15cSPeter Brune     /* Since the full step didn't work and the step is tiny, quit */
178c21ba15cSPeter Brune     if (stol * xnorm > ynorm) {
1799566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm));
1809566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
181c61ba1e2SPeter Brune       if (monitor) {
1829566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
18363a3b9bcSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n", (double)ynorm, (double)(stol * xnorm)));
1849566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
185c61ba1e2SPeter Brune       }
186*3ba16761SJacob Faibussowitsch       PetscFunctionReturn(PETSC_SUCCESS);
187c21ba15cSPeter Brune     }
18871b4ebd2SPeter Brune     /* Fit points with quadratic */
1898a430afbSPeter Brune     lambdatemp = -initslope / (g - f - 2.0 * lambda * initslope);
19071b4ebd2SPeter Brune     lambdaprev = lambda;
1918a430afbSPeter Brune     gprev      = g;
19271b4ebd2SPeter Brune     if (lambdatemp > .5 * lambda) lambdatemp = .5 * lambda;
19371b4ebd2SPeter Brune     if (lambdatemp <= .1 * lambda) lambda = .1 * lambda;
19471b4ebd2SPeter Brune     else lambda = lambdatemp;
19571b4ebd2SPeter Brune 
1969566063dSJacob Faibussowitsch     PetscCall(VecWAXPY(W, -lambda, Y, X));
1971baa6e33SBarry Smith     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
198e71169deSBarry Smith     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
19963a3b9bcSJacob Faibussowitsch       PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while attempting quadratic backtracking! %" PetscInt_FMT " \n", snes->nfuncs));
20071b4ebd2SPeter Brune       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
2019566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
202*3ba16761SJacob Faibussowitsch       PetscFunctionReturn(PETSC_SUCCESS);
20371b4ebd2SPeter Brune     }
204bd4a8a71SBarry Smith     if (objective) {
2059566063dSJacob Faibussowitsch       PetscCall(SNESComputeObjective(snes, W, &g));
2068a430afbSPeter Brune     } else {
2079566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
2089bd66eb0SPeter Brune       if (linesearch->ops->vinorm) {
2099bd66eb0SPeter Brune         gnorm = fnorm;
2109566063dSJacob Faibussowitsch         PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
2119bd66eb0SPeter Brune       } else {
2129566063dSJacob Faibussowitsch         PetscCall(VecNorm(G, NORM_2, &gnorm));
2139bd66eb0SPeter Brune       }
214f1ab69e3SPeter Brune       g = PetscSqr(gnorm);
2158a430afbSPeter Brune     }
216c98378a5SBarry Smith     if (PetscIsInfOrNanReal(g)) {
2179566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
2189566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n"));
219*3ba16761SJacob Faibussowitsch       PetscFunctionReturn(PETSC_SUCCESS);
220c98378a5SBarry Smith     }
22171b4ebd2SPeter Brune     if (monitor) {
2229566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
223bd4a8a71SBarry Smith       if (!objective) {
2249566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: gnorm after quadratic fit %14.12e\n", (double)gnorm));
2258a430afbSPeter Brune       } else {
2269566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: obj after quadratic fit %14.12e\n", (double)g));
2278a430afbSPeter Brune       }
2289566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
22971b4ebd2SPeter Brune     }
2308a430afbSPeter Brune     if (.5 * g < .5 * f + lambda * alpha * initslope) { /* sufficient reduction */
23171b4ebd2SPeter Brune       if (monitor) {
2329566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
2339566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Quadratically determined step, lambda=%18.16e\n", (double)lambda));
2349566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
23571b4ebd2SPeter Brune       }
23671b4ebd2SPeter Brune     } else {
23771b4ebd2SPeter Brune       /* Fit points with cubic */
23871b4ebd2SPeter Brune       for (count = 0; count < max_its; count++) {
23971b4ebd2SPeter Brune         if (lambda <= minlambda) {
24071b4ebd2SPeter Brune           if (monitor) {
2419566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
24263a3b9bcSJacob Faibussowitsch             PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: unable to find good step length! After %" PetscInt_FMT " tries \n", count));
243bd4a8a71SBarry Smith             if (!objective) {
2449371c9d4SSatish 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));
2458a430afbSPeter Brune             } else {
2469371c9d4SSatish 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));
2478a430afbSPeter Brune             }
2489566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
24971b4ebd2SPeter Brune           }
2509566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
251*3ba16761SJacob Faibussowitsch           PetscFunctionReturn(PETSC_SUCCESS);
25271b4ebd2SPeter Brune         }
253b000cd8dSPeter Brune         if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
2548a430afbSPeter Brune           t1 = .5 * (g - f) - lambda * initslope;
2558a430afbSPeter Brune           t2 = .5 * (gprev - f) - lambdaprev * initslope;
25671b4ebd2SPeter Brune           a  = (t1 / (lambda * lambda) - t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev);
25771b4ebd2SPeter Brune           b  = (-lambdaprev * t1 / (lambda * lambda) + lambda * t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev);
25871b4ebd2SPeter Brune           d  = b * b - 3 * a * initslope;
25971b4ebd2SPeter Brune           if (d < 0.0) d = 0.0;
260f5af7f23SKarl Rupp           if (a == 0.0) lambdatemp = -initslope / (2.0 * b);
261f5af7f23SKarl Rupp           else lambdatemp = (-b + PetscSqrtReal(d)) / (3.0 * a);
262f5af7f23SKarl Rupp 
263b000cd8dSPeter Brune         } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
2648a430afbSPeter Brune           lambdatemp = -initslope / (g - f - 2.0 * initslope);
265ce94432eSBarry Smith         } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt");
26671b4ebd2SPeter Brune         lambdaprev = lambda;
2678a430afbSPeter Brune         gprev      = g;
26871b4ebd2SPeter Brune         if (lambdatemp > .5 * lambda) lambdatemp = .5 * lambda;
26971b4ebd2SPeter Brune         if (lambdatemp <= .1 * lambda) lambda = .1 * lambda;
27071b4ebd2SPeter Brune         else lambda = lambdatemp;
2719566063dSJacob Faibussowitsch         PetscCall(VecWAXPY(W, -lambda, Y, X));
2721baa6e33SBarry Smith         if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
273e71169deSBarry Smith         if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
27463a3b9bcSJacob Faibussowitsch           PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while looking for good step length! %" PetscInt_FMT " \n", count));
27548a46eb9SPierre 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));
2769566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
27771b4ebd2SPeter Brune           snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
278*3ba16761SJacob Faibussowitsch           PetscFunctionReturn(PETSC_SUCCESS);
27971b4ebd2SPeter Brune         }
280bd4a8a71SBarry Smith         if (objective) {
2819566063dSJacob Faibussowitsch           PetscCall(SNESComputeObjective(snes, W, &g));
2828a430afbSPeter Brune         } else {
2839566063dSJacob Faibussowitsch           PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
2849bd66eb0SPeter Brune           if (linesearch->ops->vinorm) {
2859bd66eb0SPeter Brune             gnorm = fnorm;
2869566063dSJacob Faibussowitsch             PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
2879bd66eb0SPeter Brune           } else {
2889566063dSJacob Faibussowitsch             PetscCall(VecNorm(G, NORM_2, &gnorm));
2898a430afbSPeter Brune           }
290f1ab69e3SPeter Brune           g = PetscSqr(gnorm);
2919bd66eb0SPeter Brune         }
2924a2f8832SBarry Smith         if (PetscIsInfOrNanReal(g)) {
2939566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
2949566063dSJacob Faibussowitsch           PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n"));
295*3ba16761SJacob Faibussowitsch           PetscFunctionReturn(PETSC_SUCCESS);
296c98378a5SBarry Smith         }
2978a430afbSPeter Brune         if (.5 * g < .5 * f + lambda * alpha * initslope) { /* is reduction enough? */
29871b4ebd2SPeter Brune           if (monitor) {
2999566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
300bd4a8a71SBarry Smith             if (!objective) {
301b000cd8dSPeter Brune               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3029566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n", (double)gnorm, (double)lambda));
30371b4ebd2SPeter Brune               } else {
3049566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n", (double)gnorm, (double)lambda));
30571b4ebd2SPeter Brune               }
3069566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3078a430afbSPeter Brune             } else {
3088a430afbSPeter Brune               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3099566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n", (double)g, (double)lambda));
3108a430afbSPeter Brune               } else {
3119566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n", (double)g, (double)lambda));
3128a430afbSPeter Brune               }
3139566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3148a430afbSPeter Brune             }
31571b4ebd2SPeter Brune           }
31671b4ebd2SPeter Brune           break;
317f5af7f23SKarl Rupp         } else if (monitor) {
3189566063dSJacob Faibussowitsch           PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
319bd4a8a71SBarry Smith           if (!objective) {
320b000cd8dSPeter Brune             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3219566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n", (double)gnorm, (double)lambda));
32271b4ebd2SPeter Brune             } else {
3239566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n", (double)gnorm, (double)lambda));
32471b4ebd2SPeter Brune             }
3259566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3268a430afbSPeter Brune           } else {
3278a430afbSPeter Brune             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3289566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n", (double)g, (double)lambda));
3298a430afbSPeter Brune             } else {
3309566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n", (double)g, (double)lambda));
3318a430afbSPeter Brune             }
3329566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3338a430afbSPeter Brune           }
33471b4ebd2SPeter Brune         }
33571b4ebd2SPeter Brune       }
33671b4ebd2SPeter Brune     }
33771b4ebd2SPeter Brune   }
33871b4ebd2SPeter Brune 
33971b4ebd2SPeter Brune   /* postcheck */
3400c1ac483SGlenn Hammond   /* update Y to lambda*Y so that W is consistent with  X - lambda*Y */
3419566063dSJacob Faibussowitsch   PetscCall(VecScale(Y, lambda));
3429566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
34371b4ebd2SPeter Brune   if (changed_y) {
3449566063dSJacob Faibussowitsch     PetscCall(VecWAXPY(W, -1.0, Y, X));
3451baa6e33SBarry Smith     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
34671b4ebd2SPeter Brune   }
347bd4a8a71SBarry Smith   if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
3489566063dSJacob Faibussowitsch     PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
3499bd66eb0SPeter Brune     if (linesearch->ops->vinorm) {
3509bd66eb0SPeter Brune       gnorm = fnorm;
3519566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
3529bd66eb0SPeter Brune     } else {
3539566063dSJacob Faibussowitsch       PetscCall(VecNorm(G, NORM_2, &gnorm));
3549bd66eb0SPeter Brune     }
3559566063dSJacob Faibussowitsch     PetscCall(VecNorm(Y, NORM_2, &ynorm));
356c98378a5SBarry Smith     if (PetscIsInfOrNanReal(gnorm)) {
3579566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
3589566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n"));
359*3ba16761SJacob Faibussowitsch       PetscFunctionReturn(PETSC_SUCCESS);
360c98378a5SBarry Smith     }
361c427506bSPeter Brune   }
36271b4ebd2SPeter Brune 
36371b4ebd2SPeter Brune   /* copy the solution over */
3649566063dSJacob Faibussowitsch   PetscCall(VecCopy(W, X));
3659566063dSJacob Faibussowitsch   PetscCall(VecCopy(G, F));
3669566063dSJacob Faibussowitsch   PetscCall(VecNorm(X, NORM_2, &xnorm));
3679566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
3689566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm));
369*3ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
37071b4ebd2SPeter Brune }
37171b4ebd2SPeter Brune 
372d71ae5a4SJacob Faibussowitsch PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
373d71ae5a4SJacob Faibussowitsch {
3747f1410a3SPeter Brune   PetscBool          iascii;
375d7ca06c0SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
376075cc632SBarry Smith 
3777f1410a3SPeter Brune   PetscFunctionBegin;
3789566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
3797f1410a3SPeter Brune   if (iascii) {
380b000cd8dSPeter Brune     if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3819566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  interpolation: cubic\n"));
382b000cd8dSPeter Brune     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
3839566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  interpolation: quadratic\n"));
3847f1410a3SPeter Brune     }
3859566063dSJacob Faibussowitsch     PetscCall(PetscViewerASCIIPrintf(viewer, "  alpha=%e\n", (double)bt->alpha));
3867f1410a3SPeter Brune   }
387*3ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
3887f1410a3SPeter Brune }
3897f1410a3SPeter Brune 
390d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
391d71ae5a4SJacob Faibussowitsch {
39271b4ebd2SPeter Brune   PetscFunctionBegin;
3939566063dSJacob Faibussowitsch   PetscCall(PetscFree(linesearch->data));
394*3ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
39571b4ebd2SPeter Brune }
39671b4ebd2SPeter Brune 
397d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch, PetscOptionItems *PetscOptionsObject)
398d71ae5a4SJacob Faibussowitsch {
399eb23ba39SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
40071b4ebd2SPeter Brune 
4016e111a19SKarl Rupp   PetscFunctionBegin;
402d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject, "SNESLineSearch BT options");
4039566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL));
404d0609cedSBarry Smith   PetscOptionsHeadEnd();
405*3ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
40671b4ebd2SPeter Brune }
40771b4ebd2SPeter Brune 
408954494b2SJed Brown /*MC
409db609ea7SPeter Brune    SNESLINESEARCHBT - Backtracking line search.
410954494b2SJed Brown 
411db609ea7SPeter Brune    This line search finds the minimum of a polynomial fitting of the L2 norm of the
412f6dfbefdSBarry 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
413db609ea7SPeter Brune    and the fit is reattempted at most max_it times or until lambda is below minlambda.
414954494b2SJed Brown 
415cd7522eaSPeter Brune    Options Database Keys:
4163eec2f6aSPierre Jolivet +  -snes_linesearch_alpha <1e\-4> - slope descent parameter
417db609ea7SPeter Brune .  -snes_linesearch_damping <1.0> - initial step length
418eb23ba39SBarry Smith .  -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the
419e1bc860dSBarry Smith                                        step is scaled back to be of this length at the beginning of the line search
420db609ea7SPeter Brune .  -snes_linesearch_max_it <40> - maximum number of shrinking step
4213eec2f6aSPierre Jolivet .  -snes_linesearch_minlambda <1e\-12> - minimum step length allowed
422cd7522eaSPeter Brune -  -snes_linesearch_order <cubic,quadratic> - order of the approximation
423cd7522eaSPeter Brune 
424954494b2SJed Brown    Level: advanced
425954494b2SJed Brown 
426f6dfbefdSBarry Smith    Note:
427e55accfeSBarryFSmith    This line search will always produce a step that is less than or equal to, in length, the full step size.
428db609ea7SPeter Brune 
429f6dfbefdSBarry Smith    Reference:
430f6dfbefdSBarry Smith .  - * - "Numerical Methods for Unconstrained Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
431f6dfbefdSBarry Smith 
432f6dfbefdSBarry Smith .seealso: `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
433954494b2SJed Brown M*/
434d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
435d71ae5a4SJacob Faibussowitsch {
436f1c6b773SPeter Brune   SNESLineSearch_BT *bt;
43771b4ebd2SPeter Brune 
43871b4ebd2SPeter Brune   PetscFunctionBegin;
439f1c6b773SPeter Brune   linesearch->ops->apply          = SNESLineSearchApply_BT;
440f1c6b773SPeter Brune   linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
441f1c6b773SPeter Brune   linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
4420298fd71SBarry Smith   linesearch->ops->reset          = NULL;
4437f1410a3SPeter Brune   linesearch->ops->view           = SNESLineSearchView_BT;
4440298fd71SBarry Smith   linesearch->ops->setup          = NULL;
44571b4ebd2SPeter Brune 
4464dfa11a4SJacob Faibussowitsch   PetscCall(PetscNew(&bt));
447f5af7f23SKarl Rupp 
44871b4ebd2SPeter Brune   linesearch->data    = (void *)bt;
44971b4ebd2SPeter Brune   linesearch->max_its = 40;
450b000cd8dSPeter Brune   linesearch->order   = SNES_LINESEARCH_ORDER_CUBIC;
45171b4ebd2SPeter Brune   bt->alpha           = 1e-4;
452*3ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
45371b4ebd2SPeter Brune }
454