xref: /petsc/src/snes/linesearch/impls/bt/linesearchbt.c (revision c21ba15ce39d311a1cad44a785b260c339416f2e)
13d443117SPeter Brune #include <petsc-private/linesearchimpl.h> /*I  "petscsnes.h"  I*/
2b45d2f2cSJed Brown #include <petsc-private/snesimpl.h>
371b4ebd2SPeter Brune 
471b4ebd2SPeter Brune typedef struct {
571b4ebd2SPeter Brune   PetscReal        alpha; /* sufficient decrease parameter */
6f1c6b773SPeter Brune } SNESLineSearch_BT;
771b4ebd2SPeter Brune 
871b4ebd2SPeter Brune #undef __FUNCT__
92f4102e2SPeter Brune #define __FUNCT__ "SNESLineSearchBTSetAlpha"
102f4102e2SPeter Brune /*@
112f4102e2SPeter Brune    SNESLineSearchBTSetAlpha - Sets the descent parameter, alpha, in the BT linesearch variant.
122f4102e2SPeter Brune 
132f4102e2SPeter Brune    Input Parameters:
142f4102e2SPeter Brune +  linesearch - linesearch context
152f4102e2SPeter Brune -  alpha - The descent parameter
162f4102e2SPeter Brune 
172f4102e2SPeter Brune    Level: intermediate
182f4102e2SPeter Brune 
191a4f838cSPeter Brune .seealso: SNESLineSearchSetLambda(), SNESLineSearchGetTolerances() SNESLINESEARCHBT
202f4102e2SPeter Brune @*/
212f4102e2SPeter Brune PetscErrorCode SNESLineSearchBTSetAlpha(SNESLineSearch linesearch, PetscReal alpha)
222f4102e2SPeter Brune {
232f4102e2SPeter Brune   SNESLineSearch_BT  *bt;
242f4102e2SPeter Brune   PetscFunctionBegin;
252f4102e2SPeter Brune   PetscValidHeaderSpecific(linesearch,SNESLINESEARCH_CLASSID,1);
262f4102e2SPeter Brune   bt = (SNESLineSearch_BT *)linesearch->data;
272f4102e2SPeter Brune   bt->alpha = alpha;
282f4102e2SPeter Brune   PetscFunctionReturn(0);
292f4102e2SPeter Brune }
302f4102e2SPeter Brune 
312f4102e2SPeter Brune 
322f4102e2SPeter Brune #undef __FUNCT__
332f4102e2SPeter Brune #define __FUNCT__ "SNESLineSearchBTGetAlpha"
342f4102e2SPeter Brune /*@
352f4102e2SPeter Brune    SNESLineSearchBTGetAlpha - Gets the descent parameter, alpha, in the BT linesearch variant.
362f4102e2SPeter Brune 
372f4102e2SPeter Brune    Input Parameters:
382f4102e2SPeter Brune .  linesearch - linesearch context
398e557f58SPeter Brune 
408e557f58SPeter Brune    Output Parameters:
412f4102e2SPeter Brune .  alpha - The descent parameter
422f4102e2SPeter Brune 
432f4102e2SPeter Brune    Level: intermediate
442f4102e2SPeter Brune 
451a4f838cSPeter Brune .seealso: SNESLineSearchGetLambda(), SNESLineSearchGetTolerances() SNESLINESEARCHBT
462f4102e2SPeter Brune @*/
472f4102e2SPeter Brune PetscErrorCode SNESLineSearchBTGetAlpha(SNESLineSearch linesearch, PetscReal *alpha)
482f4102e2SPeter Brune {
492f4102e2SPeter Brune   SNESLineSearch_BT  *bt;
502f4102e2SPeter Brune   PetscFunctionBegin;
512f4102e2SPeter Brune   PetscValidHeaderSpecific(linesearch,SNESLINESEARCH_CLASSID,1);
522f4102e2SPeter Brune   bt = (SNESLineSearch_BT *)linesearch->data;
532f4102e2SPeter Brune   *alpha = bt->alpha;
542f4102e2SPeter Brune   PetscFunctionReturn(0);
552f4102e2SPeter Brune }
562f4102e2SPeter Brune 
572f4102e2SPeter Brune #undef __FUNCT__
58f1c6b773SPeter Brune #define __FUNCT__ "SNESLineSearchApply_BT"
59f1c6b773SPeter Brune static PetscErrorCode  SNESLineSearchApply_BT(SNESLineSearch linesearch)
6071b4ebd2SPeter Brune {
6171b4ebd2SPeter Brune   PetscBool      changed_y,changed_w;
6271b4ebd2SPeter Brune   PetscErrorCode ierr;
6371b4ebd2SPeter Brune   Vec            X,F,Y,W,G;
6471b4ebd2SPeter Brune   SNES           snes;
6571b4ebd2SPeter Brune   PetscReal      fnorm, xnorm, ynorm, gnorm, gnormprev;
663b288129SPeter Brune   PetscReal      lambda,lambdatemp,lambdaprev,minlambda,maxstep,initslope,alpha,stol;
67dfcd3828SPeter Brune   PetscReal      t1,t2,a,b,d;
6871b4ebd2SPeter Brune #if defined(PETSC_USE_COMPLEX)
6971b4ebd2SPeter Brune   PetscScalar    cinitslope;
7071b4ebd2SPeter Brune #endif
7171b4ebd2SPeter Brune   PetscBool      domainerror;
7271b4ebd2SPeter Brune   PetscViewer    monitor;
7371b4ebd2SPeter Brune   PetscInt       max_its,count;
74f1c6b773SPeter Brune   SNESLineSearch_BT  *bt;
7571b4ebd2SPeter Brune   Mat            jac;
7671b4ebd2SPeter Brune 
7771b4ebd2SPeter Brune 
7871b4ebd2SPeter Brune   PetscFunctionBegin;
7971b4ebd2SPeter Brune 
80f1c6b773SPeter Brune   ierr = SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G);CHKERRQ(ierr);
81f1c6b773SPeter Brune   ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr);
82f1c6b773SPeter Brune   ierr = SNESLineSearchGetLambda(linesearch, &lambda);CHKERRQ(ierr);
83f1c6b773SPeter Brune   ierr = SNESLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr);
84f1c6b773SPeter Brune   ierr = SNESLineSearchGetMonitor(linesearch, &monitor);CHKERRQ(ierr);
85dfcd3828SPeter Brune   ierr = SNESLineSearchGetTolerances(linesearch,&minlambda,&maxstep,PETSC_NULL,PETSC_NULL,PETSC_NULL,&max_its);
863b288129SPeter Brune   ierr = SNESGetTolerances(snes,PETSC_NULL,PETSC_NULL,&stol,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
87f1c6b773SPeter Brune   bt = (SNESLineSearch_BT *)linesearch->data;
8871b4ebd2SPeter Brune 
8971b4ebd2SPeter Brune   alpha = bt->alpha;
9071b4ebd2SPeter Brune 
9171b4ebd2SPeter Brune   ierr = SNESGetJacobian(snes, &jac, PETSC_NULL, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr);
92c69d1a72SBarry Smith   if (!jac) SETERRQ(((PetscObject)linesearch)->comm, PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix");
93c69d1a72SBarry Smith 
9471b4ebd2SPeter Brune   /* precheck */
957b1df9c1SPeter Brune   ierr = SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);CHKERRQ(ierr);
96f1c6b773SPeter Brune   ierr = SNESLineSearchSetSuccess(linesearch, PETSC_TRUE);CHKERRQ(ierr);
9771b4ebd2SPeter Brune 
983b288129SPeter Brune   ierr = VecNormBegin(Y, NORM_2, &ynorm);CHKERRQ(ierr);
993b288129SPeter Brune   ierr = VecNormBegin(X, NORM_2, &xnorm);CHKERRQ(ierr);
1003b288129SPeter Brune   ierr = VecNormEnd(Y, NORM_2, &ynorm);CHKERRQ(ierr);
1013b288129SPeter Brune   ierr = VecNormEnd(X, NORM_2, &xnorm);CHKERRQ(ierr);
1023b288129SPeter Brune 
10371b4ebd2SPeter Brune   if (ynorm == 0.0) {
10471b4ebd2SPeter Brune     if (monitor) {
105ea5d4fccSPeter Brune       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
10671b4ebd2SPeter Brune       ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Initial direction and size is 0\n");CHKERRQ(ierr);
107ea5d4fccSPeter Brune       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
10871b4ebd2SPeter Brune     }
10971b4ebd2SPeter Brune     ierr   = VecCopy(X,W);CHKERRQ(ierr);
11071b4ebd2SPeter Brune     ierr   = VecCopy(F,G);CHKERRQ(ierr);
111f1c6b773SPeter Brune     ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
11271b4ebd2SPeter Brune     PetscFunctionReturn(0);
11371b4ebd2SPeter Brune   }
11471b4ebd2SPeter Brune   if (ynorm > maxstep) {	/* Step too big, so scale back */
11571b4ebd2SPeter Brune     if (monitor) {
116ea5d4fccSPeter Brune       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
117c69d1a72SBarry Smith       ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Scaling step by %14.12e old ynorm %14.12e\n", (double)(maxstep/ynorm),(double)ynorm);CHKERRQ(ierr);
118ea5d4fccSPeter Brune       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
11971b4ebd2SPeter Brune     }
12071b4ebd2SPeter Brune     ierr = VecScale(Y,maxstep/(ynorm));CHKERRQ(ierr);
12171b4ebd2SPeter Brune     ynorm = maxstep;
12271b4ebd2SPeter Brune   }
12371b4ebd2SPeter Brune   ierr      = MatMult(jac,Y,W);CHKERRQ(ierr);
12471b4ebd2SPeter Brune #if defined(PETSC_USE_COMPLEX)
12571b4ebd2SPeter Brune   ierr      = VecDot(F,W,&cinitslope);CHKERRQ(ierr);
12671b4ebd2SPeter Brune   initslope = PetscRealPart(cinitslope);
12771b4ebd2SPeter Brune #else
12871b4ebd2SPeter Brune   ierr      = VecDot(F,W,&initslope);CHKERRQ(ierr);
12971b4ebd2SPeter Brune #endif
13071b4ebd2SPeter Brune   if (initslope > 0.0)  initslope = -initslope;
13171b4ebd2SPeter Brune   if (initslope == 0.0) initslope = -1.0;
13271b4ebd2SPeter Brune 
133c427506bSPeter Brune   ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
1349bd66eb0SPeter Brune   if (linesearch->ops->viproject) {
1359bd66eb0SPeter Brune     ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
1369bd66eb0SPeter Brune   }
13771b4ebd2SPeter Brune   if (snes->nfuncs >= snes->max_funcs) {
13871b4ebd2SPeter Brune     ierr  = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr);
13971b4ebd2SPeter Brune     snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
140f1c6b773SPeter Brune     ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
14171b4ebd2SPeter Brune     PetscFunctionReturn(0);
14271b4ebd2SPeter Brune   }
14371b4ebd2SPeter Brune   ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr);
144c427506bSPeter Brune   ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
145c427506bSPeter Brune   if (domainerror) {
146f1c6b773SPeter Brune     ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
14771b4ebd2SPeter Brune     PetscFunctionReturn(0);
14871b4ebd2SPeter Brune   }
1499bd66eb0SPeter Brune   if (linesearch->ops->vinorm) {
1509bd66eb0SPeter Brune     gnorm = fnorm;
1519bd66eb0SPeter Brune     ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
1529bd66eb0SPeter Brune   } else {
15371b4ebd2SPeter Brune     ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
1549bd66eb0SPeter Brune   }
1559bd66eb0SPeter Brune 
15671b4ebd2SPeter Brune   if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
157c69d1a72SBarry Smith   ierr = PetscInfo2(snes,"Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);CHKERRQ(ierr);
158782c0b7fSBarry Smith   if (.5*gnorm*gnorm <= .5*fnorm*fnorm + lambda*alpha*initslope) { /* Sufficient reduction or step tolerance convergence */
15971b4ebd2SPeter Brune     if (monitor) {
160ea5d4fccSPeter Brune       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
161c69d1a72SBarry Smith       ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);CHKERRQ(ierr);
162ea5d4fccSPeter Brune       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
16371b4ebd2SPeter Brune     }
16471b4ebd2SPeter Brune   } else {
165*c21ba15cSPeter Brune     /* Since the full step didn't work and the step is tiny, quit */
166*c21ba15cSPeter Brune     if (stol*xnorm > ynorm) {
167*c21ba15cSPeter Brune       ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
168*c21ba15cSPeter Brune       ierr = PetscInfo2(monitor,"Aborted due to ynorm < stol*xnorm (%14.12e < %14.12e) and inadequate full step.\n",ynorm,stol*xnorm);CHKERRQ(ierr);
169*c21ba15cSPeter Brune       PetscFunctionReturn(0);
170*c21ba15cSPeter Brune     }
17171b4ebd2SPeter Brune     /* Fit points with quadratic */
17271b4ebd2SPeter Brune     lambdatemp = -initslope/(gnorm*gnorm - fnorm*fnorm - 2.0*lambda*initslope);
17371b4ebd2SPeter Brune     lambdaprev = lambda;
17471b4ebd2SPeter Brune     gnormprev  = gnorm;
17571b4ebd2SPeter Brune     if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
17671b4ebd2SPeter Brune     if (lambdatemp <= .1*lambda) lambda = .1*lambda;
17771b4ebd2SPeter Brune     else                         lambda = lambdatemp;
17871b4ebd2SPeter Brune 
17971b4ebd2SPeter Brune     ierr  = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
1809bd66eb0SPeter Brune     if (linesearch->ops->viproject) {
1819bd66eb0SPeter Brune       ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
1829bd66eb0SPeter Brune     }
18371b4ebd2SPeter Brune     if (snes->nfuncs >= snes->max_funcs) {
18471b4ebd2SPeter Brune       ierr  = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr);
18571b4ebd2SPeter Brune       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
186f1c6b773SPeter Brune       ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
18771b4ebd2SPeter Brune       PetscFunctionReturn(0);
18871b4ebd2SPeter Brune     }
18971b4ebd2SPeter Brune     ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr);
1906188f407SPeter Brune     ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
1916188f407SPeter Brune     if (domainerror) {
19271b4ebd2SPeter Brune       PetscFunctionReturn(0);
19371b4ebd2SPeter Brune     }
1949bd66eb0SPeter Brune     if (linesearch->ops->vinorm) {
1959bd66eb0SPeter Brune       gnorm = fnorm;
1969bd66eb0SPeter Brune       ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
1979bd66eb0SPeter Brune     } else {
19871b4ebd2SPeter Brune       ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
1999bd66eb0SPeter Brune     }
20071b4ebd2SPeter Brune     if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
20171b4ebd2SPeter Brune     if (monitor) {
202ea5d4fccSPeter Brune       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
203c69d1a72SBarry Smith       ierr = PetscViewerASCIIPrintf(monitor,"    Line search: gnorm after quadratic fit %14.12e\n",(double)gnorm);CHKERRQ(ierr);
204ea5d4fccSPeter Brune       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
20571b4ebd2SPeter Brune     }
20671b4ebd2SPeter Brune     if (.5*gnorm*gnorm < .5*fnorm*fnorm + lambda*alpha*initslope) { /* sufficient reduction */
20771b4ebd2SPeter Brune       if (monitor) {
208ea5d4fccSPeter Brune         ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
20971b4ebd2SPeter Brune         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda);CHKERRQ(ierr);
210ea5d4fccSPeter Brune         ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
21171b4ebd2SPeter Brune       }
21271b4ebd2SPeter Brune     } else {
21371b4ebd2SPeter Brune       /* Fit points with cubic */
21471b4ebd2SPeter Brune       for (count = 0; count < max_its; count++) {
21571b4ebd2SPeter Brune         if (lambda <= minlambda) {
21671b4ebd2SPeter Brune           if (monitor) {
217ea5d4fccSPeter Brune             ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
21871b4ebd2SPeter Brune             ierr = PetscViewerASCIIPrintf(monitor,"    Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr);
21971b4ebd2SPeter Brune             ierr = PetscViewerASCIIPrintf(monitor,
22071b4ebd2SPeter Brune                                           "    Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
221c69d1a72SBarry Smith                                           (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);CHKERRQ(ierr);
222ea5d4fccSPeter Brune             ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
22371b4ebd2SPeter Brune           }
224f1c6b773SPeter Brune           ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
22571b4ebd2SPeter Brune           PetscFunctionReturn(0);
22671b4ebd2SPeter Brune         }
227b000cd8dSPeter Brune         if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
22871b4ebd2SPeter Brune           t1 = .5*(gnorm*gnorm - fnorm*fnorm) - lambda*initslope;
22971b4ebd2SPeter Brune           t2 = .5*(gnormprev*gnormprev  - fnorm*fnorm) - lambdaprev*initslope;
23071b4ebd2SPeter Brune           a  = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
23171b4ebd2SPeter Brune           b  = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
23271b4ebd2SPeter Brune           d  = b*b - 3*a*initslope;
23371b4ebd2SPeter Brune           if (d < 0.0) d = 0.0;
23471b4ebd2SPeter Brune           if (a == 0.0) {
23571b4ebd2SPeter Brune             lambdatemp = -initslope/(2.0*b);
23671b4ebd2SPeter Brune           } else {
23771b4ebd2SPeter Brune             lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a);
23871b4ebd2SPeter Brune           }
239b000cd8dSPeter Brune         } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
24071b4ebd2SPeter Brune           lambdatemp = -initslope/(gnorm*gnorm - fnorm*fnorm - 2.0*initslope);
24159405d5eSPeter Brune         } else {
24259405d5eSPeter Brune           SETERRQ(((PetscObject)linesearch)->comm, PETSC_ERR_SUP, "unsupported line search order for type bt");
24371b4ebd2SPeter Brune         }
24471b4ebd2SPeter Brune           lambdaprev = lambda;
24571b4ebd2SPeter Brune           gnormprev  = gnorm;
24671b4ebd2SPeter Brune         if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
24771b4ebd2SPeter Brune         if (lambdatemp <= .1*lambda) lambda     = .1*lambda;
24871b4ebd2SPeter Brune         else                         lambda     = lambdatemp;
24971b4ebd2SPeter Brune         ierr  = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
250b7b2e573SPeter Brune         if (linesearch->ops->viproject) {
251b7b2e573SPeter Brune           ierr = (*linesearch->ops->viproject)(snes,W);CHKERRQ(ierr);
252b7b2e573SPeter Brune         }
25371b4ebd2SPeter Brune         if (snes->nfuncs >= snes->max_funcs) {
25471b4ebd2SPeter Brune           ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr);
25571b4ebd2SPeter Brune           ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
256c69d1a72SBarry Smith                             (double)fnorm,(double)gnorm,(double)ynorm,(double)lambda,(double)initslope);CHKERRQ(ierr);
257f1c6b773SPeter Brune           ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
25871b4ebd2SPeter Brune           snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
25971b4ebd2SPeter Brune           PetscFunctionReturn(0);
26071b4ebd2SPeter Brune         }
26171b4ebd2SPeter Brune         ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr);
262c427506bSPeter Brune         ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
263c427506bSPeter Brune         if (domainerror) {
264f1c6b773SPeter Brune           ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
26571b4ebd2SPeter Brune           PetscFunctionReturn(0);
26671b4ebd2SPeter Brune         }
2679bd66eb0SPeter Brune         if (linesearch->ops->vinorm) {
2689bd66eb0SPeter Brune           gnorm = fnorm;
2699bd66eb0SPeter Brune           ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
2709bd66eb0SPeter Brune         } else {
27171b4ebd2SPeter Brune           ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
2729bd66eb0SPeter Brune         }
27371b4ebd2SPeter Brune         if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
27471b4ebd2SPeter Brune         if (.5*gnorm*gnorm < .5*fnorm*fnorm + lambda*alpha*initslope) { /* is reduction enough? */
27571b4ebd2SPeter Brune           if (monitor) {
276ea5d4fccSPeter Brune             ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
277b000cd8dSPeter Brune             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
278c69d1a72SBarry Smith               ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);CHKERRQ(ierr);
27971b4ebd2SPeter Brune             } else {
280c69d1a72SBarry Smith               ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);CHKERRQ(ierr);
28171b4ebd2SPeter Brune             }
282ea5d4fccSPeter Brune             ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
28371b4ebd2SPeter Brune           }
28471b4ebd2SPeter Brune           break;
28571b4ebd2SPeter Brune         } else {
28671b4ebd2SPeter Brune           if (monitor) {
287ea5d4fccSPeter Brune             ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
288b000cd8dSPeter Brune             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
289c69d1a72SBarry Smith               ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);CHKERRQ(ierr);
29071b4ebd2SPeter Brune             } else {
291c69d1a72SBarry Smith               ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);CHKERRQ(ierr);
29271b4ebd2SPeter Brune             }
293ea5d4fccSPeter Brune             ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
29471b4ebd2SPeter Brune           }
29571b4ebd2SPeter Brune         }
29671b4ebd2SPeter Brune       }
29771b4ebd2SPeter Brune     }
29871b4ebd2SPeter Brune   }
29971b4ebd2SPeter Brune 
30071b4ebd2SPeter Brune   /* postcheck */
3017b1df9c1SPeter Brune   ierr = SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);CHKERRQ(ierr);
30271b4ebd2SPeter Brune   if (changed_y) {
30371b4ebd2SPeter Brune     ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
3049bd66eb0SPeter Brune     if (linesearch->ops->viproject) {
3059bd66eb0SPeter Brune       ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
3069bd66eb0SPeter Brune     }
30771b4ebd2SPeter Brune   }
308c427506bSPeter Brune   if (changed_y || changed_w) { /* recompute the function if the step has changed */
309c427506bSPeter Brune     ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr);
31071b4ebd2SPeter Brune     ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
31171b4ebd2SPeter Brune     if (domainerror) {
312f1c6b773SPeter Brune       ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
31371b4ebd2SPeter Brune       PetscFunctionReturn(0);
31471b4ebd2SPeter Brune     }
3159bd66eb0SPeter Brune     if (linesearch->ops->vinorm) {
3169bd66eb0SPeter Brune       gnorm = fnorm;
3179bd66eb0SPeter Brune       ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
3189bd66eb0SPeter Brune     } else {
3199bd66eb0SPeter Brune       ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
3209bd66eb0SPeter Brune     }
3219bd66eb0SPeter Brune     ierr = VecNorm(Y,NORM_2,&ynorm);CHKERRQ(ierr);
322c427506bSPeter Brune     if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
3239bd66eb0SPeter Brune 
324c427506bSPeter Brune   }
32571b4ebd2SPeter Brune 
32671b4ebd2SPeter Brune   /* copy the solution over */
32771b4ebd2SPeter Brune   ierr = VecCopy(W, X);CHKERRQ(ierr);
328c427506bSPeter Brune   ierr = VecCopy(G, F);CHKERRQ(ierr);
329c427506bSPeter Brune   ierr = VecNorm(X, NORM_2, &xnorm);CHKERRQ(ierr);
330f1c6b773SPeter Brune   ierr = SNESLineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr);
331f1c6b773SPeter Brune   ierr = SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);CHKERRQ(ierr);
33271b4ebd2SPeter Brune   PetscFunctionReturn(0);
33371b4ebd2SPeter Brune }
33471b4ebd2SPeter Brune 
33571b4ebd2SPeter Brune #undef __FUNCT__
3367f1410a3SPeter Brune #define __FUNCT__ "SNESLineSearchView_BT"
3377f1410a3SPeter Brune PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
3387f1410a3SPeter Brune {
3397f1410a3SPeter Brune   PetscErrorCode    ierr;
3407f1410a3SPeter Brune   PetscBool         iascii;
3417f1410a3SPeter Brune   SNESLineSearch_BT *bt;
3427f1410a3SPeter Brune   PetscFunctionBegin;
343251f4c67SDmitry Karpeev   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
3447f1410a3SPeter Brune   bt = (SNESLineSearch_BT*)linesearch->data;
3457f1410a3SPeter Brune   if (iascii) {
346b000cd8dSPeter Brune     if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
3477f1410a3SPeter Brune     ierr = PetscViewerASCIIPrintf(viewer, "  interpolation: cubic\n");CHKERRQ(ierr);
348b000cd8dSPeter Brune     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
3497f1410a3SPeter Brune     ierr = PetscViewerASCIIPrintf(viewer, "  interpolation: quadratic\n");CHKERRQ(ierr);
3507f1410a3SPeter Brune     }
351007b6e36SPeter Brune     ierr = PetscViewerASCIIPrintf(viewer, "  alpha=%e\n", bt->alpha);CHKERRQ(ierr);
3527f1410a3SPeter Brune   }
3537f1410a3SPeter Brune   PetscFunctionReturn(0);
3547f1410a3SPeter Brune }
3557f1410a3SPeter Brune 
3567f1410a3SPeter Brune 
3577f1410a3SPeter Brune #undef __FUNCT__
358f1c6b773SPeter Brune #define __FUNCT__ "SNESLineSearchDestroy_BT"
359f1c6b773SPeter Brune static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
360731c067cSJed Brown {
361731c067cSJed Brown   PetscErrorCode ierr;
36271b4ebd2SPeter Brune 
36371b4ebd2SPeter Brune   PetscFunctionBegin;
364731c067cSJed Brown   ierr = PetscFree(linesearch->data);CHKERRQ(ierr);
36571b4ebd2SPeter Brune   PetscFunctionReturn(0);
36671b4ebd2SPeter Brune }
36771b4ebd2SPeter Brune 
36871b4ebd2SPeter Brune 
36971b4ebd2SPeter Brune #undef __FUNCT__
370f1c6b773SPeter Brune #define __FUNCT__ "SNESLineSearchSetFromOptions_BT"
371f1c6b773SPeter Brune static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch)
37271b4ebd2SPeter Brune {
37371b4ebd2SPeter Brune 
37471b4ebd2SPeter Brune   PetscErrorCode       ierr;
375f1c6b773SPeter Brune   SNESLineSearch_BT    *bt;
37671b4ebd2SPeter Brune   PetscFunctionBegin;
37771b4ebd2SPeter Brune 
378f1c6b773SPeter Brune   bt = (SNESLineSearch_BT*)linesearch->data;
37971b4ebd2SPeter Brune 
380f1c6b773SPeter Brune   ierr = PetscOptionsHead("SNESLineSearch BT options");CHKERRQ(ierr);
3817a35526eSPeter Brune   ierr = PetscOptionsReal("-snes_linesearch_alpha",   "Descent tolerance",        "SNESLineSearchBT", bt->alpha, &bt->alpha, PETSC_NULL);CHKERRQ(ierr);
38271b4ebd2SPeter Brune 
38371b4ebd2SPeter Brune   ierr = PetscOptionsTail();CHKERRQ(ierr);
38471b4ebd2SPeter Brune   PetscFunctionReturn(0);
38571b4ebd2SPeter Brune }
38671b4ebd2SPeter Brune 
38771b4ebd2SPeter Brune 
38871b4ebd2SPeter Brune #undef __FUNCT__
389f1c6b773SPeter Brune #define __FUNCT__ "SNESLineSearchCreate_BT"
390954494b2SJed Brown /*MC
391db609ea7SPeter Brune    SNESLINESEARCHBT - Backtracking line search.
392954494b2SJed Brown 
393db609ea7SPeter Brune    This line search finds the minimum of a polynomial fitting of the L2 norm of the
394db609ea7SPeter Brune    function. If this fit does not satisfy the conditions for progress, the interval shrinks
395db609ea7SPeter Brune    and the fit is reattempted at most max_it times or until lambda is below minlambda.
396954494b2SJed Brown 
397cd7522eaSPeter Brune    Options Database Keys:
398cd7522eaSPeter Brune +  -snes_linesearch_alpha<1e-4> - slope descent parameter
399db609ea7SPeter Brune .  -snes_linesearch_damping<1.0> - initial step length
400db609ea7SPeter Brune .  -snes_linesearch_max_it<40> - maximum number of shrinking step
401db609ea7SPeter Brune .  -snes_linesearch_minlambda<1e-12> - minimum step length allowed
402cd7522eaSPeter Brune -  -snes_linesearch_order<cubic,quadratic> - order of the approximation
403cd7522eaSPeter Brune 
404954494b2SJed Brown    Level: advanced
405954494b2SJed Brown 
406db609ea7SPeter Brune    Notes:
407db609ea7SPeter Brune    This line search is taken from "Numerical Methods for Unconstrained
408db609ea7SPeter Brune    Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
409db609ea7SPeter Brune 
410f1c6b773SPeter Brune .keywords: SNES, SNESLineSearch, damping
411954494b2SJed Brown 
412f1c6b773SPeter Brune .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
413954494b2SJed Brown M*/
414f1c6b773SPeter Brune PETSC_EXTERN_C PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
41571b4ebd2SPeter Brune {
41671b4ebd2SPeter Brune 
417f1c6b773SPeter Brune   SNESLineSearch_BT  *bt;
41871b4ebd2SPeter Brune   PetscErrorCode ierr;
41971b4ebd2SPeter Brune 
42071b4ebd2SPeter Brune   PetscFunctionBegin;
421f1c6b773SPeter Brune   linesearch->ops->apply          = SNESLineSearchApply_BT;
422f1c6b773SPeter Brune   linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
423f1c6b773SPeter Brune   linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
42471b4ebd2SPeter Brune   linesearch->ops->reset          = PETSC_NULL;
4257f1410a3SPeter Brune   linesearch->ops->view           = SNESLineSearchView_BT;
42671b4ebd2SPeter Brune   linesearch->ops->setup          = PETSC_NULL;
42771b4ebd2SPeter Brune 
428f1c6b773SPeter Brune   ierr = PetscNewLog(linesearch, SNESLineSearch_BT, &bt);CHKERRQ(ierr);
42971b4ebd2SPeter Brune   linesearch->data = (void *)bt;
43071b4ebd2SPeter Brune   linesearch->max_its = 40;
431b000cd8dSPeter Brune   linesearch->order = SNES_LINESEARCH_ORDER_CUBIC;
43271b4ebd2SPeter Brune   bt->alpha = 1e-4;
43371b4ebd2SPeter Brune 
43471b4ebd2SPeter Brune   PetscFunctionReturn(0);
43571b4ebd2SPeter Brune }
436