Actual source code: linesearchbt.c

petsc-3.4.4 2014-03-13
  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(monitor,"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:       SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
197:       PetscInfo2(monitor,"Aborted due to ynorm < stol*xnorm (%14.12e < %14.12e) and inadequate full step.\n",ynorm,stol*xnorm);
198:       return(0);
199:     }
200:     /* Fit points with quadratic */
201:     lambdatemp = -initslope/(g - f - 2.0*lambda*initslope);
202:     lambdaprev = lambda;
203:     gprev      = g;
204:     if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
205:     if (lambdatemp <= .1*lambda) lambda = .1*lambda;
206:     else                         lambda = lambdatemp;

208:     VecWAXPY(W,-lambda,Y,X);
209:     if (linesearch->ops->viproject) {
210:       (*linesearch->ops->viproject)(snes, W);
211:     }
212:     if (snes->nfuncs >= snes->max_funcs) {
213:       PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);
214:       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
215:       SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
216:       return(0);
217:     }
218:     if (objective) {
219:       SNESComputeObjective(snes,W,&g);
220:     } else {
221:       SNESComputeFunction(snes,W,G);
222:       SNESGetFunctionDomainError(snes, &domainerror);
223:       if (domainerror) return(0);

225:       if (linesearch->ops->vinorm) {
226:         gnorm = fnorm;
227:         (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
228:       } else {
229:         VecNorm(G,NORM_2,&gnorm);
230:       }
231:       g = PetscSqr(gnorm);
232:     }
233:     if (PetscIsInfOrNanReal(g)) {
234:       SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
235:       PetscInfo(monitor,"Aborted due to Nan or Inf in function evaluation\n");
236:       return(0);
237:     }
238:     if (monitor) {
239:       PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
240:       if (!objective) {
241:         PetscViewerASCIIPrintf(monitor,"    Line search: gnorm after quadratic fit %14.12e\n",(double)gnorm);
242:       } else {
243:         PetscViewerASCIIPrintf(monitor,"    Line search: obj after quadratic fit %14.12e\n",(double)g);
244:       }
245:       PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
246:     }
247:     if (.5*g < .5*f + lambda*alpha*initslope) { /* sufficient reduction */
248:       if (monitor) {
249:         PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
250:         PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda);
251:         PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
252:       }
253:     } else {
254:       /* Fit points with cubic */
255:       for (count = 0; count < max_its; count++) {
256:         if (lambda <= minlambda) {
257:           if (monitor) {
258:             PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
259:             PetscViewerASCIIPrintf(monitor,"    Line search: unable to find good step length! After %D tries \n",count);
260:             if (!objective) {
261:               PetscViewerASCIIPrintf(monitor,
262:                                             "    Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
263:                                             (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);
264:             } else {
265:               PetscViewerASCIIPrintf(monitor,
266:                                             "    Line search: obj(0)=%18.16e, obj=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
267:                                             (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);
268:             }
269:             PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
270:           }
271:           SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
272:           return(0);
273:         }
274:         if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
275:           t1 = .5*(g - f) - lambda*initslope;
276:           t2 = .5*(gprev  - f) - lambdaprev*initslope;
277:           a  = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
278:           b  = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
279:           d  = b*b - 3*a*initslope;
280:           if (d < 0.0) d = 0.0;
281:           if (a == 0.0) lambdatemp = -initslope/(2.0*b);
282:           else lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a);

284:         } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
285:           lambdatemp = -initslope/(g - f - 2.0*initslope);
286:         } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt");
287:         lambdaprev = lambda;
288:         gprev      = g;
289:         if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
290:         if (lambdatemp <= .1*lambda) lambda     = .1*lambda;
291:         else                         lambda     = lambdatemp;
292:         VecWAXPY(W,-lambda,Y,X);
293:         if (linesearch->ops->viproject) {
294:           (*linesearch->ops->viproject)(snes,W);
295:         }
296:         if (snes->nfuncs >= snes->max_funcs) {
297:           PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);
298:           if (!objective) {
299:             PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
300:                               (double)fnorm,(double)gnorm,(double)ynorm,(double)lambda,(double)initslope);
301:           }
302:           SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
303:           snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
304:           return(0);
305:         }
306:         if (objective) {
307:           SNESComputeObjective(snes,W,&g);
308:         } else {
309:           SNESComputeFunction(snes,W,G);
310:           SNESGetFunctionDomainError(snes, &domainerror);
311:           if (domainerror) {
312:             SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
313:             return(0);
314:           }
315:           if (linesearch->ops->vinorm) {
316:             gnorm = fnorm;
317:             (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
318:           } else {
319:             VecNorm(G,NORM_2,&gnorm);
320:           }
321:           g = PetscSqr(gnorm);
322:         }
323:         if (PetscIsInfOrNanReal(gnorm)) {
324:           SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
325:           PetscInfo(monitor,"Aborted due to Nan or Inf in function evaluation\n");
326:           return(0);
327:         }
328:         if (.5*g < .5*f + lambda*alpha*initslope) { /* is reduction enough? */
329:           if (monitor) {
330:             PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
331:             if (!objective) {
332:               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
333:                 PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
334:               } else {
335:                 PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
336:               }
337:               PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
338:             } else {
339:               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
340:                 PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);
341:               } else {
342:                 PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);
343:               }
344:               PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
345:             }
346:           }
347:           break;
348:         } else if (monitor) {
349:           PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
350:           if (!objective) {
351:             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
352:               PetscViewerASCIIPrintf(monitor,"    Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
353:             } else {
354:               PetscViewerASCIIPrintf(monitor,"    Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
355:             }
356:             PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
357:           } else {
358:             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
359:               PetscViewerASCIIPrintf(monitor,"    Line search: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);
360:             } else {
361:               PetscViewerASCIIPrintf(monitor,"    Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);
362:             }
363:             PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
364:           }
365:         }
366:       }
367:     }
368:   }

370:   /* postcheck */
371:   SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
372:   if (changed_y) {
373:     VecWAXPY(W,-lambda,Y,X);
374:     if (linesearch->ops->viproject) {
375:       (*linesearch->ops->viproject)(snes, W);
376:     }
377:   }
378:   if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
379:     SNESComputeFunction(snes,W,G);
380:     SNESGetFunctionDomainError(snes, &domainerror);
381:     if (domainerror) {
382:       SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
383:       return(0);
384:     }
385:     if (linesearch->ops->vinorm) {
386:       gnorm = fnorm;
387:       (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
388:     } else {
389:       VecNorm(G,NORM_2,&gnorm);
390:     }
391:     VecNorm(Y,NORM_2,&ynorm);
392:     if (PetscIsInfOrNanReal(gnorm)) {
393:       SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
394:       PetscInfo(monitor,"Aborted due to Nan or Inf in function evaluation\n");
395:       return(0);
396:     }
397:   }

399:   /* copy the solution over */
400:   VecCopy(W, X);
401:   VecCopy(G, F);
402:   VecNorm(X, NORM_2, &xnorm);
403:   SNESLineSearchSetLambda(linesearch, lambda);
404:   SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);
405:   return(0);
406: }

410: PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
411: {
412:   PetscErrorCode    ierr;
413:   PetscBool         iascii;
414:   SNESLineSearch_BT *bt;

417:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
418:   bt   = (SNESLineSearch_BT*)linesearch->data;
419:   if (iascii) {
420:     if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
421:       PetscViewerASCIIPrintf(viewer, "  interpolation: cubic\n");
422:     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
423:       PetscViewerASCIIPrintf(viewer, "  interpolation: quadratic\n");
424:     }
425:     PetscViewerASCIIPrintf(viewer, "  alpha=%e\n", bt->alpha);
426:   }
427:   return(0);
428: }


433: static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
434: {

438:   PetscFree(linesearch->data);
439:   return(0);
440: }


445: static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch)
446: {

448:   PetscErrorCode    ierr;
449:   SNESLineSearch_BT *bt;

452:   bt = (SNESLineSearch_BT*)linesearch->data;

454:   PetscOptionsHead("SNESLineSearch BT options");
455:   PetscOptionsReal("-snes_linesearch_alpha",   "Descent tolerance",        "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL);

457:   PetscOptionsTail();
458:   return(0);
459: }


464: /*MC
465:    SNESLINESEARCHBT - Backtracking line search.

467:    This line search finds the minimum of a polynomial fitting of the L2 norm of the
468:    function. If this fit does not satisfy the conditions for progress, the interval shrinks
469:    and the fit is reattempted at most max_it times or until lambda is below minlambda.

471:    Options Database Keys:
472: +  -snes_linesearch_alpha<1e-4> - slope descent parameter
473: .  -snes_linesearch_damping<1.0> - initial step length
474: .  -snes_linesearch_max_it<40> - maximum number of shrinking step
475: .  -snes_linesearch_minlambda<1e-12> - minimum step length allowed
476: -  -snes_linesearch_order<cubic,quadratic> - order of the approximation

478:    Level: advanced

480:    Notes:
481:    This line search is taken from "Numerical Methods for Unconstrained
482:    Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.

484: .keywords: SNES, SNESLineSearch, damping

486: .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
487: M*/
488: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
489: {

491:   SNESLineSearch_BT *bt;
492:   PetscErrorCode    ierr;

495:   linesearch->ops->apply          = SNESLineSearchApply_BT;
496:   linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
497:   linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
498:   linesearch->ops->reset          = NULL;
499:   linesearch->ops->view           = SNESLineSearchView_BT;
500:   linesearch->ops->setup          = NULL;

502:   PetscNewLog(linesearch, SNESLineSearch_BT, &bt);

504:   linesearch->data    = (void*)bt;
505:   linesearch->max_its = 40;
506:   linesearch->order   = SNES_LINESEARCH_ORDER_CUBIC;
507:   bt->alpha           = 1e-4;
508:   return(0);
509: }