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