xref: /petsc/src/snes/linesearch/impls/bt/linesearchbt.c (revision d0609ced746bc51b019815ca91d747429db24893)
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 
171a4f838cSPeter Brune .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 
401a4f838cSPeter Brune .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));
1339bd66eb0SPeter Brune     if (linesearch->ops->viproject) {
1349566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->viproject)(snes, W));
1359bd66eb0SPeter Brune     }
136e71169deSBarry Smith     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
1379566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n"));
13871b4ebd2SPeter Brune       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
1399566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
14071b4ebd2SPeter Brune       PetscFunctionReturn(0);
14171b4ebd2SPeter Brune     }
1428a430afbSPeter Brune 
143bd4a8a71SBarry Smith     if (objective) {
1449566063dSJacob Faibussowitsch       PetscCall(SNESComputeObjective(snes,W,&g));
1458a430afbSPeter Brune     } else {
1469566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->snesfunc)(snes,W,G));
1479bd66eb0SPeter Brune       if (linesearch->ops->vinorm) {
1489bd66eb0SPeter Brune         gnorm = fnorm;
1499566063dSJacob Faibussowitsch         PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
1509bd66eb0SPeter Brune       } else {
1519566063dSJacob Faibussowitsch         PetscCall(VecNorm(G,NORM_2,&gnorm));
1528a430afbSPeter Brune       }
15379e7e81bSJed Brown       g = PetscSqr(gnorm);
1549bd66eb0SPeter Brune     }
1559566063dSJacob Faibussowitsch     PetscCall(SNESLineSearchMonitor(linesearch));
1569bd66eb0SPeter Brune 
157df019d78SBarry Smith     if (!PetscIsInfOrNanReal(g)) break;
158df019d78SBarry Smith     if (monitor) {
1599566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
1609566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n",(double)lambda));
1619566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
162df019d78SBarry Smith     }
163df019d78SBarry Smith     if (lambda <= minlambda) {
1647e49c775SBarry Smith       SNESCheckFunctionNorm(snes,g);
165c98378a5SBarry Smith     }
166df019d78SBarry Smith     lambda = .5*lambda;
167df019d78SBarry Smith   }
168df019d78SBarry Smith 
169bd4a8a71SBarry Smith   if (!objective) {
1709566063dSJacob Faibussowitsch     PetscCall(PetscInfo(snes,"Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm));
1718a430afbSPeter Brune   }
1728a430afbSPeter Brune   if (.5*g <= .5*f + lambda*alpha*initslope) { /* Sufficient reduction or step tolerance convergence */
17371b4ebd2SPeter Brune     if (monitor) {
1749566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
175bd4a8a71SBarry Smith       if (!objective) {
1769566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm));
1778a430afbSPeter Brune       } else {
1789566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Using full step: obj %14.12e obj %14.12e\n", (double)f, (double)g));
1798a430afbSPeter Brune       }
1809566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
18171b4ebd2SPeter Brune     }
18271b4ebd2SPeter Brune   } else {
183c21ba15cSPeter Brune     /* Since the full step didn't work and the step is tiny, quit */
184c21ba15cSPeter Brune     if (stol*xnorm > ynorm) {
1859566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm));
1869566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
187c61ba1e2SPeter Brune       if (monitor) {
1889566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
1899566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n",(double)ynorm,(double)stol*xnorm));
1909566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
191c61ba1e2SPeter Brune       }
192c21ba15cSPeter Brune       PetscFunctionReturn(0);
193c21ba15cSPeter Brune     }
19471b4ebd2SPeter Brune     /* Fit points with quadratic */
1958a430afbSPeter Brune     lambdatemp = -initslope/(g - f - 2.0*lambda*initslope);
19671b4ebd2SPeter Brune     lambdaprev = lambda;
1978a430afbSPeter Brune     gprev      = g;
19871b4ebd2SPeter Brune     if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
19971b4ebd2SPeter Brune     if (lambdatemp <= .1*lambda) lambda = .1*lambda;
20071b4ebd2SPeter Brune     else                         lambda = lambdatemp;
20171b4ebd2SPeter Brune 
2029566063dSJacob Faibussowitsch     PetscCall(VecWAXPY(W,-lambda,Y,X));
2039bd66eb0SPeter Brune     if (linesearch->ops->viproject) {
2049566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->viproject)(snes, W));
2059bd66eb0SPeter Brune     }
206e71169deSBarry Smith     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
2079566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs));
20871b4ebd2SPeter Brune       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
2099566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
21071b4ebd2SPeter Brune       PetscFunctionReturn(0);
21171b4ebd2SPeter Brune     }
212bd4a8a71SBarry Smith     if (objective) {
2139566063dSJacob Faibussowitsch       PetscCall(SNESComputeObjective(snes,W,&g));
2148a430afbSPeter Brune     } else {
2159566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->snesfunc)(snes,W,G));
2169bd66eb0SPeter Brune       if (linesearch->ops->vinorm) {
2179bd66eb0SPeter Brune         gnorm = fnorm;
2189566063dSJacob Faibussowitsch         PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
2199bd66eb0SPeter Brune       } else {
2209566063dSJacob Faibussowitsch         PetscCall(VecNorm(G,NORM_2,&gnorm));
2219bd66eb0SPeter Brune       }
222f1ab69e3SPeter Brune       g = PetscSqr(gnorm);
2238a430afbSPeter Brune     }
224c98378a5SBarry Smith     if (PetscIsInfOrNanReal(g)) {
2259566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
2269566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n"));
227c98378a5SBarry Smith       PetscFunctionReturn(0);
228c98378a5SBarry Smith     }
22971b4ebd2SPeter Brune     if (monitor) {
2309566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
231bd4a8a71SBarry Smith       if (!objective) {
2329566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: gnorm after quadratic fit %14.12e\n",(double)gnorm));
2338a430afbSPeter Brune       } else {
2349566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: obj after quadratic fit %14.12e\n",(double)g));
2358a430afbSPeter Brune       }
2369566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
23771b4ebd2SPeter Brune     }
2388a430afbSPeter Brune     if (.5*g < .5*f + lambda*alpha*initslope) { /* sufficient reduction */
23971b4ebd2SPeter Brune       if (monitor) {
2409566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
2419566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda));
2429566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
24371b4ebd2SPeter Brune       }
24471b4ebd2SPeter Brune     } else {
24571b4ebd2SPeter Brune       /* Fit points with cubic */
24671b4ebd2SPeter Brune       for (count = 0; count < max_its; count++) {
24771b4ebd2SPeter Brune         if (lambda <= minlambda) {
24871b4ebd2SPeter Brune           if (monitor) {
2499566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
2509566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: unable to find good step length! After %D tries \n",count));
251bd4a8a71SBarry Smith             if (!objective) {
252*d0609cedSBarry 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",
253*d0609cedSBarry Smith                                                (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope));
2548a430afbSPeter Brune             } else {
255*d0609cedSBarry 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",
256*d0609cedSBarry Smith                                                (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope));
2578a430afbSPeter Brune             }
2589566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
25971b4ebd2SPeter Brune           }
2609566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
26171b4ebd2SPeter Brune           PetscFunctionReturn(0);
26271b4ebd2SPeter Brune         }
263b000cd8dSPeter Brune         if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
2648a430afbSPeter Brune           t1 = .5*(g - f) - lambda*initslope;
2658a430afbSPeter Brune           t2 = .5*(gprev  - f) - lambdaprev*initslope;
26671b4ebd2SPeter Brune           a  = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
26771b4ebd2SPeter Brune           b  = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
26871b4ebd2SPeter Brune           d  = b*b - 3*a*initslope;
26971b4ebd2SPeter Brune           if (d < 0.0) d = 0.0;
270f5af7f23SKarl Rupp           if (a == 0.0) lambdatemp = -initslope/(2.0*b);
271f5af7f23SKarl Rupp           else lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a);
272f5af7f23SKarl Rupp 
273b000cd8dSPeter Brune         } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
2748a430afbSPeter Brune           lambdatemp = -initslope/(g - f - 2.0*initslope);
275ce94432eSBarry Smith         } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt");
27671b4ebd2SPeter Brune         lambdaprev = lambda;
2778a430afbSPeter Brune         gprev      = g;
27871b4ebd2SPeter Brune         if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
27971b4ebd2SPeter Brune         if (lambdatemp <= .1*lambda) lambda     = .1*lambda;
28071b4ebd2SPeter Brune         else                         lambda     = lambdatemp;
2819566063dSJacob Faibussowitsch         PetscCall(VecWAXPY(W,-lambda,Y,X));
282b7b2e573SPeter Brune         if (linesearch->ops->viproject) {
2839566063dSJacob Faibussowitsch           PetscCall((*linesearch->ops->viproject)(snes,W));
284b7b2e573SPeter Brune         }
285e71169deSBarry Smith         if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
2869566063dSJacob Faibussowitsch           PetscCall(PetscInfo(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count));
287bd4a8a71SBarry Smith           if (!objective) {
288*d0609cedSBarry 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));
2898a430afbSPeter Brune           }
2909566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
29171b4ebd2SPeter Brune           snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
29271b4ebd2SPeter Brune           PetscFunctionReturn(0);
29371b4ebd2SPeter Brune         }
294bd4a8a71SBarry Smith         if (objective) {
2959566063dSJacob Faibussowitsch           PetscCall(SNESComputeObjective(snes,W,&g));
2968a430afbSPeter Brune         } else {
2979566063dSJacob Faibussowitsch           PetscCall((*linesearch->ops->snesfunc)(snes,W,G));
2989bd66eb0SPeter Brune           if (linesearch->ops->vinorm) {
2999bd66eb0SPeter Brune             gnorm = fnorm;
3009566063dSJacob Faibussowitsch             PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
3019bd66eb0SPeter Brune           } else {
3029566063dSJacob Faibussowitsch             PetscCall(VecNorm(G,NORM_2,&gnorm));
3038a430afbSPeter Brune           }
304f1ab69e3SPeter Brune           g = PetscSqr(gnorm);
3059bd66eb0SPeter Brune         }
3064a2f8832SBarry Smith         if (PetscIsInfOrNanReal(g)) {
3079566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
3089566063dSJacob Faibussowitsch           PetscCall(PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n"));
309c98378a5SBarry Smith           PetscFunctionReturn(0);
310c98378a5SBarry Smith         }
3118a430afbSPeter Brune         if (.5*g < .5*f + lambda*alpha*initslope) { /* is reduction enough? */
31271b4ebd2SPeter Brune           if (monitor) {
3139566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
314bd4a8a71SBarry Smith             if (!objective) {
315b000cd8dSPeter Brune               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3169566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda));
31771b4ebd2SPeter Brune               } else {
3189566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda));
31971b4ebd2SPeter Brune               }
3209566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
3218a430afbSPeter Brune             } else {
3228a430afbSPeter Brune               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3239566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda));
3248a430afbSPeter Brune               } else {
3259566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda));
3268a430afbSPeter Brune               }
3279566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
3288a430afbSPeter Brune             }
32971b4ebd2SPeter Brune           }
33071b4ebd2SPeter Brune           break;
331f5af7f23SKarl Rupp         } else if (monitor) {
3329566063dSJacob Faibussowitsch           PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
333bd4a8a71SBarry Smith           if (!objective) {
334b000cd8dSPeter Brune             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3359566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda));
33671b4ebd2SPeter Brune             } else {
3379566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda));
33871b4ebd2SPeter Brune             }
3399566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
3408a430afbSPeter Brune           } else {
3418a430afbSPeter Brune             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3429566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda));
3438a430afbSPeter Brune             } else {
3449566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda));
3458a430afbSPeter Brune             }
3469566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
3478a430afbSPeter Brune           }
34871b4ebd2SPeter Brune         }
34971b4ebd2SPeter Brune       }
35071b4ebd2SPeter Brune     }
35171b4ebd2SPeter Brune   }
35271b4ebd2SPeter Brune 
35371b4ebd2SPeter Brune   /* postcheck */
3540c1ac483SGlenn Hammond   /* update Y to lambda*Y so that W is consistent with  X - lambda*Y */
3559566063dSJacob Faibussowitsch   PetscCall(VecScale(Y,lambda));
3569566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w));
35771b4ebd2SPeter Brune   if (changed_y) {
3589566063dSJacob Faibussowitsch     PetscCall(VecWAXPY(W,-1.0,Y,X));
3599bd66eb0SPeter Brune     if (linesearch->ops->viproject) {
3609566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->viproject)(snes, W));
3619bd66eb0SPeter Brune     }
36271b4ebd2SPeter Brune   }
363bd4a8a71SBarry Smith   if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
3649566063dSJacob Faibussowitsch     PetscCall((*linesearch->ops->snesfunc)(snes,W,G));
3659bd66eb0SPeter Brune     if (linesearch->ops->vinorm) {
3669bd66eb0SPeter Brune       gnorm = fnorm;
3679566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
3689bd66eb0SPeter Brune     } else {
3699566063dSJacob Faibussowitsch       PetscCall(VecNorm(G,NORM_2,&gnorm));
3709bd66eb0SPeter Brune     }
3719566063dSJacob Faibussowitsch     PetscCall(VecNorm(Y,NORM_2,&ynorm));
372c98378a5SBarry Smith     if (PetscIsInfOrNanReal(gnorm)) {
3739566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_FAILED_NANORINF));
3749566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n"));
375c98378a5SBarry Smith       PetscFunctionReturn(0);
376c98378a5SBarry Smith     }
377c427506bSPeter Brune   }
37871b4ebd2SPeter Brune 
37971b4ebd2SPeter Brune   /* copy the solution over */
3809566063dSJacob Faibussowitsch   PetscCall(VecCopy(W, X));
3819566063dSJacob Faibussowitsch   PetscCall(VecCopy(G, F));
3829566063dSJacob Faibussowitsch   PetscCall(VecNorm(X, NORM_2, &xnorm));
3839566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
3849566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm));
38571b4ebd2SPeter Brune   PetscFunctionReturn(0);
38671b4ebd2SPeter Brune }
38771b4ebd2SPeter Brune 
3887f1410a3SPeter Brune PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
3897f1410a3SPeter Brune {
3907f1410a3SPeter Brune   PetscBool         iascii;
391d7ca06c0SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
392075cc632SBarry Smith 
3937f1410a3SPeter Brune   PetscFunctionBegin;
3949566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii));
3957f1410a3SPeter Brune   if (iascii) {
396b000cd8dSPeter Brune     if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3979566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  interpolation: cubic\n"));
398b000cd8dSPeter Brune     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
3999566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  interpolation: quadratic\n"));
4007f1410a3SPeter Brune     }
4019566063dSJacob Faibussowitsch     PetscCall(PetscViewerASCIIPrintf(viewer, "  alpha=%e\n", (double)bt->alpha));
4027f1410a3SPeter Brune   }
4037f1410a3SPeter Brune   PetscFunctionReturn(0);
4047f1410a3SPeter Brune }
4057f1410a3SPeter Brune 
406f1c6b773SPeter Brune static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
407731c067cSJed Brown {
40871b4ebd2SPeter Brune   PetscFunctionBegin;
4099566063dSJacob Faibussowitsch   PetscCall(PetscFree(linesearch->data));
41071b4ebd2SPeter Brune   PetscFunctionReturn(0);
41171b4ebd2SPeter Brune }
41271b4ebd2SPeter Brune 
4134416b707SBarry Smith static PetscErrorCode SNESLineSearchSetFromOptions_BT(PetscOptionItems *PetscOptionsObject,SNESLineSearch linesearch)
41471b4ebd2SPeter Brune {
415eb23ba39SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
41671b4ebd2SPeter Brune 
4176e111a19SKarl Rupp   PetscFunctionBegin;
418*d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject,"SNESLineSearch BT options");
4199566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_linesearch_alpha",   "Descent tolerance",        "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL));
420*d0609cedSBarry Smith   PetscOptionsHeadEnd();
42171b4ebd2SPeter Brune   PetscFunctionReturn(0);
42271b4ebd2SPeter Brune }
42371b4ebd2SPeter Brune 
424954494b2SJed Brown /*MC
425db609ea7SPeter Brune    SNESLINESEARCHBT - Backtracking line search.
426954494b2SJed Brown 
427db609ea7SPeter Brune    This line search finds the minimum of a polynomial fitting of the L2 norm of the
428eb23ba39SBarry 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
429db609ea7SPeter Brune    and the fit is reattempted at most max_it times or until lambda is below minlambda.
430954494b2SJed Brown 
431cd7522eaSPeter Brune    Options Database Keys:
4323eec2f6aSPierre Jolivet +  -snes_linesearch_alpha <1e\-4> - slope descent parameter
433db609ea7SPeter Brune .  -snes_linesearch_damping <1.0> - initial step length
434eb23ba39SBarry Smith .  -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the
435e1bc860dSBarry Smith                                        step is scaled back to be of this length at the beginning of the line search
436db609ea7SPeter Brune .  -snes_linesearch_max_it <40> - maximum number of shrinking step
4373eec2f6aSPierre Jolivet .  -snes_linesearch_minlambda <1e\-12> - minimum step length allowed
438cd7522eaSPeter Brune -  -snes_linesearch_order <cubic,quadratic> - order of the approximation
439cd7522eaSPeter Brune 
440954494b2SJed Brown    Level: advanced
441954494b2SJed Brown 
442db609ea7SPeter Brune    Notes:
443db609ea7SPeter Brune    This line search is taken from "Numerical Methods for Unconstrained
444db609ea7SPeter Brune    Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
445e55accfeSBarryFSmith 
446e55accfeSBarryFSmith    This line search will always produce a step that is less than or equal to, in length, the full step size.
447db609ea7SPeter Brune 
448f1c6b773SPeter Brune .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
449954494b2SJed Brown M*/
4508cc058d9SJed Brown PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
45171b4ebd2SPeter Brune {
45271b4ebd2SPeter Brune 
453f1c6b773SPeter Brune   SNESLineSearch_BT *bt;
45471b4ebd2SPeter Brune 
45571b4ebd2SPeter Brune   PetscFunctionBegin;
456f1c6b773SPeter Brune   linesearch->ops->apply          = SNESLineSearchApply_BT;
457f1c6b773SPeter Brune   linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
458f1c6b773SPeter Brune   linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
4590298fd71SBarry Smith   linesearch->ops->reset          = NULL;
4607f1410a3SPeter Brune   linesearch->ops->view           = SNESLineSearchView_BT;
4610298fd71SBarry Smith   linesearch->ops->setup          = NULL;
46271b4ebd2SPeter Brune 
4639566063dSJacob Faibussowitsch   PetscCall(PetscNewLog(linesearch,&bt));
464f5af7f23SKarl Rupp 
46571b4ebd2SPeter Brune   linesearch->data    = (void*)bt;
46671b4ebd2SPeter Brune   linesearch->max_its = 40;
467b000cd8dSPeter Brune   linesearch->order   = SNES_LINESEARCH_ORDER_CUBIC;
46871b4ebd2SPeter Brune   bt->alpha           = 1e-4;
46971b4ebd2SPeter Brune   PetscFunctionReturn(0);
47071b4ebd2SPeter Brune }
471