xref: /petsc/src/snes/linesearch/impls/bt/linesearchbt.c (revision ccb4e88a40f0b86eaeca07ff64c64e4de2fae686)
1 #include <petsc/private/linesearchimpl.h> /*I  "petscsnes.h"  I*/
2 #include <petsc/private/snesimpl.h>
3 
4 typedef struct {
5   PetscReal alpha;        /* sufficient decrease parameter */
6 } SNESLineSearch_BT;
7 
8 /*@
9    SNESLineSearchBTSetAlpha - Sets the descent parameter, alpha, in the BT linesearch variant.
10 
11    Input Parameters:
12 +  linesearch - linesearch context
13 -  alpha - The descent parameter
14 
15    Level: intermediate
16 
17 .seealso: SNESLineSearchSetLambda(), SNESLineSearchGetTolerances() SNESLINESEARCHBT
18 @*/
19 PetscErrorCode SNESLineSearchBTSetAlpha(SNESLineSearch linesearch, PetscReal alpha)
20 {
21   SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
22 
23   PetscFunctionBegin;
24   PetscValidHeaderSpecific(linesearch,SNESLINESEARCH_CLASSID,1);
25   bt->alpha = alpha;
26   PetscFunctionReturn(0);
27 }
28 
29 /*@
30    SNESLineSearchBTGetAlpha - Gets the descent parameter, alpha, in the BT linesearch variant.
31 
32    Input Parameters:
33 .  linesearch - linesearch context
34 
35    Output Parameters:
36 .  alpha - The descent parameter
37 
38    Level: intermediate
39 
40 .seealso: SNESLineSearchGetLambda(), SNESLineSearchGetTolerances() SNESLINESEARCHBT
41 @*/
42 PetscErrorCode SNESLineSearchBTGetAlpha(SNESLineSearch linesearch, PetscReal *alpha)
43 {
44   SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
45 
46   PetscFunctionBegin;
47   PetscValidHeaderSpecific(linesearch,SNESLINESEARCH_CLASSID,1);
48   *alpha = bt->alpha;
49   PetscFunctionReturn(0);
50 }
51 
52 static PetscErrorCode  SNESLineSearchApply_BT(SNESLineSearch linesearch)
53 {
54   PetscBool         changed_y,changed_w;
55   PetscErrorCode    ierr;
56   Vec               X,F,Y,W,G;
57   SNES              snes;
58   PetscReal         fnorm, xnorm, ynorm, gnorm;
59   PetscReal         lambda,lambdatemp,lambdaprev,minlambda,maxstep,initslope,alpha,stol;
60   PetscReal         t1,t2,a,b,d;
61   PetscReal         f;
62   PetscReal         g,gprev;
63   PetscViewer       monitor;
64   PetscInt          max_its,count;
65   SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
66   Mat               jac;
67   PetscErrorCode    (*objective)(SNES,Vec,PetscReal*,void*);
68 
69   PetscFunctionBegin;
70   ierr = SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G);CHKERRQ(ierr);
71   ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr);
72   ierr = SNESLineSearchGetLambda(linesearch, &lambda);CHKERRQ(ierr);
73   ierr = SNESLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr);
74   ierr = SNESLineSearchGetDefaultMonitor(linesearch, &monitor);CHKERRQ(ierr);
75   ierr = SNESLineSearchGetTolerances(linesearch,&minlambda,&maxstep,NULL,NULL,NULL,&max_its);CHKERRQ(ierr);
76   ierr = SNESGetTolerances(snes,NULL,NULL,&stol,NULL,NULL);CHKERRQ(ierr);
77   ierr = SNESGetObjective(snes,&objective,NULL);CHKERRQ(ierr);
78   alpha = bt->alpha;
79 
80   ierr = SNESGetJacobian(snes, &jac, NULL, NULL, NULL);CHKERRQ(ierr);
81   if (!jac && !objective) SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix");
82 
83   ierr = SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);CHKERRQ(ierr);
84   ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);CHKERRQ(ierr);
85 
86   ierr = VecNormBegin(Y, NORM_2, &ynorm);CHKERRQ(ierr);
87   ierr = VecNormBegin(X, NORM_2, &xnorm);CHKERRQ(ierr);
88   ierr = VecNormEnd(Y, NORM_2, &ynorm);CHKERRQ(ierr);
89   ierr = VecNormEnd(X, NORM_2, &xnorm);CHKERRQ(ierr);
90 
91   if (ynorm == 0.0) {
92     if (monitor) {
93       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
94       ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Initial direction and size is 0\n");CHKERRQ(ierr);
95       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
96     }
97     ierr = VecCopy(X,W);CHKERRQ(ierr);
98     ierr = VecCopy(F,G);CHKERRQ(ierr);
99     ierr = SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm);CHKERRQ(ierr);
100     ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);CHKERRQ(ierr);
101     PetscFunctionReturn(0);
102   }
103   if (ynorm > maxstep) {        /* Step too big, so scale back */
104     if (monitor) {
105       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
106       ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Scaling step by %14.12e old ynorm %14.12e\n", (double)(maxstep/ynorm),(double)ynorm);CHKERRQ(ierr);
107       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
108     }
109     ierr  = VecScale(Y,maxstep/(ynorm));CHKERRQ(ierr);
110     ynorm = maxstep;
111   }
112 
113   /* if the SNES has an objective set, use that instead of the function value */
114   if (objective) {
115     ierr = SNESComputeObjective(snes,X,&f);CHKERRQ(ierr);
116   } else {
117     f = fnorm*fnorm;
118   }
119 
120   /* compute the initial slope */
121   if (objective) {
122     /* slope comes from the function (assumed to be the gradient of the objective */
123     ierr = VecDotRealPart(Y,F,&initslope);CHKERRQ(ierr);
124   } else {
125     /* slope comes from the normal equations */
126     ierr = MatMult(jac,Y,W);CHKERRQ(ierr);
127     ierr = VecDotRealPart(F,W,&initslope);CHKERRQ(ierr);
128     if (initslope > 0.0)  initslope = -initslope;
129     if (initslope == 0.0) initslope = -1.0;
130   }
131 
132   while (PETSC_TRUE) {
133     ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
134     if (linesearch->ops->viproject) {
135       ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
136     }
137     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
138       ierr         = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr);
139       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
140       ierr         = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);CHKERRQ(ierr);
141       PetscFunctionReturn(0);
142     }
143 
144     if (objective) {
145       ierr = SNESComputeObjective(snes,W,&g);CHKERRQ(ierr);
146     } else {
147       ierr = (*linesearch->ops->snesfunc)(snes,W,G);CHKERRQ(ierr);
148       if (linesearch->ops->vinorm) {
149         gnorm = fnorm;
150         ierr  = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
151       } else {
152         ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
153       }
154       g = PetscSqr(gnorm);
155     }
156     ierr = SNESLineSearchMonitor(linesearch);CHKERRQ(ierr);
157 
158     if (!PetscIsInfOrNanReal(g)) break;
159     if (monitor) {
160       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
161       ierr = PetscViewerASCIIPrintf(monitor,"    Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n",(double)lambda);CHKERRQ(ierr);
162       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
163     }
164     if (lambda <= minlambda) {
165       SNESCheckFunctionNorm(snes,g);
166     }
167     lambda = .5*lambda;
168   }
169 
170   if (!objective) {
171     ierr = PetscInfo2(snes,"Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);CHKERRQ(ierr);
172   }
173   if (.5*g <= .5*f + lambda*alpha*initslope) { /* Sufficient reduction or step tolerance convergence */
174     if (monitor) {
175       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
176       if (!objective) {
177         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);CHKERRQ(ierr);
178       } else {
179         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Using full step: obj %14.12e obj %14.12e\n", (double)f, (double)g);CHKERRQ(ierr);
180       }
181       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
182     }
183   } else {
184     /* Since the full step didn't work and the step is tiny, quit */
185     if (stol*xnorm > ynorm) {
186       ierr = SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm);CHKERRQ(ierr);
187       ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);CHKERRQ(ierr);
188       if (monitor) {
189         ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
190         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n",(double)ynorm,(double)stol*xnorm);CHKERRQ(ierr);
191         ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
192       }
193       PetscFunctionReturn(0);
194     }
195     /* Fit points with quadratic */
196     lambdatemp = -initslope/(g - f - 2.0*lambda*initslope);
197     lambdaprev = lambda;
198     gprev      = g;
199     if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
200     if (lambdatemp <= .1*lambda) lambda = .1*lambda;
201     else                         lambda = lambdatemp;
202 
203     ierr  = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
204     if (linesearch->ops->viproject) {
205       ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
206     }
207     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
208       ierr         = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr);
209       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
210       ierr         = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);CHKERRQ(ierr);
211       PetscFunctionReturn(0);
212     }
213     if (objective) {
214       ierr = SNESComputeObjective(snes,W,&g);CHKERRQ(ierr);
215     } else {
216       ierr = (*linesearch->ops->snesfunc)(snes,W,G);CHKERRQ(ierr);
217       if (linesearch->ops->vinorm) {
218         gnorm = fnorm;
219         ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
220       } else {
221         ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
222       }
223       g = PetscSqr(gnorm);
224     }
225     if (PetscIsInfOrNanReal(g)) {
226       ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF);CHKERRQ(ierr);
227       ierr = PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");CHKERRQ(ierr);
228       PetscFunctionReturn(0);
229     }
230     if (monitor) {
231       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
232       if (!objective) {
233         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: gnorm after quadratic fit %14.12e\n",(double)gnorm);CHKERRQ(ierr);
234       } else {
235         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: obj after quadratic fit %14.12e\n",(double)g);CHKERRQ(ierr);
236       }
237       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
238     }
239     if (.5*g < .5*f + lambda*alpha*initslope) { /* sufficient reduction */
240       if (monitor) {
241         ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
242         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda);CHKERRQ(ierr);
243         ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
244       }
245     } else {
246       /* Fit points with cubic */
247       for (count = 0; count < max_its; count++) {
248         if (lambda <= minlambda) {
249           if (monitor) {
250             ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
251             ierr = PetscViewerASCIIPrintf(monitor,"    Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr);
252             if (!objective) {
253               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",
254                                                          (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);CHKERRQ(ierr);
255             } else {
256               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",
257                                                          (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);CHKERRQ(ierr);
258             }
259             ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
260           }
261           ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);CHKERRQ(ierr);
262           PetscFunctionReturn(0);
263         }
264         if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
265           t1 = .5*(g - f) - lambda*initslope;
266           t2 = .5*(gprev  - f) - lambdaprev*initslope;
267           a  = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
268           b  = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
269           d  = b*b - 3*a*initslope;
270           if (d < 0.0) d = 0.0;
271           if (a == 0.0) lambdatemp = -initslope/(2.0*b);
272           else lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a);
273 
274         } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
275           lambdatemp = -initslope/(g - f - 2.0*initslope);
276         } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt");
277         lambdaprev = lambda;
278         gprev      = g;
279         if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
280         if (lambdatemp <= .1*lambda) lambda     = .1*lambda;
281         else                         lambda     = lambdatemp;
282         ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
283         if (linesearch->ops->viproject) {
284           ierr = (*linesearch->ops->viproject)(snes,W);CHKERRQ(ierr);
285         }
286         if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
287           ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr);
288           if (!objective) {
289             ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
290                               (double)fnorm,(double)gnorm,(double)ynorm,(double)lambda,(double)initslope);CHKERRQ(ierr);
291           }
292           ierr         = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);CHKERRQ(ierr);
293           snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
294           PetscFunctionReturn(0);
295         }
296         if (objective) {
297           ierr = SNESComputeObjective(snes,W,&g);CHKERRQ(ierr);
298         } else {
299           ierr = (*linesearch->ops->snesfunc)(snes,W,G);CHKERRQ(ierr);
300           if (linesearch->ops->vinorm) {
301             gnorm = fnorm;
302             ierr  = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
303           } else {
304             ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
305           }
306           g = PetscSqr(gnorm);
307         }
308         if (PetscIsInfOrNanReal(g)) {
309           ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF);CHKERRQ(ierr);
310           ierr = PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");CHKERRQ(ierr);
311           PetscFunctionReturn(0);
312         }
313         if (.5*g < .5*f + lambda*alpha*initslope) { /* is reduction enough? */
314           if (monitor) {
315             ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
316             if (!objective) {
317               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
318                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);CHKERRQ(ierr);
319               } else {
320                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);CHKERRQ(ierr);
321               }
322               ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
323             } else {
324               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
325                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
326               } else {
327                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
328               }
329               ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
330             }
331           }
332           break;
333         } else if (monitor) {
334           ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
335           if (!objective) {
336             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
337               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);
338             } else {
339               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);
340             }
341             ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
342           } else {
343             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
344               ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
345             } else {
346               ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
347             }
348             ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
349           }
350         }
351       }
352     }
353   }
354 
355   /* postcheck */
356   /* update Y to lambda*Y so that W is consistent with  X - lambda*Y */
357   ierr = VecScale(Y,lambda);CHKERRQ(ierr);
358   ierr = SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);CHKERRQ(ierr);
359   if (changed_y) {
360     ierr = VecWAXPY(W,-1.0,Y,X);CHKERRQ(ierr);
361     if (linesearch->ops->viproject) {
362       ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
363     }
364   }
365   if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
366     ierr = (*linesearch->ops->snesfunc)(snes,W,G);CHKERRQ(ierr);
367     if (linesearch->ops->vinorm) {
368       gnorm = fnorm;
369       ierr  = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
370     } else {
371       ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
372     }
373     ierr = VecNorm(Y,NORM_2,&ynorm);CHKERRQ(ierr);
374     if (PetscIsInfOrNanReal(gnorm)) {
375       ierr = SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_FAILED_NANORINF);CHKERRQ(ierr);
376       ierr = PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");CHKERRQ(ierr);
377       PetscFunctionReturn(0);
378     }
379   }
380 
381   /* copy the solution over */
382   ierr = VecCopy(W, X);CHKERRQ(ierr);
383   ierr = VecCopy(G, F);CHKERRQ(ierr);
384   ierr = VecNorm(X, NORM_2, &xnorm);CHKERRQ(ierr);
385   ierr = SNESLineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr);
386   ierr = SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);CHKERRQ(ierr);
387   PetscFunctionReturn(0);
388 }
389 
390 PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
391 {
392   PetscErrorCode    ierr;
393   PetscBool         iascii;
394   SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
395 
396   PetscFunctionBegin;
397   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
398   if (iascii) {
399     if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
400       ierr = PetscViewerASCIIPrintf(viewer, "  interpolation: cubic\n");CHKERRQ(ierr);
401     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
402       ierr = PetscViewerASCIIPrintf(viewer, "  interpolation: quadratic\n");CHKERRQ(ierr);
403     }
404     ierr = PetscViewerASCIIPrintf(viewer, "  alpha=%e\n", (double)bt->alpha);CHKERRQ(ierr);
405   }
406   PetscFunctionReturn(0);
407 }
408 
409 static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
410 {
411   PetscErrorCode ierr;
412 
413   PetscFunctionBegin;
414   ierr = PetscFree(linesearch->data);CHKERRQ(ierr);
415   PetscFunctionReturn(0);
416 }
417 
418 static PetscErrorCode SNESLineSearchSetFromOptions_BT(PetscOptionItems *PetscOptionsObject,SNESLineSearch linesearch)
419 {
420   PetscErrorCode    ierr;
421   SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
422 
423   PetscFunctionBegin;
424   ierr = PetscOptionsHead(PetscOptionsObject,"SNESLineSearch BT options");CHKERRQ(ierr);
425   ierr = PetscOptionsReal("-snes_linesearch_alpha",   "Descent tolerance",        "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL);CHKERRQ(ierr);
426   ierr = PetscOptionsTail();CHKERRQ(ierr);
427   PetscFunctionReturn(0);
428 }
429 
430 /*MC
431    SNESLINESEARCHBT - Backtracking line search.
432 
433    This line search finds the minimum of a polynomial fitting of the L2 norm of the
434    function or the objective function if it is provided with SNESSetObjective(). If this fit does not satisfy the conditions for progress, the interval shrinks
435    and the fit is reattempted at most max_it times or until lambda is below minlambda.
436 
437    Options Database Keys:
438 +  -snes_linesearch_alpha <1e\-4> - slope descent parameter
439 .  -snes_linesearch_damping <1.0> - initial step length
440 .  -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the
441                                        step is scaled back to be of this length at the beginning of the line search
442 .  -snes_linesearch_max_it <40> - maximum number of shrinking step
443 .  -snes_linesearch_minlambda <1e\-12> - minimum step length allowed
444 -  -snes_linesearch_order <cubic,quadratic> - order of the approximation
445 
446    Level: advanced
447 
448    Notes:
449    This line search is taken from "Numerical Methods for Unconstrained
450    Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
451 
452    This line search will always produce a step that is less than or equal to, in length, the full step size.
453 
454 .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
455 M*/
456 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
457 {
458 
459   SNESLineSearch_BT *bt;
460   PetscErrorCode    ierr;
461 
462   PetscFunctionBegin;
463   linesearch->ops->apply          = SNESLineSearchApply_BT;
464   linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
465   linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
466   linesearch->ops->reset          = NULL;
467   linesearch->ops->view           = SNESLineSearchView_BT;
468   linesearch->ops->setup          = NULL;
469 
470   ierr = PetscNewLog(linesearch,&bt);CHKERRQ(ierr);
471 
472   linesearch->data    = (void*)bt;
473   linesearch->max_its = 40;
474   linesearch->order   = SNES_LINESEARCH_ORDER_CUBIC;
475   bt->alpha           = 1e-4;
476   PetscFunctionReturn(0);
477 }
478