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