Actual source code: linesearchbt.c

petsc-master 2017-02-19
Report Typos and Errors
  1:  #include <petsc/private/linesearchimpl.h>
  2:  #include <petsc/private/snesimpl.h>

  4: typedef struct {
  5:   PetscReal alpha;        /* sufficient decrease parameter */
  6: } SNESLineSearch_BT;

  8: /*@
  9:    SNESLineSearchBTSetAlpha - Sets the descent parameter, alpha, in the BT linesearch variant.

 11:    Input Parameters:
 12: +  linesearch - linesearch context
 13: -  alpha - The descent parameter

 15:    Level: intermediate

 17: .seealso: SNESLineSearchSetLambda(), SNESLineSearchGetTolerances() SNESLINESEARCHBT
 18: @*/
 19: PetscErrorCode SNESLineSearchBTSetAlpha(SNESLineSearch linesearch, PetscReal alpha)
 20: {
 21:   SNESLineSearch_BT *bt;

 25:   bt        = (SNESLineSearch_BT*)linesearch->data;
 26:   bt->alpha = alpha;
 27:   return(0);
 28: }


 31: /*@
 32:    SNESLineSearchBTGetAlpha - Gets the descent parameter, alpha, in the BT linesearch variant.

 34:    Input Parameters:
 35: .  linesearch - linesearch context

 37:    Output Parameters:
 38: .  alpha - The descent parameter

 40:    Level: intermediate

 42: .seealso: SNESLineSearchGetLambda(), SNESLineSearchGetTolerances() SNESLINESEARCHBT
 43: @*/
 44: PetscErrorCode SNESLineSearchBTGetAlpha(SNESLineSearch linesearch, PetscReal *alpha)
 45: {
 46:   SNESLineSearch_BT *bt;

 50:   bt     = (SNESLineSearch_BT*)linesearch->data;
 51:   *alpha = bt->alpha;
 52:   return(0);
 53: }

 55: static PetscErrorCode  SNESLineSearchApply_BT(SNESLineSearch linesearch)
 56: {
 57:   PetscBool         changed_y,changed_w;
 58:   PetscErrorCode    ierr;
 59:   Vec               X,F,Y,W,G;
 60:   SNES              snes;
 61:   PetscReal         fnorm, xnorm, ynorm, gnorm;
 62:   PetscReal         lambda,lambdatemp,lambdaprev,minlambda,maxstep,initslope,alpha,stol;
 63:   PetscReal         t1,t2,a,b,d;
 64:   PetscReal         f;
 65:   PetscReal         g,gprev;
 66:   PetscViewer       monitor;
 67:   PetscInt          max_its,count;
 68:   SNESLineSearch_BT *bt;
 69:   Mat               jac;
 70:   PetscErrorCode    (*objective)(SNES,Vec,PetscReal*,void*);

 73:   SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G);
 74:   SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);
 75:   SNESLineSearchGetLambda(linesearch, &lambda);
 76:   SNESLineSearchGetSNES(linesearch, &snes);
 77:   SNESLineSearchGetDefaultMonitor(linesearch, &monitor);
 78:   SNESLineSearchGetTolerances(linesearch,&minlambda,&maxstep,NULL,NULL,NULL,&max_its);
 79:   SNESGetTolerances(snes,NULL,NULL,&stol,NULL,NULL);
 80:   SNESGetObjective(snes,&objective,NULL);
 81:   bt   = (SNESLineSearch_BT*)linesearch->data;
 82:   alpha = bt->alpha;

 84:   SNESGetJacobian(snes, &jac, NULL, NULL, NULL);
 85:   if (!jac && !objective) SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix");

 87:   SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);
 88:   SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);

 90:   VecNormBegin(Y, NORM_2, &ynorm);
 91:   VecNormBegin(X, NORM_2, &xnorm);
 92:   VecNormEnd(Y, NORM_2, &ynorm);
 93:   VecNormEnd(X, NORM_2, &xnorm);

 95:   if (ynorm == 0.0) {
 96:     if (monitor) {
 97:       PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
 98:       PetscViewerASCIIPrintf(monitor,"    Line search: Initial direction and size is 0\n");
 99:       PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
100:     }
101:     VecCopy(X,W);
102:     VecCopy(F,G);
103:     SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm);
104:     SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
105:     return(0);
106:   }
107:   if (ynorm > maxstep) {        /* Step too big, so scale back */
108:     if (monitor) {
109:       PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
110:       PetscViewerASCIIPrintf(monitor,"    Line search: Scaling step by %14.12e old ynorm %14.12e\n", (double)(maxstep/ynorm),(double)ynorm);
111:       PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
112:     }
113:     VecScale(Y,maxstep/(ynorm));
114:     ynorm = maxstep;
115:   }

117:   /* if the SNES has an objective set, use that instead of the function value */
118:   if (objective) {
119:     SNESComputeObjective(snes,X,&f);
120:   } else {
121:     f = fnorm*fnorm;
122:   }

124:   /* compute the initial slope */
125:   if (objective) {
126:     /* slope comes from the function (assumed to be the gradient of the objective */
127:     VecDotRealPart(Y,F,&initslope);
128:   } else {
129:     /* slope comes from the normal equations */
130:     MatMult(jac,Y,W);
131:     VecDotRealPart(F,W,&initslope);
132:     if (initslope > 0.0)  initslope = -initslope;
133:     if (initslope == 0.0) initslope = -1.0;
134:   }

136:   VecWAXPY(W,-lambda,Y,X);
137:   if (linesearch->ops->viproject) {
138:     (*linesearch->ops->viproject)(snes, W);
139:   }
140:   if (snes->nfuncs >= snes->max_funcs) {
141:     PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");
142:     snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
143:     SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);
144:     return(0);
145:   }

147:   if (objective) {
148:     SNESComputeObjective(snes,W,&g);
149:   } else {
150:     (*linesearch->ops->snesfunc)(snes,W,G);
151:     if (linesearch->ops->vinorm) {
152:       gnorm = fnorm;
153:       (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
154:     } else {
155:       VecNorm(G,NORM_2,&gnorm);
156:     }
157:     g = PetscSqr(gnorm);
158:   }
159:   SNESLineSearchMonitor(linesearch);

161:   if (PetscIsInfOrNanReal(g)) {
162:     SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF);
163:     PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");
164:     return(0);
165:   }
166:   if (!objective) {
167:     PetscInfo2(snes,"Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);
168:   }
169:   if (.5*g <= .5*f + lambda*alpha*initslope) { /* Sufficient reduction or step tolerance convergence */
170:     if (monitor) {
171:       PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
172:       if (!objective) {
173:         PetscViewerASCIIPrintf(monitor,"    Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);
174:       } else {
175:         PetscViewerASCIIPrintf(monitor,"    Line search: Using full step: obj0 %14.12e obj %14.12e\n", (double)f, (double)g);
176:       }
177:       PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
178:     }
179:   } else {
180:     /* Since the full step didn't work and the step is tiny, quit */
181:     if (stol*xnorm > ynorm) {
182:       SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm);
183:       SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
184:       if (monitor) {
185:         PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
186:         PetscViewerASCIIPrintf(monitor,"    Line search: Aborted due to ynorm < stol*xnorm (%14.12e < %14.12e) and inadequate full step.\n",(double)ynorm,(double)stol*xnorm);
187:         PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
188:       }
189:       return(0);
190:     }
191:     /* Fit points with quadratic */
192:     lambdatemp = -initslope/(g - f - 2.0*lambda*initslope);
193:     lambdaprev = lambda;
194:     gprev      = g;
195:     if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
196:     if (lambdatemp <= .1*lambda) lambda = .1*lambda;
197:     else                         lambda = lambdatemp;

199:     VecWAXPY(W,-lambda,Y,X);
200:     if (linesearch->ops->viproject) {
201:       (*linesearch->ops->viproject)(snes, W);
202:     }
203:     if (snes->nfuncs >= snes->max_funcs) {
204:       PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);
205:       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
206:       SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);
207:       return(0);
208:     }
209:     if (objective) {
210:       SNESComputeObjective(snes,W,&g);
211:     } else {
212:       (*linesearch->ops->snesfunc)(snes,W,G);
213:       if (linesearch->ops->vinorm) {
214:         gnorm = fnorm;
215:         (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
216:       } else {
217:         VecNorm(G,NORM_2,&gnorm);
218:       }
219:       g = PetscSqr(gnorm);
220:     }
221:     if (PetscIsInfOrNanReal(g)) {
222:       SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF);
223:       PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");
224:       return(0);
225:     }
226:     if (monitor) {
227:       PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
228:       if (!objective) {
229:         PetscViewerASCIIPrintf(monitor,"    Line search: gnorm after quadratic fit %14.12e\n",(double)gnorm);
230:       } else {
231:         PetscViewerASCIIPrintf(monitor,"    Line search: obj after quadratic fit %14.12e\n",(double)g);
232:       }
233:       PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
234:     }
235:     if (.5*g < .5*f + lambda*alpha*initslope) { /* sufficient reduction */
236:       if (monitor) {
237:         PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
238:         PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda);
239:         PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
240:       }
241:     } else {
242:       /* Fit points with cubic */
243:       for (count = 0; count < max_its; count++) {
244:         if (lambda <= minlambda) {
245:           if (monitor) {
246:             PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
247:             PetscViewerASCIIPrintf(monitor,"    Line search: unable to find good step length! After %D tries \n",count);
248:             if (!objective) {
249:               PetscViewerASCIIPrintf(monitor,"    Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
250:                                                          (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);
251:             } else {
252:               PetscViewerASCIIPrintf(monitor,"    Line search: obj(0)=%18.16e, obj=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
253:                                                          (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);
254:             }
255:             PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
256:           }
257:           SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
258:           return(0);
259:         }
260:         if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
261:           t1 = .5*(g - f) - lambda*initslope;
262:           t2 = .5*(gprev  - f) - lambdaprev*initslope;
263:           a  = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
264:           b  = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
265:           d  = b*b - 3*a*initslope;
266:           if (d < 0.0) d = 0.0;
267:           if (a == 0.0) lambdatemp = -initslope/(2.0*b);
268:           else lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a);

270:         } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
271:           lambdatemp = -initslope/(g - f - 2.0*initslope);
272:         } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt");
273:         lambdaprev = lambda;
274:         gprev      = g;
275:         if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
276:         if (lambdatemp <= .1*lambda) lambda     = .1*lambda;
277:         else                         lambda     = lambdatemp;
278:         VecWAXPY(W,-lambda,Y,X);
279:         if (linesearch->ops->viproject) {
280:           (*linesearch->ops->viproject)(snes,W);
281:         }
282:         if (snes->nfuncs >= snes->max_funcs) {
283:           PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);
284:           if (!objective) {
285:             PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",
286:                               (double)fnorm,(double)gnorm,(double)ynorm,(double)lambda,(double)initslope);
287:           }
288:           SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION);
289:           snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
290:           return(0);
291:         }
292:         if (objective) {
293:           SNESComputeObjective(snes,W,&g);
294:         } else {
295:           (*linesearch->ops->snesfunc)(snes,W,G);
296:           if (linesearch->ops->vinorm) {
297:             gnorm = fnorm;
298:             (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
299:           } else {
300:             VecNorm(G,NORM_2,&gnorm);
301:           }
302:           g = PetscSqr(gnorm);
303:         }
304:         if (PetscIsInfOrNanReal(g)) {
305:           SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF);
306:           PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");
307:           return(0);
308:         }
309:         if (.5*g < .5*f + lambda*alpha*initslope) { /* is reduction enough? */
310:           if (monitor) {
311:             PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
312:             if (!objective) {
313:               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
314:                 PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
315:               } else {
316:                 PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
317:               }
318:               PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
319:             } else {
320:               if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
321:                 PetscViewerASCIIPrintf(monitor,"    Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);
322:               } else {
323:                 PetscViewerASCIIPrintf(monitor,"    Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda);
324:               }
325:               PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
326:             }
327:           }
328:           break;
329:         } else if (monitor) {
330:           PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
331:           if (!objective) {
332:             if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
333:               PetscViewerASCIIPrintf(monitor,"    Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);
334:             } else {
335:               PetscViewerASCIIPrintf(monitor,"    Line search: Quadratic step no good, shrinking lambda, current gnorm %12.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: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);
341:             } else {
342:               PetscViewerASCIIPrintf(monitor,"    Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda);
343:             }
344:             PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
345:           }
346:         }
347:       }
348:     }
349:   }

351:   /* postcheck */
352:   SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
353:   if (changed_y) {
354:     VecWAXPY(W,-lambda,Y,X);
355:     if (linesearch->ops->viproject) {
356:       (*linesearch->ops->viproject)(snes, W);
357:     }
358:   }
359:   if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */
360:     (*linesearch->ops->snesfunc)(snes,W,G);
361:     if (linesearch->ops->vinorm) {
362:       gnorm = fnorm;
363:       (*linesearch->ops->vinorm)(snes, G, W, &gnorm);
364:     } else {
365:       VecNorm(G,NORM_2,&gnorm);
366:     }
367:     VecNorm(Y,NORM_2,&ynorm);
368:     if (PetscIsInfOrNanReal(gnorm)) {
369:       SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_FAILED_NANORINF);
370:       PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n");
371:       return(0);
372:     }
373:   }

375:   /* copy the solution over */
376:   VecCopy(W, X);
377:   VecCopy(G, F);
378:   VecNorm(X, NORM_2, &xnorm);
379:   SNESLineSearchSetLambda(linesearch, lambda);
380:   SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);
381:   return(0);
382: }

384: PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer)
385: {
386:   PetscErrorCode    ierr;
387:   PetscBool         iascii;
388:   SNESLineSearch_BT *bt;

391:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
392:   bt   = (SNESLineSearch_BT*)linesearch->data;
393:   if (iascii) {
394:     if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) {
395:       PetscViewerASCIIPrintf(viewer, "  interpolation: cubic\n");
396:     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
397:       PetscViewerASCIIPrintf(viewer, "  interpolation: quadratic\n");
398:     }
399:     PetscViewerASCIIPrintf(viewer, "  alpha=%e\n", (double)bt->alpha);
400:   }
401:   return(0);
402: }


405: static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch)
406: {

410:   PetscFree(linesearch->data);
411:   return(0);
412: }


415: static PetscErrorCode SNESLineSearchSetFromOptions_BT(PetscOptionItems *PetscOptionsObject,SNESLineSearch linesearch)
416: {

418:   PetscErrorCode    ierr;
419:   SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data;

422:   PetscOptionsHead(PetscOptionsObject,"SNESLineSearch BT options");
423:   PetscOptionsReal("-snes_linesearch_alpha",   "Descent tolerance",        "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL);
424:   PetscOptionsTail();
425:   return(0);
426: }


429: /*MC
430:    SNESLINESEARCHBT - Backtracking line search.

432:    This line search finds the minimum of a polynomial fitting of the L2 norm of the
433:    function or the objective function if it is provided with SNESSetObjective(). If this fit does not satisfy the conditions for progress, the interval shrinks
434:    and the fit is reattempted at most max_it times or until lambda is below minlambda.

436:    Options Database Keys:
437: +  -snes_linesearch_alpha<1e-4> - slope descent parameter
438: .  -snes_linesearch_damping<1.0> - initial step length
439: .  -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the
440:                                        step is scaled back to be of this length at the beginning of the line search
441: .  -snes_linesearch_max_it<40> - maximum number of shrinking step
442: .  -snes_linesearch_minlambda<1e-12> - minimum step length allowed
443: -  -snes_linesearch_order<cubic,quadratic> - order of the approximation

445:    Level: advanced

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

451: .keywords: SNES, SNESLineSearch, damping

453: .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
454: M*/
455: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
456: {

458:   SNESLineSearch_BT *bt;
459:   PetscErrorCode    ierr;

462:   linesearch->ops->apply          = SNESLineSearchApply_BT;
463:   linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
464:   linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
465:   linesearch->ops->reset          = NULL;
466:   linesearch->ops->view           = SNESLineSearchView_BT;
467:   linesearch->ops->setup          = NULL;

469:   PetscNewLog(linesearch,&bt);

471:   linesearch->data    = (void*)bt;
472:   linesearch->max_its = 40;
473:   linesearch->order   = SNES_LINESEARCH_ORDER_CUBIC;
474:   bt->alpha           = 1e-4;
475:   return(0);
476: }