Actual source code: linesearchbt.c
petsc-3.3-p7 2013-05-11
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;
26: bt = (SNESLineSearch_BT *)linesearch->data;
27: bt->alpha = alpha;
28: return(0);
29: }
34: /*@
35: SNESLineSearchBTGetAlpha - Gets the descent parameter, alpha, in the BT linesearch variant.
37: Input Parameters:
38: . linesearch - linesearch context
40: Output Parameters:
41: . alpha - The descent parameter
43: Level: intermediate
45: .seealso: SNESLineSearchGetLambda(), SNESLineSearchGetTolerances() SNESLINESEARCHBT
46: @*/
47: PetscErrorCode SNESLineSearchBTGetAlpha(SNESLineSearch linesearch, PetscReal *alpha)
48: {
49: SNESLineSearch_BT *bt;
52: bt = (SNESLineSearch_BT *)linesearch->data;
53: *alpha = bt->alpha;
54: return(0);
55: }
59: static PetscErrorCode SNESLineSearchApply_BT(SNESLineSearch linesearch)
60: {
61: PetscBool changed_y,changed_w;
63: Vec X,F,Y,W,G;
64: SNES snes;
65: PetscReal fnorm, xnorm, ynorm, gnorm, gnormprev;
66: PetscReal lambda,lambdatemp,lambdaprev,minlambda,maxstep,initslope,alpha,stol;
67: PetscReal t1,t2,a,b,d;
68: #if defined(PETSC_USE_COMPLEX)
69: PetscScalar cinitslope;
70: #endif
71: PetscBool domainerror;
72: PetscViewer monitor;
73: PetscInt max_its,count;
74: SNESLineSearch_BT *bt;
75: Mat jac;
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,PETSC_NULL,PETSC_NULL,PETSC_NULL,&max_its);
86: SNESGetTolerances(snes,PETSC_NULL,PETSC_NULL,&stol,PETSC_NULL,PETSC_NULL);
87: bt = (SNESLineSearch_BT *)linesearch->data;
89: alpha = bt->alpha;
91: SNESGetJacobian(snes, &jac, PETSC_NULL, PETSC_NULL, PETSC_NULL);
92: if (!jac) {
93: SETERRQ(((PetscObject)linesearch)->comm, PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix");
94: }
95: /* precheck */
96: SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);
97: SNESLineSearchSetSuccess(linesearch, PETSC_TRUE);
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: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
113: return(0);
114: }
115: if (ynorm > maxstep) { /* Step too big, so scale back */
116: if (monitor) {
117: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
118: PetscViewerASCIIPrintf(monitor," Line search: Scaling step by %14.12e old ynorm %14.12e\n", (maxstep/ynorm),ynorm);
119: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
120: }
121: VecScale(Y,maxstep/(ynorm));
122: ynorm = maxstep;
123: }
124: MatMult(jac,Y,W);
125: #if defined(PETSC_USE_COMPLEX)
126: VecDot(F,W,&cinitslope);
127: initslope = PetscRealPart(cinitslope);
128: #else
129: VecDot(F,W,&initslope);
130: #endif
131: if (initslope > 0.0) initslope = -initslope;
132: if (initslope == 0.0) initslope = -1.0;
134: VecWAXPY(W,-lambda,Y,X);
135: if (linesearch->ops->viproject) {
136: (*linesearch->ops->viproject)(snes, W);
137: }
138: if (snes->nfuncs >= snes->max_funcs) {
139: PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");
140: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
141: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
142: return(0);
143: }
144: SNESComputeFunction(snes,W,G);
145: SNESGetFunctionDomainError(snes, &domainerror);
146: if (domainerror) {
147: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
148: return(0);
149: }
150: if (linesearch->ops->vinorm) {
151: gnorm = fnorm;
152: (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
153: } else {
154: VecNorm(G,NORM_2,&gnorm);
155: }
157: if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
158: PetscInfo2(snes,"Initial fnorm %14.12e gnorm %14.12e\n", fnorm, gnorm);
159: if (.5*gnorm*gnorm <= .5*fnorm*fnorm + lambda*alpha*initslope || ynorm < stol*xnorm) { /* Sufficient reduction or step tolerance convergence */
160: if (monitor) {
161: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
162: PetscViewerASCIIPrintf(monitor," Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", fnorm, gnorm);
163: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
164: }
165: } else {
166: /* Fit points with quadratic */
167: lambdatemp = -initslope/(gnorm*gnorm - fnorm*fnorm - 2.0*lambda*initslope);
168: lambdaprev = lambda;
169: gnormprev = gnorm;
170: if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
171: if (lambdatemp <= .1*lambda) lambda = .1*lambda;
172: else lambda = lambdatemp;
174: VecWAXPY(W,-lambda,Y,X);
175: if (linesearch->ops->viproject) {
176: (*linesearch->ops->viproject)(snes, W);
177: }
178: if (snes->nfuncs >= snes->max_funcs) {
179: PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);
180: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
181: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
182: return(0);
183: }
184: SNESComputeFunction(snes,W,G);
185: SNESGetFunctionDomainError(snes, &domainerror);
186: if (domainerror) {
187: return(0);
188: }
189: if (linesearch->ops->vinorm) {
190: gnorm = fnorm;
191: (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
192: } else {
193: VecNorm(G,NORM_2,&gnorm);
194: }
195: if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
196: if (monitor) {
197: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
198: PetscViewerASCIIPrintf(monitor," Line search: gnorm after quadratic fit %14.12e\n",gnorm);
199: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
200: }
201: if (.5*gnorm*gnorm < .5*fnorm*fnorm + lambda*alpha*initslope) { /* sufficient reduction */
202: if (monitor) {
203: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
204: PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda);
205: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
206: }
207: } else {
208: /* Fit points with cubic */
209: for (count = 0; count < max_its; count++) {
210: if (lambda <= minlambda) {
211: if (monitor) {
212: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
213: PetscViewerASCIIPrintf(monitor," Line search: unable to find good step length! After %D tries \n",count);
214: PetscViewerASCIIPrintf(monitor,
215: " Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
216: fnorm, gnorm, ynorm, minlambda, lambda, initslope);
217: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
218: }
219: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
220: return(0);
221: }
222: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
223: t1 = .5*(gnorm*gnorm - fnorm*fnorm) - lambda*initslope;
224: t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope;
225: a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
226: b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
227: d = b*b - 3*a*initslope;
228: if (d < 0.0) d = 0.0;
229: if (a == 0.0) {
230: lambdatemp = -initslope/(2.0*b);
231: } else {
232: lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a);
233: }
234: } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
235: lambdatemp = -initslope/(gnorm*gnorm - fnorm*fnorm - 2.0*initslope);
236: } else {
237: SETERRQ(((PetscObject)linesearch)->comm, PETSC_ERR_SUP, "unsupported line search order for type bt");
238: }
239: lambdaprev = lambda;
240: gnormprev = gnorm;
241: if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
242: if (lambdatemp <= .1*lambda) lambda = .1*lambda;
243: else lambda = lambdatemp;
244: VecWAXPY(W,-lambda,Y,X);
245: if (linesearch->ops->viproject) {
246: (*linesearch->ops->viproject)(snes,W);
247: }
248: if (snes->nfuncs >= snes->max_funcs) {
249: PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);
250: PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
251: fnorm,gnorm,ynorm,lambda,initslope);
252: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
253: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
254: return(0);
255: }
256: SNESComputeFunction(snes,W,G);
257: SNESGetFunctionDomainError(snes, &domainerror);
258: if (domainerror) {
259: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
260: return(0);
261: }
262: if (linesearch->ops->vinorm) {
263: gnorm = fnorm;
264: (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
265: } else {
266: VecNorm(G,NORM_2,&gnorm);
267: }
268: if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
269: if (.5*gnorm*gnorm < .5*fnorm*fnorm + lambda*alpha*initslope) { /* is reduction enough? */
270: if (monitor) {
271: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
272: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
273: PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",gnorm,lambda);
274: } else {
275: PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",gnorm,lambda);
276: }
277: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
278: }
279: break;
280: } else {
281: if (monitor) {
282: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
283: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
284: PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",gnorm,lambda);
285: } else {
286: PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",gnorm,lambda);
287: }
288: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
289: }
290: }
291: }
292: }
293: }
295: /* postcheck */
296: SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
297: if (changed_y) {
298: VecWAXPY(W,-lambda,Y,X);
299: if (linesearch->ops->viproject) {
300: (*linesearch->ops->viproject)(snes, W);
301: }
302: }
303: if (changed_y || changed_w) { /* recompute the function if the step has changed */
304: SNESComputeFunction(snes,W,G);
305: SNESGetFunctionDomainError(snes, &domainerror);
306: if (domainerror) {
307: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
308: return(0);
309: }
310: if (linesearch->ops->vinorm) {
311: gnorm = fnorm;
312: (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
313: } else {
314: VecNorm(G,NORM_2,&gnorm);
315: }
316: VecNorm(Y,NORM_2,&ynorm);
317: if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
319: }
321: /* copy the solution over */
322: VecCopy(W, X);
323: VecCopy(G, F);
324: VecNorm(X, NORM_2, &xnorm);
325: SNESLineSearchSetLambda(linesearch, lambda);
326: SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);
327: return(0);
328: }
332: PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
333: {
334: PetscErrorCode ierr;
335: PetscBool iascii;
336: SNESLineSearch_BT *bt;
338: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
339: bt = (SNESLineSearch_BT*)linesearch->data;
340: if (iascii) {
341: if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
342: PetscViewerASCIIPrintf(viewer, " interpolation: cubic\n");
343: } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
344: PetscViewerASCIIPrintf(viewer, " interpolation: quadratic\n");
345: }
346: PetscViewerASCIIPrintf(viewer, " alpha=%e\n", bt->alpha);
347: }
348: return(0);
349: }
354: static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
355: {
359: PetscFree(linesearch->data);
360: return(0);
361: }
366: static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch)
367: {
369: PetscErrorCode ierr;
370: SNESLineSearch_BT *bt;
373: bt = (SNESLineSearch_BT*)linesearch->data;
375: PetscOptionsHead("SNESLineSearch BT options");
376: PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, PETSC_NULL);
378: PetscOptionsTail();
379: return(0);
380: }
385: /*MC
386: SNESLINESEARCHBT - Backtracking line search.
388: This line search finds the minimum of a polynomial fitting of the L2 norm of the
389: function. If this fit does not satisfy the conditions for progress, the interval shrinks
390: and the fit is reattempted at most max_it times or until lambda is below minlambda.
392: Options Database Keys:
393: + -snes_linesearch_alpha<1e-4> - slope descent parameter
394: . -snes_linesearch_damping<1.0> - initial step length
395: . -snes_linesearch_max_it<40> - maximum number of shrinking step
396: . -snes_linesearch_minlambda<1e-12> - minimum step length allowed
397: - -snes_linesearch_order<cubic,quadratic> - order of the approximation
399: Level: advanced
401: Notes:
402: This line search is taken from "Numerical Methods for Unconstrained
403: Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.
405: .keywords: SNES, SNESLineSearch, damping
407: .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
408: M*/
410: {
412: SNESLineSearch_BT *bt;
416: linesearch->ops->apply = SNESLineSearchApply_BT;
417: linesearch->ops->destroy = SNESLineSearchDestroy_BT;
418: linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
419: linesearch->ops->reset = PETSC_NULL;
420: linesearch->ops->view = SNESLineSearchView_BT;
421: linesearch->ops->setup = PETSC_NULL;
423: PetscNewLog(linesearch, SNESLineSearch_BT, &bt);
424: linesearch->data = (void *)bt;
425: linesearch->max_its = 40;
426: linesearch->order = SNES_LINESEARCH_ORDER_CUBIC;
427: bt->alpha = 1e-4;
429: return(0);
430: }