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