Actual source code: linesearchcp.c

petsc-3.8.3 2017-12-09
Report Typos and Errors
  1:  #include <petsc/private/linesearchimpl.h>
  2:  #include <petscsnes.h>

  4: static PetscErrorCode SNESLineSearchApply_CP(SNESLineSearch linesearch)
  5: {
  6:   PetscBool      changed_y, changed_w;
  8:   Vec            X, Y, F, W;
  9:   SNES           snes;
 10:   PetscReal      xnorm, ynorm, gnorm, steptol, atol, rtol, ltol, maxstep;

 12:   PetscReal   lambda, lambda_old, lambda_update, delLambda;
 13:   PetscScalar fty, fty_init, fty_old, fty_mid1, fty_mid2, s;
 14:   PetscInt    i, max_its;

 16:   PetscViewer monitor;

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

 27:   /* precheck */
 28:   SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);
 29:   lambda_old = 0.0;

 31:   VecDot(F,Y,&fty_old);
 32:   if (PetscAbsScalar(fty_old) < atol) {
 33:     if (monitor) {
 34:       PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
 35:       PetscViewerASCIIPrintf(monitor,"    Line search terminated ended at initial point because dot(F,Y) = %g < atol = %g\n",(double)PetscAbsScalar(fty_old), (double)atol);
 36:       PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
 37:     }
 38:     return(0);
 39:   }

 41:   fty_init = fty_old;

 43:   for (i = 0; i < max_its; i++) {
 44:     /* compute the norm at lambda */
 45:     VecCopy(X, W);
 46:     VecAXPY(W, -lambda, Y);
 47:     if (linesearch->ops->viproject) {
 48:       (*linesearch->ops->viproject)(snes, W);
 49:     }
 50:     (*linesearch->ops->snesfunc)(snes,W,F);
 51:     VecDot(F,Y,&fty);

 53:     delLambda = lambda - lambda_old;

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

 65:     /* compute the search direction */
 66:     if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) {
 67:       s = (fty - fty_old) / delLambda;
 68:     } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
 69:       VecCopy(X, W);
 70:       VecAXPY(W, -0.5*(lambda + lambda_old), Y);
 71:       if (linesearch->ops->viproject) {
 72:         (*linesearch->ops->viproject)(snes, W);
 73:       }
 74:       (*linesearch->ops->snesfunc)(snes,W,F);
 75:       VecDot(F, Y, &fty_mid1);
 76:       s    = (3.*fty - 4.*fty_mid1 + fty_old) / delLambda;
 77:     } else {
 78:       VecCopy(X, W);
 79:       VecAXPY(W, -0.5*(lambda + lambda_old), Y);
 80:       if (linesearch->ops->viproject) {
 81:         (*linesearch->ops->viproject)(snes, W);
 82:       }
 83:       (*linesearch->ops->snesfunc)(snes,W,F);
 84:       VecDot(F, Y, &fty_mid1);
 85:       VecCopy(X, W);
 86:       VecAXPY(W, -(lambda + 0.5*(lambda - lambda_old)), Y);
 87:       if (linesearch->ops->viproject) {
 88:         (*linesearch->ops->viproject)(snes, W);
 89:       }
 90:       (*linesearch->ops->snesfunc)(snes, W, F);
 91:       VecDot(F, Y, &fty_mid2);
 92:       s    = (2.*fty_mid2 + 3.*fty - 6.*fty_mid1 + fty_old) / (3.*delLambda);
 93:     }
 94:     /* if the solve is going in the wrong direction, fix it */
 95:     if (PetscRealPart(s) > 0.) s = -s;
 96:     lambda_update =  lambda - PetscRealPart(fty / s);

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

101:     if (PetscIsInfOrNanReal(lambda_update)) break;
102:     if (lambda_update > maxstep) break;

104:     /* compute the new state of the line search */
105:     lambda_old = lambda;
106:     lambda     = lambda_update;
107:     fty_old    = fty;
108:   }
109:   /* construct the solution */
110:   VecCopy(X, W);
111:   VecAXPY(W, -lambda, Y);
112:   if (linesearch->ops->viproject) {
113:     (*linesearch->ops->viproject)(snes, W);
114:   }
115:   /* postcheck */
116:   SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
117:   if (changed_y) {
118:     VecAXPY(X, -lambda, Y);
119:     if (linesearch->ops->viproject) {
120:       (*linesearch->ops->viproject)(snes, X);
121:     }
122:   } else {
123:     VecCopy(W, X);
124:   }
125:   (*linesearch->ops->snesfunc)(snes,X,F);

127:   SNESLineSearchComputeNorms(linesearch);
128:   SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);

130:   SNESLineSearchSetLambda(linesearch, lambda);

132:   if (monitor) {
133:     PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
134:     PetscViewerASCIIPrintf(monitor,"    Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm);
135:     PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
136:   }
137:   if (lambda <= steptol) {
138:     SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
139:   }
140:   return(0);
141: }

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

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

154:    Notes:
155:    This method does NOT use the objective function if it is provided with SNESSetObjective().

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

159:    Level: advanced

161: .keywords: SNES, SNESLineSearch, damping

163: .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
164: M*/
165: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_CP(SNESLineSearch linesearch)
166: {
168:   linesearch->ops->apply          = SNESLineSearchApply_CP;
169:   linesearch->ops->destroy        = NULL;
170:   linesearch->ops->setfromoptions = NULL;
171:   linesearch->ops->reset          = NULL;
172:   linesearch->ops->view           = NULL;
173:   linesearch->ops->setup          = NULL;
174:   linesearch->order               = SNES_LINESEARCH_ORDER_LINEAR;

176:   linesearch->max_its = 1;
177:   return(0);
178: }