Actual source code: linesearchbt.c

petsc-3.6.4 2016-04-12
Report Typos and Errors
  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: }