xref: /petsc/src/snes/linesearch/impls/bt/linesearchbt.c (revision bd89dbf26d8a5efecb980364933175da61864cd7)
1af0996ceSBarry Smith #include <petsc/private/linesearchimpl.h> /*I  "petscsnes.h"  I*/
2af0996ceSBarry Smith #include <petsc/private/snesimpl.h>
371b4ebd2SPeter Brune 
471b4ebd2SPeter Brune typedef struct {
571b4ebd2SPeter Brune   PetscReal alpha; /* sufficient decrease parameter */
6f1c6b773SPeter Brune } SNESLineSearch_BT;
771b4ebd2SPeter Brune 
82f4102e2SPeter Brune /*@
9420bcc1bSBarry Smith   SNESLineSearchBTSetAlpha - Sets the descent parameter, `alpha`, in the `SNESLINESEARCHBT` `SNESLineSearch` variant.
102f4102e2SPeter Brune 
112f4102e2SPeter Brune   Input Parameters:
122f4102e2SPeter Brune + linesearch - linesearch context
132f4102e2SPeter Brune - alpha      - The descent parameter
142f4102e2SPeter Brune 
152f4102e2SPeter Brune   Level: intermediate
162f4102e2SPeter Brune 
17420bcc1bSBarry Smith .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchSetLambda()`, `SNESLineSearchGetTolerances()`, `SNESLINESEARCHBT`, `SNESLineSearchBTGetAlpha()`
182f4102e2SPeter Brune @*/
SNESLineSearchBTSetAlpha(SNESLineSearch linesearch,PetscReal alpha)19d71ae5a4SJacob Faibussowitsch PetscErrorCode SNESLineSearchBTSetAlpha(SNESLineSearch linesearch, PetscReal alpha)
20d71ae5a4SJacob Faibussowitsch {
21d7ca06c0SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
22d1b38757SStefano Zampini   PetscBool          isbt;
236e111a19SKarl Rupp 
242f4102e2SPeter Brune   PetscFunctionBegin;
252f4102e2SPeter Brune   PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 1);
26d1b38757SStefano Zampini   PetscCall(PetscObjectTypeCompare((PetscObject)linesearch, SNESLINESEARCHBT, &isbt));
27d1b38757SStefano Zampini   if (isbt) bt->alpha = alpha;
283ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
292f4102e2SPeter Brune }
302f4102e2SPeter Brune 
312f4102e2SPeter Brune /*@
32420bcc1bSBarry Smith   SNESLineSearchBTGetAlpha - Gets the descent parameter, `alpha`, in the `SNESLINESEARCHBT` variant that was set with `SNESLineSearchBTSetAlpha()`
332f4102e2SPeter Brune 
34f6dfbefdSBarry Smith   Input Parameter:
352f4102e2SPeter Brune . linesearch - linesearch context
368e557f58SPeter Brune 
37f6dfbefdSBarry Smith   Output Parameter:
382f4102e2SPeter Brune . alpha - The descent parameter
392f4102e2SPeter Brune 
402f4102e2SPeter Brune   Level: intermediate
412f4102e2SPeter Brune 
420241b274SPierre Jolivet .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchGetLambda()`, `SNESLineSearchGetTolerances()`, `SNESLINESEARCHBT`, `SNESLineSearchBTSetAlpha()`
432f4102e2SPeter Brune @*/
SNESLineSearchBTGetAlpha(SNESLineSearch linesearch,PetscReal * alpha)44d71ae5a4SJacob Faibussowitsch PetscErrorCode SNESLineSearchBTGetAlpha(SNESLineSearch linesearch, PetscReal *alpha)
45d71ae5a4SJacob Faibussowitsch {
46d7ca06c0SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
47d1b38757SStefano Zampini   PetscBool          isbt;
486e111a19SKarl Rupp 
492f4102e2SPeter Brune   PetscFunctionBegin;
502f4102e2SPeter Brune   PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 1);
51d1b38757SStefano Zampini   PetscCall(PetscObjectTypeCompare((PetscObject)linesearch, SNESLINESEARCHBT, &isbt));
52d1b38757SStefano Zampini   PetscCheck(isbt, PetscObjectComm((PetscObject)linesearch), PETSC_ERR_USER, "Not for type %s", ((PetscObject)linesearch)->type_name);
532f4102e2SPeter Brune   *alpha = bt->alpha;
543ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
552f4102e2SPeter Brune }
562f4102e2SPeter Brune 
SNESLineSearchApply_BT(SNESLineSearch linesearch)57d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchApply_BT(SNESLineSearch linesearch)
58d71ae5a4SJacob Faibussowitsch {
596b72add0SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
6071b4ebd2SPeter Brune   PetscBool          changed_y, changed_w;
6171b4ebd2SPeter Brune   Vec                X, F, Y, W, G;
6271b4ebd2SPeter Brune   SNES               snes;
638a430afbSPeter Brune   PetscReal          fnorm, xnorm, ynorm, gnorm;
64a99ef635SJonas Heinzmann   PetscReal          lambda, lambdatemp, lambdaprev, minlambda, initslope, alpha, stol;
65dfcd3828SPeter Brune   PetscReal          t1, t2, a, b, d;
668a430afbSPeter Brune   PetscReal          f;
678a430afbSPeter Brune   PetscReal          g, gprev;
6871b4ebd2SPeter Brune   PetscViewer        monitor;
69a99ef635SJonas Heinzmann   PetscInt           max_it, count;
7071b4ebd2SPeter Brune   Mat                jac;
716b72add0SBarry Smith   SNESObjectiveFn   *objective;
7232661dccSStefano Zampini   const char *const  ordStr[] = {"Linear", "Quadratic", "Cubic"};
7371b4ebd2SPeter Brune 
7471b4ebd2SPeter Brune   PetscFunctionBegin;
759566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G));
76a68bbae5SBarry Smith   PetscCall(SNESLineSearchGetNorms(linesearch, NULL, &fnorm, NULL));
779566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetLambda(linesearch, &lambda));
789566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
799566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor));
80a99ef635SJonas Heinzmann   PetscCall(SNESLineSearchGetTolerances(linesearch, &minlambda, NULL, NULL, NULL, NULL, &max_it));
819566063dSJacob Faibussowitsch   PetscCall(SNESGetTolerances(snes, NULL, NULL, &stol, NULL, NULL));
829566063dSJacob Faibussowitsch   PetscCall(SNESGetObjective(snes, &objective, NULL));
8371b4ebd2SPeter Brune   alpha = bt->alpha;
8471b4ebd2SPeter Brune 
859566063dSJacob Faibussowitsch   PetscCall(SNESGetJacobian(snes, &jac, NULL, NULL, NULL));
8608401ef6SPierre Jolivet   PetscCheck(jac || objective, PetscObjectComm((PetscObject)linesearch), PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix");
87c69d1a72SBarry Smith 
889566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y));
899566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
9071b4ebd2SPeter Brune 
919566063dSJacob Faibussowitsch   PetscCall(VecNormBegin(Y, NORM_2, &ynorm));
929566063dSJacob Faibussowitsch   PetscCall(VecNormBegin(X, NORM_2, &xnorm));
939566063dSJacob Faibussowitsch   PetscCall(VecNormEnd(Y, NORM_2, &ynorm));
949566063dSJacob Faibussowitsch   PetscCall(VecNormEnd(X, NORM_2, &xnorm));
953b288129SPeter Brune 
9671b4ebd2SPeter Brune   if (ynorm == 0.0) {
9771b4ebd2SPeter Brune     if (monitor) {
989566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
999566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Initial direction and size is 0\n"));
1009566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
10171b4ebd2SPeter Brune     }
1029566063dSJacob Faibussowitsch     PetscCall(VecCopy(X, W));
1039566063dSJacob Faibussowitsch     PetscCall(VecCopy(F, G));
1049566063dSJacob Faibussowitsch     PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm));
1059566063dSJacob Faibussowitsch     PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
1063ba16761SJacob Faibussowitsch     PetscFunctionReturn(PETSC_SUCCESS);
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));
112*76c63389SBarry Smith     SNESLineSearchCheckObjectiveDomainError(snes, f);
1138a430afbSPeter Brune   } else {
114e3ce9940SStefano Zampini     f = 0.5 * PetscSqr(fnorm);
1158a430afbSPeter Brune   }
1168a430afbSPeter Brune 
1178a430afbSPeter Brune   /* compute the initial slope */
118bd4a8a71SBarry Smith   if (objective) {
11932661dccSStefano Zampini     /* slope comes from the function (assumed to be the gradient of the objective) */
1209566063dSJacob Faibussowitsch     PetscCall(VecDotRealPart(Y, F, &initslope));
1218a430afbSPeter Brune   } else {
1228a430afbSPeter Brune     /* slope comes from the normal equations */
1239566063dSJacob Faibussowitsch     PetscCall(MatMult(jac, Y, W));
1249566063dSJacob Faibussowitsch     PetscCall(VecDotRealPart(F, W, &initslope));
12571b4ebd2SPeter Brune     if (initslope > 0.0) initslope = -initslope;
12671b4ebd2SPeter Brune     if (initslope == 0.0) initslope = -1.0;
1278a430afbSPeter Brune   }
12871b4ebd2SPeter Brune 
129*76c63389SBarry Smith   /* repeatedly cut lambda until the norm or objective function is not infinity or NaN or lambda is too small */
130df019d78SBarry Smith   while (PETSC_TRUE) {
1319566063dSJacob Faibussowitsch     PetscCall(VecWAXPY(W, -lambda, Y, X));
1321baa6e33SBarry Smith     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
133e71169deSBarry Smith     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
1349566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while checking full step length!\n"));
13571b4ebd2SPeter Brune       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
1369566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
1373ba16761SJacob Faibussowitsch       PetscFunctionReturn(PETSC_SUCCESS);
13871b4ebd2SPeter Brune     }
1398a430afbSPeter Brune 
140bd4a8a71SBarry Smith     if (objective) {
1419566063dSJacob Faibussowitsch       PetscCall(SNESComputeObjective(snes, W, &g));
1428a430afbSPeter Brune     } else {
1439566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
1449bd66eb0SPeter Brune       if (linesearch->ops->vinorm) {
1459bd66eb0SPeter Brune         gnorm = fnorm;
1469566063dSJacob Faibussowitsch         PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
1479bd66eb0SPeter Brune       } else {
1489566063dSJacob Faibussowitsch         PetscCall(VecNorm(G, NORM_2, &gnorm));
1498a430afbSPeter Brune       }
150e3ce9940SStefano Zampini       g = 0.5 * PetscSqr(gnorm);
1519bd66eb0SPeter Brune     }
1529566063dSJacob Faibussowitsch     PetscCall(SNESLineSearchMonitor(linesearch));
1539bd66eb0SPeter Brune 
154df019d78SBarry Smith     if (!PetscIsInfOrNanReal(g)) break;
155df019d78SBarry Smith     if (monitor) {
1569566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
157*76c63389SBarry Smith       PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: objective function at lambdas = %g is infinity or NaN, cutting lambda\n", (double)lambda));
1589566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
159df019d78SBarry Smith     }
160*76c63389SBarry Smith     if (lambda <= minlambda) SNESCheckFunctionDomainError(snes, g);
16132661dccSStefano Zampini     lambda *= .5;
162df019d78SBarry Smith   }
163df019d78SBarry Smith 
16448a46eb9SPierre Jolivet   if (!objective) PetscCall(PetscInfo(snes, "Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm));
165e3ce9940SStefano Zampini   if (g <= f + lambda * alpha * initslope) { /* Sufficient reduction or step tolerance convergence */
16671b4ebd2SPeter Brune     if (monitor) {
1679566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
168bd4a8a71SBarry Smith       if (!objective) {
1699566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm));
1708a430afbSPeter Brune       } else {
171e3ce9940SStefano Zampini         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Using full step: old obj %14.12e new obj %14.12e\n", (double)f, (double)g));
1728a430afbSPeter Brune       }
1739566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
17471b4ebd2SPeter Brune     }
17571b4ebd2SPeter Brune   } else {
176c21ba15cSPeter Brune     if (stol * xnorm > ynorm) {
177a68bbae5SBarry Smith       /* Since the full step didn't give sufficient decrease and the step is tiny, exit */
1789566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm));
1799566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
180c61ba1e2SPeter Brune       if (monitor) {
1819566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
18263a3b9bcSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n", (double)ynorm, (double)(stol * xnorm)));
1839566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
184c61ba1e2SPeter Brune       }
1853ba16761SJacob Faibussowitsch       PetscFunctionReturn(PETSC_SUCCESS);
186c21ba15cSPeter Brune     }
18732661dccSStefano Zampini     /* Here to avoid -Wmaybe-uninitiliazed warnings */
18871b4ebd2SPeter Brune     lambdaprev = lambda;
1898a430afbSPeter Brune     gprev      = g;
19032661dccSStefano Zampini     if (linesearch->order != SNES_LINESEARCH_ORDER_LINEAR) {
19132661dccSStefano Zampini       /* Fit points with quadratic */
192e3ce9940SStefano Zampini       lambdatemp = -initslope * PetscSqr(lambda) / (2.0 * (g - f - lambda * initslope));
193e3ce9940SStefano Zampini       lambda     = PetscClipInterval(lambdatemp, .1 * lambda, .5 * lambda);
19471b4ebd2SPeter Brune 
1959566063dSJacob Faibussowitsch       PetscCall(VecWAXPY(W, -lambda, Y, X));
1961baa6e33SBarry Smith       if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
197e71169deSBarry Smith       if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
19863a3b9bcSJacob Faibussowitsch         PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while attempting quadratic backtracking! %" PetscInt_FMT " \n", snes->nfuncs));
19971b4ebd2SPeter Brune         snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
2009566063dSJacob Faibussowitsch         PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
2013ba16761SJacob Faibussowitsch         PetscFunctionReturn(PETSC_SUCCESS);
20271b4ebd2SPeter Brune       }
203bd4a8a71SBarry Smith       if (objective) {
2049566063dSJacob Faibussowitsch         PetscCall(SNESComputeObjective(snes, W, &g));
205*76c63389SBarry Smith         SNESLineSearchCheckObjectiveDomainError(snes, 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         }
214*76c63389SBarry Smith         SNESLineSearchCheckFunctionDomainError(snes, linesearch, gnorm);
215e3ce9940SStefano Zampini         g = 0.5 * PetscSqr(gnorm);
2168a430afbSPeter Brune       }
217c98378a5SBarry Smith       if (PetscIsInfOrNanReal(g)) {
2189566063dSJacob Faibussowitsch         PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
219*76c63389SBarry Smith         PetscCall(PetscInfo(snes, "Aborted due to infinity or NaN in function evaluation\n"));
2203ba16761SJacob Faibussowitsch         PetscFunctionReturn(PETSC_SUCCESS);
221c98378a5SBarry Smith       }
22271b4ebd2SPeter Brune       if (monitor) {
2239566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
224bd4a8a71SBarry Smith         if (!objective) {
2259566063dSJacob Faibussowitsch           PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: gnorm after quadratic fit %14.12e\n", (double)gnorm));
2268a430afbSPeter Brune         } else {
2279566063dSJacob Faibussowitsch           PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: obj after quadratic fit %14.12e\n", (double)g));
2288a430afbSPeter Brune         }
2299566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
23071b4ebd2SPeter Brune       }
23132661dccSStefano Zampini     }
232e3ce9940SStefano Zampini     if (linesearch->order != SNES_LINESEARCH_ORDER_LINEAR && g <= f + lambda * alpha * initslope) { /* sufficient reduction */
23371b4ebd2SPeter Brune       if (monitor) {
2349566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
2359566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: Quadratically determined step, lambda=%18.16e\n", (double)lambda));
2369566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
23771b4ebd2SPeter Brune       }
23871b4ebd2SPeter Brune     } else {
239a99ef635SJonas Heinzmann       for (count = 0; count < max_it; count++) {
24071b4ebd2SPeter Brune         if (lambda <= minlambda) {
24171b4ebd2SPeter Brune           if (monitor) {
2429566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
24363a3b9bcSJacob Faibussowitsch             PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: unable to find good step length! After %" PetscInt_FMT " tries \n", count));
244bd4a8a71SBarry Smith             if (!objective) {
2459371c9d4SSatish Balay               PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope));
2468a430afbSPeter Brune             } else {
247e3ce9940SStefano Zampini               PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: obj(0)=%18.16e, obj=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope));
2488a430afbSPeter Brune             }
2499566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
25071b4ebd2SPeter Brune           }
2519566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
2523ba16761SJacob Faibussowitsch           PetscFunctionReturn(PETSC_SUCCESS);
25371b4ebd2SPeter Brune         }
254b000cd8dSPeter Brune         if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
255e3ce9940SStefano Zampini           /* Fit points with cubic */
256e3ce9940SStefano Zampini           t1 = g - f - lambda * initslope;
257e3ce9940SStefano Zampini           t2 = gprev - f - lambdaprev * initslope;
25871b4ebd2SPeter Brune           a  = (t1 / (lambda * lambda) - t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev);
25971b4ebd2SPeter Brune           b  = (-lambdaprev * t1 / (lambda * lambda) + lambda * t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev);
26071b4ebd2SPeter Brune           d  = b * b - 3 * a * initslope;
26171b4ebd2SPeter Brune           if (d < 0.0) d = 0.0;
262f5af7f23SKarl Rupp           if (a == 0.0) lambdatemp = -initslope / (2.0 * b);
263f5af7f23SKarl Rupp           else lambdatemp = (-b + PetscSqrtReal(d)) / (3.0 * a);
264b000cd8dSPeter Brune         } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
265e3ce9940SStefano Zampini           lambdatemp = -initslope * PetscSqr(lambda) / (2.0 * (g - f - lambda * initslope));
26632661dccSStefano Zampini         } else if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) { /* Just backtrack */
26732661dccSStefano Zampini           lambdatemp = .5 * lambda;
26832661dccSStefano Zampini         } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "Line search order %" PetscInt_FMT " for type bt", linesearch->order);
26971b4ebd2SPeter Brune         lambdaprev = lambda;
2708a430afbSPeter Brune         gprev      = g;
271e3ce9940SStefano Zampini 
272e3ce9940SStefano Zampini         lambda = PetscClipInterval(lambdatemp, .1 * lambda, .5 * lambda);
2739566063dSJacob Faibussowitsch         PetscCall(VecWAXPY(W, -lambda, Y, X));
2741baa6e33SBarry Smith         if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
275e71169deSBarry Smith         if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
27663a3b9bcSJacob Faibussowitsch           PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while looking for good step length! %" PetscInt_FMT " \n", count));
27748a46eb9SPierre Jolivet           if (!objective) PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", (double)fnorm, (double)gnorm, (double)ynorm, (double)lambda, (double)initslope));
2789566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
27971b4ebd2SPeter Brune           snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
2803ba16761SJacob Faibussowitsch           PetscFunctionReturn(PETSC_SUCCESS);
28171b4ebd2SPeter Brune         }
282bd4a8a71SBarry Smith         if (objective) {
2839566063dSJacob Faibussowitsch           PetscCall(SNESComputeObjective(snes, W, &g));
284*76c63389SBarry Smith           SNESLineSearchCheckObjectiveDomainError(snes, g);
2858a430afbSPeter Brune         } else {
2869566063dSJacob Faibussowitsch           PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
2879bd66eb0SPeter Brune           if (linesearch->ops->vinorm) {
2889bd66eb0SPeter Brune             gnorm = fnorm;
2899566063dSJacob Faibussowitsch             PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
2909bd66eb0SPeter Brune           } else {
2919566063dSJacob Faibussowitsch             PetscCall(VecNorm(G, NORM_2, &gnorm));
2928a430afbSPeter Brune           }
293*76c63389SBarry Smith           SNESLineSearchCheckFunctionDomainError(snes, linesearch, gnorm);
294e3ce9940SStefano Zampini           g = 0.5 * PetscSqr(gnorm);
2959bd66eb0SPeter Brune         }
296e3ce9940SStefano Zampini         if (g <= f + lambda * alpha * initslope) { /* is reduction enough? */
29771b4ebd2SPeter Brune           if (monitor) {
2989566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
299bd4a8a71SBarry Smith             if (!objective) {
30032661dccSStefano Zampini               PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: %s step, current gnorm %14.12e lambda=%18.16e\n", ordStr[linesearch->order - 1], (double)gnorm, (double)lambda));
3019566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3028a430afbSPeter Brune             } else {
30332661dccSStefano Zampini               PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: %s step, obj %14.12e lambda=%18.16e\n", ordStr[linesearch->order - 1], (double)g, (double)lambda));
3049566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3058a430afbSPeter Brune             }
30671b4ebd2SPeter Brune           }
30771b4ebd2SPeter Brune           break;
308f5af7f23SKarl Rupp         } else if (monitor) {
3099566063dSJacob Faibussowitsch           PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
310bd4a8a71SBarry Smith           if (!objective) {
31132661dccSStefano 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));
3129566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3138a430afbSPeter Brune           } else {
31432661dccSStefano 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));
3159566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
3168a430afbSPeter Brune           }
31771b4ebd2SPeter Brune         }
31871b4ebd2SPeter Brune       }
31971b4ebd2SPeter Brune     }
32071b4ebd2SPeter Brune   }
32171b4ebd2SPeter Brune 
32271b4ebd2SPeter Brune   /* postcheck */
3236b095a85SStefano Zampini   PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
3249566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
32571b4ebd2SPeter Brune   if (changed_y) {
3266b095a85SStefano Zampini     if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X));
3271baa6e33SBarry Smith     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
32871b4ebd2SPeter Brune   }
329bd4a8a71SBarry Smith   if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
3309566063dSJacob Faibussowitsch     PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
3319bd66eb0SPeter Brune     if (linesearch->ops->vinorm) {
3329bd66eb0SPeter Brune       gnorm = fnorm;
3339566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
3349bd66eb0SPeter Brune     } else {
3359566063dSJacob Faibussowitsch       PetscCall(VecNorm(G, NORM_2, &gnorm));
3369bd66eb0SPeter Brune     }
337*76c63389SBarry Smith     SNESLineSearchCheckFunctionDomainError(snes, linesearch, gnorm);
3389566063dSJacob Faibussowitsch     PetscCall(VecNorm(Y, NORM_2, &ynorm));
339c427506bSPeter Brune   }
34071b4ebd2SPeter Brune 
34171b4ebd2SPeter Brune   /* copy the solution over */
3429566063dSJacob Faibussowitsch   PetscCall(VecCopy(W, X));
3439566063dSJacob Faibussowitsch   PetscCall(VecCopy(G, F));
3449566063dSJacob Faibussowitsch   PetscCall(VecNorm(X, NORM_2, &xnorm));
3459566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm));
3463ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
34771b4ebd2SPeter Brune }
34871b4ebd2SPeter Brune 
SNESLineSearchView_BT(SNESLineSearch linesearch,PetscViewer viewer)34966976f2fSJacob Faibussowitsch static PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
350d71ae5a4SJacob Faibussowitsch {
3519f196a02SMartin Diehl   PetscBool          isascii;
352d7ca06c0SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
353075cc632SBarry Smith 
3547f1410a3SPeter Brune   PetscFunctionBegin;
3559f196a02SMartin Diehl   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
3569f196a02SMartin Diehl   if (isascii) {
357b000cd8dSPeter Brune     if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3589566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  interpolation: cubic\n"));
359b000cd8dSPeter Brune     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
3609566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  interpolation: quadratic\n"));
3617f1410a3SPeter Brune     }
3629566063dSJacob Faibussowitsch     PetscCall(PetscViewerASCIIPrintf(viewer, "  alpha=%e\n", (double)bt->alpha));
3637f1410a3SPeter Brune   }
3643ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
3657f1410a3SPeter Brune }
3667f1410a3SPeter Brune 
SNESLineSearchDestroy_BT(SNESLineSearch linesearch)367d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
368d71ae5a4SJacob Faibussowitsch {
36971b4ebd2SPeter Brune   PetscFunctionBegin;
3709566063dSJacob Faibussowitsch   PetscCall(PetscFree(linesearch->data));
3713ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
37271b4ebd2SPeter Brune }
37371b4ebd2SPeter Brune 
SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch,PetscOptionItems PetscOptionsObject)374ce78bad3SBarry Smith static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch, PetscOptionItems PetscOptionsObject)
375d71ae5a4SJacob Faibussowitsch {
376eb23ba39SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data;
37771b4ebd2SPeter Brune 
3786e111a19SKarl Rupp   PetscFunctionBegin;
379d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject, "SNESLineSearch BT options");
3809566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL));
381d0609cedSBarry Smith   PetscOptionsHeadEnd();
3823ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
38371b4ebd2SPeter Brune }
38471b4ebd2SPeter Brune 
385954494b2SJed Brown /*MC
3861d27aa22SBarry Smith    SNESLINESEARCHBT - Backtracking line search {cite}`dennis:83`.
387954494b2SJed Brown 
388a99ef635SJonas Heinzmann    This line search finds the minimum of a polynomial fitting either $1/2 ||F(x_k + \lambda Y_k)||_2^2$,
389a99ef635SJonas Heinzmann    or the objective function $G(x_k + \lambda Y_k)$ if it is provided with `SNESSetObjective()`.
390a99ef635SJonas Heinzmann    If this fit does not satisfy the sufficient decrease conditions, the interval shrinks
39189dfbbd5SBarry Smith    and the fit is reattempted at most `max_it` times or until $\lambda$ is below `minlambda`.
392954494b2SJed Brown 
393cd7522eaSPeter Brune    Options Database Keys:
3943eec2f6aSPierre Jolivet +  -snes_linesearch_alpha <1e\-4>      - slope descent parameter
395a99ef635SJonas Heinzmann .  -snes_linesearch_damping <1.0>      - initial `lambda` on entry to the line search
396a99ef635SJonas Heinzmann .  -snes_linesearch_max_it <40>        - maximum number of shrinking iterations in the line search
397a99ef635SJonas Heinzmann .  -snes_linesearch_minlambda <1e\-12> - minimum `lambda` (scaling of solution update) allowed
398a99ef635SJonas Heinzmann -  -snes_linesearch_order <3>          - order of the polynomial fit, must be 1, 2, or 3. With order 1, it performs a simple backtracking without any curve fitting
399cd7522eaSPeter Brune 
400954494b2SJed Brown    Level: advanced
401954494b2SJed Brown 
402f6dfbefdSBarry Smith    Note:
403e55accfeSBarryFSmith    This line search will always produce a step that is less than or equal to, in length, the full step size.
404db609ea7SPeter Brune 
405420bcc1bSBarry Smith .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
406954494b2SJed Brown M*/
SNESLineSearchCreate_BT(SNESLineSearch linesearch)407d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
408d71ae5a4SJacob Faibussowitsch {
409f1c6b773SPeter Brune   SNESLineSearch_BT *bt;
41071b4ebd2SPeter Brune 
41171b4ebd2SPeter Brune   PetscFunctionBegin;
412f1c6b773SPeter Brune   linesearch->ops->apply          = SNESLineSearchApply_BT;
413f1c6b773SPeter Brune   linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
414f1c6b773SPeter Brune   linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
4150298fd71SBarry Smith   linesearch->ops->reset          = NULL;
4167f1410a3SPeter Brune   linesearch->ops->view           = SNESLineSearchView_BT;
4170298fd71SBarry Smith   linesearch->ops->setup          = NULL;
41871b4ebd2SPeter Brune 
4194dfa11a4SJacob Faibussowitsch   PetscCall(PetscNew(&bt));
420f5af7f23SKarl Rupp 
42171b4ebd2SPeter Brune   linesearch->data   = (void *)bt;
422a99ef635SJonas Heinzmann   linesearch->max_it = 40;
423b000cd8dSPeter Brune   linesearch->order  = SNES_LINESEARCH_ORDER_CUBIC;
42471b4ebd2SPeter Brune   bt->alpha          = 1e-4;
4253ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
42671b4ebd2SPeter Brune }
427