Actual source code: linesearchl2.c

petsc-3.9.2 2018-05-20
Report Typos and Errors
  1: #include <petsc/private/linesearchimpl.h>
  2: #include <petscsnes.h>

  4: static PetscErrorCode  SNESLineSearchApply_L2(SNESLineSearch linesearch)
  5: {

  7:   PetscBool      changed_y, changed_w;
  9:   Vec            X;
 10:   Vec            F;
 11:   Vec            Y;
 12:   Vec            W;
 13:   SNES           snes;
 14:   PetscReal      gnorm;
 15:   PetscReal      ynorm;
 16:   PetscReal      xnorm;
 17:   PetscReal      steptol, maxstep, rtol, atol, ltol;

 19:   PetscViewer monitor;
 20:   PetscReal   lambda, lambda_old, lambda_mid, lambda_update, delLambda;
 21:   PetscReal   fnrm, fnrm_old, fnrm_mid;
 22:   PetscReal   delFnrm, delFnrm_old, del2Fnrm;
 23:   PetscInt    i, max_its;

 25:   PetscErrorCode (*objective)(SNES,Vec,PetscReal*,void*);

 28:   SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL);
 29:   SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
 30:   SNESLineSearchGetLambda(linesearch, &lambda);
 31:   SNESLineSearchGetSNES(linesearch, &snes);
 32:   SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);
 33:   SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, &ltol, &max_its);
 34:   SNESLineSearchGetDefaultMonitor(linesearch, &monitor);

 36:   SNESGetObjective(snes,&objective,NULL);

 38:   /* precheck */
 39:   SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);
 40:   lambda_old = 0.0;
 41:   if (!objective) {
 42:     fnrm_old = gnorm*gnorm;
 43:   } else {
 44:     SNESComputeObjective(snes,X,&fnrm_old);
 45:   }
 46:   lambda_mid = 0.5*(lambda + lambda_old);

 48:   for (i = 0; i < max_its; i++) {

 50:     while (PETSC_TRUE) {
 51:       VecCopy(X, W);
 52:       VecAXPY(W, -lambda_mid, Y);
 53:       if (linesearch->ops->viproject) {
 54:         (*linesearch->ops->viproject)(snes, W);
 55:       }
 56:       if (!objective) {
 57:         /* compute the norm at the midpoint */
 58:         (*linesearch->ops->snesfunc)(snes, W, F);
 59:         if (linesearch->ops->vinorm) {
 60:           fnrm_mid = gnorm;
 61:           (*linesearch->ops->vinorm)(snes, F, W, &fnrm_mid);
 62:         } else {
 63:           VecNorm(F,NORM_2,&fnrm_mid);
 64:         }

 66:         /* compute the norm at the new endpoit */
 67:         VecCopy(X, W);
 68:         VecAXPY(W, -lambda, Y);
 69:         if (linesearch->ops->viproject) {
 70:           (*linesearch->ops->viproject)(snes, W);
 71:         }
 72:         (*linesearch->ops->snesfunc)(snes, W, F);
 73:         if (linesearch->ops->vinorm) {
 74:           fnrm = gnorm;
 75:           (*linesearch->ops->vinorm)(snes, F, W, &fnrm);
 76:         } else {
 77:           VecNorm(F,NORM_2,&fnrm);
 78:         }
 79:         fnrm_mid = fnrm_mid*fnrm_mid;
 80:         fnrm = fnrm*fnrm;
 81:       } else {
 82:         /* compute the objective at the midpoint */
 83:         VecCopy(X, W);
 84:         VecAXPY(W, -lambda_mid, Y);
 85:         SNESComputeObjective(snes,W,&fnrm_mid);

 87:         /* compute the objective at the new endpoint */
 88:         VecCopy(X, W);
 89:         VecAXPY(W, -lambda, Y);
 90:         SNESComputeObjective(snes,W,&fnrm);
 91:       }
 92:       if (!PetscIsInfOrNanReal(fnrm)) break;
 93:       if (monitor) {
 94:         PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
 95:         PetscViewerASCIIPrintf(monitor,"    Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n",(double)lambda);
 96:         PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
 97:       }
 98:       if (lambda <= steptol) {
 99:         SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
100:         return(0);
101:       }
102:       maxstep = .95*lambda; /* forbid the search from ever going back to the "failed" length that generates Nan or Inf */
103:       lambda  = .5*(lambda + lambda_old);
104:       lambda_mid = .5*(lambda + lambda_old);
105:     }

107:     delLambda   = lambda - lambda_old;
108:     /* compute f'() at the end points using second order one sided differencing */
109:     delFnrm     = (3.*fnrm - 4.*fnrm_mid + 1.*fnrm_old) / delLambda;
110:     delFnrm_old = (-3.*fnrm_old + 4.*fnrm_mid -1.*fnrm) / delLambda;
111:     /* compute f''() at the midpoint using centered differencing */
112:     del2Fnrm    = (delFnrm - delFnrm_old) / delLambda;

114:     if (monitor) {
115:       PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
116:       if (!objective) {
117:         PetscViewerASCIIPrintf(monitor,"    Line search: lambdas = [%g, %g, %g], fnorms = [%g, %g, %g]\n",(double)lambda, (double)lambda_mid, (double)lambda_old, (double)PetscSqrtReal(fnrm), (double)PetscSqrtReal(fnrm_mid), (double)PetscSqrtReal(fnrm_old));
118:       } else {
119:         PetscViewerASCIIPrintf(monitor,"    Line search: lambdas = [%g, %g, %g], obj = [%g, %g, %g]\n",(double)lambda, (double)lambda_mid, (double)lambda_old, (double)fnrm, (double)fnrm_mid, (double)fnrm_old);
120:       }
121:       PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
122:     }

124:     /* compute the secant (Newton) update -- always go downhill */
125:     if (del2Fnrm > 0.) lambda_update = lambda - delFnrm / del2Fnrm;
126:     else lambda_update = lambda + delFnrm / del2Fnrm;

128:     if (lambda_update < steptol) lambda_update = 0.5*(lambda + lambda_old);

130:     if (PetscIsInfOrNanReal(lambda_update)) break;

132:     if (lambda_update > maxstep) break;

134:     /* update the endpoints and the midpoint of the bracketed secant region */
135:     lambda_old = lambda;
136:     lambda     = lambda_update;
137:     fnrm_old   = fnrm;
138:     lambda_mid = 0.5*(lambda + lambda_old);
139:   }
140:   /* construct the solution */
141:   VecCopy(X, W);
142:   VecAXPY(W, -lambda, Y);
143:   if (linesearch->ops->viproject) {
144:     (*linesearch->ops->viproject)(snes, W);
145:   }

147:   /* postcheck */
148:   SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
149:   if (changed_y) {
150:     VecAXPY(X, -lambda, Y);
151:     if (linesearch->ops->viproject) {
152:       (*linesearch->ops->viproject)(snes, X);
153:     }
154:   } else {
155:     VecCopy(W, X);
156:   }
157:   (*linesearch->ops->snesfunc)(snes,X,F);

159:   SNESLineSearchSetLambda(linesearch, lambda);
160:   SNESLineSearchComputeNorms(linesearch);
161:   SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);

163:   if (monitor) {
164:     PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
165:     PetscViewerASCIIPrintf(monitor,"    Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm);
166:     PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
167:   }
168:   if (lambda <= steptol) {
169:     SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
170:   }
171:   return(0);
172: }

174: /*MC
175:    SNESLINESEARCHL2 - Secant search in the L2 norm of the function or the objective function, if it is provided with SNESSetObjective().

177:    Attempts to solve min_lambda f(x + lambda y) using the secant method with the initial bracketing of lambda between [0,damping]. Differences of f()
178:    are used to approximate the first and second derivative of f() with respect to lambda, f'() and f''(). The secant method is run for maxit iterations.

180:    When an objective function is provided f(w) is the objective function otherwise f(w) = ||F(w)||^2. x is the current step and y is the search direction.

182:    This has no checks on whether the secant method is actually converging.

184:    Options Database Keys:
185: +  -snes_linesearch_max_it <maxit> - maximum number of iterations, default is 1
186: .  -snes_linesearch_maxstep <length> - the algorithm insures that a step length is never longer than this value
187: .  -snes_linesearch_damping <damping> - initial step is scaled back by this factor, default is 1.0
188: -  -snes_linesearch_minlambda <minlambda> - minimum allowable lambda

190:    Level: advanced

192:    Developer Notes: A better name for this method might be SNESLINESEARCHSECANT, L2 is not descriptive

194: .keywords: SNES, nonlinear, line search, norm, secant

196: .seealso: SNESLINESEARCHBT, SNESLINESEARCHCP, SNESLineSearch, SNESLineSearchCreate(), SNESLineSearchSetType()
197: M*/
198: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_L2(SNESLineSearch linesearch)
199: {
201:   linesearch->ops->apply          = SNESLineSearchApply_L2;
202:   linesearch->ops->destroy        = NULL;
203:   linesearch->ops->setfromoptions = NULL;
204:   linesearch->ops->reset          = NULL;
205:   linesearch->ops->view           = NULL;
206:   linesearch->ops->setup          = NULL;

208:   linesearch->max_its = 1;
209:   return(0);
210: }