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