xref: /petsc/src/snes/linesearch/impls/bt/linesearchbt.c (revision 1baa6e3354bfe224b33a0c158f482508792a8a6e)
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 @*/
192f4102e2SPeter Brune PetscErrorCode SNESLineSearchBTSetAlpha(SNESLineSearch linesearch, PetscReal alpha)
202f4102e2SPeter Brune {
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;
262f4102e2SPeter Brune   PetscFunctionReturn(0);
272f4102e2SPeter Brune }
282f4102e2SPeter Brune 
292f4102e2SPeter Brune /*@
302f4102e2SPeter Brune    SNESLineSearchBTGetAlpha - Gets the descent parameter, alpha, in the BT linesearch variant.
312f4102e2SPeter Brune 
322f4102e2SPeter Brune    Input Parameters:
332f4102e2SPeter Brune .  linesearch - linesearch context
348e557f58SPeter Brune 
358e557f58SPeter Brune    Output Parameters:
362f4102e2SPeter Brune .  alpha - The descent parameter
372f4102e2SPeter Brune 
382f4102e2SPeter Brune    Level: intermediate
392f4102e2SPeter Brune 
40db781477SPatrick Sanan .seealso: `SNESLineSearchGetLambda()`, `SNESLineSearchGetTolerances()` `SNESLINESEARCHBT`
412f4102e2SPeter Brune @*/
422f4102e2SPeter Brune PetscErrorCode SNESLineSearchBTGetAlpha(SNESLineSearch linesearch, PetscReal *alpha)
432f4102e2SPeter Brune {
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;
492f4102e2SPeter Brune   PetscFunctionReturn(0);
502f4102e2SPeter Brune }
512f4102e2SPeter Brune 
52f1c6b773SPeter Brune static PetscErrorCode  SNESLineSearchApply_BT(SNESLineSearch linesearch)
5371b4ebd2SPeter Brune {
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));
10071b4ebd2SPeter Brune     PetscFunctionReturn(0);
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));
133*1baa6e33SBarry 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));
13871b4ebd2SPeter Brune       PetscFunctionReturn(0);
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     }
161df019d78SBarry Smith     if (lambda <= minlambda) {
1627e49c775SBarry Smith       SNESCheckFunctionNorm(snes,g);
163c98378a5SBarry Smith     }
164df019d78SBarry Smith     lambda = .5*lambda;
165df019d78SBarry Smith   }
166df019d78SBarry Smith 
167bd4a8a71SBarry Smith   if (!objective) {
1689566063dSJacob Faibussowitsch     PetscCall(PetscInfo(snes,"Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm));
1698a430afbSPeter Brune   }
1708a430afbSPeter Brune   if (.5*g <= .5*f + lambda*alpha*initslope) { /* Sufficient reduction or step tolerance convergence */
17171b4ebd2SPeter Brune     if (monitor) {
1729566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
173bd4a8a71SBarry Smith       if (!objective) {
1749566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm));
1758a430afbSPeter Brune       } else {
1769566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Using full step: obj %14.12e obj %14.12e\n", (double)f, (double)g));
1778a430afbSPeter Brune       }
1789566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
17971b4ebd2SPeter Brune     }
18071b4ebd2SPeter Brune   } else {
181c21ba15cSPeter Brune     /* Since the full step didn't work and the step is tiny, quit */
182c21ba15cSPeter Brune     if (stol*xnorm > ynorm) {
1839566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm));
1849566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
185c61ba1e2SPeter Brune       if (monitor) {
1869566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
18763a3b9bcSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n",(double)ynorm,(double)(stol*xnorm)));
1889566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
189c61ba1e2SPeter Brune       }
190c21ba15cSPeter Brune       PetscFunctionReturn(0);
191c21ba15cSPeter Brune     }
19271b4ebd2SPeter Brune     /* Fit points with quadratic */
1938a430afbSPeter Brune     lambdatemp = -initslope/(g - f - 2.0*lambda*initslope);
19471b4ebd2SPeter Brune     lambdaprev = lambda;
1958a430afbSPeter Brune     gprev      = g;
19671b4ebd2SPeter Brune     if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
19771b4ebd2SPeter Brune     if (lambdatemp <= .1*lambda) lambda = .1*lambda;
19871b4ebd2SPeter Brune     else                         lambda = lambdatemp;
19971b4ebd2SPeter Brune 
2009566063dSJacob Faibussowitsch     PetscCall(VecWAXPY(W,-lambda,Y,X));
201*1baa6e33SBarry Smith     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
202e71169deSBarry Smith     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
20363a3b9bcSJacob Faibussowitsch       PetscCall(PetscInfo(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %" PetscInt_FMT " \n",snes->nfuncs));
20471b4ebd2SPeter Brune       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
2059566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
20671b4ebd2SPeter Brune       PetscFunctionReturn(0);
20771b4ebd2SPeter Brune     }
208bd4a8a71SBarry Smith     if (objective) {
2099566063dSJacob Faibussowitsch       PetscCall(SNESComputeObjective(snes,W,&g));
2108a430afbSPeter Brune     } else {
2119566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->snesfunc)(snes,W,G));
2129bd66eb0SPeter Brune       if (linesearch->ops->vinorm) {
2139bd66eb0SPeter Brune         gnorm = fnorm;
2149566063dSJacob Faibussowitsch         PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
2159bd66eb0SPeter Brune       } else {
2169566063dSJacob Faibussowitsch         PetscCall(VecNorm(G,NORM_2,&gnorm));
2179bd66eb0SPeter Brune       }
218f1ab69e3SPeter Brune       g = PetscSqr(gnorm);
2198a430afbSPeter Brune     }
220c98378a5SBarry Smith     if (PetscIsInfOrNanReal(g)) {
2219566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
2229566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n"));
223c98378a5SBarry Smith       PetscFunctionReturn(0);
224c98378a5SBarry Smith     }
22571b4ebd2SPeter Brune     if (monitor) {
2269566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
227bd4a8a71SBarry Smith       if (!objective) {
2289566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: gnorm after quadratic fit %14.12e\n",(double)gnorm));
2298a430afbSPeter Brune       } else {
2309566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: obj after quadratic fit %14.12e\n",(double)g));
2318a430afbSPeter Brune       }
2329566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
23371b4ebd2SPeter Brune     }
2348a430afbSPeter Brune     if (.5*g < .5*f + lambda*alpha*initslope) { /* sufficient reduction */
23571b4ebd2SPeter Brune       if (monitor) {
2369566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
2379566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda));
2389566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
23971b4ebd2SPeter Brune       }
24071b4ebd2SPeter Brune     } else {
24171b4ebd2SPeter Brune       /* Fit points with cubic */
24271b4ebd2SPeter Brune       for (count = 0; count < max_its; count++) {
24371b4ebd2SPeter Brune         if (lambda <= minlambda) {
24471b4ebd2SPeter Brune           if (monitor) {
2459566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
24663a3b9bcSJacob Faibussowitsch             PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: unable to find good step length! After %" PetscInt_FMT " tries \n",count));
247bd4a8a71SBarry Smith             if (!objective) {
248d0609cedSBarry Smith               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",
249d0609cedSBarry Smith                                                (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope));
2508a430afbSPeter Brune             } else {
251d0609cedSBarry Smith               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",
252d0609cedSBarry Smith                                                (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope));
2538a430afbSPeter Brune             }
2549566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
25571b4ebd2SPeter Brune           }
2569566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
25771b4ebd2SPeter Brune           PetscFunctionReturn(0);
25871b4ebd2SPeter Brune         }
259b000cd8dSPeter Brune         if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
2608a430afbSPeter Brune           t1 = .5*(g - f) - lambda*initslope;
2618a430afbSPeter Brune           t2 = .5*(gprev  - f) - lambdaprev*initslope;
26271b4ebd2SPeter Brune           a  = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
26371b4ebd2SPeter Brune           b  = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
26471b4ebd2SPeter Brune           d  = b*b - 3*a*initslope;
26571b4ebd2SPeter Brune           if (d < 0.0) d = 0.0;
266f5af7f23SKarl Rupp           if (a == 0.0) lambdatemp = -initslope/(2.0*b);
267f5af7f23SKarl Rupp           else lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a);
268f5af7f23SKarl Rupp 
269b000cd8dSPeter Brune         } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
2708a430afbSPeter Brune           lambdatemp = -initslope/(g - f - 2.0*initslope);
271ce94432eSBarry Smith         } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt");
27271b4ebd2SPeter Brune         lambdaprev = lambda;
2738a430afbSPeter Brune         gprev      = g;
27471b4ebd2SPeter Brune         if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
27571b4ebd2SPeter Brune         if (lambdatemp <= .1*lambda) lambda     = .1*lambda;
27671b4ebd2SPeter Brune         else                         lambda     = lambdatemp;
2779566063dSJacob Faibussowitsch         PetscCall(VecWAXPY(W,-lambda,Y,X));
278*1baa6e33SBarry Smith         if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes,W));
279e71169deSBarry Smith         if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
28063a3b9bcSJacob Faibussowitsch           PetscCall(PetscInfo(snes,"Exceeded maximum function evaluations, while looking for good step length! %" PetscInt_FMT " \n",count));
281bd4a8a71SBarry Smith           if (!objective) {
282d0609cedSBarry Smith             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));
2838a430afbSPeter Brune           }
2849566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
28571b4ebd2SPeter Brune           snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
28671b4ebd2SPeter Brune           PetscFunctionReturn(0);
28771b4ebd2SPeter Brune         }
288bd4a8a71SBarry Smith         if (objective) {
2899566063dSJacob Faibussowitsch           PetscCall(SNESComputeObjective(snes,W,&g));
2908a430afbSPeter Brune         } else {
2919566063dSJacob Faibussowitsch           PetscCall((*linesearch->ops->snesfunc)(snes,W,G));
2929bd66eb0SPeter Brune           if (linesearch->ops->vinorm) {
2939bd66eb0SPeter Brune             gnorm = fnorm;
2949566063dSJacob Faibussowitsch             PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
2959bd66eb0SPeter Brune           } else {
2969566063dSJacob Faibussowitsch             PetscCall(VecNorm(G,NORM_2,&gnorm));
2978a430afbSPeter Brune           }
298f1ab69e3SPeter Brune           g = PetscSqr(gnorm);
2999bd66eb0SPeter Brune         }
3004a2f8832SBarry Smith         if (PetscIsInfOrNanReal(g)) {
3019566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
3029566063dSJacob Faibussowitsch           PetscCall(PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n"));
303c98378a5SBarry Smith           PetscFunctionReturn(0);
304c98378a5SBarry Smith         }
3058a430afbSPeter Brune         if (.5*g < .5*f + lambda*alpha*initslope) { /* is reduction enough? */
30671b4ebd2SPeter Brune           if (monitor) {
3079566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
308bd4a8a71SBarry Smith             if (!objective) {
309b000cd8dSPeter Brune               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3109566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda));
31171b4ebd2SPeter Brune               } else {
3129566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda));
31371b4ebd2SPeter Brune               }
3149566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
3158a430afbSPeter Brune             } else {
3168a430afbSPeter Brune               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3179566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda));
3188a430afbSPeter Brune               } else {
3199566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda));
3208a430afbSPeter Brune               }
3219566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
3228a430afbSPeter Brune             }
32371b4ebd2SPeter Brune           }
32471b4ebd2SPeter Brune           break;
325f5af7f23SKarl Rupp         } else if (monitor) {
3269566063dSJacob Faibussowitsch           PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
327bd4a8a71SBarry Smith           if (!objective) {
328b000cd8dSPeter Brune             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3299566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda));
33071b4ebd2SPeter Brune             } else {
3319566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda));
33271b4ebd2SPeter Brune             }
3339566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
3348a430afbSPeter Brune           } else {
3358a430afbSPeter Brune             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3369566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda));
3378a430afbSPeter Brune             } else {
3389566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda));
3398a430afbSPeter Brune             }
3409566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
3418a430afbSPeter Brune           }
34271b4ebd2SPeter Brune         }
34371b4ebd2SPeter Brune       }
34471b4ebd2SPeter Brune     }
34571b4ebd2SPeter Brune   }
34671b4ebd2SPeter Brune 
34771b4ebd2SPeter Brune   /* postcheck */
3480c1ac483SGlenn Hammond   /* update Y to lambda*Y so that W is consistent with  X - lambda*Y */
3499566063dSJacob Faibussowitsch   PetscCall(VecScale(Y,lambda));
3509566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w));
35171b4ebd2SPeter Brune   if (changed_y) {
3529566063dSJacob Faibussowitsch     PetscCall(VecWAXPY(W,-1.0,Y,X));
353*1baa6e33SBarry Smith     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
35471b4ebd2SPeter Brune   }
355bd4a8a71SBarry Smith   if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
3569566063dSJacob Faibussowitsch     PetscCall((*linesearch->ops->snesfunc)(snes,W,G));
3579bd66eb0SPeter Brune     if (linesearch->ops->vinorm) {
3589bd66eb0SPeter Brune       gnorm = fnorm;
3599566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
3609bd66eb0SPeter Brune     } else {
3619566063dSJacob Faibussowitsch       PetscCall(VecNorm(G,NORM_2,&gnorm));
3629bd66eb0SPeter Brune     }
3639566063dSJacob Faibussowitsch     PetscCall(VecNorm(Y,NORM_2,&ynorm));
364c98378a5SBarry Smith     if (PetscIsInfOrNanReal(gnorm)) {
3659566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_FAILED_NANORINF));
3669566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n"));
367c98378a5SBarry Smith       PetscFunctionReturn(0);
368c98378a5SBarry Smith     }
369c427506bSPeter Brune   }
37071b4ebd2SPeter Brune 
37171b4ebd2SPeter Brune   /* copy the solution over */
3729566063dSJacob Faibussowitsch   PetscCall(VecCopy(W, X));
3739566063dSJacob Faibussowitsch   PetscCall(VecCopy(G, F));
3749566063dSJacob Faibussowitsch   PetscCall(VecNorm(X, NORM_2, &xnorm));
3759566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
3769566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm));
37771b4ebd2SPeter Brune   PetscFunctionReturn(0);
37871b4ebd2SPeter Brune }
37971b4ebd2SPeter Brune 
3807f1410a3SPeter Brune PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
3817f1410a3SPeter Brune {
3827f1410a3SPeter Brune   PetscBool         iascii;
383d7ca06c0SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
384075cc632SBarry Smith 
3857f1410a3SPeter Brune   PetscFunctionBegin;
3869566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii));
3877f1410a3SPeter Brune   if (iascii) {
388b000cd8dSPeter Brune     if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3899566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  interpolation: cubic\n"));
390b000cd8dSPeter Brune     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
3919566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  interpolation: quadratic\n"));
3927f1410a3SPeter Brune     }
3939566063dSJacob Faibussowitsch     PetscCall(PetscViewerASCIIPrintf(viewer, "  alpha=%e\n", (double)bt->alpha));
3947f1410a3SPeter Brune   }
3957f1410a3SPeter Brune   PetscFunctionReturn(0);
3967f1410a3SPeter Brune }
3977f1410a3SPeter Brune 
398f1c6b773SPeter Brune static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
399731c067cSJed Brown {
40071b4ebd2SPeter Brune   PetscFunctionBegin;
4019566063dSJacob Faibussowitsch   PetscCall(PetscFree(linesearch->data));
40271b4ebd2SPeter Brune   PetscFunctionReturn(0);
40371b4ebd2SPeter Brune }
40471b4ebd2SPeter Brune 
4054416b707SBarry Smith static PetscErrorCode SNESLineSearchSetFromOptions_BT(PetscOptionItems *PetscOptionsObject,SNESLineSearch linesearch)
40671b4ebd2SPeter Brune {
407eb23ba39SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
40871b4ebd2SPeter Brune 
4096e111a19SKarl Rupp   PetscFunctionBegin;
410d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject,"SNESLineSearch BT options");
4119566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_linesearch_alpha",   "Descent tolerance",        "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL));
412d0609cedSBarry Smith   PetscOptionsHeadEnd();
41371b4ebd2SPeter Brune   PetscFunctionReturn(0);
41471b4ebd2SPeter Brune }
41571b4ebd2SPeter Brune 
416954494b2SJed Brown /*MC
417db609ea7SPeter Brune    SNESLINESEARCHBT - Backtracking line search.
418954494b2SJed Brown 
419db609ea7SPeter Brune    This line search finds the minimum of a polynomial fitting of the L2 norm of the
420eb23ba39SBarry 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
421db609ea7SPeter Brune    and the fit is reattempted at most max_it times or until lambda is below minlambda.
422954494b2SJed Brown 
423cd7522eaSPeter Brune    Options Database Keys:
4243eec2f6aSPierre Jolivet +  -snes_linesearch_alpha <1e\-4> - slope descent parameter
425db609ea7SPeter Brune .  -snes_linesearch_damping <1.0> - initial step length
426eb23ba39SBarry Smith .  -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the
427e1bc860dSBarry Smith                                        step is scaled back to be of this length at the beginning of the line search
428db609ea7SPeter Brune .  -snes_linesearch_max_it <40> - maximum number of shrinking step
4293eec2f6aSPierre Jolivet .  -snes_linesearch_minlambda <1e\-12> - minimum step length allowed
430cd7522eaSPeter Brune -  -snes_linesearch_order <cubic,quadratic> - order of the approximation
431cd7522eaSPeter Brune 
432954494b2SJed Brown    Level: advanced
433954494b2SJed Brown 
434db609ea7SPeter Brune    Notes:
435db609ea7SPeter Brune    This line search is taken from "Numerical Methods for Unconstrained
436db609ea7SPeter Brune    Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
437e55accfeSBarryFSmith 
438e55accfeSBarryFSmith    This line search will always produce a step that is less than or equal to, in length, the full step size.
439db609ea7SPeter Brune 
440db781477SPatrick Sanan .seealso: `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
441954494b2SJed Brown M*/
4428cc058d9SJed Brown PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
44371b4ebd2SPeter Brune {
44471b4ebd2SPeter Brune 
445f1c6b773SPeter Brune   SNESLineSearch_BT *bt;
44671b4ebd2SPeter Brune 
44771b4ebd2SPeter Brune   PetscFunctionBegin;
448f1c6b773SPeter Brune   linesearch->ops->apply          = SNESLineSearchApply_BT;
449f1c6b773SPeter Brune   linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
450f1c6b773SPeter Brune   linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
4510298fd71SBarry Smith   linesearch->ops->reset          = NULL;
4527f1410a3SPeter Brune   linesearch->ops->view           = SNESLineSearchView_BT;
4530298fd71SBarry Smith   linesearch->ops->setup          = NULL;
45471b4ebd2SPeter Brune 
4559566063dSJacob Faibussowitsch   PetscCall(PetscNewLog(linesearch,&bt));
456f5af7f23SKarl Rupp 
45771b4ebd2SPeter Brune   linesearch->data    = (void*)bt;
45871b4ebd2SPeter Brune   linesearch->max_its = 40;
459b000cd8dSPeter Brune   linesearch->order   = SNES_LINESEARCH_ORDER_CUBIC;
46071b4ebd2SPeter Brune   bt->alpha           = 1e-4;
46171b4ebd2SPeter Brune   PetscFunctionReturn(0);
46271b4ebd2SPeter Brune }
463