Actual source code: linesearchbt.c
petsc-3.8.4 2018-03-24
1: #include <petsc/private/linesearchimpl.h>
2: #include <petsc/private/snesimpl.h>
4: typedef struct {
5: PetscReal alpha; /* sufficient decrease parameter */
6: } SNESLineSearch_BT;
8: /*@
9: SNESLineSearchBTSetAlpha - Sets the descent parameter, alpha, in the BT linesearch variant.
11: Input Parameters:
12: + linesearch - linesearch context
13: - alpha - The descent parameter
15: Level: intermediate
17: .seealso: SNESLineSearchSetLambda(), SNESLineSearchGetTolerances() SNESLINESEARCHBT
18: @*/
19: PetscErrorCode SNESLineSearchBTSetAlpha(SNESLineSearch linesearch, PetscReal alpha)
20: {
21: SNESLineSearch_BT *bt;
25: bt = (SNESLineSearch_BT*)linesearch->data;
26: bt->alpha = alpha;
27: return(0);
28: }
31: /*@
32: SNESLineSearchBTGetAlpha - Gets the descent parameter, alpha, in the BT linesearch variant.
34: Input Parameters:
35: . linesearch - linesearch context
37: Output Parameters:
38: . alpha - The descent parameter
40: Level: intermediate
42: .seealso: SNESLineSearchGetLambda(), SNESLineSearchGetTolerances() SNESLINESEARCHBT
43: @*/
44: PetscErrorCode SNESLineSearchBTGetAlpha(SNESLineSearch linesearch, PetscReal *alpha)
45: {
46: SNESLineSearch_BT *bt;
50: bt = (SNESLineSearch_BT*)linesearch->data;
51: *alpha = bt->alpha;
52: return(0);
53: }
55: static PetscErrorCode SNESLineSearchApply_BT(SNESLineSearch linesearch)
56: {
57: PetscBool changed_y,changed_w;
58: PetscErrorCode ierr;
59: Vec X,F,Y,W,G;
60: SNES snes;
61: PetscReal fnorm, xnorm, ynorm, gnorm;
62: PetscReal lambda,lambdatemp,lambdaprev,minlambda,maxstep,initslope,alpha,stol;
63: PetscReal t1,t2,a,b,d;
64: PetscReal f;
65: PetscReal g,gprev;
66: PetscViewer monitor;
67: PetscInt max_its,count;
68: SNESLineSearch_BT *bt;
69: Mat jac;
70: PetscErrorCode (*objective)(SNES,Vec,PetscReal*,void*);
73: SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G);
74: SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);
75: SNESLineSearchGetLambda(linesearch, &lambda);
76: SNESLineSearchGetSNES(linesearch, &snes);
77: SNESLineSearchGetDefaultMonitor(linesearch, &monitor);
78: SNESLineSearchGetTolerances(linesearch,&minlambda,&maxstep,NULL,NULL,NULL,&max_its);
79: SNESGetTolerances(snes,NULL,NULL,&stol,NULL,NULL);
80: SNESGetObjective(snes,&objective,NULL);
81: bt = (SNESLineSearch_BT*)linesearch->data;
82: alpha = bt->alpha;
84: SNESGetJacobian(snes, &jac, NULL, NULL, NULL);
85: if (!jac && !objective) SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix");
87: SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);
88: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);
90: VecNormBegin(Y, NORM_2, &ynorm);
91: VecNormBegin(X, NORM_2, &xnorm);
92: VecNormEnd(Y, NORM_2, &ynorm);
93: VecNormEnd(X, NORM_2, &xnorm);
95: if (ynorm == 0.0) {
96: if (monitor) {
97: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
98: PetscViewerASCIIPrintf(monitor," Line search: Initial direction and size is 0\n");
99: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
100: }
101: VecCopy(X,W);
102: VecCopy(F,G);
103: SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm);
104: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
105: return(0);
106: }
107: if (ynorm > maxstep) { /* Step too big, so scale back */
108: if (monitor) {
109: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
110: PetscViewerASCIIPrintf(monitor," Line search: Scaling step by %14.12e old ynorm %14.12e\n", (double)(maxstep/ynorm),(double)ynorm);
111: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
112: }
113: VecScale(Y,maxstep/(ynorm));
114: ynorm = maxstep;
115: }
117: /* if the SNES has an objective set, use that instead of the function value */
118: if (objective) {
119: SNESComputeObjective(snes,X,&f);
120: } else {
121: f = fnorm*fnorm;
122: }
124: /* compute the initial slope */
125: if (objective) {
126: /* slope comes from the function (assumed to be the gradient of the objective */
127: VecDotRealPart(Y,F,&initslope);
128: } else {
129: /* slope comes from the normal equations */
130: MatMult(jac,Y,W);
131: VecDotRealPart(F,W,&initslope);
132: if (initslope > 0.0) initslope = -initslope;
133: if (initslope == 0.0) initslope = -1.0;
134: }
136: while (PETSC_TRUE) {
137: VecWAXPY(W,-lambda,Y,X);
138: if (linesearch->ops->viproject) {
139: (*linesearch->ops->viproject)(snes, W);
140: }
141: if (snes->nfuncs >= snes->max_funcs) {
142: PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");
143: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
144: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);
145: return(0);
146: }
148: if (objective) {
149: SNESComputeObjective(snes,W,&g);
150: } else {
151: (*linesearch->ops->snesfunc)(snes,W,G);
152: if (linesearch->ops->vinorm) {
153: gnorm = fnorm;
154: (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
155: } else {
156: VecNorm(G,NORM_2,&gnorm);
157: }
158: g = PetscSqr(gnorm);
159: }
160: SNESLineSearchMonitor(linesearch);
162: if (!PetscIsInfOrNanReal(g)) break;
163: if (monitor) {
164: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
165: PetscViewerASCIIPrintf(monitor," Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n",(double)lambda);
166: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
167: }
168: if (lambda <= minlambda) {
169: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
170: return(0);
171: }
172: lambda = .5*lambda;
173: }
175: if (!objective) {
176: PetscInfo2(snes,"Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);
177: }
178: if (.5*g <= .5*f + lambda*alpha*initslope) { /* Sufficient reduction or step tolerance convergence */
179: if (monitor) {
180: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
181: if (!objective) {
182: PetscViewerASCIIPrintf(monitor," Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);
183: } else {
184: PetscViewerASCIIPrintf(monitor," Line search: Using full step: obj %14.12e obj %14.12e\n", (double)f, (double)g);
185: }
186: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
187: }
188: } else {
189: /* Since the full step didn't work and the step is tiny, quit */
190: if (stol*xnorm > ynorm) {
191: SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm);
192: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);
193: if (monitor) {
194: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
195: PetscViewerASCIIPrintf(monitor," Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n",(double)ynorm,(double)stol*xnorm);
196: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
197: }
198: return(0);
199: }
200: /* Fit points with quadratic */
201: lambdatemp = -initslope/(g - f - 2.0*lambda*initslope);
202: lambdaprev = lambda;
203: gprev = g;
204: if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
205: if (lambdatemp <= .1*lambda) lambda = .1*lambda;
206: else lambda = lambdatemp;
208: VecWAXPY(W,-lambda,Y,X);
209: if (linesearch->ops->viproject) {
210: (*linesearch->ops->viproject)(snes, W);
211: }
212: if (snes->nfuncs >= snes->max_funcs) {
213: PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);
214: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
215: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);
216: return(0);
217: }
218: if (objective) {
219: SNESComputeObjective(snes,W,&g);
220: } else {
221: (*linesearch->ops->snesfunc)(snes,W,G);
222: if (linesearch->ops->vinorm) {
223: gnorm = fnorm;
224: (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
225: } else {
226: VecNorm(G,NORM_2,&gnorm);
227: }
228: g = PetscSqr(gnorm);
229: }
230: if (PetscIsInfOrNanReal(g)) {
231: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF);
232: PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");
233: return(0);
234: }
235: if (monitor) {
236: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
237: if (!objective) {
238: PetscViewerASCIIPrintf(monitor," Line search: gnorm after quadratic fit %14.12e\n",(double)gnorm);
239: } else {
240: PetscViewerASCIIPrintf(monitor," Line search: obj after quadratic fit %14.12e\n",(double)g);
241: }
242: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
243: }
244: if (.5*g < .5*f + lambda*alpha*initslope) { /* sufficient reduction */
245: if (monitor) {
246: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
247: PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda);
248: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
249: }
250: } else {
251: /* Fit points with cubic */
252: for (count = 0; count < max_its; count++) {
253: if (lambda <= minlambda) {
254: if (monitor) {
255: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
256: PetscViewerASCIIPrintf(monitor," Line search: unable to find good step length! After %D tries \n",count);
257: if (!objective) {
258: PetscViewerASCIIPrintf(monitor," 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);
260: } else {
261: PetscViewerASCIIPrintf(monitor," Line search: obj(0)=%18.16e, obj=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
262: (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);
263: }
264: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
265: }
266: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
267: return(0);
268: }
269: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
270: t1 = .5*(g - f) - lambda*initslope;
271: t2 = .5*(gprev - f) - lambdaprev*initslope;
272: a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
273: b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
274: d = b*b - 3*a*initslope;
275: if (d < 0.0) d = 0.0;
276: if (a == 0.0) lambdatemp = -initslope/(2.0*b);
277: else lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a);
279: } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
280: lambdatemp = -initslope/(g - f - 2.0*initslope);
281: } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt");
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: VecWAXPY(W,-lambda,Y,X);
288: if (linesearch->ops->viproject) {
289: (*linesearch->ops->viproject)(snes,W);
290: }
291: if (snes->nfuncs >= snes->max_funcs) {
292: PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);
293: if (!objective) {
294: 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);
296: }
297: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);
298: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
299: return(0);
300: }
301: if (objective) {
302: SNESComputeObjective(snes,W,&g);
303: } else {
304: (*linesearch->ops->snesfunc)(snes,W,G);
305: if (linesearch->ops->vinorm) {
306: gnorm = fnorm;
307: (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
308: } else {
309: VecNorm(G,NORM_2,&gnorm);
310: }
311: g = PetscSqr(gnorm);
312: }
313: if (PetscIsInfOrNanReal(g)) {
314: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF);
315: PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");
316: return(0);
317: }
318: if (.5*g < .5*f + lambda*alpha*initslope) { /* is reduction enough? */
319: if (monitor) {
320: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
321: if (!objective) {
322: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
323: PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
324: } else {
325: PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
326: }
327: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
328: } else {
329: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
330: PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);
331: } else {
332: PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);
333: }
334: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
335: }
336: }
337: break;
338: } else if (monitor) {
339: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
340: if (!objective) {
341: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
342: PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
343: } else {
344: PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
345: }
346: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
347: } else {
348: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
349: PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);
350: } else {
351: PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);
352: }
353: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
354: }
355: }
356: }
357: }
358: }
360: /* postcheck */
361: SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
362: if (changed_y) {
363: VecWAXPY(W,-lambda,Y,X);
364: if (linesearch->ops->viproject) {
365: (*linesearch->ops->viproject)(snes, W);
366: }
367: }
368: if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
369: (*linesearch->ops->snesfunc)(snes,W,G);
370: if (linesearch->ops->vinorm) {
371: gnorm = fnorm;
372: (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
373: } else {
374: VecNorm(G,NORM_2,&gnorm);
375: }
376: VecNorm(Y,NORM_2,&ynorm);
377: if (PetscIsInfOrNanReal(gnorm)) {
378: SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_FAILED_NANORINF);
379: PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");
380: return(0);
381: }
382: }
384: /* copy the solution over */
385: VecCopy(W, X);
386: VecCopy(G, F);
387: VecNorm(X, NORM_2, &xnorm);
388: SNESLineSearchSetLambda(linesearch, lambda);
389: SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);
390: return(0);
391: }
393: PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
394: {
395: PetscErrorCode ierr;
396: PetscBool iascii;
397: SNESLineSearch_BT *bt;
400: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
401: bt = (SNESLineSearch_BT*)linesearch->data;
402: if (iascii) {
403: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
404: PetscViewerASCIIPrintf(viewer, " interpolation: cubic\n");
405: } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
406: PetscViewerASCIIPrintf(viewer, " interpolation: quadratic\n");
407: }
408: PetscViewerASCIIPrintf(viewer, " alpha=%e\n", (double)bt->alpha);
409: }
410: return(0);
411: }
414: static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
415: {
419: PetscFree(linesearch->data);
420: return(0);
421: }
424: static PetscErrorCode SNESLineSearchSetFromOptions_BT(PetscOptionItems *PetscOptionsObject,SNESLineSearch linesearch)
425: {
427: PetscErrorCode ierr;
428: SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
431: PetscOptionsHead(PetscOptionsObject,"SNESLineSearch BT options");
432: PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL);
433: PetscOptionsTail();
434: return(0);
435: }
438: /*MC
439: SNESLINESEARCHBT - Backtracking line search.
441: This line search finds the minimum of a polynomial fitting of the L2 norm of the
442: function or the objective function if it is provided with SNESSetObjective(). If this fit does not satisfy the conditions for progress, the interval shrinks
443: and the fit is reattempted at most max_it times or until lambda is below minlambda.
445: Options Database Keys:
446: + -snes_linesearch_alpha<1e-4> - slope descent parameter
447: . -snes_linesearch_damping<1.0> - initial step length
448: . -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the
449: step is scaled back to be of this length at the beginning of the line search
450: . -snes_linesearch_max_it<40> - maximum number of shrinking step
451: . -snes_linesearch_minlambda<1e-12> - minimum step length allowed
452: - -snes_linesearch_order<cubic,quadratic> - order of the approximation
454: Level: advanced
456: Notes:
457: This line search is taken from "Numerical Methods for Unconstrained
458: Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
460: .keywords: SNES, SNESLineSearch, damping
462: .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
463: M*/
464: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
465: {
467: SNESLineSearch_BT *bt;
468: PetscErrorCode ierr;
471: linesearch->ops->apply = SNESLineSearchApply_BT;
472: linesearch->ops->destroy = SNESLineSearchDestroy_BT;
473: linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
474: linesearch->ops->reset = NULL;
475: linesearch->ops->view = SNESLineSearchView_BT;
476: linesearch->ops->setup = NULL;
478: PetscNewLog(linesearch,&bt);
480: linesearch->data = (void*)bt;
481: linesearch->max_its = 40;
482: linesearch->order = SNES_LINESEARCH_ORDER_CUBIC;
483: bt->alpha = 1e-4;
484: return(0);
485: }