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