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