xref: /petsc/src/snes/linesearch/impls/bt/linesearchbt.c (revision d54338ec4fdf48b56b5efdd9bf04ca9d50d6ab04)
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 {
281           SETERRQ(((PetscObject)linesearch)->comm, PETSC_ERR_SUP, "unsupported line search order for type bt");
282         }
283           lambdaprev = lambda;
284           gprev  = g;
285         if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
286         if (lambdatemp <= .1*lambda) lambda     = .1*lambda;
287         else                         lambda     = lambdatemp;
288         ierr  = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
289         if (linesearch->ops->viproject) {
290           ierr = (*linesearch->ops->viproject)(snes,W);CHKERRQ(ierr);
291         }
292         if (snes->nfuncs >= snes->max_funcs) {
293           ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr);
294           if (!objective) {
295             ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
296                               (double)fnorm,(double)gnorm,(double)ynorm,(double)lambda,(double)initslope);CHKERRQ(ierr);
297           }
298           ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
299           snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
300           PetscFunctionReturn(0);
301         }
302         if (objective) {
303           ierr = SNESComputeObjective(snes,W,&g);CHKERRQ(ierr);
304         } else {
305           ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr);
306           ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
307           if (domainerror) {
308             ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
309             PetscFunctionReturn(0);
310           }
311           if (linesearch->ops->vinorm) {
312             gnorm = fnorm;
313             ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
314           } else {
315             ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
316             g = gnorm*gnorm;
317           }
318         }
319         if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
320         if (.5*g < .5*f + lambda*alpha*initslope) { /* is reduction enough? */
321           if (monitor) {
322             ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
323             if (!objective) {
324               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
325                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);CHKERRQ(ierr);
326               } else {
327                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);CHKERRQ(ierr);
328               }
329               ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
330             } else {
331               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
332                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
333               } else {
334                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
335               }
336               ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
337             }
338           }
339           break;
340         } else {
341           if (monitor) {
342             ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
343             if (!objective) {
344               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
345                 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);
346               } else {
347                 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);
348               }
349               ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
350             } else {
351               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
352                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
353               } else {
354                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
355               }
356               ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
357             }
358           }
359         }
360       }
361     }
362   }
363 
364   /* postcheck */
365   ierr = SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);CHKERRQ(ierr);
366   if (changed_y) {
367     ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
368     if (linesearch->ops->viproject) {
369       ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
370     }
371   }
372   if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
373     ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr);
374     ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
375     if (domainerror) {
376       ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
377       PetscFunctionReturn(0);
378     }
379     if (linesearch->ops->vinorm) {
380       gnorm = fnorm;
381       ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
382     } else {
383       ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
384     }
385     ierr = VecNorm(Y,NORM_2,&ynorm);CHKERRQ(ierr);
386     if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
387 
388   }
389 
390   /* copy the solution over */
391   ierr = VecCopy(W, X);CHKERRQ(ierr);
392   ierr = VecCopy(G, F);CHKERRQ(ierr);
393   ierr = VecNorm(X, NORM_2, &xnorm);CHKERRQ(ierr);
394   ierr = SNESLineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr);
395   ierr = SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);CHKERRQ(ierr);
396   PetscFunctionReturn(0);
397 }
398 
399 #undef __FUNCT__
400 #define __FUNCT__ "SNESLineSearchView_BT"
401 PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
402 {
403   PetscErrorCode    ierr;
404   PetscBool         iascii;
405   SNESLineSearch_BT *bt;
406 
407   PetscFunctionBegin;
408   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
409   bt = (SNESLineSearch_BT*)linesearch->data;
410   if (iascii) {
411     if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
412     ierr = PetscViewerASCIIPrintf(viewer, "  interpolation: cubic\n");CHKERRQ(ierr);
413     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
414     ierr = PetscViewerASCIIPrintf(viewer, "  interpolation: quadratic\n");CHKERRQ(ierr);
415     }
416     ierr = PetscViewerASCIIPrintf(viewer, "  alpha=%e\n", bt->alpha);CHKERRQ(ierr);
417   }
418   PetscFunctionReturn(0);
419 }
420 
421 
422 #undef __FUNCT__
423 #define __FUNCT__ "SNESLineSearchDestroy_BT"
424 static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
425 {
426   PetscErrorCode ierr;
427 
428   PetscFunctionBegin;
429   ierr = PetscFree(linesearch->data);CHKERRQ(ierr);
430   PetscFunctionReturn(0);
431 }
432 
433 
434 #undef __FUNCT__
435 #define __FUNCT__ "SNESLineSearchSetFromOptions_BT"
436 static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch)
437 {
438 
439   PetscErrorCode       ierr;
440   SNESLineSearch_BT    *bt;
441 
442   PetscFunctionBegin;
443   bt = (SNESLineSearch_BT*)linesearch->data;
444 
445   ierr = PetscOptionsHead("SNESLineSearch BT options");CHKERRQ(ierr);
446   ierr = PetscOptionsReal("-snes_linesearch_alpha",   "Descent tolerance",        "SNESLineSearchBT", bt->alpha, &bt->alpha, PETSC_NULL);CHKERRQ(ierr);
447 
448   ierr = PetscOptionsTail();CHKERRQ(ierr);
449   PetscFunctionReturn(0);
450 }
451 
452 
453 #undef __FUNCT__
454 #define __FUNCT__ "SNESLineSearchCreate_BT"
455 /*MC
456    SNESLINESEARCHBT - Backtracking line search.
457 
458    This line search finds the minimum of a polynomial fitting of the L2 norm of the
459    function. If this fit does not satisfy the conditions for progress, the interval shrinks
460    and the fit is reattempted at most max_it times or until lambda is below minlambda.
461 
462    Options Database Keys:
463 +  -snes_linesearch_alpha<1e-4> - slope descent parameter
464 .  -snes_linesearch_damping<1.0> - initial step length
465 .  -snes_linesearch_max_it<40> - maximum number of shrinking step
466 .  -snes_linesearch_minlambda<1e-12> - minimum step length allowed
467 -  -snes_linesearch_order<cubic,quadratic> - order of the approximation
468 
469    Level: advanced
470 
471    Notes:
472    This line search is taken from "Numerical Methods for Unconstrained
473    Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
474 
475 .keywords: SNES, SNESLineSearch, damping
476 
477 .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
478 M*/
479 PETSC_EXTERN_C PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
480 {
481 
482   SNESLineSearch_BT  *bt;
483   PetscErrorCode ierr;
484 
485   PetscFunctionBegin;
486   linesearch->ops->apply          = SNESLineSearchApply_BT;
487   linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
488   linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
489   linesearch->ops->reset          = PETSC_NULL;
490   linesearch->ops->view           = SNESLineSearchView_BT;
491   linesearch->ops->setup          = PETSC_NULL;
492 
493   ierr = PetscNewLog(linesearch, SNESLineSearch_BT, &bt);CHKERRQ(ierr);
494   linesearch->data = (void *)bt;
495   linesearch->max_its = 40;
496   linesearch->order = SNES_LINESEARCH_ORDER_CUBIC;
497   bt->alpha = 1e-4;
498 
499   PetscFunctionReturn(0);
500 }
501