xref: /petsc/src/snes/linesearch/impls/bt/linesearchbt.c (revision 670f3ff9a12f3667ff7d4cebfc50ebc8579c9e33) 
1 #include <private/linesearchimpl.h>
2 #include <private/snesimpl.h>
3 
4 typedef struct {
5   PetscReal        alpha; /* sufficient decrease parameter */
6 } SNESLineSearch_BT;
7 
8 #undef __FUNCT__
9 #define __FUNCT__ "SNESLineSearchApply_BT"
10 static PetscErrorCode  SNESLineSearchApply_BT(SNESLineSearch linesearch)
11 {
12   PetscBool      changed_y, changed_w;
13   PetscErrorCode ierr;
14   Vec            X, F, Y, W, G;
15   SNES           snes;
16   PetscReal      fnorm, xnorm, ynorm, gnorm, gnormprev;
17   PetscReal      lambda, lambdatemp, lambdaprev, minlambda, maxstep, rellength, initslope, alpha;
18   PetscReal      t1, t2, a, b, d, steptol;
19 #if defined(PETSC_USE_COMPLEX)
20   PetscScalar    cinitslope;
21 #endif
22   PetscBool      domainerror;
23   PetscViewer    monitor;
24   PetscInt       max_its, count;
25   SNESLineSearch_BT  *bt;
26   Mat            jac;
27 
28 
29   PetscFunctionBegin;
30 
31   ierr = SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G);CHKERRQ(ierr);
32   ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr);
33   ierr = SNESLineSearchGetLambda(linesearch, &lambda);CHKERRQ(ierr);
34   ierr = SNESLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr);
35   ierr = SNESLineSearchGetMonitor(linesearch, &monitor);CHKERRQ(ierr);
36   ierr = SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, PETSC_NULL, PETSC_NULL, PETSC_NULL, &max_its);
37   bt = (SNESLineSearch_BT *)linesearch->data;
38 
39   alpha = bt->alpha;
40 
41   ierr = SNESGetJacobian(snes, &jac, PETSC_NULL, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr);
42   if (!jac) {
43     SETERRQ(((PetscObject)linesearch)->comm, PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix");
44   }
45   /* precheck */
46   ierr = SNESLineSearchPreCheck(linesearch, &changed_y);CHKERRQ(ierr);
47   ierr = SNESLineSearchSetSuccess(linesearch, PETSC_TRUE);CHKERRQ(ierr);
48 
49   ierr = VecNorm(Y, NORM_2, &ynorm);CHKERRQ(ierr);
50   if (ynorm == 0.0) {
51     if (monitor) {
52       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
53       ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Initial direction and size is 0\n");CHKERRQ(ierr);
54       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
55     }
56     ierr   = VecCopy(X,W);CHKERRQ(ierr);
57     ierr   = VecCopy(F,G);CHKERRQ(ierr);
58     ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
59     PetscFunctionReturn(0);
60   }
61   if (ynorm > maxstep) {	/* Step too big, so scale back */
62     if (monitor) {
63       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
64       ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Scaling step by %14.12e old ynorm %14.12e\n", (maxstep/ynorm),ynorm);CHKERRQ(ierr);
65       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
66     }
67     ierr = VecScale(Y,maxstep/(ynorm));CHKERRQ(ierr);
68     ynorm = maxstep;
69   }
70   ierr      = VecMaxPointwiseDivide(Y,X,&rellength);CHKERRQ(ierr);
71   minlambda = steptol/rellength;
72   ierr      = MatMult(jac,Y,W);CHKERRQ(ierr);
73 #if defined(PETSC_USE_COMPLEX)
74   ierr      = VecDot(F,W,&cinitslope);CHKERRQ(ierr);
75   initslope = PetscRealPart(cinitslope);
76 #else
77   ierr      = VecDot(F,W,&initslope);CHKERRQ(ierr);
78 #endif
79   if (initslope > 0.0)  initslope = -initslope;
80   if (initslope == 0.0) initslope = -1.0;
81 
82   ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
83   if (linesearch->ops->viproject) {
84     ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
85   }
86   if (snes->nfuncs >= snes->max_funcs) {
87     ierr  = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr);
88     snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
89     ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
90     PetscFunctionReturn(0);
91   }
92   ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr);
93   ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
94   if (domainerror) {
95     ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
96     PetscFunctionReturn(0);
97   }
98   if (linesearch->ops->vinorm) {
99     gnorm = fnorm;
100     ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
101   } else {
102     ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
103   }
104 
105   if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
106   ierr = PetscInfo2(snes,"Initial fnorm %14.12e gnorm %14.12e\n", fnorm, gnorm);CHKERRQ(ierr);
107   if (.5*gnorm*gnorm <= .5*fnorm*fnorm + lambda*alpha*initslope) { /* Sufficient reduction */
108     if (monitor) {
109       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
110       ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", fnorm, gnorm);CHKERRQ(ierr);
111       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
112     }
113   } else {
114     /* Fit points with quadratic */
115     lambdatemp = -initslope/(gnorm*gnorm - fnorm*fnorm - 2.0*lambda*initslope);
116     lambdaprev = lambda;
117     gnormprev  = gnorm;
118     if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
119     if (lambdatemp <= .1*lambda) lambda = .1*lambda;
120     else                         lambda = lambdatemp;
121 
122     ierr  = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
123     if (linesearch->ops->viproject) {
124       ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
125     }
126     if (snes->nfuncs >= snes->max_funcs) {
127       ierr  = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr);
128       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
129       ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
130       PetscFunctionReturn(0);
131     }
132     ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr);
133     ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
134     if (domainerror) {
135       PetscFunctionReturn(0);
136     }
137     if (linesearch->ops->vinorm) {
138       gnorm = fnorm;
139       ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
140     } else {
141       ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
142     }
143     if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
144     if (monitor) {
145       ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
146       ierr = PetscViewerASCIIPrintf(monitor,"    Line search: gnorm after quadratic fit %14.12e\n",gnorm);CHKERRQ(ierr);
147       ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
148     }
149     if (.5*gnorm*gnorm < .5*fnorm*fnorm + lambda*alpha*initslope) { /* sufficient reduction */
150       if (monitor) {
151         ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
152         ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda);CHKERRQ(ierr);
153         ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
154       }
155     } else {
156       /* Fit points with cubic */
157       for (count = 0; count < max_its; count++) {
158         if (lambda <= minlambda) {
159           if (monitor) {
160             ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
161             ierr = PetscViewerASCIIPrintf(monitor,"    Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr);
162             ierr = PetscViewerASCIIPrintf(monitor,
163                                           "    Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
164                                           fnorm, gnorm, ynorm, minlambda, lambda, initslope);CHKERRQ(ierr);
165             ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
166           }
167           ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
168           PetscFunctionReturn(0);
169         }
170         if (linesearch->order == SNES_LINESEARCH_CUBIC) {
171           t1 = .5*(gnorm*gnorm - fnorm*fnorm) - lambda*initslope;
172           t2 = .5*(gnormprev*gnormprev  - fnorm*fnorm) - lambdaprev*initslope;
173           a  = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
174           b  = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
175           d  = b*b - 3*a*initslope;
176           if (d < 0.0) d = 0.0;
177           if (a == 0.0) {
178             lambdatemp = -initslope/(2.0*b);
179           } else {
180             lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a);
181           }
182         } else if (linesearch->order == SNES_LINESEARCH_QUADRATIC) {
183           lambdatemp = -initslope/(gnorm*gnorm - fnorm*fnorm - 2.0*initslope);
184         } else {
185           SETERRQ(((PetscObject)linesearch)->comm, PETSC_ERR_SUP, "unsupported line search order for type bt");
186         }
187           lambdaprev = lambda;
188           gnormprev  = gnorm;
189         if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
190         if (lambdatemp <= .1*lambda) lambda     = .1*lambda;
191         else                         lambda     = lambdatemp;
192         ierr  = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
193         if (snes->nfuncs >= snes->max_funcs) {
194           ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr);
195           ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
196                             fnorm,gnorm,ynorm,lambda,initslope);CHKERRQ(ierr);
197           ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
198           snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
199           PetscFunctionReturn(0);
200         }
201         ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr);
202         ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
203         if (domainerror) {
204           ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
205           PetscFunctionReturn(0);
206         }
207         if (linesearch->ops->vinorm) {
208           gnorm = fnorm;
209           ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
210         } else {
211           ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
212         }
213         if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
214         if (.5*gnorm*gnorm < .5*fnorm*fnorm + lambda*alpha*initslope) { /* is reduction enough? */
215           if (monitor) {
216             ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
217             if (linesearch->order == SNES_LINESEARCH_CUBIC) {
218               ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr);
219             } else {
220               ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr);
221             }
222             ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
223           }
224           break;
225         } else {
226           if (monitor) {
227             ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
228             if (linesearch->order == SNES_LINESEARCH_CUBIC) {
229               ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr);
230             } else {
231               ierr = PetscViewerASCIIPrintf(monitor,"    Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr);
232             }
233             ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
234           }
235         }
236       }
237     }
238   }
239 
240   /* postcheck */
241   ierr = SNESLineSearchPostCheck(linesearch, &changed_y, &changed_w);CHKERRQ(ierr);
242   if (changed_y) {
243     ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
244     if (linesearch->ops->viproject) {
245       ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
246     }
247   }
248   if (changed_y || changed_w) { /* recompute the function if the step has changed */
249     ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr);
250     ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
251     if (domainerror) {
252       ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr);
253       PetscFunctionReturn(0);
254     }
255     if (linesearch->ops->vinorm) {
256       gnorm = fnorm;
257       ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr);
258     } else {
259       ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr);
260     }
261     ierr = VecNorm(Y,NORM_2,&ynorm);CHKERRQ(ierr);
262     if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
263 
264   }
265 
266   /* copy the solution over */
267   ierr = VecCopy(W, X);CHKERRQ(ierr);
268   ierr = VecCopy(G, F);CHKERRQ(ierr);
269   ierr = VecNorm(X, NORM_2, &xnorm);CHKERRQ(ierr);
270   ierr = SNESLineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr);
271   ierr = SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);CHKERRQ(ierr);
272   PetscFunctionReturn(0);
273 }
274 
275 #undef __FUNCT__
276 #define __FUNCT__ "SNESLineSearchView_BT"
277 PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
278 {
279   PetscErrorCode    ierr;
280   PetscBool         iascii;
281   SNESLineSearch_BT *bt;
282   PetscFunctionBegin;
283   ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
284   bt = (SNESLineSearch_BT*)linesearch->data;
285   if (iascii) {
286     if (linesearch->order == SNES_LINESEARCH_CUBIC) {
287     ierr = PetscViewerASCIIPrintf(viewer, "  interpolation: cubic\n");CHKERRQ(ierr);
288     } else if (linesearch->order == SNES_LINESEARCH_QUADRATIC) {
289     ierr = PetscViewerASCIIPrintf(viewer, "  interpolation: quadratic\n");CHKERRQ(ierr);
290     }
291     ierr = PetscViewerASCIIPrintf(viewer, "  alpha=%G\n", bt->alpha);CHKERRQ(ierr);
292   }
293   PetscFunctionReturn(0);
294 }
295 
296 
297 #undef __FUNCT__
298 #define __FUNCT__ "SNESLineSearchDestroy_BT"
299 static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
300 {
301   PetscErrorCode ierr;
302 
303   PetscFunctionBegin;
304   ierr = PetscFree(linesearch->data);CHKERRQ(ierr);
305   PetscFunctionReturn(0);
306 }
307 
308 
309 #undef __FUNCT__
310 #define __FUNCT__ "SNESLineSearchSetFromOptions_BT"
311 static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch)
312 {
313 
314   PetscErrorCode       ierr;
315   SNESLineSearch_BT    *bt;
316   PetscFunctionBegin;
317 
318   bt = (SNESLineSearch_BT*)linesearch->data;
319 
320   ierr = PetscOptionsHead("SNESLineSearch BT options");CHKERRQ(ierr);
321   ierr = PetscOptionsReal("-snes_linesearch_alpha",   "Descent tolerance",        "SNESLineSearchBT", bt->alpha, &bt->alpha, PETSC_NULL);CHKERRQ(ierr);
322 
323   ierr = PetscOptionsTail();CHKERRQ(ierr);
324   PetscFunctionReturn(0);
325 }
326 
327 
328 #undef __FUNCT__
329 #define __FUNCT__ "SNESLineSearchCreate_BT"
330 /*MC
331    SNES_LINESEARCH_BT - Backtracking line searches.
332 
333    These linesearches try a polynomial fit for the L2 norm of the error
334    using the gradient.  Failing that, they step back and try again.
335 
336    Options Database Keys:
337 +  -snes_linesearch_alpha<1e-4> - slope descent parameter
338 .  -snes_linesearch_damping<1.0> - full steplength
339 -  -snes_linesearch_order<cubic, quadratic> - order of the approximation
340 
341    Level: advanced
342 
343 .keywords: SNES, SNESLineSearch, damping
344 
345 .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
346 M*/
347 PETSC_EXTERN_C PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
348 {
349 
350   SNESLineSearch_BT  *bt;
351   PetscErrorCode ierr;
352 
353   PetscFunctionBegin;
354   linesearch->ops->apply          = SNESLineSearchApply_BT;
355   linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
356   linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
357   linesearch->ops->reset          = PETSC_NULL;
358   linesearch->ops->view           = SNESLineSearchView_BT;
359   linesearch->ops->setup          = PETSC_NULL;
360 
361   ierr = PetscNewLog(linesearch, SNESLineSearch_BT, &bt);CHKERRQ(ierr);
362   linesearch->data = (void *)bt;
363   linesearch->max_its = 40;
364   linesearch->order = SNES_LINESEARCH_CUBIC;
365   bt->alpha = 1e-4;
366 
367   PetscFunctionReturn(0);
368 }
369