xref: /petsc/src/snes/linesearch/impls/bt/linesearchbt.c (revision bebe2cf65d55febe21a5af8db2bd2e168caaa2e7) ! 
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   PetscViewer       monitor;
73   PetscInt          max_its,count;
74   SNESLineSearch_BT *bt;
75   Mat               jac;
76   PetscErrorCode    (*objective)(SNES,Vec,PetscReal*,void*);
77 
78   PetscFunctionBegin;
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,NULL,NULL,NULL,&max_its);CHKERRQ(ierr);
85   ierr = SNESGetTolerances(snes,NULL,NULL,&stol,NULL,NULL);CHKERRQ(ierr);
86   ierr = SNESGetObjective(snes,&objective,NULL);CHKERRQ(ierr);
87   bt   = (SNESLineSearch_BT*)linesearch->data;
88 
89   alpha = bt->alpha;
90 
91   ierr = SNESGetJacobian(snes, &jac, NULL, NULL, NULL);CHKERRQ(ierr);
92 
93   if (!jac && !objective) SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix");
94 
95   /* precheck */
96   ierr = SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);CHKERRQ(ierr);
97   ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);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 = SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm);CHKERRQ(ierr);
113     ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);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         = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);CHKERRQ(ierr);
153     PetscFunctionReturn(0);
154   }
155 
156   if (objective) {
157     ierr = SNESComputeObjective(snes,W,&g);CHKERRQ(ierr);
158   } else {
159     ierr = (*linesearch->ops->snesfunc)(snes,W,G);CHKERRQ(ierr);
160     if (linesearch->ops->vinorm) {
161       gnorm = fnorm;
162       ierr  = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
163     } else {
164       ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
165     }
166     g = PetscSqr(gnorm);
167   }
168 
169   if (PetscIsInfOrNanReal(g)) {
170     ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF);CHKERRQ(ierr);
171     ierr = PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");CHKERRQ(ierr);
172     PetscFunctionReturn(0);
173   }
174   if (!objective) {
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 (!objective) {
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 = SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm);CHKERRQ(ierr);
191       ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);CHKERRQ(ierr);
192       if (monitor) {
193         ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
194         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Aborted due to ynorm < stol*xnorm (%14.12e < %14.12e) and inadequate full step.\n",(double)ynorm,(double)stol*xnorm);CHKERRQ(ierr);
195         ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
196       }
197       PetscFunctionReturn(0);
198     }
199     /* Fit points with quadratic */
200     lambdatemp = -initslope/(g - f - 2.0*lambda*initslope);
201     lambdaprev = lambda;
202     gprev      = g;
203     if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
204     if (lambdatemp <= .1*lambda) lambda = .1*lambda;
205     else                         lambda = lambdatemp;
206 
207     ierr  = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
208     if (linesearch->ops->viproject) {
209       ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
210     }
211     if (snes->nfuncs >= snes->max_funcs) {
212       ierr         = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr);
213       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
214       ierr         = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);CHKERRQ(ierr);
215       PetscFunctionReturn(0);
216     }
217     if (objective) {
218       ierr = SNESComputeObjective(snes,W,&g);CHKERRQ(ierr);
219     } else {
220       ierr = (*linesearch->ops->snesfunc)(snes,W,G);CHKERRQ(ierr);
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       }
227       g = PetscSqr(gnorm);
228     }
229     if (PetscIsInfOrNanReal(g)) {
230       ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF);CHKERRQ(ierr);
231       ierr = PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");CHKERRQ(ierr);
232       PetscFunctionReturn(0);
233     }
234     if (monitor) {
235       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
236       if (!objective) {
237         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: gnorm after quadratic fit %14.12e\n",(double)gnorm);CHKERRQ(ierr);
238       } else {
239         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: obj after quadratic fit %14.12e\n",(double)g);CHKERRQ(ierr);
240       }
241       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
242     }
243     if (.5*g < .5*f + lambda*alpha*initslope) { /* sufficient reduction */
244       if (monitor) {
245         ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
246         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda);CHKERRQ(ierr);
247         ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
248       }
249     } else {
250       /* Fit points with cubic */
251       for (count = 0; count < max_its; count++) {
252         if (lambda <= minlambda) {
253           if (monitor) {
254             ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
255             ierr = PetscViewerASCIIPrintf(monitor,"    Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr);
256             if (!objective) {
257               ierr = PetscViewerASCIIPrintf(monitor,
258                                             "    Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
259                                             (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);CHKERRQ(ierr);
260             } else {
261               ierr = PetscViewerASCIIPrintf(monitor,
262                                             "    Line search: obj(0)=%18.16e, obj=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
263                                             (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);CHKERRQ(ierr);
264             }
265             ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
266           }
267           ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);CHKERRQ(ierr);
268           PetscFunctionReturn(0);
269         }
270         if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
271           t1 = .5*(g - f) - lambda*initslope;
272           t2 = .5*(gprev  - f) - lambdaprev*initslope;
273           a  = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
274           b  = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
275           d  = b*b - 3*a*initslope;
276           if (d < 0.0) d = 0.0;
277           if (a == 0.0) lambdatemp = -initslope/(2.0*b);
278           else lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a);
279 
280         } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
281           lambdatemp = -initslope/(g - f - 2.0*initslope);
282         } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt");
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         = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);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 = (*linesearch->ops->snesfunc)(snes,W,G);CHKERRQ(ierr);
306           if (linesearch->ops->vinorm) {
307             gnorm = fnorm;
308             ierr  = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
309           } else {
310             ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
311           }
312           g = PetscSqr(gnorm);
313         }
314         if (PetscIsInfOrNanReal(gnorm)) {
315           ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF);CHKERRQ(ierr);
316           ierr = PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");CHKERRQ(ierr);
317           PetscFunctionReturn(0);
318         }
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 (!objective) {
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 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   /* postcheck */
362   ierr = SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);CHKERRQ(ierr);
363   if (changed_y) {
364     ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
365     if (linesearch->ops->viproject) {
366       ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
367     }
368   }
369   if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
370     ierr = (*linesearch->ops->snesfunc)(snes,W,G);CHKERRQ(ierr);
371     if (linesearch->ops->vinorm) {
372       gnorm = fnorm;
373       ierr  = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
374     } else {
375       ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
376     }
377     ierr = VecNorm(Y,NORM_2,&ynorm);CHKERRQ(ierr);
378     if (PetscIsInfOrNanReal(gnorm)) {
379       ierr = SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_FAILED_NANORINF);CHKERRQ(ierr);
380       ierr = PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");CHKERRQ(ierr);
381       PetscFunctionReturn(0);
382     }
383   }
384 
385   /* copy the solution over */
386   ierr = VecCopy(W, X);CHKERRQ(ierr);
387   ierr = VecCopy(G, F);CHKERRQ(ierr);
388   ierr = VecNorm(X, NORM_2, &xnorm);CHKERRQ(ierr);
389   ierr = SNESLineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr);
390   ierr = SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);CHKERRQ(ierr);
391   PetscFunctionReturn(0);
392 }
393 
394 #undef __FUNCT__
395 #define __FUNCT__ "SNESLineSearchView_BT"
396 PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
397 {
398   PetscErrorCode    ierr;
399   PetscBool         iascii;
400   SNESLineSearch_BT *bt;
401 
402   PetscFunctionBegin;
403   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
404   bt   = (SNESLineSearch_BT*)linesearch->data;
405   if (iascii) {
406     if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
407       ierr = PetscViewerASCIIPrintf(viewer, "  interpolation: cubic\n");CHKERRQ(ierr);
408     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
409       ierr = PetscViewerASCIIPrintf(viewer, "  interpolation: quadratic\n");CHKERRQ(ierr);
410     }
411     ierr = PetscViewerASCIIPrintf(viewer, "  alpha=%e\n", (double)bt->alpha);CHKERRQ(ierr);
412   }
413   PetscFunctionReturn(0);
414 }
415 
416 
417 #undef __FUNCT__
418 #define __FUNCT__ "SNESLineSearchDestroy_BT"
419 static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
420 {
421   PetscErrorCode ierr;
422 
423   PetscFunctionBegin;
424   ierr = PetscFree(linesearch->data);CHKERRQ(ierr);
425   PetscFunctionReturn(0);
426 }
427 
428 
429 #undef __FUNCT__
430 #define __FUNCT__ "SNESLineSearchSetFromOptions_BT"
431 static PetscErrorCode SNESLineSearchSetFromOptions_BT(PetscOptions *PetscOptionsObject,SNESLineSearch linesearch)
432 {
433 
434   PetscErrorCode    ierr;
435   SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
436 
437   PetscFunctionBegin;
438   ierr = PetscOptionsHead(PetscOptionsObject,"SNESLineSearch BT options");CHKERRQ(ierr);
439   ierr = PetscOptionsReal("-snes_linesearch_alpha",   "Descent tolerance",        "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL);CHKERRQ(ierr);
440   ierr = PetscOptionsTail();CHKERRQ(ierr);
441   PetscFunctionReturn(0);
442 }
443 
444 
445 #undef __FUNCT__
446 #define __FUNCT__ "SNESLineSearchCreate_BT"
447 /*MC
448    SNESLINESEARCHBT - Backtracking line search.
449 
450    This line search finds the minimum of a polynomial fitting of the L2 norm of the
451    function or the objective function if it is provided with SNESSetObjective(). If this fit does not satisfy the conditions for progress, the interval shrinks
452    and the fit is reattempted at most max_it times or until lambda is below minlambda.
453 
454    Options Database Keys:
455 +  -snes_linesearch_alpha<1e-4> - slope descent parameter
456 .  -snes_linesearch_damping<1.0> - initial step length
457 .  -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the
458                                        step is scaled back to be of this length at the beginning of the line search
459 .  -snes_linesearch_max_it<40> - maximum number of shrinking step
460 .  -snes_linesearch_minlambda<1e-12> - minimum step length allowed
461 -  -snes_linesearch_order<cubic,quadratic> - order of the approximation
462 
463    Level: advanced
464 
465    Notes:
466    This line search is taken from "Numerical Methods for Unconstrained
467    Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
468 
469 .keywords: SNES, SNESLineSearch, damping
470 
471 .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
472 M*/
473 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
474 {
475 
476   SNESLineSearch_BT *bt;
477   PetscErrorCode    ierr;
478 
479   PetscFunctionBegin;
480   linesearch->ops->apply          = SNESLineSearchApply_BT;
481   linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
482   linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
483   linesearch->ops->reset          = NULL;
484   linesearch->ops->view           = SNESLineSearchView_BT;
485   linesearch->ops->setup          = NULL;
486 
487   ierr = PetscNewLog(linesearch,&bt);CHKERRQ(ierr);
488 
489   linesearch->data    = (void*)bt;
490   linesearch->max_its = 40;
491   linesearch->order   = SNES_LINESEARCH_ORDER_CUBIC;
492   bt->alpha           = 1e-4;
493   PetscFunctionReturn(0);
494 }
495