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