xref: /petsc/src/snes/linesearch/impls/bt/linesearchbt.c (revision bf388a1f1fb7cbb324cb4eba80ecc0a331d4a847) 
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   PetscBool      domainerror;
71   PetscViewer    monitor;
72   PetscInt       max_its,count;
73   SNESLineSearch_BT  *bt;
74   Mat            jac;
75   SNESObjective  obj;
76 
77   PetscFunctionBegin;
78 
79   ierr = SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G);CHKERRQ(ierr);
80   ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr);
81   ierr = SNESLineSearchGetLambda(linesearch, &lambda);CHKERRQ(ierr);
82   ierr = SNESLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr);
83   ierr = SNESLineSearchGetMonitor(linesearch, &monitor);CHKERRQ(ierr);
84   ierr = SNESLineSearchGetTolerances(linesearch,&minlambda,&maxstep,PETSC_NULL,PETSC_NULL,PETSC_NULL,&max_its);CHKERRQ(ierr);
85   ierr = SNESGetTolerances(snes,PETSC_NULL,PETSC_NULL,&stol,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
86   ierr = SNESGetObjective(snes,&obj,PETSC_NULL);CHKERRQ(ierr);
87   bt = (SNESLineSearch_BT *)linesearch->data;
88 
89   alpha = bt->alpha;
90 
91   ierr = SNESGetJacobian(snes, &jac, PETSC_NULL, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr);
92 
93   if (!jac && !obj) SETERRQ(((PetscObject)linesearch)->comm, PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix");
94 
95   /* precheck */
96   ierr = SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);CHKERRQ(ierr);
97   ierr = SNESLineSearchSetSuccess(linesearch, PETSC_TRUE);CHKERRQ(ierr);
98 
99   ierr = VecNormBegin(Y, NORM_2, &ynorm);CHKERRQ(ierr);
100   ierr = VecNormBegin(X, NORM_2, &xnorm);CHKERRQ(ierr);
101   ierr = VecNormEnd(Y, NORM_2, &ynorm);CHKERRQ(ierr);
102   ierr = VecNormEnd(X, NORM_2, &xnorm);CHKERRQ(ierr);
103 
104   if (ynorm == 0.0) {
105     if (monitor) {
106       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
107       ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Initial direction and size is 0\n");CHKERRQ(ierr);
108       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
109     }
110     ierr   = VecCopy(X,W);CHKERRQ(ierr);
111     ierr   = VecCopy(F,G);CHKERRQ(ierr);
112     ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
113     PetscFunctionReturn(0);
114   }
115   if (ynorm > maxstep) {	/* Step too big, so scale back */
116     if (monitor) {
117       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
118       ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Scaling step by %14.12e old ynorm %14.12e\n", (double)(maxstep/ynorm),(double)ynorm);CHKERRQ(ierr);
119       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
120     }
121     ierr = VecScale(Y,maxstep/(ynorm));CHKERRQ(ierr);
122     ynorm = maxstep;
123   }
124 
125   /* if the SNES has an objective set, use that instead of the function value */
126   if (obj) {
127     ierr = SNESComputeObjective(snes,X,&f);CHKERRQ(ierr);
128   } else {
129     f = fnorm*fnorm;
130   }
131 
132   /* compute the initial slope */
133   if (obj) {
134     /* slope comes from the function (assumed to be the gradient of the objective */
135     ierr = VecDotRealPart(Y,F,&initslope);CHKERRQ(ierr);
136   } else {
137     /* slope comes from the normal equations */
138     ierr      = MatMult(jac,Y,W);CHKERRQ(ierr);
139     ierr      = VecDotRealPart(F,W,&initslope);CHKERRQ(ierr);
140     if (initslope > 0.0)  initslope = -initslope;
141     if (initslope == 0.0) initslope = -1.0;
142   }
143 
144   ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
145   if (linesearch->ops->viproject) {
146     ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
147   }
148   if (snes->nfuncs >= snes->max_funcs) {
149     ierr  = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr);
150     snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
151     ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
152     PetscFunctionReturn(0);
153   }
154 
155   if (obj) {
156     ierr = SNESComputeObjective(snes,W,&g);CHKERRQ(ierr);
157   } else {
158     ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr);
159     ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
160     if (domainerror) {
161       ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
162       PetscFunctionReturn(0);
163     }
164     if (linesearch->ops->vinorm) {
165       gnorm = fnorm;
166       ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
167     } else {
168       ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
169     }
170     g = PetscSqr(gnorm);
171   }
172 
173   if (PetscIsInfOrNanReal(g)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
174   if (!obj) {
175     ierr = PetscInfo2(snes,"Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);CHKERRQ(ierr);
176   }
177   if (.5*g <= .5*f + lambda*alpha*initslope) { /* Sufficient reduction or step tolerance convergence */
178     if (monitor) {
179       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
180       if (!obj) {
181         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);CHKERRQ(ierr);
182       } else {
183         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Using full step: obj0 %14.12e obj %14.12e\n", (double)f, (double)g);CHKERRQ(ierr);
184       }
185       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
186     }
187   } else {
188     /* Since the full step didn't work and the step is tiny, quit */
189     if (stol*xnorm > ynorm) {
190       ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
191       ierr = PetscInfo2(monitor,"Aborted due to ynorm < stol*xnorm (%14.12e < %14.12e) and inadequate full step.\n",ynorm,stol*xnorm);CHKERRQ(ierr);
192       PetscFunctionReturn(0);
193     }
194     /* Fit points with quadratic */
195     lambdatemp = -initslope/(g - f - 2.0*lambda*initslope);
196     lambdaprev = lambda;
197     gprev  = g;
198     if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
199     if (lambdatemp <= .1*lambda) lambda = .1*lambda;
200     else                         lambda = lambdatemp;
201 
202     ierr  = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
203     if (linesearch->ops->viproject) {
204       ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
205     }
206     if (snes->nfuncs >= snes->max_funcs) {
207       ierr  = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr);
208       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
209       ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
210       PetscFunctionReturn(0);
211     }
212     if (obj) {
213       ierr = SNESComputeObjective(snes,W,&g);CHKERRQ(ierr);
214     } else {
215       ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr);
216       ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
217       if (domainerror) {
218         PetscFunctionReturn(0);
219       }
220       if (linesearch->ops->vinorm) {
221         gnorm = fnorm;
222         ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
223       } else {
224         ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
225         g = gnorm*gnorm;
226       }
227     }
228     if (PetscIsInfOrNanReal(g)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
229     if (monitor) {
230       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
231       if (!obj) {
232         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: gnorm after quadratic fit %14.12e\n",(double)gnorm);CHKERRQ(ierr);
233       } else {
234         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: obj after quadratic fit %14.12e\n",(double)g);CHKERRQ(ierr);
235       }
236       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
237     }
238     if (.5*g < .5*f + lambda*alpha*initslope) { /* sufficient reduction */
239       if (monitor) {
240         ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
241         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda);CHKERRQ(ierr);
242         ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
243       }
244     } else {
245       /* Fit points with cubic */
246       for (count = 0; count < max_its; count++) {
247         if (lambda <= minlambda) {
248           if (monitor) {
249             ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
250             ierr = PetscViewerASCIIPrintf(monitor,"    Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr);
251             if (!obj) {
252             ierr = PetscViewerASCIIPrintf(monitor,
253                                           "    Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
254                                           (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);CHKERRQ(ierr);
255             } else {
256               ierr = PetscViewerASCIIPrintf(monitor,
257                                             "    Line search: obj(0)=%18.16e, obj=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
258                                             (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);CHKERRQ(ierr);
259             }
260             ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
261           }
262           ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
263           PetscFunctionReturn(0);
264         }
265         if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
266           t1 = .5*(g - f) - lambda*initslope;
267           t2 = .5*(gprev  - f) - lambdaprev*initslope;
268           a  = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
269           b  = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
270           d  = b*b - 3*a*initslope;
271           if (d < 0.0) d = 0.0;
272           if (a == 0.0) {
273             lambdatemp = -initslope/(2.0*b);
274           } else {
275             lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a);
276           }
277         } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
278           lambdatemp = -initslope/(g - f - 2.0*initslope);
279         } else {
280           SETERRQ(((PetscObject)linesearch)->comm, PETSC_ERR_SUP, "unsupported line search order for type bt");
281         }
282           lambdaprev = lambda;
283           gprev  = g;
284         if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
285         if (lambdatemp <= .1*lambda) lambda     = .1*lambda;
286         else                         lambda     = lambdatemp;
287         ierr  = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
288         if (linesearch->ops->viproject) {
289           ierr = (*linesearch->ops->viproject)(snes,W);CHKERRQ(ierr);
290         }
291         if (snes->nfuncs >= snes->max_funcs) {
292           ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr);
293           if (!obj) {
294             ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
295                               (double)fnorm,(double)gnorm,(double)ynorm,(double)lambda,(double)initslope);CHKERRQ(ierr);
296           }
297           ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
298           snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
299           PetscFunctionReturn(0);
300         }
301         if (obj) {
302           ierr = SNESComputeObjective(snes,W,&g);CHKERRQ(ierr);
303         } else {
304           ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr);
305           ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
306           if (domainerror) {
307             ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
308             PetscFunctionReturn(0);
309           }
310           if (linesearch->ops->vinorm) {
311             gnorm = fnorm;
312             ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
313           } else {
314             ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
315             g = gnorm*gnorm;
316           }
317         }
318         if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
319         if (.5*g < .5*f + lambda*alpha*initslope) { /* is reduction enough? */
320           if (monitor) {
321             ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
322             if (!obj) {
323               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
324                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);CHKERRQ(ierr);
325               } else {
326                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);CHKERRQ(ierr);
327               }
328               ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
329             } else {
330               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
331                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
332               } else {
333                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
334               }
335               ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
336             }
337           }
338           break;
339         } else {
340           if (monitor) {
341             ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
342             if (!obj) {
343               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
344                 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);
345               } else {
346                 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);
347               }
348               ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
349             } else {
350               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
351                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
352               } else {
353                 ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);CHKERRQ(ierr);
354               }
355               ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
356             }
357           }
358         }
359       }
360     }
361   }
362 
363   /* postcheck */
364   ierr = SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);CHKERRQ(ierr);
365   if (changed_y) {
366     ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
367     if (linesearch->ops->viproject) {
368       ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
369     }
370   }
371   if (changed_y || changed_w || obj) { /* recompute the function norm if the step has changed or the objective isn't the norm */
372     ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr);
373     ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
374     if (domainerror) {
375       ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
376       PetscFunctionReturn(0);
377     }
378     if (linesearch->ops->vinorm) {
379       gnorm = fnorm;
380       ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
381     } else {
382       ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
383     }
384     ierr = VecNorm(Y,NORM_2,&ynorm);CHKERRQ(ierr);
385     if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
386 
387   }
388 
389   /* copy the solution over */
390   ierr = VecCopy(W, X);CHKERRQ(ierr);
391   ierr = VecCopy(G, F);CHKERRQ(ierr);
392   ierr = VecNorm(X, NORM_2, &xnorm);CHKERRQ(ierr);
393   ierr = SNESLineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr);
394   ierr = SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);CHKERRQ(ierr);
395   PetscFunctionReturn(0);
396 }
397 
398 #undef __FUNCT__
399 #define __FUNCT__ "SNESLineSearchView_BT"
400 PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
401 {
402   PetscErrorCode    ierr;
403   PetscBool         iascii;
404   SNESLineSearch_BT *bt;
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   PetscFunctionBegin;
440 
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