Actual source code: linesearchl2.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  1: #include <petsc-private/linesearchimpl.h>
  2: #include <petscsnes.h>

  6: static PetscErrorCode  SNESLineSearchApply_L2(SNESLineSearch linesearch)
  7: {

  9:   PetscBool      changed_y, changed_w;
 11:   Vec            X;
 12:   Vec            F;
 13:   Vec            Y;
 14:   Vec            W;
 15:   SNES           snes;
 16:   PetscReal      gnorm;
 17:   PetscReal      ynorm;
 18:   PetscReal      xnorm;
 19:   PetscReal      steptol, maxstep, rtol, atol, ltol;

 21:   PetscViewer monitor;
 22:   PetscBool   domainerror;
 23:   PetscReal   lambda, lambda_old, lambda_mid, lambda_update, delLambda;
 24:   PetscReal   fnrm, fnrm_old, fnrm_mid;
 25:   PetscReal   delFnrm, delFnrm_old, del2Fnrm;
 26:   PetscInt    i, max_its;

 28:   PetscErrorCode (*objective)(SNES,Vec,PetscReal*,void*);

 31:   SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL);
 32:   SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
 33:   SNESLineSearchGetLambda(linesearch, &lambda);
 34:   SNESLineSearchGetSNES(linesearch, &snes);
 35:   SNESLineSearchSetSuccess(linesearch, PETSC_TRUE);
 36:   SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, &ltol, &max_its);
 37:   SNESLineSearchGetMonitor(linesearch, &monitor);

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

 41:   /* precheck */
 42:   SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);
 43:   lambda_old = 0.0;
 44:   if (!objective) {
 45:     fnrm_old = gnorm*gnorm;
 46:   } else {
 47:     SNESComputeObjective(snes,X,&fnrm_old);
 48:   }
 49:   lambda_mid = 0.5*(lambda + lambda_old);

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

 53:     /* compute the norm at the midpoint */
 54:     VecCopy(X, W);
 55:     VecAXPY(W, -lambda_mid, Y);
 56:     if (linesearch->ops->viproject) {
 57:       (*linesearch->ops->viproject)(snes, W);
 58:     }
 59:     if (!objective) {
 60:       (*linesearch->ops->snesfunc)(snes, W, F);
 61:       if (linesearch->ops->vinorm) {
 62:         fnrm_mid = gnorm;
 63:         (*linesearch->ops->vinorm)(snes, F, W, &fnrm_mid);
 64:       } else {
 65:         VecNorm(F,NORM_2,&fnrm_mid);
 66:       }

 68:       /* compute the norm at lambda */
 69:       VecCopy(X, W);
 70:       VecAXPY(W, -lambda, Y);
 71:       if (linesearch->ops->viproject) {
 72:         (*linesearch->ops->viproject)(snes, W);
 73:       }
 74:       (*linesearch->ops->snesfunc)(snes, W, F);
 75:       if (linesearch->ops->vinorm) {
 76:         fnrm = gnorm;
 77:         (*linesearch->ops->vinorm)(snes, F, W, &fnrm);
 78:       } else {
 79:         VecNorm(F,NORM_2,&fnrm);
 80:       }
 81:       fnrm_mid = fnrm_mid*fnrm_mid;
 82:       fnrm = fnrm*fnrm;
 83:     } else {
 84:       /* compute the objective at the midpoint */
 85:       VecCopy(X, W);
 86:       VecAXPY(W, -lambda_mid, Y);
 87:       SNESComputeObjective(snes,W,&fnrm_mid);

 89:       /* compute the objective at the midpoint */
 90:       VecCopy(X, W);
 91:       VecAXPY(W, -lambda, Y);
 92:       SNESComputeObjective(snes,W,&fnrm);
 93:     }
 94:     /* this gives us the derivatives at the endpoints -- compute them from here

 96:      a = x - a0

 98:      p_0(x) = (x / dA - 1)(2x / dA - 1)
 99:      p_1(x) = 4(x / dA)(1 - x / dA)
100:      p_2(x) = (x / dA)(2x / dA - 1)

102:      dp_0[0] / dx = 3 / dA
103:      dp_1[0] / dx = -4 / dA
104:      dp_2[0] / dx = 1 / dA

106:      dp_0[dA] / dx = -1 / dA
107:      dp_1[dA] / dx = 4 / dA
108:      dp_2[dA] / dx = -3 / dA

110:      d^2p_0[0] / dx^2 =  4 / dA^2
111:      d^2p_1[0] / dx^2 = -8 / dA^2
112:      d^2p_2[0] / dx^2 =  4 / dA^2
113:      */

115:     delLambda   = lambda - lambda_old;
116:     delFnrm     = (3.*fnrm - 4.*fnrm_mid + 1.*fnrm_old) / delLambda;
117:     delFnrm_old = (-3.*fnrm_old + 4.*fnrm_mid -1.*fnrm) / delLambda;
118:     del2Fnrm    = (delFnrm - delFnrm_old) / delLambda;

120:     if (monitor) {
121:       PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
122:       if (!objective) {
123:         PetscViewerASCIIPrintf(monitor,"    Line search: lambdas = [%g, %g, %g], fnorms = [%g, %g, %g]\n",
124:                                       (double)lambda, (double)lambda_mid, (double)lambda_old, (double)PetscSqrtReal(fnrm), (double)PetscSqrtReal(fnrm_mid), (double)PetscSqrtReal(fnrm_old));
125:       } else {
126:         PetscViewerASCIIPrintf(monitor,"    Line search: lambdas = [%g, %g, %g], obj = [%g, %g, %g]\n",
127:                                       (double)lambda, (double)lambda_mid, (double)lambda_old, (double)fnrm, (double)fnrm_mid, (double)fnrm_old);

129:       }
130:       PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
131:     }

133:     /* compute the search direction -- always go downhill */
134:     if (del2Fnrm > 0.) lambda_update = lambda - delFnrm / del2Fnrm;
135:     else lambda_update = lambda + delFnrm / del2Fnrm;

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

139:     if (PetscIsInfOrNanScalar(lambda_update)) break;

141:     if (lambda_update > maxstep) break;

143:     /* compute the new state of the line search */
144:     lambda_old = lambda;
145:     lambda     = lambda_update;
146:     fnrm_old   = fnrm;
147:     lambda_mid = 0.5*(lambda + lambda_old);
148:   }
149:   /* construct the solution */
150:   VecCopy(X, W);
151:   VecAXPY(W, -lambda, Y);
152:   if (linesearch->ops->viproject) {
153:     (*linesearch->ops->viproject)(snes, W);
154:   }

156:   /* postcheck */
157:   SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
158:   if (changed_y) {
159:     VecAXPY(X, -lambda, Y);
160:     if (linesearch->ops->viproject) {
161:       (*linesearch->ops->viproject)(snes, X);
162:     }
163:   } else {
164:     VecCopy(W, X);
165:   }
166:   (*linesearch->ops->snesfunc)(snes,X,F);
167:   SNESGetFunctionDomainError(snes, &domainerror);
168:   if (domainerror) {
169:     SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
170:     return(0);
171:   }

173:   SNESLineSearchSetLambda(linesearch, lambda);
174:   SNESLineSearchComputeNorms(linesearch);
175:   SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);

177:   if (monitor) {
178:     PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
179:     PetscViewerASCIIPrintf(monitor,"    Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm);
180:     PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
181:   }
182:   if (lambda <= steptol) {
183:     SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
184:   }
185:   return(0);
186: }

190: /*MC
191:    SNESLINESEARCHL2 - Secant search in the L2 norm of the function.

193:    The function norm is evaluated at points in [0, damping] to construct
194:    a polynomial fitting.  This fitting is used to construct a new lambda
195:    based upon secant descent.  The process is repeated on the new
196:    interval, [lambda, lambda_old], max_it - 1 times.

198:    Options Database Keys:
199: +  -snes_linesearch_max_it<1> - maximum number of iterations
200: .  -snes_linesearch_damping<1.0> - initial steplength
201: -  -snes_linesearch_minlambda - minimum allowable lambda

203:    Level: advanced

205: .keywords: SNES, nonlinear, line search, norm, secant

207: .seealso: SNESLineSearchBT, SNESLineSearchCP, SNESLineSearch
208: M*/
209: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_L2(SNESLineSearch linesearch)
210: {
212:   linesearch->ops->apply          = SNESLineSearchApply_L2;
213:   linesearch->ops->destroy        = NULL;
214:   linesearch->ops->setfromoptions = NULL;
215:   linesearch->ops->reset          = NULL;
216:   linesearch->ops->view           = NULL;
217:   linesearch->ops->setup          = NULL;

219:   linesearch->max_its = 1;
220:   return(0);
221: }