Actual source code: mlSNES.c

  1: #ifdef PETSC_RCS_HEADER
  2: static char vcid[] = "$Id: mlSNES.c,v 1.18 2001/01/27 21:32:36 knepley Exp $";
  3: #endif

 5:  #include petscsnes.h
 6:  #include gvec.h

  8: #undef __FUNCT__
 10: /*@C
 11:   PCMultiLevelCubicLineSearch - This function performs a cubic line search, but
 12:   the function evaluations are P^T f so that the residual is computed on the
 13:   constrained manifold.

 15:   Input Parameters:
 16: + snes  - The SNES
 17: . ctx   - The optional user context
 18: . x     - The current iterate
 19: . f     - The residual evaluated at x (contains projected residual on output)
 20: . y     - The search direction (contains new iterate on output)
 21: . w     - A work vector
 22: - fnorm - The 2-norm of f

 24:   Output Parameters:
 25: + f     - The projected residual P^T g evaluated at new iterate y
 26: . g     - The residual evaluated at new iterate y
 27: . y     - The new iterate (contains search direction on input)
 28: . gnorm - The 2-norm of P^T g
 29: . ynorm - The 2-norm of y, the search length
 30: - flag  - 0 if line search succeeds; -1 on failure.

 32:   Level: advanced

 34:   Notes:
 35:   This line search is taken from "Numerical Methods for Unconstrained
 36:   Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.

 38: .keywords: SNES, nonlinear, line search, cubic

 40: .seealso: SNESNoLineSearch(), SNESNoLineSearch(), SNESSetLineSearch()
 41: @*/
 42: int PCMultiLevelCubicLineSearch(SNES snes, void *ctx, Vec x, Vec f, Vec g, Vec y, Vec w, PetscReal fnorm, PetscReal *ynorm, PetscReal *gnorm, int *flag)
 43: {
 44:   /*
 45:     Note that for line search purposes we work with with the related
 46:     minimization problem:
 47:       min  z(x):  R^n -> R,
 48:     where z(x) = .5 * fnorm*fnorm, and fnorm = || f ||_2.
 49:    */

 51:   SLES        sles;
 52:   PC          pc;
 53:   Mat         jac;
 54:   PetscReal   abstol, steptol, initslope, lambdaprev, gnormprev, a, b, d, t1, t2;
 55:   PetscReal   maxstep, minlambda, alpha, lambda, lambdatemp, lambdaneg;
 56: #if defined(PETSC_USE_COMPLEX)
 57:   PetscScalar cinitslope, clambda;
 58: #endif
 59:   PetscScalar scale;
 60:   PetscScalar mone = -1.0;
 61:   int         count;
 62:   int         ierr;

 65:   PetscLogEventBegin(SNES_LineSearch, snes, x, f, g);
 66:   SNESGetTolerances(snes, &abstol, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL);
 67:   SNESGetJacobian(snes, &jac, PETSC_NULL, PETSC_NULL, PETSC_NULL);
 68:   SNESGetLineSearchParams(snes, &alpha, &maxstep, &steptol);
 69:   SNESGetSLES(snes, &sles);
 70:   SLESGetPC(sles, &pc);
 71:   *flag = 0;

 73:   /* First we must recover the true solution */
 74:   PCMultiLevelBuildSolution(pc, y);

 76:   VecNorm(y, NORM_2, ynorm);
 77:   if (*ynorm < abstol) {
 78:     PetscLogInfo(snes, "PCMultiLevelCubicLineSearch: Search direction and size are nearly 0n");
 79:     *gnorm = fnorm;
 80:     VecCopy(x, y);
 81:     VecCopy(f, g);
 82:     PCApplySymmetricLeft(pc, g, f);
 83:     goto theend1;
 84:   }
 85:   if (*ynorm > maxstep) {
 86:     /* Step too big, so scale back */
 87:     scale  = maxstep/(*ynorm);
 88: #if defined(PETSC_USE_COMPLEX)
 89:     PetscLogInfo(snes, "PCMultiLevelCubicLineSearch: Scaling step by %gn", real(scale));
 90: #else
 91:     PetscLogInfo(snes, "PCMultiLevelCubicLineSearch: Scaling step by %gn", scale);
 92: #endif
 93:     VecScale(&scale, y);
 94:     *ynorm = maxstep;
 95:   }
 96:   minlambda = steptol/(*ynorm);
 97:   MatMult(jac, y, g);
 98:   PCApplySymmetricLeft(pc, g, w);
 99:   PCApplySymmetricLeft(pc, f, g);
100: #if defined(USE_PETSC_COMPLEX)
101:   VecDot(g, w, &cinitslope);
102:   initslope = real(cinitslope);
103: #else
104:   VecDot(g, w, &initslope);
105: #endif
106:   if (initslope > 0.0)  initslope = -initslope;
107:   if (initslope == 0.0) initslope = -1.0;

109:   VecCopy(y, w);
110:   VecAYPX(&mone, x, w);
111:   SNESComputeFunction(snes, w, g);
112:   PCApplySymmetricLeft(pc, g, f);
113:   VecNorm(f, NORM_2, gnorm);
114:   if ((*gnorm)*(*gnorm)*0.5 <= fnorm*fnorm*0.5 + alpha*initslope) {
115:     /* Sufficient reduction */
116:     VecCopy(w, y);
117:     PetscLogInfo(snes, "PCMultiLevelCubicLineSearch: Using full stepn");
118:     goto theend1;
119:   }

121:   /* Fit points with quadratic */
122:   count      = 0;
123:   lambda     = 1.0;
124:   lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope);
125:   lambdaprev = lambda;
126:   gnormprev  = *gnorm;
127:   if (lambdatemp <= 0.1*lambda) {
128:     lambda   = 0.1*lambda;
129:   } else {
130:     lambda   = lambdatemp;
131:   }
132:   VecCopy(x, w);
133:   lambdaneg  = -lambda;
134: #if defined(PETSC_USE_COMPLEX)
135:   clambda    = lambdaneg;
136:   VecAXPY(&clambda, y, w);
137: #else
138:   VecAXPY(&lambdaneg, y, w);
139: #endif
140:   SNESComputeFunction(snes, w, g);
141:   PCApplySymmetricLeft(pc, g, f);
142:   VecNorm(f, NORM_2, gnorm);
143:   if ((*gnorm)*(*gnorm)*0.5 <= fnorm*fnorm*0.5 + alpha*initslope) {
144:     /* Sufficient reduction */
145:     VecCopy(w, y);
146:     PetscLogInfo(snes, "PCMultiLevelCubicLineSearch: Quadratically determined step, lambda=%gn", lambda);
147:     goto theend1;
148:   }

150:   /* Fit points with cubic */
151:   count = 1;
152:   while (1) {
153:     if (lambda <= minlambda) {
154:       /* Bad luck; use full step */
155:       PetscLogInfo(snes, "PCMultiLevelCubicLineSearch:Unable to find good step length! %d n",count);
156:       PetscLogInfo(snes, "PCMultiLevelCubicLineSearch:fnorm=%g, gnorm=%g, ynorm=%g, lambda=%g, initial slope=%gn",
157:                    fnorm, *gnorm, *ynorm, lambda, initslope);
158:       ierr  = VecCopy(w, y);
159:       *flag = -1;
160:       break;
161:     }
162:     t1 = ((*gnorm)*(*gnorm)    - fnorm*fnorm)*0.5 - lambda*initslope;
163:     t2 = (gnormprev*gnormprev  - fnorm*fnorm)*0.5 - lambdaprev*initslope;
164:     a  = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
165:     b  = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
166:     d  = b*b - 3*a*initslope;
167:     if (d < 0.0) d = 0.0;
168:     if (a == 0.0) {
169:       lambdatemp = -initslope/(2.0*b);
170:     } else {
171:       lambdatemp = (-b + sqrt(d))/(3.0*a);
172:     }
173:     if (lambdatemp > .5*lambda) lambdatemp = .5*lambda;
174:     lambdaprev = lambda;
175:     gnormprev  = *gnorm;
176:     if (lambdatemp <= .1*lambda) {
177:       lambda   = .1*lambda;
178:     } else {
179:       lambda   = lambdatemp;
180:     }
181:     VecCopy(x, w);
182:     lambdaneg  = -lambda;
183: #if defined(USE_PETSC_COMPLEX)
184:     clambda    = lambdaneg;
185:     VecAXPY(&clambda, y, w);
186: #else
187:     VecAXPY(&lambdaneg, y, w);
188: #endif
189:     SNESComputeFunction(snes, w, g);
190:     PCApplySymmetricLeft(pc, g, f);
191:     VecNorm(f, NORM_2, gnorm);
192:     if ((*gnorm)*(*gnorm)*0.5 <= fnorm*fnorm*0.5 + alpha*initslope) {
193:       /* Sufficient reduction */
194:       VecCopy(w, y);
195:       PetscLogInfo(snes, "PCMultiLevelCubicLineSearch: Cubically determined step, lambda=%gn", lambda);
196:       goto theend1;
197:     } else {
198:       PetscLogInfo(snes, "PCMultiLevelCubicLineSearch: Cubic step no good, shrinking lambda, lambda=%gn", lambda);
199:     }
200:     count++;
201:   }

203:   theend1:
204:   PetscLogEventEnd(SNES_LineSearch, snes, x, f, g);
205:   return(0);
206: }

208: #undef __FUNCT__
210: /*@C
211:   PCMultiLevelUpdateSNES - This function updates the nonlinear iterate for the PC
212:   context at each Newton step so that the appropriate system reduction may be
213:   performed.

215:   Collective on SNES

217:   Input Parameters:
218: + snes - The SNES
219: - step - The Newton step

221:   Level: advanced

223:   Note:
224:   This is usually used internally to keep track of the current
225:   solution in an nonlinear iteration.

227: .keywords SNES, multilevel, update
228: .seealso PCMultiLevelCubicLineSearch
229: @*/
230: int PCMultiLevelUpdateSNES(SNES snes, int step)
231: {
232:   SLES  sles;
233:   PC    pc;
234:   GVec  sol;
235:   int   ierr;

239:   SNESGetSLES(snes, &sles);
240:   SLESGetPC(sles, &pc);
241:   SNESGetSolution(snes, &sol);
242:   PCMultiLevelSetNonlinearIterate(pc, sol);
243:   return(0);
244: }