xref: /petsc/src/snes/linesearch/impls/bt/linesearchbt.c (revision 08401ef684002a709c6d3db98a0c9f54a8bcf1ec)
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   PetscErrorCode    ierr;
5671b4ebd2SPeter Brune   Vec               X,F,Y,W,G;
5771b4ebd2SPeter Brune   SNES              snes;
588a430afbSPeter Brune   PetscReal         fnorm, xnorm, ynorm, gnorm;
593b288129SPeter Brune   PetscReal         lambda,lambdatemp,lambdaprev,minlambda,maxstep,initslope,alpha,stol;
60dfcd3828SPeter Brune   PetscReal         t1,t2,a,b,d;
618a430afbSPeter Brune   PetscReal         f;
628a430afbSPeter Brune   PetscReal         g,gprev;
6371b4ebd2SPeter Brune   PetscViewer       monitor;
6471b4ebd2SPeter Brune   PetscInt          max_its,count;
65d7ca06c0SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
6671b4ebd2SPeter Brune   Mat               jac;
67bd4a8a71SBarry Smith   PetscErrorCode    (*objective)(SNES,Vec,PetscReal*,void*);
6871b4ebd2SPeter Brune 
6971b4ebd2SPeter Brune   PetscFunctionBegin;
709566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G));
719566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm));
729566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetLambda(linesearch, &lambda));
739566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
749566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor));
759566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchGetTolerances(linesearch,&minlambda,&maxstep,NULL,NULL,NULL,&max_its));
769566063dSJacob Faibussowitsch   PetscCall(SNESGetTolerances(snes,NULL,NULL,&stol,NULL,NULL));
779566063dSJacob Faibussowitsch   PetscCall(SNESGetObjective(snes,&objective,NULL));
7871b4ebd2SPeter Brune   alpha = bt->alpha;
7971b4ebd2SPeter Brune 
809566063dSJacob Faibussowitsch   PetscCall(SNESGetJacobian(snes, &jac, NULL, NULL, NULL));
81*08401ef6SPierre Jolivet   PetscCheck(jac || objective,PetscObjectComm((PetscObject)linesearch), PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix");
82c69d1a72SBarry Smith 
839566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchPreCheck(linesearch,X,Y,&changed_y));
849566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
8571b4ebd2SPeter Brune 
869566063dSJacob Faibussowitsch   PetscCall(VecNormBegin(Y, NORM_2, &ynorm));
879566063dSJacob Faibussowitsch   PetscCall(VecNormBegin(X, NORM_2, &xnorm));
889566063dSJacob Faibussowitsch   PetscCall(VecNormEnd(Y, NORM_2, &ynorm));
899566063dSJacob Faibussowitsch   PetscCall(VecNormEnd(X, NORM_2, &xnorm));
903b288129SPeter Brune 
9171b4ebd2SPeter Brune   if (ynorm == 0.0) {
9271b4ebd2SPeter Brune     if (monitor) {
939566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
949566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Initial direction and size is 0\n"));
959566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
9671b4ebd2SPeter Brune     }
979566063dSJacob Faibussowitsch     PetscCall(VecCopy(X,W));
989566063dSJacob Faibussowitsch     PetscCall(VecCopy(F,G));
999566063dSJacob Faibussowitsch     PetscCall(SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm));
1009566063dSJacob Faibussowitsch     PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
10171b4ebd2SPeter Brune     PetscFunctionReturn(0);
10271b4ebd2SPeter Brune   }
10371b4ebd2SPeter Brune   if (ynorm > maxstep) {        /* Step too big, so scale back */
10471b4ebd2SPeter Brune     if (monitor) {
1059566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
1069566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Scaling step by %14.12e old ynorm %14.12e\n", (double)(maxstep/ynorm),(double)ynorm));
1079566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
10871b4ebd2SPeter Brune     }
1099566063dSJacob Faibussowitsch     PetscCall(VecScale(Y,maxstep/(ynorm)));
11071b4ebd2SPeter Brune     ynorm = maxstep;
11171b4ebd2SPeter Brune   }
1128a430afbSPeter Brune 
1138a430afbSPeter Brune   /* if the SNES has an objective set, use that instead of the function value */
114bd4a8a71SBarry Smith   if (objective) {
1159566063dSJacob Faibussowitsch     PetscCall(SNESComputeObjective(snes,X,&f));
1168a430afbSPeter Brune   } else {
1178a430afbSPeter Brune     f = fnorm*fnorm;
1188a430afbSPeter Brune   }
1198a430afbSPeter Brune 
1208a430afbSPeter Brune   /* compute the initial slope */
121bd4a8a71SBarry Smith   if (objective) {
1228a430afbSPeter Brune     /* slope comes from the function (assumed to be the gradient of the objective */
1239566063dSJacob Faibussowitsch     PetscCall(VecDotRealPart(Y,F,&initslope));
1248a430afbSPeter Brune   } else {
1258a430afbSPeter Brune     /* slope comes from the normal equations */
1269566063dSJacob Faibussowitsch     PetscCall(MatMult(jac,Y,W));
1279566063dSJacob Faibussowitsch     PetscCall(VecDotRealPart(F,W,&initslope));
12871b4ebd2SPeter Brune     if (initslope > 0.0)  initslope = -initslope;
12971b4ebd2SPeter Brune     if (initslope == 0.0) initslope = -1.0;
1308a430afbSPeter Brune   }
13171b4ebd2SPeter Brune 
132df019d78SBarry Smith   while (PETSC_TRUE) {
1339566063dSJacob Faibussowitsch     PetscCall(VecWAXPY(W,-lambda,Y,X));
1349bd66eb0SPeter Brune     if (linesearch->ops->viproject) {
1359566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->viproject)(snes, W));
1369bd66eb0SPeter Brune     }
137e71169deSBarry Smith     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
1389566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n"));
13971b4ebd2SPeter Brune       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
1409566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
14171b4ebd2SPeter Brune       PetscFunctionReturn(0);
14271b4ebd2SPeter Brune     }
1438a430afbSPeter Brune 
144bd4a8a71SBarry Smith     if (objective) {
1459566063dSJacob Faibussowitsch       PetscCall(SNESComputeObjective(snes,W,&g));
1468a430afbSPeter Brune     } else {
1479566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->snesfunc)(snes,W,G));
1489bd66eb0SPeter Brune       if (linesearch->ops->vinorm) {
1499bd66eb0SPeter Brune         gnorm = fnorm;
1509566063dSJacob Faibussowitsch         PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
1519bd66eb0SPeter Brune       } else {
1529566063dSJacob Faibussowitsch         PetscCall(VecNorm(G,NORM_2,&gnorm));
1538a430afbSPeter Brune       }
15479e7e81bSJed Brown       g = PetscSqr(gnorm);
1559bd66eb0SPeter Brune     }
1569566063dSJacob Faibussowitsch     PetscCall(SNESLineSearchMonitor(linesearch));
1579bd66eb0SPeter Brune 
158df019d78SBarry Smith     if (!PetscIsInfOrNanReal(g)) break;
159df019d78SBarry Smith     if (monitor) {
1609566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
1619566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n",(double)lambda));
1629566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
163df019d78SBarry Smith     }
164df019d78SBarry Smith     if (lambda <= minlambda) {
1657e49c775SBarry Smith       SNESCheckFunctionNorm(snes,g);
166c98378a5SBarry Smith     }
167df019d78SBarry Smith     lambda = .5*lambda;
168df019d78SBarry Smith   }
169df019d78SBarry Smith 
170bd4a8a71SBarry Smith   if (!objective) {
1719566063dSJacob Faibussowitsch     PetscCall(PetscInfo(snes,"Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm));
1728a430afbSPeter Brune   }
1738a430afbSPeter Brune   if (.5*g <= .5*f + lambda*alpha*initslope) { /* Sufficient reduction or step tolerance convergence */
17471b4ebd2SPeter Brune     if (monitor) {
1759566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
176bd4a8a71SBarry Smith       if (!objective) {
1779566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm));
1788a430afbSPeter Brune       } else {
1799566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Using full step: obj %14.12e obj %14.12e\n", (double)f, (double)g));
1808a430afbSPeter Brune       }
1819566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
18271b4ebd2SPeter Brune     }
18371b4ebd2SPeter Brune   } else {
184c21ba15cSPeter Brune     /* Since the full step didn't work and the step is tiny, quit */
185c21ba15cSPeter Brune     if (stol*xnorm > ynorm) {
1869566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm));
1879566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
188c61ba1e2SPeter Brune       if (monitor) {
1899566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
1909566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n",(double)ynorm,(double)stol*xnorm));
1919566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
192c61ba1e2SPeter Brune       }
193c21ba15cSPeter Brune       PetscFunctionReturn(0);
194c21ba15cSPeter Brune     }
19571b4ebd2SPeter Brune     /* Fit points with quadratic */
1968a430afbSPeter Brune     lambdatemp = -initslope/(g - f - 2.0*lambda*initslope);
19771b4ebd2SPeter Brune     lambdaprev = lambda;
1988a430afbSPeter Brune     gprev      = g;
19971b4ebd2SPeter Brune     if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
20071b4ebd2SPeter Brune     if (lambdatemp <= .1*lambda) lambda = .1*lambda;
20171b4ebd2SPeter Brune     else                         lambda = lambdatemp;
20271b4ebd2SPeter Brune 
2039566063dSJacob Faibussowitsch     PetscCall(VecWAXPY(W,-lambda,Y,X));
2049bd66eb0SPeter Brune     if (linesearch->ops->viproject) {
2059566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->viproject)(snes, W));
2069bd66eb0SPeter Brune     }
207e71169deSBarry Smith     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
2089566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs));
20971b4ebd2SPeter Brune       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
2109566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
21171b4ebd2SPeter Brune       PetscFunctionReturn(0);
21271b4ebd2SPeter Brune     }
213bd4a8a71SBarry Smith     if (objective) {
2149566063dSJacob Faibussowitsch       PetscCall(SNESComputeObjective(snes,W,&g));
2158a430afbSPeter Brune     } else {
2169566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->snesfunc)(snes,W,G));
2179bd66eb0SPeter Brune       if (linesearch->ops->vinorm) {
2189bd66eb0SPeter Brune         gnorm = fnorm;
2199566063dSJacob Faibussowitsch         PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
2209bd66eb0SPeter Brune       } else {
2219566063dSJacob Faibussowitsch         PetscCall(VecNorm(G,NORM_2,&gnorm));
2229bd66eb0SPeter Brune       }
223f1ab69e3SPeter Brune       g = PetscSqr(gnorm);
2248a430afbSPeter Brune     }
225c98378a5SBarry Smith     if (PetscIsInfOrNanReal(g)) {
2269566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
2279566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n"));
228c98378a5SBarry Smith       PetscFunctionReturn(0);
229c98378a5SBarry Smith     }
23071b4ebd2SPeter Brune     if (monitor) {
2319566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
232bd4a8a71SBarry Smith       if (!objective) {
2339566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: gnorm after quadratic fit %14.12e\n",(double)gnorm));
2348a430afbSPeter Brune       } else {
2359566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: obj after quadratic fit %14.12e\n",(double)g));
2368a430afbSPeter Brune       }
2379566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
23871b4ebd2SPeter Brune     }
2398a430afbSPeter Brune     if (.5*g < .5*f + lambda*alpha*initslope) { /* sufficient reduction */
24071b4ebd2SPeter Brune       if (monitor) {
2419566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
2429566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda));
2439566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
24471b4ebd2SPeter Brune       }
24571b4ebd2SPeter Brune     } else {
24671b4ebd2SPeter Brune       /* Fit points with cubic */
24771b4ebd2SPeter Brune       for (count = 0; count < max_its; count++) {
24871b4ebd2SPeter Brune         if (lambda <= minlambda) {
24971b4ebd2SPeter Brune           if (monitor) {
2509566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
2519566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: unable to find good step length! After %D tries \n",count));
252bd4a8a71SBarry Smith             if (!objective) {
2534a2f8832SBarry Smith               ierr = PetscViewerASCIIPrintf(monitor,"    Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
2549566063dSJacob Faibussowitsch                                                          (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);PetscCall(ierr);
2558a430afbSPeter Brune             } else {
2564a2f8832SBarry Smith               ierr = 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",
2579566063dSJacob Faibussowitsch                                                          (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);PetscCall(ierr);
2588a430afbSPeter Brune             }
2599566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
26071b4ebd2SPeter Brune           }
2619566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
26271b4ebd2SPeter Brune           PetscFunctionReturn(0);
26371b4ebd2SPeter Brune         }
264b000cd8dSPeter Brune         if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
2658a430afbSPeter Brune           t1 = .5*(g - f) - lambda*initslope;
2668a430afbSPeter Brune           t2 = .5*(gprev  - f) - lambdaprev*initslope;
26771b4ebd2SPeter Brune           a  = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
26871b4ebd2SPeter Brune           b  = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
26971b4ebd2SPeter Brune           d  = b*b - 3*a*initslope;
27071b4ebd2SPeter Brune           if (d < 0.0) d = 0.0;
271f5af7f23SKarl Rupp           if (a == 0.0) lambdatemp = -initslope/(2.0*b);
272f5af7f23SKarl Rupp           else lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a);
273f5af7f23SKarl Rupp 
274b000cd8dSPeter Brune         } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
2758a430afbSPeter Brune           lambdatemp = -initslope/(g - f - 2.0*initslope);
276ce94432eSBarry Smith         } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt");
27771b4ebd2SPeter Brune         lambdaprev = lambda;
2788a430afbSPeter Brune         gprev      = g;
27971b4ebd2SPeter Brune         if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
28071b4ebd2SPeter Brune         if (lambdatemp <= .1*lambda) lambda     = .1*lambda;
28171b4ebd2SPeter Brune         else                         lambda     = lambdatemp;
2829566063dSJacob Faibussowitsch         PetscCall(VecWAXPY(W,-lambda,Y,X));
283b7b2e573SPeter Brune         if (linesearch->ops->viproject) {
2849566063dSJacob Faibussowitsch           PetscCall((*linesearch->ops->viproject)(snes,W));
285b7b2e573SPeter Brune         }
286e71169deSBarry Smith         if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
2879566063dSJacob Faibussowitsch           PetscCall(PetscInfo(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count));
288bd4a8a71SBarry Smith           if (!objective) {
2897d3de750SJacob Faibussowitsch             ierr = PetscInfo(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
2909566063dSJacob Faibussowitsch                               (double)fnorm,(double)gnorm,(double)ynorm,(double)lambda,(double)initslope);PetscCall(ierr);
2918a430afbSPeter Brune           }
2929566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
29371b4ebd2SPeter Brune           snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
29471b4ebd2SPeter Brune           PetscFunctionReturn(0);
29571b4ebd2SPeter Brune         }
296bd4a8a71SBarry Smith         if (objective) {
2979566063dSJacob Faibussowitsch           PetscCall(SNESComputeObjective(snes,W,&g));
2988a430afbSPeter Brune         } else {
2999566063dSJacob Faibussowitsch           PetscCall((*linesearch->ops->snesfunc)(snes,W,G));
3009bd66eb0SPeter Brune           if (linesearch->ops->vinorm) {
3019bd66eb0SPeter Brune             gnorm = fnorm;
3029566063dSJacob Faibussowitsch             PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
3039bd66eb0SPeter Brune           } else {
3049566063dSJacob Faibussowitsch             PetscCall(VecNorm(G,NORM_2,&gnorm));
3058a430afbSPeter Brune           }
306f1ab69e3SPeter Brune           g = PetscSqr(gnorm);
3079bd66eb0SPeter Brune         }
3084a2f8832SBarry Smith         if (PetscIsInfOrNanReal(g)) {
3099566063dSJacob Faibussowitsch           PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
3109566063dSJacob Faibussowitsch           PetscCall(PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n"));
311c98378a5SBarry Smith           PetscFunctionReturn(0);
312c98378a5SBarry Smith         }
3138a430afbSPeter Brune         if (.5*g < .5*f + lambda*alpha*initslope) { /* is reduction enough? */
31471b4ebd2SPeter Brune           if (monitor) {
3159566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
316bd4a8a71SBarry Smith             if (!objective) {
317b000cd8dSPeter Brune               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3189566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda));
31971b4ebd2SPeter Brune               } else {
3209566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda));
32171b4ebd2SPeter Brune               }
3229566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
3238a430afbSPeter Brune             } else {
3248a430afbSPeter Brune               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3259566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda));
3268a430afbSPeter Brune               } else {
3279566063dSJacob Faibussowitsch                 PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda));
3288a430afbSPeter Brune               }
3299566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
3308a430afbSPeter Brune             }
33171b4ebd2SPeter Brune           }
33271b4ebd2SPeter Brune           break;
333f5af7f23SKarl Rupp         } else if (monitor) {
3349566063dSJacob Faibussowitsch           PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel));
335bd4a8a71SBarry Smith           if (!objective) {
336b000cd8dSPeter Brune             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3379566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda));
33871b4ebd2SPeter Brune             } else {
3399566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda));
34071b4ebd2SPeter Brune             }
3419566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
3428a430afbSPeter Brune           } else {
3438a430afbSPeter Brune             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3449566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda));
3458a430afbSPeter Brune             } else {
3469566063dSJacob Faibussowitsch               PetscCall(PetscViewerASCIIPrintf(monitor,"    Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda));
3478a430afbSPeter Brune             }
3489566063dSJacob Faibussowitsch             PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel));
3498a430afbSPeter Brune           }
35071b4ebd2SPeter Brune         }
35171b4ebd2SPeter Brune       }
35271b4ebd2SPeter Brune     }
35371b4ebd2SPeter Brune   }
35471b4ebd2SPeter Brune 
35571b4ebd2SPeter Brune   /* postcheck */
3560c1ac483SGlenn Hammond   /* update Y to lambda*Y so that W is consistent with  X - lambda*Y */
3579566063dSJacob Faibussowitsch   PetscCall(VecScale(Y,lambda));
3589566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w));
35971b4ebd2SPeter Brune   if (changed_y) {
3609566063dSJacob Faibussowitsch     PetscCall(VecWAXPY(W,-1.0,Y,X));
3619bd66eb0SPeter Brune     if (linesearch->ops->viproject) {
3629566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->viproject)(snes, W));
3639bd66eb0SPeter Brune     }
36471b4ebd2SPeter Brune   }
365bd4a8a71SBarry Smith   if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
3669566063dSJacob Faibussowitsch     PetscCall((*linesearch->ops->snesfunc)(snes,W,G));
3679bd66eb0SPeter Brune     if (linesearch->ops->vinorm) {
3689bd66eb0SPeter Brune       gnorm = fnorm;
3699566063dSJacob Faibussowitsch       PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm));
3709bd66eb0SPeter Brune     } else {
3719566063dSJacob Faibussowitsch       PetscCall(VecNorm(G,NORM_2,&gnorm));
3729bd66eb0SPeter Brune     }
3739566063dSJacob Faibussowitsch     PetscCall(VecNorm(Y,NORM_2,&ynorm));
374c98378a5SBarry Smith     if (PetscIsInfOrNanReal(gnorm)) {
3759566063dSJacob Faibussowitsch       PetscCall(SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_FAILED_NANORINF));
3769566063dSJacob Faibussowitsch       PetscCall(PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n"));
377c98378a5SBarry Smith       PetscFunctionReturn(0);
378c98378a5SBarry Smith     }
379c427506bSPeter Brune   }
38071b4ebd2SPeter Brune 
38171b4ebd2SPeter Brune   /* copy the solution over */
3829566063dSJacob Faibussowitsch   PetscCall(VecCopy(W, X));
3839566063dSJacob Faibussowitsch   PetscCall(VecCopy(G, F));
3849566063dSJacob Faibussowitsch   PetscCall(VecNorm(X, NORM_2, &xnorm));
3859566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
3869566063dSJacob Faibussowitsch   PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm));
38771b4ebd2SPeter Brune   PetscFunctionReturn(0);
38871b4ebd2SPeter Brune }
38971b4ebd2SPeter Brune 
3907f1410a3SPeter Brune PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
3917f1410a3SPeter Brune {
3927f1410a3SPeter Brune   PetscBool         iascii;
393d7ca06c0SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
394075cc632SBarry Smith 
3957f1410a3SPeter Brune   PetscFunctionBegin;
3969566063dSJacob Faibussowitsch   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii));
3977f1410a3SPeter Brune   if (iascii) {
398b000cd8dSPeter Brune     if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3999566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  interpolation: cubic\n"));
400b000cd8dSPeter Brune     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
4019566063dSJacob Faibussowitsch       PetscCall(PetscViewerASCIIPrintf(viewer, "  interpolation: quadratic\n"));
4027f1410a3SPeter Brune     }
4039566063dSJacob Faibussowitsch     PetscCall(PetscViewerASCIIPrintf(viewer, "  alpha=%e\n", (double)bt->alpha));
4047f1410a3SPeter Brune   }
4057f1410a3SPeter Brune   PetscFunctionReturn(0);
4067f1410a3SPeter Brune }
4077f1410a3SPeter Brune 
408f1c6b773SPeter Brune static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
409731c067cSJed Brown {
41071b4ebd2SPeter Brune   PetscFunctionBegin;
4119566063dSJacob Faibussowitsch   PetscCall(PetscFree(linesearch->data));
41271b4ebd2SPeter Brune   PetscFunctionReturn(0);
41371b4ebd2SPeter Brune }
41471b4ebd2SPeter Brune 
4154416b707SBarry Smith static PetscErrorCode SNESLineSearchSetFromOptions_BT(PetscOptionItems *PetscOptionsObject,SNESLineSearch linesearch)
41671b4ebd2SPeter Brune {
417eb23ba39SBarry Smith   SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
41871b4ebd2SPeter Brune 
4196e111a19SKarl Rupp   PetscFunctionBegin;
4209566063dSJacob Faibussowitsch   PetscCall(PetscOptionsHead(PetscOptionsObject,"SNESLineSearch BT options"));
4219566063dSJacob Faibussowitsch   PetscCall(PetscOptionsReal("-snes_linesearch_alpha",   "Descent tolerance",        "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL));
4229566063dSJacob Faibussowitsch   PetscCall(PetscOptionsTail());
42371b4ebd2SPeter Brune   PetscFunctionReturn(0);
42471b4ebd2SPeter Brune }
42571b4ebd2SPeter Brune 
426954494b2SJed Brown /*MC
427db609ea7SPeter Brune    SNESLINESEARCHBT - Backtracking line search.
428954494b2SJed Brown 
429db609ea7SPeter Brune    This line search finds the minimum of a polynomial fitting of the L2 norm of the
430eb23ba39SBarry 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
431db609ea7SPeter Brune    and the fit is reattempted at most max_it times or until lambda is below minlambda.
432954494b2SJed Brown 
433cd7522eaSPeter Brune    Options Database Keys:
4343eec2f6aSPierre Jolivet +  -snes_linesearch_alpha <1e\-4> - slope descent parameter
435db609ea7SPeter Brune .  -snes_linesearch_damping <1.0> - initial step length
436eb23ba39SBarry Smith .  -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the
437e1bc860dSBarry Smith                                        step is scaled back to be of this length at the beginning of the line search
438db609ea7SPeter Brune .  -snes_linesearch_max_it <40> - maximum number of shrinking step
4393eec2f6aSPierre Jolivet .  -snes_linesearch_minlambda <1e\-12> - minimum step length allowed
440cd7522eaSPeter Brune -  -snes_linesearch_order <cubic,quadratic> - order of the approximation
441cd7522eaSPeter Brune 
442954494b2SJed Brown    Level: advanced
443954494b2SJed Brown 
444db609ea7SPeter Brune    Notes:
445db609ea7SPeter Brune    This line search is taken from "Numerical Methods for Unconstrained
446db609ea7SPeter Brune    Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
447e55accfeSBarryFSmith 
448e55accfeSBarryFSmith    This line search will always produce a step that is less than or equal to, in length, the full step size.
449db609ea7SPeter Brune 
450f1c6b773SPeter Brune .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
451954494b2SJed Brown M*/
4528cc058d9SJed Brown PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
45371b4ebd2SPeter Brune {
45471b4ebd2SPeter Brune 
455f1c6b773SPeter Brune   SNESLineSearch_BT *bt;
45671b4ebd2SPeter Brune 
45771b4ebd2SPeter Brune   PetscFunctionBegin;
458f1c6b773SPeter Brune   linesearch->ops->apply          = SNESLineSearchApply_BT;
459f1c6b773SPeter Brune   linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
460f1c6b773SPeter Brune   linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
4610298fd71SBarry Smith   linesearch->ops->reset          = NULL;
4627f1410a3SPeter Brune   linesearch->ops->view           = SNESLineSearchView_BT;
4630298fd71SBarry Smith   linesearch->ops->setup          = NULL;
46471b4ebd2SPeter Brune 
4659566063dSJacob Faibussowitsch   PetscCall(PetscNewLog(linesearch,&bt));
466f5af7f23SKarl Rupp 
46771b4ebd2SPeter Brune   linesearch->data    = (void*)bt;
46871b4ebd2SPeter Brune   linesearch->max_its = 40;
469b000cd8dSPeter Brune   linesearch->order   = SNES_LINESEARCH_ORDER_CUBIC;
47071b4ebd2SPeter Brune   bt->alpha           = 1e-4;
47171b4ebd2SPeter Brune   PetscFunctionReturn(0);
47271b4ebd2SPeter Brune }
473