Actual source code: linesearchcp.c

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

  4: static PetscErrorCode SNESLineSearchApply_CP(SNESLineSearch linesearch)
  5: {
  6:   PetscBool   changed_y, changed_w;
  7:   Vec         X, Y, F, W;
  8:   SNES        snes;
  9:   PetscReal   xnorm, ynorm, gnorm, steptol, atol, rtol, ltol, maxstep;
 10:   PetscReal   lambda, lambda_old, lambda_update, delLambda;
 11:   PetscScalar fty, fty_init, fty_old, fty_mid1, fty_mid2, s;
 12:   PetscInt    i, max_its;
 13:   PetscViewer monitor;

 15:   SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL);
 16:   SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
 17:   SNESLineSearchGetSNES(linesearch, &snes);
 18:   SNESLineSearchGetLambda(linesearch, &lambda);
 19:   SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, &ltol, &max_its);
 20:   SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);
 21:   SNESLineSearchGetDefaultMonitor(linesearch, &monitor);

 23:   /* precheck */
 24:   SNESLineSearchPreCheck(linesearch, X, Y, &changed_y);
 25:   lambda_old = 0.0;

 27:   VecDot(F, Y, &fty_old);
 28:   if (PetscAbsScalar(fty_old) < atol * ynorm) {
 29:     if (monitor) {
 30:       PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel);
 31:       PetscViewerASCIIPrintf(monitor, "    Line search terminated at initial point because dot(F,Y) = %g < atol*||y|| = %g\n", (double)PetscAbsScalar(fty_old), (double)(atol * ynorm));
 32:       PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel);
 33:     }
 34:     SNESSetConvergedReason(linesearch->snes, SNES_CONVERGED_FNORM_ABS);
 35:     return 0;
 36:   }

 38:   fty_init = fty_old;

 40:   for (i = 0; i < max_its; i++) {
 41:     /* compute the norm at lambda */
 42:     VecWAXPY(W, -lambda, Y, X);
 43:     if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, W);
 44:     (*linesearch->ops->snesfunc)(snes, W, F);
 45:     VecDot(F, Y, &fty);

 47:     delLambda = lambda - lambda_old;

 49:     /* check for convergence */
 50:     if (PetscAbsReal(delLambda) < steptol * lambda) break;
 51:     if (PetscAbsScalar(fty) / PetscAbsScalar(fty_init) < rtol) break;
 52:     if (PetscAbsScalar(fty) < atol * ynorm && i > 0) break;
 53:     if (monitor) {
 54:       PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel);
 55:       PetscViewerASCIIPrintf(monitor, "    Line search: lambdas = [%g, %g], ftys = [%g, %g]\n", (double)lambda, (double)lambda_old, (double)PetscRealPart(fty), (double)PetscRealPart(fty_old));
 56:       PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel);
 57:     }

 59:     /* compute the search direction */
 60:     if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) {
 61:       s = (fty - fty_old) / delLambda;
 62:     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
 63:       VecWAXPY(W, -0.5 * (lambda + lambda_old), Y, X);
 64:       if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, W);
 65:       (*linesearch->ops->snesfunc)(snes, W, F);
 66:       VecDot(F, Y, &fty_mid1);
 67:       s = (3. * fty - 4. * fty_mid1 + fty_old) / delLambda;
 68:     } else {
 69:       VecWAXPY(W, -0.5 * (lambda + lambda_old), Y, X);
 70:       if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, W);
 71:       (*linesearch->ops->snesfunc)(snes, W, F);
 72:       VecDot(F, Y, &fty_mid1);
 73:       VecWAXPY(W, -(lambda + 0.5 * (lambda - lambda_old)), Y, X);
 74:       if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, W);
 75:       (*linesearch->ops->snesfunc)(snes, W, F);
 76:       VecDot(F, Y, &fty_mid2);
 77:       s = (2. * fty_mid2 + 3. * fty - 6. * fty_mid1 + fty_old) / (3. * delLambda);
 78:     }
 79:     /* if the solve is going in the wrong direction, fix it */
 80:     if (PetscRealPart(s) > 0.) s = -s;
 81:     if (s == 0.0) break;
 82:     lambda_update = lambda - PetscRealPart(fty / s);

 84:     /* switch directions if we stepped out of bounds */
 85:     if (lambda_update < steptol) lambda_update = lambda + PetscRealPart(fty / s);

 87:     if (PetscIsInfOrNanReal(lambda_update)) break;
 88:     if (lambda_update > maxstep) break;

 90:     /* compute the new state of the line search */
 91:     lambda_old = lambda;
 92:     lambda     = lambda_update;
 93:     fty_old    = fty;
 94:   }
 95:   /* construct the solution */
 96:   VecWAXPY(W, -lambda, Y, X);
 97:   if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, W);
 98:   /* postcheck */
 99:   SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w);
100:   if (changed_y) {
101:     VecAXPY(X, -lambda, Y);
102:     if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, X);
103:   } else {
104:     VecCopy(W, X);
105:   }
106:   (*linesearch->ops->snesfunc)(snes, X, F);

108:   SNESLineSearchComputeNorms(linesearch);
109:   SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);

111:   SNESLineSearchSetLambda(linesearch, lambda);

113:   if (monitor) {
114:     PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel);
115:     PetscViewerASCIIPrintf(monitor, "    Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm);
116:     PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel);
117:   }
118:   if (lambda <= steptol) SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
119:   return 0;
120: }

122: /*MC
123:    SNESLINESEARCHCP - Critical point line search. This line search assumes that there exists some
124:    artificial G(x) for which the `SNESFunction` F(x) = grad G(x).  Therefore, this line search seeks
125:    to find roots of dot(F, Y) via a secant method.

127:    Options Database Keys:
128: +  -snes_linesearch_minlambda <minlambda> - the minimum acceptable lambda
129: .  -snes_linesearch_maxstep <length> - the algorithm insures that a step length is never longer than this value
130: .  -snes_linesearch_damping <damping> - initial trial step length is scaled by this factor, default is 1.0
131: -  -snes_linesearch_max_it <max_it> - the maximum number of secant steps performed.

133:    Notes:
134:    This method does NOT use the objective function if it is provided with `SNESSetObjective()`.

136:    This method is the preferred line search for `SNESQN` and `SNESNCG`.

138:    Level: advanced

140: .seealso: `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
141: M*/
142: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_CP(SNESLineSearch linesearch)
143: {
144:   linesearch->ops->apply          = SNESLineSearchApply_CP;
145:   linesearch->ops->destroy        = NULL;
146:   linesearch->ops->setfromoptions = NULL;
147:   linesearch->ops->reset          = NULL;
148:   linesearch->ops->view           = NULL;
149:   linesearch->ops->setup          = NULL;
150:   linesearch->order               = SNES_LINESEARCH_ORDER_LINEAR;

152:   linesearch->max_its = 1;
153:   return 0;
154: }