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