xref: /petsc/src/snes/linesearch/impls/bt/linesearchbt.c (revision 65f8aed5f7eaa1e2ef2ddeffe666264e0669c876)
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       ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);CHKERRQ(ierr);
166       PetscFunctionReturn(0);
167     }
168     lambda = .5*lambda;
169   }
170 
171   if (!objective) {
172     ierr = PetscInfo2(snes,"Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);CHKERRQ(ierr);
173   }
174   if (.5*g <= .5*f + lambda*alpha*initslope) { /* Sufficient reduction or step tolerance convergence */
175     if (monitor) {
176       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
177       if (!objective) {
178         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);CHKERRQ(ierr);
179       } else {
180         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Using full step: obj %14.12e obj %14.12e\n", (double)f, (double)g);CHKERRQ(ierr);
181       }
182       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
183     }
184   } else {
185     /* Since the full step didn't work and the step is tiny, quit */
186     if (stol*xnorm > ynorm) {
187       ierr = SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm);CHKERRQ(ierr);
188       ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);CHKERRQ(ierr);
189       if (monitor) {
190         ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
191         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n",(double)ynorm,(double)stol*xnorm);CHKERRQ(ierr);
192         ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
193       }
194       PetscFunctionReturn(0);
195     }
196     /* Fit points with quadratic */
197     lambdatemp = -initslope/(g - f - 2.0*lambda*initslope);
198     lambdaprev = lambda;
199     gprev      = g;
200     if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
201     if (lambdatemp <= .1*lambda) lambda = .1*lambda;
202     else                         lambda = lambdatemp;
203 
204     ierr  = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
205     if (linesearch->ops->viproject) {
206       ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
207     }
208     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
209       ierr         = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr);
210       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
211       ierr         = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);CHKERRQ(ierr);
212       PetscFunctionReturn(0);
213     }
214     if (objective) {
215       ierr = SNESComputeObjective(snes,W,&g);CHKERRQ(ierr);
216     } else {
217       ierr = (*linesearch->ops->snesfunc)(snes,W,G);CHKERRQ(ierr);
218       if (linesearch->ops->vinorm) {
219         gnorm = fnorm;
220         ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
221       } else {
222         ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
223       }
224       g = PetscSqr(gnorm);
225     }
226     if (PetscIsInfOrNanReal(g)) {
227       ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF);CHKERRQ(ierr);
228       ierr = PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");CHKERRQ(ierr);
229       PetscFunctionReturn(0);
230     }
231     if (monitor) {
232       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
233       if (!objective) {
234         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: gnorm after quadratic fit %14.12e\n",(double)gnorm);CHKERRQ(ierr);
235       } else {
236         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: obj after quadratic fit %14.12e\n",(double)g);CHKERRQ(ierr);
237       }
238       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
239     }
240     if (.5*g < .5*f + lambda*alpha*initslope) { /* sufficient reduction */
241       if (monitor) {
242         ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
243         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda);CHKERRQ(ierr);
244         ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
245       }
246     } else {
247       /* Fit points with cubic */
248       for (count = 0; count < max_its; count++) {
249         if (lambda <= minlambda) {
250           if (monitor) {
251             ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
252             ierr = PetscViewerASCIIPrintf(monitor,"    Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr);
253             if (!objective) {
254               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",
255                                                          (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);CHKERRQ(ierr);
256             } else {
257               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",
258                                                          (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);CHKERRQ(ierr);
259             }
260             ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
261           }
262           ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);CHKERRQ(ierr);
263           PetscFunctionReturn(0);
264         }
265         if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
266           t1 = .5*(g - f) - lambda*initslope;
267           t2 = .5*(gprev  - f) - lambdaprev*initslope;
268           a  = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
269           b  = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
270           d  = b*b - 3*a*initslope;
271           if (d < 0.0) d = 0.0;
272           if (a == 0.0) lambdatemp = -initslope/(2.0*b);
273           else lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a);
274 
275         } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
276           lambdatemp = -initslope/(g - f - 2.0*initslope);
277         } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt");
278         lambdaprev = lambda;
279         gprev      = g;
280         if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
281         if (lambdatemp <= .1*lambda) lambda     = .1*lambda;
282         else                         lambda     = lambdatemp;
283         ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
284         if (linesearch->ops->viproject) {
285           ierr = (*linesearch->ops->viproject)(snes,W);CHKERRQ(ierr);
286         }
287         if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
288           ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr);
289           if (!objective) {
290             ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
291                               (double)fnorm,(double)gnorm,(double)ynorm,(double)lambda,(double)initslope);CHKERRQ(ierr);
292           }
293           ierr         = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);CHKERRQ(ierr);
294           snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
295           PetscFunctionReturn(0);
296         }
297         if (objective) {
298           ierr = SNESComputeObjective(snes,W,&g);CHKERRQ(ierr);
299         } else {
300           ierr = (*linesearch->ops->snesfunc)(snes,W,G);CHKERRQ(ierr);
301           if (linesearch->ops->vinorm) {
302             gnorm = fnorm;
303             ierr  = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
304           } else {
305             ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
306           }
307           g = PetscSqr(gnorm);
308         }
309         if (PetscIsInfOrNanReal(g)) {
310           ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF);CHKERRQ(ierr);
311           ierr = PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");CHKERRQ(ierr);
312           PetscFunctionReturn(0);
313         }
314         if (.5*g < .5*f + lambda*alpha*initslope) { /* is reduction enough? */
315           if (monitor) {
316             ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
317             if (!objective) {
318               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
319                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);CHKERRQ(ierr);
320               } else {
321                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);CHKERRQ(ierr);
322               }
323               ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
324             } else {
325               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
326                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
327               } else {
328                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
329               }
330               ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
331             }
332           }
333           break;
334         } else if (monitor) {
335           ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
336           if (!objective) {
337             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
338               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);
339             } else {
340               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);
341             }
342             ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
343           } else {
344             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
345               ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
346             } else {
347               ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
348             }
349             ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
350           }
351         }
352       }
353     }
354   }
355 
356   /* postcheck */
357   ierr = SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);CHKERRQ(ierr);
358   if (changed_y) {
359     ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
360     if (linesearch->ops->viproject) {
361       ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
362     }
363   }
364   if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
365     ierr = (*linesearch->ops->snesfunc)(snes,W,G);CHKERRQ(ierr);
366     if (linesearch->ops->vinorm) {
367       gnorm = fnorm;
368       ierr  = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
369     } else {
370       ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
371     }
372     ierr = VecNorm(Y,NORM_2,&ynorm);CHKERRQ(ierr);
373     if (PetscIsInfOrNanReal(gnorm)) {
374       ierr = SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_FAILED_NANORINF);CHKERRQ(ierr);
375       ierr = PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");CHKERRQ(ierr);
376       PetscFunctionReturn(0);
377     }
378   }
379 
380   /* copy the solution over */
381   ierr = VecCopy(W, X);CHKERRQ(ierr);
382   ierr = VecCopy(G, F);CHKERRQ(ierr);
383   ierr = VecNorm(X, NORM_2, &xnorm);CHKERRQ(ierr);
384   ierr = SNESLineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr);
385   ierr = SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);CHKERRQ(ierr);
386   PetscFunctionReturn(0);
387 }
388 
389 PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
390 {
391   PetscErrorCode    ierr;
392   PetscBool         iascii;
393   SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
394 
395   PetscFunctionBegin;
396   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
397   if (iascii) {
398     if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
399       ierr = PetscViewerASCIIPrintf(viewer, "  interpolation: cubic\n");CHKERRQ(ierr);
400     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
401       ierr = PetscViewerASCIIPrintf(viewer, "  interpolation: quadratic\n");CHKERRQ(ierr);
402     }
403     ierr = PetscViewerASCIIPrintf(viewer, "  alpha=%e\n", (double)bt->alpha);CHKERRQ(ierr);
404   }
405   PetscFunctionReturn(0);
406 }
407 
408 static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
409 {
410   PetscErrorCode ierr;
411 
412   PetscFunctionBegin;
413   ierr = PetscFree(linesearch->data);CHKERRQ(ierr);
414   PetscFunctionReturn(0);
415 }
416 
417 static PetscErrorCode SNESLineSearchSetFromOptions_BT(PetscOptionItems *PetscOptionsObject,SNESLineSearch linesearch)
418 {
419   PetscErrorCode    ierr;
420   SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
421 
422   PetscFunctionBegin;
423   ierr = PetscOptionsHead(PetscOptionsObject,"SNESLineSearch BT options");CHKERRQ(ierr);
424   ierr = PetscOptionsReal("-snes_linesearch_alpha",   "Descent tolerance",        "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL);CHKERRQ(ierr);
425   ierr = PetscOptionsTail();CHKERRQ(ierr);
426   PetscFunctionReturn(0);
427 }
428 
429 /*MC
430    SNESLINESEARCHBT - Backtracking line search.
431 
432    This line search finds the minimum of a polynomial fitting of the L2 norm of the
433    function or the objective function if it is provided with SNESSetObjective(). If this fit does not satisfy the conditions for progress, the interval shrinks
434    and the fit is reattempted at most max_it times or until lambda is below minlambda.
435 
436    Options Database Keys:
437 +  -snes_linesearch_alpha<1e-4> - slope descent parameter
438 .  -snes_linesearch_damping<1.0> - initial step length
439 .  -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the
440                                        step is scaled back to be of this length at the beginning of the line search
441 .  -snes_linesearch_max_it<40> - maximum number of shrinking step
442 .  -snes_linesearch_minlambda<1e-12> - minimum step length allowed
443 -  -snes_linesearch_order<cubic,quadratic> - order of the approximation
444 
445    Level: advanced
446 
447    Notes:
448    This line search is taken from "Numerical Methods for Unconstrained
449    Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
450 
451 .keywords: SNES, SNESLineSearch, damping
452 
453 .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
454 M*/
455 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
456 {
457 
458   SNESLineSearch_BT *bt;
459   PetscErrorCode    ierr;
460 
461   PetscFunctionBegin;
462   linesearch->ops->apply          = SNESLineSearchApply_BT;
463   linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
464   linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
465   linesearch->ops->reset          = NULL;
466   linesearch->ops->view           = SNESLineSearchView_BT;
467   linesearch->ops->setup          = NULL;
468 
469   ierr = PetscNewLog(linesearch,&bt);CHKERRQ(ierr);
470 
471   linesearch->data    = (void*)bt;
472   linesearch->max_its = 40;
473   linesearch->order   = SNES_LINESEARCH_ORDER_CUBIC;
474   bt->alpha           = 1e-4;
475   PetscFunctionReturn(0);
476 }
477