Actual source code: linesearchbt.c
petsc-3.4.5 2014-06-29
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: PetscBool domainerror;
73: PetscViewer monitor;
74: PetscInt max_its,count;
75: SNESLineSearch_BT *bt;
76: Mat jac;
77: PetscErrorCode (*objective)(SNES,Vec,PetscReal*,void*);
80: SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G);
81: SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);
82: SNESLineSearchGetLambda(linesearch, &lambda);
83: SNESLineSearchGetSNES(linesearch, &snes);
84: SNESLineSearchGetMonitor(linesearch, &monitor);
85: SNESLineSearchGetTolerances(linesearch,&minlambda,&maxstep,NULL,NULL,NULL,&max_its);
86: SNESGetTolerances(snes,NULL,NULL,&stol,NULL,NULL);
87: SNESGetObjective(snes,&objective,NULL);
88: bt = (SNESLineSearch_BT*)linesearch->data;
90: alpha = bt->alpha;
92: SNESGetJacobian(snes, &jac, NULL, NULL, NULL);
94: if (!jac && !objective) SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix");
96: /* precheck */
97: SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);
98: SNESLineSearchSetSuccess(linesearch, PETSC_TRUE);
100: VecNormBegin(Y, NORM_2, &ynorm);
101: VecNormBegin(X, NORM_2, &xnorm);
102: VecNormEnd(Y, NORM_2, &ynorm);
103: VecNormEnd(X, NORM_2, &xnorm);
105: if (ynorm == 0.0) {
106: if (monitor) {
107: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
108: PetscViewerASCIIPrintf(monitor," Line search: Initial direction and size is 0\n");
109: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
110: }
111: VecCopy(X,W);
112: VecCopy(F,G);
113: SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm);
114: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
115: return(0);
116: }
117: if (ynorm > maxstep) { /* Step too big, so scale back */
118: if (monitor) {
119: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
120: PetscViewerASCIIPrintf(monitor," Line search: Scaling step by %14.12e old ynorm %14.12e\n", (double)(maxstep/ynorm),(double)ynorm);
121: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
122: }
123: VecScale(Y,maxstep/(ynorm));
124: ynorm = maxstep;
125: }
127: /* if the SNES has an objective set, use that instead of the function value */
128: if (objective) {
129: SNESComputeObjective(snes,X,&f);
130: } else {
131: f = fnorm*fnorm;
132: }
134: /* compute the initial slope */
135: if (objective) {
136: /* slope comes from the function (assumed to be the gradient of the objective */
137: VecDotRealPart(Y,F,&initslope);
138: } else {
139: /* slope comes from the normal equations */
140: MatMult(jac,Y,W);
141: VecDotRealPart(F,W,&initslope);
142: if (initslope > 0.0) initslope = -initslope;
143: if (initslope == 0.0) initslope = -1.0;
144: }
146: VecWAXPY(W,-lambda,Y,X);
147: if (linesearch->ops->viproject) {
148: (*linesearch->ops->viproject)(snes, W);
149: }
150: if (snes->nfuncs >= snes->max_funcs) {
151: PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");
152: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
153: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
154: return(0);
155: }
157: if (objective) {
158: SNESComputeObjective(snes,W,&g);
159: } else {
160: SNESComputeFunction(snes,W,G);
161: SNESGetFunctionDomainError(snes, &domainerror);
162: if (domainerror) {
163: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
164: return(0);
165: }
166: if (linesearch->ops->vinorm) {
167: gnorm = fnorm;
168: (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
169: } else {
170: VecNorm(G,NORM_2,&gnorm);
171: }
172: g = PetscSqr(gnorm);
173: }
175: if (PetscIsInfOrNanReal(g)) {
176: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
177: PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");
178: return(0);
179: }
180: if (!objective) {
181: PetscInfo2(snes,"Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);
182: }
183: if (.5*g <= .5*f + lambda*alpha*initslope) { /* Sufficient reduction or step tolerance convergence */
184: if (monitor) {
185: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
186: if (!objective) {
187: PetscViewerASCIIPrintf(monitor," Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);
188: } else {
189: PetscViewerASCIIPrintf(monitor," Line search: Using full step: obj0 %14.12e obj %14.12e\n", (double)f, (double)g);
190: }
191: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
192: }
193: } else {
194: /* Since the full step didn't work and the step is tiny, quit */
195: if (stol*xnorm > ynorm) {
196: SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm);
197: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
198: if (monitor) {
199: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
200: PetscViewerASCIIPrintf(monitor," Line search: Aborted due to ynorm < stol*xnorm (%14.12e < %14.12e) and inadequate full step.\n",(double)ynorm,(double)stol*xnorm);
201: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
202: }
203: return(0);
204: }
205: /* Fit points with quadratic */
206: lambdatemp = -initslope/(g - f - 2.0*lambda*initslope);
207: lambdaprev = lambda;
208: gprev = g;
209: if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
210: if (lambdatemp <= .1*lambda) lambda = .1*lambda;
211: else lambda = lambdatemp;
213: VecWAXPY(W,-lambda,Y,X);
214: if (linesearch->ops->viproject) {
215: (*linesearch->ops->viproject)(snes, W);
216: }
217: if (snes->nfuncs >= snes->max_funcs) {
218: PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);
219: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
220: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
221: return(0);
222: }
223: if (objective) {
224: SNESComputeObjective(snes,W,&g);
225: } else {
226: SNESComputeFunction(snes,W,G);
227: SNESGetFunctionDomainError(snes, &domainerror);
228: if (domainerror) return(0);
230: if (linesearch->ops->vinorm) {
231: gnorm = fnorm;
232: (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
233: } else {
234: VecNorm(G,NORM_2,&gnorm);
235: }
236: g = PetscSqr(gnorm);
237: }
238: if (PetscIsInfOrNanReal(g)) {
239: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
240: PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");
241: return(0);
242: }
243: if (monitor) {
244: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
245: if (!objective) {
246: PetscViewerASCIIPrintf(monitor," Line search: gnorm after quadratic fit %14.12e\n",(double)gnorm);
247: } else {
248: PetscViewerASCIIPrintf(monitor," Line search: obj after quadratic fit %14.12e\n",(double)g);
249: }
250: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
251: }
252: if (.5*g < .5*f + lambda*alpha*initslope) { /* sufficient reduction */
253: if (monitor) {
254: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
255: PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda);
256: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
257: }
258: } else {
259: /* Fit points with cubic */
260: for (count = 0; count < max_its; count++) {
261: if (lambda <= minlambda) {
262: if (monitor) {
263: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
264: PetscViewerASCIIPrintf(monitor," Line search: unable to find good step length! After %D tries \n",count);
265: if (!objective) {
266: PetscViewerASCIIPrintf(monitor,
267: " Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
268: (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);
269: } else {
270: PetscViewerASCIIPrintf(monitor,
271: " Line search: obj(0)=%18.16e, obj=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
272: (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);
273: }
274: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
275: }
276: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
277: return(0);
278: }
279: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
280: t1 = .5*(g - f) - lambda*initslope;
281: t2 = .5*(gprev - f) - lambdaprev*initslope;
282: a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
283: b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
284: d = b*b - 3*a*initslope;
285: if (d < 0.0) d = 0.0;
286: if (a == 0.0) lambdatemp = -initslope/(2.0*b);
287: else lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a);
289: } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
290: lambdatemp = -initslope/(g - f - 2.0*initslope);
291: } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt");
292: lambdaprev = lambda;
293: gprev = g;
294: if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
295: if (lambdatemp <= .1*lambda) lambda = .1*lambda;
296: else lambda = lambdatemp;
297: VecWAXPY(W,-lambda,Y,X);
298: if (linesearch->ops->viproject) {
299: (*linesearch->ops->viproject)(snes,W);
300: }
301: if (snes->nfuncs >= snes->max_funcs) {
302: PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);
303: if (!objective) {
304: PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
305: (double)fnorm,(double)gnorm,(double)ynorm,(double)lambda,(double)initslope);
306: }
307: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
308: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
309: return(0);
310: }
311: if (objective) {
312: SNESComputeObjective(snes,W,&g);
313: } else {
314: SNESComputeFunction(snes,W,G);
315: SNESGetFunctionDomainError(snes, &domainerror);
316: if (domainerror) {
317: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
318: return(0);
319: }
320: if (linesearch->ops->vinorm) {
321: gnorm = fnorm;
322: (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
323: } else {
324: VecNorm(G,NORM_2,&gnorm);
325: }
326: g = PetscSqr(gnorm);
327: }
328: if (PetscIsInfOrNanReal(gnorm)) {
329: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
330: PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");
331: return(0);
332: }
333: if (.5*g < .5*f + lambda*alpha*initslope) { /* is reduction enough? */
334: if (monitor) {
335: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
336: if (!objective) {
337: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
338: PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
339: } else {
340: PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, current gnorm %14.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: Cubically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);
346: } else {
347: PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);
348: }
349: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
350: }
351: }
352: break;
353: } else if (monitor) {
354: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
355: if (!objective) {
356: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
357: PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
358: } else {
359: PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
360: }
361: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
362: } else {
363: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
364: PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);
365: } else {
366: PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);
367: }
368: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
369: }
370: }
371: }
372: }
373: }
375: /* postcheck */
376: SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
377: if (changed_y) {
378: VecWAXPY(W,-lambda,Y,X);
379: if (linesearch->ops->viproject) {
380: (*linesearch->ops->viproject)(snes, W);
381: }
382: }
383: if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
384: SNESComputeFunction(snes,W,G);
385: SNESGetFunctionDomainError(snes, &domainerror);
386: if (domainerror) {
387: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
388: return(0);
389: }
390: if (linesearch->ops->vinorm) {
391: gnorm = fnorm;
392: (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
393: } else {
394: VecNorm(G,NORM_2,&gnorm);
395: }
396: VecNorm(Y,NORM_2,&ynorm);
397: if (PetscIsInfOrNanReal(gnorm)) {
398: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
399: PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");
400: return(0);
401: }
402: }
404: /* copy the solution over */
405: VecCopy(W, X);
406: VecCopy(G, F);
407: VecNorm(X, NORM_2, &xnorm);
408: SNESLineSearchSetLambda(linesearch, lambda);
409: SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);
410: return(0);
411: }
415: PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
416: {
417: PetscErrorCode ierr;
418: PetscBool iascii;
419: SNESLineSearch_BT *bt;
422: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
423: bt = (SNESLineSearch_BT*)linesearch->data;
424: if (iascii) {
425: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
426: PetscViewerASCIIPrintf(viewer, " interpolation: cubic\n");
427: } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
428: PetscViewerASCIIPrintf(viewer, " interpolation: quadratic\n");
429: }
430: PetscViewerASCIIPrintf(viewer, " alpha=%e\n", bt->alpha);
431: }
432: return(0);
433: }
438: static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
439: {
443: PetscFree(linesearch->data);
444: return(0);
445: }
450: static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch)
451: {
453: PetscErrorCode ierr;
454: SNESLineSearch_BT *bt;
457: bt = (SNESLineSearch_BT*)linesearch->data;
459: PetscOptionsHead("SNESLineSearch BT options");
460: PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL);
462: PetscOptionsTail();
463: return(0);
464: }
469: /*MC
470: SNESLINESEARCHBT - Backtracking line search.
472: This line search finds the minimum of a polynomial fitting of the L2 norm of the
473: function. If this fit does not satisfy the conditions for progress, the interval shrinks
474: and the fit is reattempted at most max_it times or until lambda is below minlambda.
476: Options Database Keys:
477: + -snes_linesearch_alpha<1e-4> - slope descent parameter
478: . -snes_linesearch_damping<1.0> - initial step length
479: . -snes_linesearch_max_it<40> - maximum number of shrinking step
480: . -snes_linesearch_minlambda<1e-12> - minimum step length allowed
481: - -snes_linesearch_order<cubic,quadratic> - order of the approximation
483: Level: advanced
485: Notes:
486: This line search is taken from "Numerical Methods for Unconstrained
487: Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
489: .keywords: SNES, SNESLineSearch, damping
491: .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
492: M*/
493: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
494: {
496: SNESLineSearch_BT *bt;
497: PetscErrorCode ierr;
500: linesearch->ops->apply = SNESLineSearchApply_BT;
501: linesearch->ops->destroy = SNESLineSearchDestroy_BT;
502: linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
503: linesearch->ops->reset = NULL;
504: linesearch->ops->view = SNESLineSearchView_BT;
505: linesearch->ops->setup = NULL;
507: PetscNewLog(linesearch, SNESLineSearch_BT, &bt);
509: linesearch->data = (void*)bt;
510: linesearch->max_its = 40;
511: linesearch->order = SNES_LINESEARCH_ORDER_CUBIC;
512: bt->alpha = 1e-4;
513: return(0);
514: }