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