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