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