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