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: }