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