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