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