Actual source code: linesearchl2.c
petsc-3.7.5 2017-01-01
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: PetscReal lambda, lambda_old, lambda_mid, lambda_update, delLambda;
23: PetscReal fnrm, fnrm_old, fnrm_mid;
24: PetscReal delFnrm, delFnrm_old, del2Fnrm;
25: PetscInt i, max_its;
27: PetscErrorCode (*objective)(SNES,Vec,PetscReal*,void*);
30: SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL);
31: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
32: SNESLineSearchGetLambda(linesearch, &lambda);
33: SNESLineSearchGetSNES(linesearch, &snes);
34: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);
35: SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, <ol, &max_its);
36: SNESLineSearchGetDefaultMonitor(linesearch, &monitor);
38: SNESGetObjective(snes,&objective,NULL);
40: /* precheck */
41: SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);
42: lambda_old = 0.0;
43: if (!objective) {
44: fnrm_old = gnorm*gnorm;
45: } else {
46: SNESComputeObjective(snes,X,&fnrm_old);
47: }
48: lambda_mid = 0.5*(lambda + lambda_old);
50: for (i = 0; i < max_its; i++) {
52: VecCopy(X, W);
53: VecAXPY(W, -lambda_mid, Y);
54: if (linesearch->ops->viproject) {
55: (*linesearch->ops->viproject)(snes, W);
56: }
57: if (!objective) {
58: /* compute the norm at the midpoint */
59: (*linesearch->ops->snesfunc)(snes, W, F);
60: if (linesearch->ops->vinorm) {
61: fnrm_mid = gnorm;
62: (*linesearch->ops->vinorm)(snes, F, W, &fnrm_mid);
63: } else {
64: VecNorm(F,NORM_2,&fnrm_mid);
65: }
67: /* compute the norm at the new endpoit */
68: VecCopy(X, W);
69: VecAXPY(W, -lambda, Y);
70: if (linesearch->ops->viproject) {
71: (*linesearch->ops->viproject)(snes, W);
72: }
73: (*linesearch->ops->snesfunc)(snes, W, F);
74: if (linesearch->ops->vinorm) {
75: fnrm = gnorm;
76: (*linesearch->ops->vinorm)(snes, F, W, &fnrm);
77: } else {
78: VecNorm(F,NORM_2,&fnrm);
79: }
80: fnrm_mid = fnrm_mid*fnrm_mid;
81: fnrm = fnrm*fnrm;
82: } else {
83: /* compute the objective at the midpoint */
84: VecCopy(X, W);
85: VecAXPY(W, -lambda_mid, Y);
86: SNESComputeObjective(snes,W,&fnrm_mid);
88: /* compute the objective at the new endpoint */
89: VecCopy(X, W);
90: VecAXPY(W, -lambda, Y);
91: SNESComputeObjective(snes,W,&fnrm);
92: }
94: delLambda = lambda - lambda_old;
95: /* compute f'() at the end points using second order one sided differencing */
96: delFnrm = (3.*fnrm - 4.*fnrm_mid + 1.*fnrm_old) / delLambda;
97: delFnrm_old = (-3.*fnrm_old + 4.*fnrm_mid -1.*fnrm) / delLambda;
98: /* compute f''() at the midpoint using centered differencing */
99: del2Fnrm = (delFnrm - delFnrm_old) / delLambda;
101: if (monitor) {
102: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
103: if (!objective) {
104: 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));
105: } else {
106: 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);
107: }
108: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
109: }
111: /* compute the secant (Newton) update -- always go downhill */
112: if (del2Fnrm > 0.) lambda_update = lambda - delFnrm / del2Fnrm;
113: else lambda_update = lambda + delFnrm / del2Fnrm;
115: if (lambda_update < steptol) lambda_update = 0.5*(lambda + lambda_old);
117: if (PetscIsInfOrNanReal(lambda_update)) break;
119: if (lambda_update > maxstep) break;
121: /* update the endpoints and the midpoint of the bracketed secant region */
122: lambda_old = lambda;
123: lambda = lambda_update;
124: fnrm_old = fnrm;
125: lambda_mid = 0.5*(lambda + lambda_old);
126: }
127: /* construct the solution */
128: VecCopy(X, W);
129: VecAXPY(W, -lambda, Y);
130: if (linesearch->ops->viproject) {
131: (*linesearch->ops->viproject)(snes, W);
132: }
134: /* postcheck */
135: SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
136: if (changed_y) {
137: VecAXPY(X, -lambda, Y);
138: if (linesearch->ops->viproject) {
139: (*linesearch->ops->viproject)(snes, X);
140: }
141: } else {
142: VecCopy(W, X);
143: }
144: (*linesearch->ops->snesfunc)(snes,X,F);
146: SNESLineSearchSetLambda(linesearch, lambda);
147: SNESLineSearchComputeNorms(linesearch);
148: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
150: if (monitor) {
151: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
152: PetscViewerASCIIPrintf(monitor," Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm);
153: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
154: }
155: if (lambda <= steptol) {
156: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
157: }
158: return(0);
159: }
163: /*MC
164: SNESLINESEARCHL2 - Secant search in the L2 norm of the function or the objective function, if it is provided with SNESSetObjective().
166: 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()
167: 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.
169: 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.
171: This has no checks on whether the secant method is actually converging.
173: Options Database Keys:
174: + -snes_linesearch_max_it <maxit> - maximum number of iterations, default is 1
175: . -snes_linesearch_maxstep <length> - the algorithm insures that a step length is never longer than this value
176: . -snes_linesearch_damping <damping> - initial step is scaled back by this factor, default is 1.0
177: - -snes_linesearch_minlambda <minlambda> - minimum allowable lambda
179: Level: advanced
181: Developer Notes: A better name for this method might be SNESLINESEARCHSECANT, L2 is not descriptive
183: .keywords: SNES, nonlinear, line search, norm, secant
185: .seealso: SNESLINESEARCHBT, SNESLINESEARCHCP, SNESLineSearch, SNESLineSearchCreate(), SNESLineSearchSetType()
186: M*/
187: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_L2(SNESLineSearch linesearch)
188: {
190: linesearch->ops->apply = SNESLineSearchApply_L2;
191: linesearch->ops->destroy = NULL;
192: linesearch->ops->setfromoptions = NULL;
193: linesearch->ops->reset = NULL;
194: linesearch->ops->view = NULL;
195: linesearch->ops->setup = NULL;
197: linesearch->max_its = 1;
198: return(0);
199: }