Actual source code: linesearchl2.c

  1: #include <petsc/private/linesearchimpl.h>
  2: #include <petscsnes.h>

  4: static PetscErrorCode SNESLineSearchApply_L2(SNESLineSearch linesearch)
  5: {
  6:   PetscBool   changed_y, changed_w;
  7:   Vec         X;
  8:   Vec         F;
  9:   Vec         Y;
 10:   Vec         W;
 11:   SNES        snes;
 12:   PetscReal   gnorm;
 13:   PetscReal   ynorm;
 14:   PetscReal   xnorm;
 15:   PetscReal   steptol, maxstep, rtol, atol, ltol;
 16:   PetscViewer monitor;
 17:   PetscReal   lambda, lambda_old, lambda_mid, lambda_update, delLambda;
 18:   PetscReal   fnrm, fnrm_old, fnrm_mid;
 19:   PetscReal   delFnrm, delFnrm_old, del2Fnrm;
 20:   PetscInt    i, max_its;
 21:   PetscErrorCode (*objective)(SNES, Vec, PetscReal *, void *);

 23:   SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL);
 24:   SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
 25:   SNESLineSearchGetLambda(linesearch, &lambda);
 26:   SNESLineSearchGetSNES(linesearch, &snes);
 27:   SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);
 28:   SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, &ltol, &max_its);
 29:   SNESLineSearchGetDefaultMonitor(linesearch, &monitor);

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

 33:   /* precheck */
 34:   SNESLineSearchPreCheck(linesearch, X, Y, &changed_y);
 35:   lambda_old = 0.0;
 36:   if (!objective) {
 37:     fnrm_old = gnorm * gnorm;
 38:   } else {
 39:     SNESComputeObjective(snes, X, &fnrm_old);
 40:   }
 41:   lambda_mid = 0.5 * (lambda + lambda_old);

 43:   for (i = 0; i < max_its; i++) {
 44:     while (PETSC_TRUE) {
 45:       VecWAXPY(W, -lambda_mid, Y, X);
 46:       if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, W);
 47:       if (!objective) {
 48:         /* compute the norm at the midpoint */
 49:         (*linesearch->ops->snesfunc)(snes, W, F);
 50:         if (linesearch->ops->vinorm) {
 51:           fnrm_mid = gnorm;
 52:           (*linesearch->ops->vinorm)(snes, F, W, &fnrm_mid);
 53:         } else {
 54:           VecNorm(F, NORM_2, &fnrm_mid);
 55:         }

 57:         /* compute the norm at the new endpoit */
 58:         VecWAXPY(W, -lambda, Y, X);
 59:         if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, W);
 60:         (*linesearch->ops->snesfunc)(snes, W, F);
 61:         if (linesearch->ops->vinorm) {
 62:           fnrm = gnorm;
 63:           (*linesearch->ops->vinorm)(snes, F, W, &fnrm);
 64:         } else {
 65:           VecNorm(F, NORM_2, &fnrm);
 66:         }
 67:         fnrm_mid = fnrm_mid * fnrm_mid;
 68:         fnrm     = fnrm * fnrm;
 69:       } else {
 70:         /* compute the objective at the midpoint */
 71:         VecWAXPY(W, -lambda_mid, Y, X);
 72:         SNESComputeObjective(snes, W, &fnrm_mid);

 74:         /* compute the objective at the new endpoint */
 75:         VecWAXPY(W, -lambda, Y, X);
 76:         SNESComputeObjective(snes, W, &fnrm);
 77:       }
 78:       if (!PetscIsInfOrNanReal(fnrm)) break;
 79:       if (monitor) {
 80:         PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel);
 81:         PetscViewerASCIIPrintf(monitor, "    Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n", (double)lambda);
 82:         PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel);
 83:       }
 84:       if (lambda <= steptol) {
 85:         SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
 86:         return 0;
 87:       }
 88:       maxstep    = .95 * lambda; /* forbid the search from ever going back to the "failed" length that generates Nan or Inf */
 89:       lambda     = .5 * (lambda + lambda_old);
 90:       lambda_mid = .5 * (lambda + lambda_old);
 91:     }

 93:     delLambda = lambda - lambda_old;
 94:     /* compute f'() at the end points using second order one sided differencing */
 95:     delFnrm     = (3. * fnrm - 4. * fnrm_mid + 1. * fnrm_old) / delLambda;
 96:     delFnrm_old = (-3. * fnrm_old + 4. * fnrm_mid - 1. * fnrm) / delLambda;
 97:     /* compute f''() at the midpoint using centered differencing */
 98:     del2Fnrm = (delFnrm - delFnrm_old) / delLambda;

100:     if (monitor) {
101:       PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel);
102:       if (!objective) {
103:         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));
104:       } else {
105:         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);
106:       }
107:       PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel);
108:     }

110:     /* compute the secant (Newton) update -- always go downhill */
111:     if (del2Fnrm > 0.) lambda_update = lambda - delFnrm / del2Fnrm;
112:     else if (del2Fnrm < 0.) lambda_update = lambda + delFnrm / del2Fnrm;
113:     else break;

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

117:     if (PetscIsInfOrNanReal(lambda_update)) break;

119:     if (lambda_update > maxstep) break;

121:     /* update the endpoints and the midpoint of the bracketed secant region */
122:     lambda_old = lambda;
123:     lambda     = lambda_update;
124:     fnrm_old   = fnrm;
125:     lambda_mid = 0.5 * (lambda + lambda_old);
126:   }
127:   /* construct the solution */
128:   VecWAXPY(W, -lambda, Y, X);
129:   if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, W);

131:   /* postcheck */
132:   SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w);
133:   if (changed_y) {
134:     VecAXPY(X, -lambda, Y);
135:     if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, X);
136:   } else {
137:     VecCopy(W, X);
138:   }
139:   (*linesearch->ops->snesfunc)(snes, X, F);

141:   SNESLineSearchSetLambda(linesearch, lambda);
142:   SNESLineSearchComputeNorms(linesearch);
143:   SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);

145:   if (monitor) {
146:     PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel);
147:     PetscViewerASCIIPrintf(monitor, "    Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm);
148:     PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel);
149:   }
150:   if (lambda <= steptol) SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
151:   return 0;
152: }

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

157:    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()
158:    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.

160:    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.

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

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

170:    Level: advanced

172:    Developer Note:
173:     A better name for this method might be `SNESLINESEARCHSECANT`, L2 is not descriptive

175: .seealso: `SNESLINESEARCHBT`, `SNESLINESEARCHCP`, `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
176: M*/
177: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_L2(SNESLineSearch linesearch)
178: {
179:   linesearch->ops->apply          = SNESLineSearchApply_L2;
180:   linesearch->ops->destroy        = NULL;
181:   linesearch->ops->setfromoptions = NULL;
182:   linesearch->ops->reset          = NULL;
183:   linesearch->ops->view           = NULL;
184:   linesearch->ops->setup          = NULL;

186:   linesearch->max_its = 1;
187:   return 0;
188: }