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:   PetscFunctionBegin;
 16:   PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL));
 17:   PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm));
 18:   PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
 19:   PetscCall(SNESLineSearchGetLambda(linesearch, &lambda));
 20:   PetscCall(SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, &ltol, &max_its));
 21:   PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
 22:   PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor));

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

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

 39:   fty_init = fty_old;

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

 48:     delLambda = lambda - lambda_old;

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

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

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

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

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

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

112:   PetscCall(SNESLineSearchSetLambda(linesearch, lambda));

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

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

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

134:    Level: advanced

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

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

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

154:   linesearch->max_its = 1;
155:   PetscFunctionReturn(PETSC_SUCCESS);
156: }