Actual source code: linesearchbt.c
petsc-3.12.5 2020-03-29
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: SNESCheckFunctionNorm(snes,g);
166: }
167: lambda = .5*lambda;
168: }
170: if (!objective) {
171: PetscInfo2(snes,"Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);
172: }
173: if (.5*g <= .5*f + lambda*alpha*initslope) { /* Sufficient reduction or step tolerance convergence */
174: if (monitor) {
175: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
176: if (!objective) {
177: PetscViewerASCIIPrintf(monitor," Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);
178: } else {
179: PetscViewerASCIIPrintf(monitor," Line search: Using full step: obj %14.12e obj %14.12e\n", (double)f, (double)g);
180: }
181: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
182: }
183: } else {
184: /* Since the full step didn't work and the step is tiny, quit */
185: if (stol*xnorm > ynorm) {
186: SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm);
187: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);
188: if (monitor) {
189: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
190: PetscViewerASCIIPrintf(monitor," Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n",(double)ynorm,(double)stol*xnorm);
191: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
192: }
193: return(0);
194: }
195: /* Fit points with quadratic */
196: lambdatemp = -initslope/(g - f - 2.0*lambda*initslope);
197: lambdaprev = lambda;
198: gprev = g;
199: if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
200: if (lambdatemp <= .1*lambda) lambda = .1*lambda;
201: else lambda = lambdatemp;
203: VecWAXPY(W,-lambda,Y,X);
204: if (linesearch->ops->viproject) {
205: (*linesearch->ops->viproject)(snes, W);
206: }
207: if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
208: PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);
209: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
210: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);
211: return(0);
212: }
213: if (objective) {
214: SNESComputeObjective(snes,W,&g);
215: } else {
216: (*linesearch->ops->snesfunc)(snes,W,G);
217: if (linesearch->ops->vinorm) {
218: gnorm = fnorm;
219: (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
220: } else {
221: VecNorm(G,NORM_2,&gnorm);
222: }
223: g = PetscSqr(gnorm);
224: }
225: if (PetscIsInfOrNanReal(g)) {
226: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF);
227: PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");
228: return(0);
229: }
230: if (monitor) {
231: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
232: if (!objective) {
233: PetscViewerASCIIPrintf(monitor," Line search: gnorm after quadratic fit %14.12e\n",(double)gnorm);
234: } else {
235: PetscViewerASCIIPrintf(monitor," Line search: obj after quadratic fit %14.12e\n",(double)g);
236: }
237: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
238: }
239: if (.5*g < .5*f + lambda*alpha*initslope) { /* sufficient reduction */
240: if (monitor) {
241: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
242: PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda);
243: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
244: }
245: } else {
246: /* Fit points with cubic */
247: for (count = 0; count < max_its; count++) {
248: if (lambda <= minlambda) {
249: if (monitor) {
250: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
251: PetscViewerASCIIPrintf(monitor," Line search: unable to find good step length! After %D tries \n",count);
252: if (!objective) {
253: PetscViewerASCIIPrintf(monitor," Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
254: (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);
255: } else {
256: 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",
257: (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);
258: }
259: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
260: }
261: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
262: return(0);
263: }
264: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
265: t1 = .5*(g - f) - lambda*initslope;
266: t2 = .5*(gprev - f) - lambdaprev*initslope;
267: a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
268: b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
269: d = b*b - 3*a*initslope;
270: if (d < 0.0) d = 0.0;
271: if (a == 0.0) lambdatemp = -initslope/(2.0*b);
272: else lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a);
274: } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
275: lambdatemp = -initslope/(g - f - 2.0*initslope);
276: } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt");
277: lambdaprev = lambda;
278: gprev = g;
279: if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
280: if (lambdatemp <= .1*lambda) lambda = .1*lambda;
281: else lambda = lambdatemp;
282: VecWAXPY(W,-lambda,Y,X);
283: if (linesearch->ops->viproject) {
284: (*linesearch->ops->viproject)(snes,W);
285: }
286: if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
287: PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);
288: if (!objective) {
289: PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
290: (double)fnorm,(double)gnorm,(double)ynorm,(double)lambda,(double)initslope);
291: }
292: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);
293: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
294: return(0);
295: }
296: if (objective) {
297: SNESComputeObjective(snes,W,&g);
298: } else {
299: (*linesearch->ops->snesfunc)(snes,W,G);
300: if (linesearch->ops->vinorm) {
301: gnorm = fnorm;
302: (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
303: } else {
304: VecNorm(G,NORM_2,&gnorm);
305: }
306: g = PetscSqr(gnorm);
307: }
308: if (PetscIsInfOrNanReal(g)) {
309: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF);
310: PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");
311: return(0);
312: }
313: if (.5*g < .5*f + lambda*alpha*initslope) { /* is reduction enough? */
314: if (monitor) {
315: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
316: if (!objective) {
317: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
318: PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
319: } else {
320: PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
321: }
322: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
323: } else {
324: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
325: PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);
326: } else {
327: PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);
328: }
329: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
330: }
331: }
332: break;
333: } else if (monitor) {
334: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
335: if (!objective) {
336: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
337: PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
338: } else {
339: PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
340: }
341: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
342: } else {
343: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
344: PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);
345: } else {
346: PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);
347: }
348: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
349: }
350: }
351: }
352: }
353: }
355: /* postcheck */
356: /* update Y to lambda*Y so that W is consistent with X - lambda*Y */
357: VecScale(Y,lambda);
358: SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
359: if (changed_y) {
360: VecWAXPY(W,-1.0,Y,X);
361: if (linesearch->ops->viproject) {
362: (*linesearch->ops->viproject)(snes, W);
363: }
364: }
365: if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
366: (*linesearch->ops->snesfunc)(snes,W,G);
367: if (linesearch->ops->vinorm) {
368: gnorm = fnorm;
369: (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
370: } else {
371: VecNorm(G,NORM_2,&gnorm);
372: }
373: VecNorm(Y,NORM_2,&ynorm);
374: if (PetscIsInfOrNanReal(gnorm)) {
375: SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_FAILED_NANORINF);
376: PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");
377: return(0);
378: }
379: }
381: /* copy the solution over */
382: VecCopy(W, X);
383: VecCopy(G, F);
384: VecNorm(X, NORM_2, &xnorm);
385: SNESLineSearchSetLambda(linesearch, lambda);
386: SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);
387: return(0);
388: }
390: PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
391: {
392: PetscErrorCode ierr;
393: PetscBool iascii;
394: SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
397: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
398: if (iascii) {
399: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
400: PetscViewerASCIIPrintf(viewer, " interpolation: cubic\n");
401: } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
402: PetscViewerASCIIPrintf(viewer, " interpolation: quadratic\n");
403: }
404: PetscViewerASCIIPrintf(viewer, " alpha=%e\n", (double)bt->alpha);
405: }
406: return(0);
407: }
409: static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
410: {
414: PetscFree(linesearch->data);
415: return(0);
416: }
418: static PetscErrorCode SNESLineSearchSetFromOptions_BT(PetscOptionItems *PetscOptionsObject,SNESLineSearch linesearch)
419: {
420: PetscErrorCode ierr;
421: SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
424: PetscOptionsHead(PetscOptionsObject,"SNESLineSearch BT options");
425: PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL);
426: PetscOptionsTail();
427: return(0);
428: }
430: /*MC
431: SNESLINESEARCHBT - Backtracking line search.
433: This line search finds the minimum of a polynomial fitting of the L2 norm of the
434: function or the objective function if it is provided with SNESSetObjective(). If this fit does not satisfy the conditions for progress, the interval shrinks
435: and the fit is reattempted at most max_it times or until lambda is below minlambda.
437: Options Database Keys:
438: + -snes_linesearch_alpha<1e-4> - slope descent parameter
439: . -snes_linesearch_damping<1.0> - initial step length
440: . -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the
441: step is scaled back to be of this length at the beginning of the line search
442: . -snes_linesearch_max_it<40> - maximum number of shrinking step
443: . -snes_linesearch_minlambda<1e-12> - minimum step length allowed
444: - -snes_linesearch_order<cubic,quadratic> - order of the approximation
446: Level: advanced
448: Notes:
449: This line search is taken from "Numerical Methods for Unconstrained
450: Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
451:
452: This line search will always produce a step that is less than or equal to, in length, the full step size.
454: .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
455: M*/
456: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
457: {
459: SNESLineSearch_BT *bt;
460: PetscErrorCode ierr;
463: linesearch->ops->apply = SNESLineSearchApply_BT;
464: linesearch->ops->destroy = SNESLineSearchDestroy_BT;
465: linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
466: linesearch->ops->reset = NULL;
467: linesearch->ops->view = SNESLineSearchView_BT;
468: linesearch->ops->setup = NULL;
470: PetscNewLog(linesearch,&bt);
472: linesearch->data = (void*)bt;
473: linesearch->max_its = 40;
474: linesearch->order = SNES_LINESEARCH_ORDER_CUBIC;
475: bt->alpha = 1e-4;
476: return(0);
477: }