Actual source code: linesearchl2.c
1: #include <petsc/private/linesearchimpl.h>
2: #include <petscsnes.h>
4: static PetscErrorCode SNESLineSearchApply_L2(SNESLineSearch linesearch)
5: {
6: PetscBool changed_y, changed_w;
7: Vec X;
8: Vec F;
9: Vec Y;
10: Vec W;
11: SNES snes;
12: PetscReal gnorm;
13: PetscReal ynorm;
14: PetscReal xnorm;
15: PetscReal steptol, maxstep, rtol, atol, ltol;
16: PetscViewer monitor;
17: PetscReal lambda, lambda_old, lambda_mid, lambda_update, delLambda;
18: PetscReal fnrm, fnrm_old, fnrm_mid;
19: PetscReal delFnrm, delFnrm_old, del2Fnrm;
20: PetscInt i, max_its;
21: PetscErrorCode (*objective)(SNES,Vec,PetscReal*,void*);
23: SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL);
24: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
25: SNESLineSearchGetLambda(linesearch, &lambda);
26: SNESLineSearchGetSNES(linesearch, &snes);
27: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);
28: SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, <ol, &max_its);
29: SNESLineSearchGetDefaultMonitor(linesearch, &monitor);
31: SNESGetObjective(snes,&objective,NULL);
33: /* precheck */
34: SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);
35: lambda_old = 0.0;
36: if (!objective) {
37: fnrm_old = gnorm*gnorm;
38: } else {
39: SNESComputeObjective(snes,X,&fnrm_old);
40: }
41: lambda_mid = 0.5*(lambda + lambda_old);
43: for (i = 0; i < max_its; i++) {
45: while (PETSC_TRUE) {
46: VecCopy(X, W);
47: VecAXPY(W, -lambda_mid, Y);
48: if (linesearch->ops->viproject) {
49: (*linesearch->ops->viproject)(snes, W);
50: }
51: if (!objective) {
52: /* compute the norm at the midpoint */
53: (*linesearch->ops->snesfunc)(snes, W, F);
54: if (linesearch->ops->vinorm) {
55: fnrm_mid = gnorm;
56: (*linesearch->ops->vinorm)(snes, F, W, &fnrm_mid);
57: } else {
58: VecNorm(F,NORM_2,&fnrm_mid);
59: }
61: /* compute the norm at the new endpoit */
62: VecCopy(X, W);
63: VecAXPY(W, -lambda, Y);
64: if (linesearch->ops->viproject) {
65: (*linesearch->ops->viproject)(snes, W);
66: }
67: (*linesearch->ops->snesfunc)(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_mid = fnrm_mid*fnrm_mid;
75: fnrm = fnrm*fnrm;
76: } else {
77: /* compute the objective at the midpoint */
78: VecCopy(X, W);
79: VecAXPY(W, -lambda_mid, Y);
80: SNESComputeObjective(snes,W,&fnrm_mid);
82: /* compute the objective at the new endpoint */
83: VecCopy(X, W);
84: VecAXPY(W, -lambda, Y);
85: SNESComputeObjective(snes,W,&fnrm);
86: }
87: if (!PetscIsInfOrNanReal(fnrm)) break;
88: if (monitor) {
89: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
90: PetscViewerASCIIPrintf(monitor," Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n",(double)lambda);
91: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
92: }
93: if (lambda <= steptol) {
94: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
95: return 0;
96: }
97: maxstep = .95*lambda; /* forbid the search from ever going back to the "failed" length that generates Nan or Inf */
98: lambda = .5*(lambda + lambda_old);
99: lambda_mid = .5*(lambda + lambda_old);
100: }
102: delLambda = lambda - lambda_old;
103: /* compute f'() at the end points using second order one sided differencing */
104: delFnrm = (3.*fnrm - 4.*fnrm_mid + 1.*fnrm_old) / delLambda;
105: delFnrm_old = (-3.*fnrm_old + 4.*fnrm_mid -1.*fnrm) / delLambda;
106: /* compute f''() at the midpoint using centered differencing */
107: del2Fnrm = (delFnrm - delFnrm_old) / delLambda;
109: if (monitor) {
110: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
111: if (!objective) {
112: PetscViewerASCIIPrintf(monitor," Line search: lambdas = [%g, %g, %g], fnorms = [%g, %g, %g]\n",(double)lambda, (double)lambda_mid, (double)lambda_old, (double)PetscSqrtReal(fnrm), (double)PetscSqrtReal(fnrm_mid), (double)PetscSqrtReal(fnrm_old));
113: } else {
114: PetscViewerASCIIPrintf(monitor," Line search: lambdas = [%g, %g, %g], obj = [%g, %g, %g]\n",(double)lambda, (double)lambda_mid, (double)lambda_old, (double)fnrm, (double)fnrm_mid, (double)fnrm_old);
115: }
116: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
117: }
119: /* compute the secant (Newton) update -- always go downhill */
120: if (del2Fnrm > 0.) lambda_update = lambda - delFnrm / del2Fnrm;
121: else if (del2Fnrm < 0.) lambda_update = lambda + delFnrm / del2Fnrm;
122: else break;
124: if (lambda_update < steptol) lambda_update = 0.5*(lambda + lambda_old);
126: if (PetscIsInfOrNanReal(lambda_update)) break;
128: if (lambda_update > maxstep) break;
130: /* update the endpoints and the midpoint of the bracketed secant region */
131: lambda_old = lambda;
132: lambda = lambda_update;
133: fnrm_old = fnrm;
134: lambda_mid = 0.5*(lambda + lambda_old);
135: }
136: /* construct the solution */
137: VecCopy(X, W);
138: VecAXPY(W, -lambda, Y);
139: if (linesearch->ops->viproject) {
140: (*linesearch->ops->viproject)(snes, W);
141: }
143: /* postcheck */
144: SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
145: if (changed_y) {
146: VecAXPY(X, -lambda, Y);
147: if (linesearch->ops->viproject) {
148: (*linesearch->ops->viproject)(snes, X);
149: }
150: } else {
151: VecCopy(W, X);
152: }
153: (*linesearch->ops->snesfunc)(snes,X,F);
155: SNESLineSearchSetLambda(linesearch, lambda);
156: SNESLineSearchComputeNorms(linesearch);
157: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
159: if (monitor) {
160: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
161: PetscViewerASCIIPrintf(monitor," Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm);
162: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
163: }
164: if (lambda <= steptol) {
165: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
166: }
167: return 0;
168: }
170: /*MC
171: SNESLINESEARCHL2 - Secant search in the L2 norm of the function or the objective function, if it is provided with SNESSetObjective().
173: Attempts to solve min_lambda f(x + lambda y) using the secant method with the initial bracketing of lambda between [0,damping]. Differences of f()
174: are used to approximate the first and second derivative of f() with respect to lambda, f'() and f''(). The secant method is run for maxit iterations.
176: When an objective function is provided f(w) is the objective function otherwise f(w) = ||F(w)||^2. x is the current step and y is the search direction.
178: This has no checks on whether the secant method is actually converging.
180: Options Database Keys:
181: + -snes_linesearch_max_it <maxit> - maximum number of iterations, default is 1
182: . -snes_linesearch_maxstep <length> - the algorithm insures that a step length is never longer than this value
183: . -snes_linesearch_damping <damping> - initial step is scaled back by this factor, default is 1.0
184: - -snes_linesearch_minlambda <minlambda> - minimum allowable lambda
186: Level: advanced
188: Developer Notes:
189: A better name for this method might be SNESLINESEARCHSECANT, L2 is not descriptive
191: .seealso: SNESLINESEARCHBT, SNESLINESEARCHCP, SNESLineSearch, SNESLineSearchCreate(), SNESLineSearchSetType()
192: M*/
193: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_L2(SNESLineSearch linesearch)
194: {
195: linesearch->ops->apply = SNESLineSearchApply_L2;
196: linesearch->ops->destroy = NULL;
197: linesearch->ops->setfromoptions = NULL;
198: linesearch->ops->reset = NULL;
199: linesearch->ops->view = NULL;
200: linesearch->ops->setup = NULL;
202: linesearch->max_its = 1;
203: return 0;
204: }