xref: /petsc/src/snes/linesearch/impls/bt/linesearchbt.c (revision 4fc747eaadbeca11629f314a99edccbc2ed7b3d3) 
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,"    Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
256                                                          (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);CHKERRQ(ierr);
257             } else {
258               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",
259                                                          (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);CHKERRQ(ierr);
260             }
261             ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
262           }
263           ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);CHKERRQ(ierr);
264           PetscFunctionReturn(0);
265         }
266         if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
267           t1 = .5*(g - f) - lambda*initslope;
268           t2 = .5*(gprev  - f) - lambdaprev*initslope;
269           a  = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
270           b  = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
271           d  = b*b - 3*a*initslope;
272           if (d < 0.0) d = 0.0;
273           if (a == 0.0) lambdatemp = -initslope/(2.0*b);
274           else lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a);
275 
276         } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
277           lambdatemp = -initslope/(g - f - 2.0*initslope);
278         } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt");
279         lambdaprev = lambda;
280         gprev      = g;
281         if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
282         if (lambdatemp <= .1*lambda) lambda     = .1*lambda;
283         else                         lambda     = lambdatemp;
284         ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
285         if (linesearch->ops->viproject) {
286           ierr = (*linesearch->ops->viproject)(snes,W);CHKERRQ(ierr);
287         }
288         if (snes->nfuncs >= snes->max_funcs) {
289           ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr);
290           if (!objective) {
291             ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
292                               (double)fnorm,(double)gnorm,(double)ynorm,(double)lambda,(double)initslope);CHKERRQ(ierr);
293           }
294           ierr         = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);CHKERRQ(ierr);
295           snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
296           PetscFunctionReturn(0);
297         }
298         if (objective) {
299           ierr = SNESComputeObjective(snes,W,&g);CHKERRQ(ierr);
300         } else {
301           ierr = (*linesearch->ops->snesfunc)(snes,W,G);CHKERRQ(ierr);
302           if (linesearch->ops->vinorm) {
303             gnorm = fnorm;
304             ierr  = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
305           } else {
306             ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
307           }
308           g = PetscSqr(gnorm);
309         }
310         if (PetscIsInfOrNanReal(g)) {
311           ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF);CHKERRQ(ierr);
312           ierr = PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");CHKERRQ(ierr);
313           PetscFunctionReturn(0);
314         }
315         if (.5*g < .5*f + lambda*alpha*initslope) { /* is reduction enough? */
316           if (monitor) {
317             ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
318             if (!objective) {
319               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
320                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);CHKERRQ(ierr);
321               } else {
322                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);CHKERRQ(ierr);
323               }
324               ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
325             } else {
326               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
327                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
328               } else {
329                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
330               }
331               ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
332             }
333           }
334           break;
335         } else if (monitor) {
336           ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
337           if (!objective) {
338             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
339               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);
340             } else {
341               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);
342             }
343             ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
344           } else {
345             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
346               ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
347             } else {
348               ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
349             }
350             ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
351           }
352         }
353       }
354     }
355   }
356 
357   /* postcheck */
358   ierr = SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);CHKERRQ(ierr);
359   if (changed_y) {
360     ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
361     if (linesearch->ops->viproject) {
362       ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
363     }
364   }
365   if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
366     ierr = (*linesearch->ops->snesfunc)(snes,W,G);CHKERRQ(ierr);
367     if (linesearch->ops->vinorm) {
368       gnorm = fnorm;
369       ierr  = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
370     } else {
371       ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
372     }
373     ierr = VecNorm(Y,NORM_2,&ynorm);CHKERRQ(ierr);
374     if (PetscIsInfOrNanReal(gnorm)) {
375       ierr = SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_FAILED_NANORINF);CHKERRQ(ierr);
376       ierr = PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");CHKERRQ(ierr);
377       PetscFunctionReturn(0);
378     }
379   }
380 
381   /* copy the solution over */
382   ierr = VecCopy(W, X);CHKERRQ(ierr);
383   ierr = VecCopy(G, F);CHKERRQ(ierr);
384   ierr = VecNorm(X, NORM_2, &xnorm);CHKERRQ(ierr);
385   ierr = SNESLineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr);
386   ierr = SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);CHKERRQ(ierr);
387   PetscFunctionReturn(0);
388 }
389 
390 #undef __FUNCT__
391 #define __FUNCT__ "SNESLineSearchView_BT"
392 PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
393 {
394   PetscErrorCode    ierr;
395   PetscBool         iascii;
396   SNESLineSearch_BT *bt;
397 
398   PetscFunctionBegin;
399   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
400   bt   = (SNESLineSearch_BT*)linesearch->data;
401   if (iascii) {
402     if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
403       ierr = PetscViewerASCIIPrintf(viewer, "  interpolation: cubic\n");CHKERRQ(ierr);
404     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
405       ierr = PetscViewerASCIIPrintf(viewer, "  interpolation: quadratic\n");CHKERRQ(ierr);
406     }
407     ierr = PetscViewerASCIIPrintf(viewer, "  alpha=%e\n", (double)bt->alpha);CHKERRQ(ierr);
408   }
409   PetscFunctionReturn(0);
410 }
411 
412 
413 #undef __FUNCT__
414 #define __FUNCT__ "SNESLineSearchDestroy_BT"
415 static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
416 {
417   PetscErrorCode ierr;
418 
419   PetscFunctionBegin;
420   ierr = PetscFree(linesearch->data);CHKERRQ(ierr);
421   PetscFunctionReturn(0);
422 }
423 
424 
425 #undef __FUNCT__
426 #define __FUNCT__ "SNESLineSearchSetFromOptions_BT"
427 static PetscErrorCode SNESLineSearchSetFromOptions_BT(PetscOptionItems *PetscOptionsObject,SNESLineSearch linesearch)
428 {
429 
430   PetscErrorCode    ierr;
431   SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
432 
433   PetscFunctionBegin;
434   ierr = PetscOptionsHead(PetscOptionsObject,"SNESLineSearch BT options");CHKERRQ(ierr);
435   ierr = PetscOptionsReal("-snes_linesearch_alpha",   "Descent tolerance",        "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL);CHKERRQ(ierr);
436   ierr = PetscOptionsTail();CHKERRQ(ierr);
437   PetscFunctionReturn(0);
438 }
439 
440 
441 #undef __FUNCT__
442 #define __FUNCT__ "SNESLineSearchCreate_BT"
443 /*MC
444    SNESLINESEARCHBT - Backtracking line search.
445 
446    This line search finds the minimum of a polynomial fitting of the L2 norm of the
447    function or the objective function if it is provided with SNESSetObjective(). If this fit does not satisfy the conditions for progress, the interval shrinks
448    and the fit is reattempted at most max_it times or until lambda is below minlambda.
449 
450    Options Database Keys:
451 +  -snes_linesearch_alpha<1e-4> - slope descent parameter
452 .  -snes_linesearch_damping<1.0> - initial step length
453 .  -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the
454                                        step is scaled back to be of this length at the beginning of the line search
455 .  -snes_linesearch_max_it<40> - maximum number of shrinking step
456 .  -snes_linesearch_minlambda<1e-12> - minimum step length allowed
457 -  -snes_linesearch_order<cubic,quadratic> - order of the approximation
458 
459    Level: advanced
460 
461    Notes:
462    This line search is taken from "Numerical Methods for Unconstrained
463    Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
464 
465 .keywords: SNES, SNESLineSearch, damping
466 
467 .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
468 M*/
469 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
470 {
471 
472   SNESLineSearch_BT *bt;
473   PetscErrorCode    ierr;
474 
475   PetscFunctionBegin;
476   linesearch->ops->apply          = SNESLineSearchApply_BT;
477   linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
478   linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
479   linesearch->ops->reset          = NULL;
480   linesearch->ops->view           = SNESLineSearchView_BT;
481   linesearch->ops->setup          = NULL;
482 
483   ierr = PetscNewLog(linesearch,&bt);CHKERRQ(ierr);
484 
485   linesearch->data    = (void*)bt;
486   linesearch->max_its = 40;
487   linesearch->order   = SNES_LINESEARCH_ORDER_CUBIC;
488   bt->alpha           = 1e-4;
489   PetscFunctionReturn(0);
490 }
491