Actual source code: linesearchl2.c
petsc-3.3-p7 2013-05-11
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;
30: SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, PETSC_NULL);
31: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
32: SNESLineSearchGetLambda(linesearch, &lambda);
33: SNESLineSearchGetSNES(linesearch, &snes);
34: SNESLineSearchSetSuccess(linesearch, PETSC_TRUE);
35: SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, <ol, &max_its);
36: SNESLineSearchGetMonitor(linesearch, &monitor);
38: /* precheck */
39: SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);
40: lambda_old = 0.0;
41: fnrm_old = gnorm*gnorm;
42: lambda_mid = 0.5*(lambda + lambda_old);
44: for (i = 0; i < max_its; i++) {
46: /* compute the norm at the midpoint */
47: VecCopy(X, W);
48: VecAXPY(W, -lambda_mid, Y);
49: if (linesearch->ops->viproject) {
50: (*linesearch->ops->viproject)(snes, W);
51: }
52: SNESComputeFunction(snes, W, F);
53: if (linesearch->ops->vinorm) {
54: fnrm_mid = gnorm;
55: (*linesearch->ops->vinorm)(snes, F, W, &fnrm_mid);
56: } else {
57: VecNorm(F, NORM_2, &fnrm_mid);
58: }
59: fnrm_mid = fnrm_mid*fnrm_mid;
61: /* compute the norm at lambda */
62: VecCopy(X, W);
63: VecAXPY(W, -lambda, Y);
64: if (linesearch->ops->viproject) {
65: (*linesearch->ops->viproject)(snes, W);
66: }
67: SNESComputeFunction(snes, W, F);
68: if (linesearch->ops->vinorm) {
69: fnrm = gnorm;
70: (*linesearch->ops->vinorm)(snes, F, W, &fnrm);
71: } else {
72: VecNorm(F, NORM_2, &fnrm);
73: }
74: fnrm = fnrm*fnrm;
76: /* this gives us the derivatives at the endpoints -- compute them from here
78: a = x - a0
80: p_0(x) = (x / dA - 1)(2x / dA - 1)
81: p_1(x) = 4(x / dA)(1 - x / dA)
82: p_2(x) = (x / dA)(2x / dA - 1)
84: dp_0[0] / dx = 3 / dA
85: dp_1[0] / dx = -4 / dA
86: dp_2[0] / dx = 1 / dA
88: dp_0[dA] / dx = -1 / dA
89: dp_1[dA] / dx = 4 / dA
90: dp_2[dA] / dx = -3 / dA
92: d^2p_0[0] / dx^2 = 4 / dA^2
93: d^2p_1[0] / dx^2 = -8 / dA^2
94: d^2p_2[0] / dx^2 = 4 / dA^2
95: */
97: delLambda = lambda - lambda_old;
98: delFnrm = (3.*fnrm - 4.*fnrm_mid + 1.*fnrm_old) / delLambda;
99: delFnrm_old = (-3.*fnrm_old + 4.*fnrm_mid -1.*fnrm) / delLambda;
100: del2Fnrm = (delFnrm - delFnrm_old) / delLambda;
102: if (monitor) {
103: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
104: PetscViewerASCIIPrintf(monitor," Line search: lambdas = [%g, %g, %g], fnorms = [%g, %g, %g]\n",
105: lambda, lambda_mid, lambda_old, PetscSqrtReal(fnrm), PetscSqrtReal(fnrm_mid), PetscSqrtReal(fnrm_old));
106: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
107: }
109: /* compute the search direction -- always go downhill */
110: if (del2Fnrm > 0.) {
111: lambda_update = lambda - delFnrm / del2Fnrm;
112: } else {
113: lambda_update = lambda + delFnrm / del2Fnrm;
114: }
116: if (lambda_update < steptol) {
117: lambda_update = 0.5*(lambda + lambda_old);
118: }
120: if (PetscIsInfOrNanScalar(lambda_update)) break;
122: if (lambda_update > maxstep) {
123: break;
124: }
126: /* compute the new state of the line search */
127: lambda_old = lambda;
128: lambda = lambda_update;
129: fnrm_old = fnrm;
130: lambda_mid = 0.5*(lambda + lambda_old);
131: }
132: /* construct the solution */
133: VecCopy(X, W);
134: VecAXPY(W, -lambda, Y);
135: if (linesearch->ops->viproject) {
136: (*linesearch->ops->viproject)(snes, W);
137: }
139: /* postcheck */
140: SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
141: if (changed_y) {
142: VecAXPY(X, -lambda, Y);
143: if (linesearch->ops->viproject) {
144: (*linesearch->ops->viproject)(snes, X);
145: }
146: } else {
147: VecCopy(W, X);
148: }
149: SNESComputeFunction(snes,X,F);
150: SNESGetFunctionDomainError(snes, &domainerror);
151: if (domainerror) {
152: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
153: return(0);
154: }
156: SNESLineSearchSetLambda(linesearch, lambda);
157: SNESLineSearchComputeNorms(linesearch);
158: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
160: if (monitor) {
161: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
162: PetscViewerASCIIPrintf(monitor," Line search terminated: lambda = %g, fnorms = %g\n", lambda, gnorm);
163: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
164: }
165: if (lambda <= steptol) {
166: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
167: }
168: return(0);
169: }
173: /*MC
174: SNESLINESEARCHL2 - Secant search in the L2 norm of the function.
176: The function norm is evaluated at points in [0, damping] to construct
177: a polynomial fitting. This fitting is used to construct a new lambda
178: based upon secant descent. The process is repeated on the new
179: interval, [lambda, lambda_old], max_it - 1 times.
181: Options Database Keys:
182: + -snes_linesearch_max_it<1> - maximum number of iterations
183: . -snes_linesearch_damping<1.0> - initial steplength
184: - -snes_linesearch_minlambda - minimum allowable lambda
186: Level: advanced
188: .keywords: SNES, nonlinear, line search, norm, secant
190: .seealso: SNESLineSearchBT, SNESLineSearchCP, SNESLineSearch
191: M*/
193: {
195: linesearch->ops->apply = SNESLineSearchApply_L2;
196: linesearch->ops->destroy = PETSC_NULL;
197: linesearch->ops->setfromoptions = PETSC_NULL;
198: linesearch->ops->reset = PETSC_NULL;
199: linesearch->ops->view = PETSC_NULL;
200: linesearch->ops->setup = PETSC_NULL;
202: linesearch->max_its = 1;
204: return(0);
205: }