Actual source code: linesearchcp.c

petsc-3.4.4 2014-03-13
  1: #include <petsc-private/linesearchimpl.h>
  2: #include <petscsnes.h>

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

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

 18:   PetscViewer monitor;

 21:   SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL);
 22:   SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
 23:   SNESLineSearchGetSNES(linesearch, &snes);
 24:   SNESLineSearchGetLambda(linesearch, &lambda);
 25:   SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, &ltol, &max_its);
 26:   SNESLineSearchSetSuccess(linesearch, PETSC_TRUE);
 27:   SNESLineSearchGetMonitor(linesearch, &monitor);

 29:   /* precheck */
 30:   SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);
 31:   lambda_old = 0.0;
 32:   VecDot(F, Y, &fty_old);
 33:   fty_init   = fty_old;

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

 37:     /* compute the norm at lambda */
 38:     VecCopy(X, W);
 39:     VecAXPY(W, -lambda, Y);
 40:     if (linesearch->ops->viproject) {
 41:       (*linesearch->ops->viproject)(snes, W);
 42:     }
 43:     SNESComputeFunction(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 && 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:       VecCopy(X, W);
 64:       VecAXPY(W, -0.5*(lambda + lambda_old), Y);
 65:       if (linesearch->ops->viproject) {
 66:         (*linesearch->ops->viproject)(snes, W);
 67:       }
 68:       SNESComputeFunction(snes, W, F);
 69:       VecDot(F, Y, &fty_mid1);
 70:       s    = (3.*fty - 4.*fty_mid1 + fty_old) / delLambda;
 71:     } else {
 72:       VecCopy(X, W);
 73:       VecAXPY(W, -0.5*(lambda + lambda_old), Y);
 74:       if (linesearch->ops->viproject) {
 75:         (*linesearch->ops->viproject)(snes, W);
 76:       }
 77:       SNESComputeFunction(snes, W, F);
 78:       VecDot(F, Y, &fty_mid1);
 79:       VecCopy(X, W);
 80:       VecAXPY(W, -(lambda + 0.5*(lambda - lambda_old)), Y);
 81:       if (linesearch->ops->viproject) {
 82:         (*linesearch->ops->viproject)(snes, W);
 83:       }
 84:       SNESComputeFunction(snes, W, F);
 85:       VecDot(F, Y, &fty_mid2);
 86:       s    = (2.*fty_mid2 + 3.*fty - 6.*fty_mid1 + fty_old) / (3.*delLambda);
 87:     }
 88:     /* if the solve is going in the wrong direction, fix it */
 89:     if (PetscRealPart(s) > 0.) s = -s;
 90:     lambda_update =  lambda - PetscRealPart(fty / s);

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

 95:     if (PetscIsInfOrNanScalar(lambda_update)) break;
 96:     if (lambda_update > maxstep) break;

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

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

129:   SNESLineSearchSetLambda(linesearch, lambda);

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

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

149:    Options Database Keys:
150: +  -snes_linesearch_minlambda - the minimum acceptable lambda
151: .  -snes_linesearch_damping - initial trial step length
152: -  -snes_linesearch_max_it  - the maximum number of secant steps performed.

154:    Notes:
155:    This method is the preferred line search for SNESQN and SNESNCG.

157:    Level: advanced

159: .keywords: SNES, SNESLineSearch, damping

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

174:   linesearch->max_its = 1;
175:   return(0);
176: }