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