Actual source code: linesearchbt.c
petsc-3.6.1 2015-08-06
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: SNESLineSearchGetMonitor(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;
89: alpha = bt->alpha;
91: SNESGetJacobian(snes, &jac, NULL, NULL, NULL);
93: if (!jac && !objective) SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix");
95: /* precheck */
96: SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);
97: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);
99: VecNormBegin(Y, NORM_2, &ynorm);
100: VecNormBegin(X, NORM_2, &xnorm);
101: VecNormEnd(Y, NORM_2, &ynorm);
102: VecNormEnd(X, NORM_2, &xnorm);
104: if (ynorm == 0.0) {
105: if (monitor) {
106: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
107: PetscViewerASCIIPrintf(monitor," Line search: Initial direction and size is 0\n");
108: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
109: }
110: VecCopy(X,W);
111: VecCopy(F,G);
112: SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm);
113: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
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: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);
153: return(0);
154: }
156: if (objective) {
157: SNESComputeObjective(snes,W,&g);
158: } else {
159: (*linesearch->ops->snesfunc)(snes,W,G);
160: if (linesearch->ops->vinorm) {
161: gnorm = fnorm;
162: (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
163: } else {
164: VecNorm(G,NORM_2,&gnorm);
165: }
166: g = PetscSqr(gnorm);
167: }
169: if (PetscIsInfOrNanReal(g)) {
170: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF);
171: PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");
172: return(0);
173: }
174: if (!objective) {
175: PetscInfo2(snes,"Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);
176: }
177: if (.5*g <= .5*f + lambda*alpha*initslope) { /* Sufficient reduction or step tolerance convergence */
178: if (monitor) {
179: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
180: if (!objective) {
181: PetscViewerASCIIPrintf(monitor," Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);
182: } else {
183: PetscViewerASCIIPrintf(monitor," Line search: Using full step: obj0 %14.12e obj %14.12e\n", (double)f, (double)g);
184: }
185: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
186: }
187: } else {
188: /* Since the full step didn't work and the step is tiny, quit */
189: if (stol*xnorm > ynorm) {
190: SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm);
191: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
192: if (monitor) {
193: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
194: PetscViewerASCIIPrintf(monitor," Line search: Aborted due to ynorm < stol*xnorm (%14.12e < %14.12e) and inadequate full step.\n",(double)ynorm,(double)stol*xnorm);
195: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
196: }
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: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);
215: return(0);
216: }
217: if (objective) {
218: SNESComputeObjective(snes,W,&g);
219: } else {
220: (*linesearch->ops->snesfunc)(snes,W,G);
221: if (linesearch->ops->vinorm) {
222: gnorm = fnorm;
223: (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
224: } else {
225: VecNorm(G,NORM_2,&gnorm);
226: }
227: g = PetscSqr(gnorm);
228: }
229: if (PetscIsInfOrNanReal(g)) {
230: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF);
231: PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");
232: return(0);
233: }
234: if (monitor) {
235: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
236: if (!objective) {
237: PetscViewerASCIIPrintf(monitor," Line search: gnorm after quadratic fit %14.12e\n",(double)gnorm);
238: } else {
239: PetscViewerASCIIPrintf(monitor," Line search: obj after quadratic fit %14.12e\n",(double)g);
240: }
241: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
242: }
243: if (.5*g < .5*f + lambda*alpha*initslope) { /* sufficient reduction */
244: if (monitor) {
245: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
246: PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda);
247: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
248: }
249: } else {
250: /* Fit points with cubic */
251: for (count = 0; count < max_its; count++) {
252: if (lambda <= minlambda) {
253: if (monitor) {
254: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
255: PetscViewerASCIIPrintf(monitor," Line search: unable to find good step length! After %D tries \n",count);
256: if (!objective) {
257: PetscViewerASCIIPrintf(monitor,
258: " Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
259: (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);
260: } else {
261: PetscViewerASCIIPrintf(monitor,
262: " Line search: obj(0)=%18.16e, obj=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
263: (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);
264: }
265: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
266: }
267: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
268: return(0);
269: }
270: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
271: t1 = .5*(g - f) - lambda*initslope;
272: t2 = .5*(gprev - f) - lambdaprev*initslope;
273: a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
274: b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
275: d = b*b - 3*a*initslope;
276: if (d < 0.0) d = 0.0;
277: if (a == 0.0) lambdatemp = -initslope/(2.0*b);
278: else lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a);
280: } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
281: lambdatemp = -initslope/(g - f - 2.0*initslope);
282: } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt");
283: lambdaprev = lambda;
284: gprev = g;
285: if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
286: if (lambdatemp <= .1*lambda) lambda = .1*lambda;
287: else lambda = lambdatemp;
288: VecWAXPY(W,-lambda,Y,X);
289: if (linesearch->ops->viproject) {
290: (*linesearch->ops->viproject)(snes,W);
291: }
292: if (snes->nfuncs >= snes->max_funcs) {
293: PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);
294: if (!objective) {
295: PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
296: (double)fnorm,(double)gnorm,(double)ynorm,(double)lambda,(double)initslope);
297: }
298: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);
299: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
300: return(0);
301: }
302: if (objective) {
303: SNESComputeObjective(snes,W,&g);
304: } else {
305: (*linesearch->ops->snesfunc)(snes,W,G);
306: if (linesearch->ops->vinorm) {
307: gnorm = fnorm;
308: (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
309: } else {
310: VecNorm(G,NORM_2,&gnorm);
311: }
312: g = PetscSqr(gnorm);
313: }
314: if (PetscIsInfOrNanReal(gnorm)) {
315: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF);
316: PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");
317: return(0);
318: }
319: if (.5*g < .5*f + lambda*alpha*initslope) { /* is reduction enough? */
320: if (monitor) {
321: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
322: if (!objective) {
323: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
324: PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
325: } else {
326: PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
327: }
328: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
329: } else {
330: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
331: PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);
332: } else {
333: PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);
334: }
335: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
336: }
337: }
338: break;
339: } else if (monitor) {
340: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
341: if (!objective) {
342: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
343: PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
344: } else {
345: PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
346: }
347: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
348: } else {
349: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
350: PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);
351: } else {
352: PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);
353: }
354: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
355: }
356: }
357: }
358: }
359: }
361: /* postcheck */
362: SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
363: if (changed_y) {
364: VecWAXPY(W,-lambda,Y,X);
365: if (linesearch->ops->viproject) {
366: (*linesearch->ops->viproject)(snes, W);
367: }
368: }
369: if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
370: (*linesearch->ops->snesfunc)(snes,W,G);
371: if (linesearch->ops->vinorm) {
372: gnorm = fnorm;
373: (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
374: } else {
375: VecNorm(G,NORM_2,&gnorm);
376: }
377: VecNorm(Y,NORM_2,&ynorm);
378: if (PetscIsInfOrNanReal(gnorm)) {
379: SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_FAILED_NANORINF);
380: PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");
381: return(0);
382: }
383: }
385: /* copy the solution over */
386: VecCopy(W, X);
387: VecCopy(G, F);
388: VecNorm(X, NORM_2, &xnorm);
389: SNESLineSearchSetLambda(linesearch, lambda);
390: SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);
391: return(0);
392: }
396: PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
397: {
398: PetscErrorCode ierr;
399: PetscBool iascii;
400: SNESLineSearch_BT *bt;
403: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
404: bt = (SNESLineSearch_BT*)linesearch->data;
405: if (iascii) {
406: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
407: PetscViewerASCIIPrintf(viewer, " interpolation: cubic\n");
408: } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
409: PetscViewerASCIIPrintf(viewer, " interpolation: quadratic\n");
410: }
411: PetscViewerASCIIPrintf(viewer, " alpha=%e\n", (double)bt->alpha);
412: }
413: return(0);
414: }
419: static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
420: {
424: PetscFree(linesearch->data);
425: return(0);
426: }
431: static PetscErrorCode SNESLineSearchSetFromOptions_BT(PetscOptions *PetscOptionsObject,SNESLineSearch linesearch)
432: {
434: PetscErrorCode ierr;
435: SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;
438: PetscOptionsHead(PetscOptionsObject,"SNESLineSearch BT options");
439: PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL);
440: PetscOptionsTail();
441: return(0);
442: }
447: /*MC
448: SNESLINESEARCHBT - Backtracking line search.
450: This line search finds the minimum of a polynomial fitting of the L2 norm of the
451: function or the objective function if it is provided with SNESSetObjective(). If this fit does not satisfy the conditions for progress, the interval shrinks
452: and the fit is reattempted at most max_it times or until lambda is below minlambda.
454: Options Database Keys:
455: + -snes_linesearch_alpha<1e-4> - slope descent parameter
456: . -snes_linesearch_damping<1.0> - initial step length
457: . -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the
458: step is scaled back to be of this length at the beginning of the line search
459: . -snes_linesearch_max_it<40> - maximum number of shrinking step
460: . -snes_linesearch_minlambda<1e-12> - minimum step length allowed
461: - -snes_linesearch_order<cubic,quadratic> - order of the approximation
463: Level: advanced
465: Notes:
466: This line search is taken from "Numerical Methods for Unconstrained
467: Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
469: .keywords: SNES, SNESLineSearch, damping
471: .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
472: M*/
473: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
474: {
476: SNESLineSearch_BT *bt;
477: PetscErrorCode ierr;
480: linesearch->ops->apply = SNESLineSearchApply_BT;
481: linesearch->ops->destroy = SNESLineSearchDestroy_BT;
482: linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
483: linesearch->ops->reset = NULL;
484: linesearch->ops->view = SNESLineSearchView_BT;
485: linesearch->ops->setup = NULL;
487: PetscNewLog(linesearch,&bt);
489: linesearch->data = (void*)bt;
490: linesearch->max_its = 40;
491: linesearch->order = SNES_LINESEARCH_ORDER_CUBIC;
492: bt->alpha = 1e-4;
493: return(0);
494: }