xref: /petsc/src/snes/linesearch/impls/bt/linesearchbt.c (revision e5a36eccef3d6b83a2c625c30d0dfd5adb4001f2)
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   ierr = SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);CHKERRQ(ierr);
357   if (changed_y) {
358     ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
359     if (linesearch->ops->viproject) {
360       ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
361     }
362   }
363   if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
364     ierr = (*linesearch->ops->snesfunc)(snes,W,G);CHKERRQ(ierr);
365     if (linesearch->ops->vinorm) {
366       gnorm = fnorm;
367       ierr  = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
368     } else {
369       ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
370     }
371     ierr = VecNorm(Y,NORM_2,&ynorm);CHKERRQ(ierr);
372     if (PetscIsInfOrNanReal(gnorm)) {
373       ierr = SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_FAILED_NANORINF);CHKERRQ(ierr);
374       ierr = PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");CHKERRQ(ierr);
375       PetscFunctionReturn(0);
376     }
377   }
378 
379   /* copy the solution over */
380   ierr = VecCopy(W, X);CHKERRQ(ierr);
381   ierr = VecCopy(G, F);CHKERRQ(ierr);
382   ierr = VecNorm(X, NORM_2, &xnorm);CHKERRQ(ierr);
383   ierr = SNESLineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr);
384   ierr = SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);CHKERRQ(ierr);
385   PetscFunctionReturn(0);
386 }
387 
388 PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
389 {
390   PetscErrorCode    ierr;
391   PetscBool         iascii;
392   SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
393 
394   PetscFunctionBegin;
395   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
396   if (iascii) {
397     if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
398       ierr = PetscViewerASCIIPrintf(viewer, "  interpolation: cubic\n");CHKERRQ(ierr);
399     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
400       ierr = PetscViewerASCIIPrintf(viewer, "  interpolation: quadratic\n");CHKERRQ(ierr);
401     }
402     ierr = PetscViewerASCIIPrintf(viewer, "  alpha=%e\n", (double)bt->alpha);CHKERRQ(ierr);
403   }
404   PetscFunctionReturn(0);
405 }
406 
407 static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
408 {
409   PetscErrorCode ierr;
410 
411   PetscFunctionBegin;
412   ierr = PetscFree(linesearch->data);CHKERRQ(ierr);
413   PetscFunctionReturn(0);
414 }
415 
416 static PetscErrorCode SNESLineSearchSetFromOptions_BT(PetscOptionItems *PetscOptionsObject,SNESLineSearch linesearch)
417 {
418   PetscErrorCode    ierr;
419   SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
420 
421   PetscFunctionBegin;
422   ierr = PetscOptionsHead(PetscOptionsObject,"SNESLineSearch BT options");CHKERRQ(ierr);
423   ierr = PetscOptionsReal("-snes_linesearch_alpha",   "Descent tolerance",        "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL);CHKERRQ(ierr);
424   ierr = PetscOptionsTail();CHKERRQ(ierr);
425   PetscFunctionReturn(0);
426 }
427 
428 /*MC
429    SNESLINESEARCHBT - Backtracking line search.
430 
431    This line search finds the minimum of a polynomial fitting of the L2 norm of the
432    function or the objective function if it is provided with SNESSetObjective(). If this fit does not satisfy the conditions for progress, the interval shrinks
433    and the fit is reattempted at most max_it times or until lambda is below minlambda.
434 
435    Options Database Keys:
436 +  -snes_linesearch_alpha<1e-4> - slope descent parameter
437 .  -snes_linesearch_damping<1.0> - initial step length
438 .  -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the
439                                        step is scaled back to be of this length at the beginning of the line search
440 .  -snes_linesearch_max_it<40> - maximum number of shrinking step
441 .  -snes_linesearch_minlambda<1e-12> - minimum step length allowed
442 -  -snes_linesearch_order<cubic,quadratic> - order of the approximation
443 
444    Level: advanced
445 
446    Notes:
447    This line search is taken from "Numerical Methods for Unconstrained
448    Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
449 
450 .keywords: SNES, SNESLineSearch, damping
451 
452 .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
453 M*/
454 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
455 {
456 
457   SNESLineSearch_BT *bt;
458   PetscErrorCode    ierr;
459 
460   PetscFunctionBegin;
461   linesearch->ops->apply          = SNESLineSearchApply_BT;
462   linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
463   linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
464   linesearch->ops->reset          = NULL;
465   linesearch->ops->view           = SNESLineSearchView_BT;
466   linesearch->ops->setup          = NULL;
467 
468   ierr = PetscNewLog(linesearch,&bt);CHKERRQ(ierr);
469 
470   linesearch->data    = (void*)bt;
471   linesearch->max_its = 40;
472   linesearch->order   = SNES_LINESEARCH_ORDER_CUBIC;
473   bt->alpha           = 1e-4;
474   PetscFunctionReturn(0);
475 }
476