Actual source code: linesearchl2.c
petsc-3.6.1 2015-08-06
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: SNESLineSearchGetMonitor(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: /* compute the norm at the midpoint */
53: VecCopy(X, W);
54: VecAXPY(W, -lambda_mid, Y);
55: if (linesearch->ops->viproject) {
56: (*linesearch->ops->viproject)(snes, W);
57: }
58: if (!objective) {
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 lambda */
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 midpoint */
89: VecCopy(X, W);
90: VecAXPY(W, -lambda, Y);
91: SNESComputeObjective(snes,W,&fnrm);
92: }
93: /* this gives us the derivatives at the endpoints -- compute them from here
95: a = x - a0
97: p_0(x) = (x / dA - 1)(2x / dA - 1)
98: p_1(x) = 4(x / dA)(1 - x / dA)
99: p_2(x) = (x / dA)(2x / dA - 1)
101: dp_0[0] / dx = 3 / dA
102: dp_1[0] / dx = -4 / dA
103: dp_2[0] / dx = 1 / dA
105: dp_0[dA] / dx = -1 / dA
106: dp_1[dA] / dx = 4 / dA
107: dp_2[dA] / dx = -3 / dA
109: d^2p_0[0] / dx^2 = 4 / dA^2
110: d^2p_1[0] / dx^2 = -8 / dA^2
111: d^2p_2[0] / dx^2 = 4 / dA^2
112: */
114: delLambda = lambda - lambda_old;
115: delFnrm = (3.*fnrm - 4.*fnrm_mid + 1.*fnrm_old) / delLambda;
116: delFnrm_old = (-3.*fnrm_old + 4.*fnrm_mid -1.*fnrm) / delLambda;
117: del2Fnrm = (delFnrm - delFnrm_old) / delLambda;
119: if (monitor) {
120: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
121: if (!objective) {
122: PetscViewerASCIIPrintf(monitor," Line search: lambdas = [%g, %g, %g], fnorms = [%g, %g, %g]\n",
123: (double)lambda, (double)lambda_mid, (double)lambda_old, (double)PetscSqrtReal(fnrm), (double)PetscSqrtReal(fnrm_mid), (double)PetscSqrtReal(fnrm_old));
124: } else {
125: PetscViewerASCIIPrintf(monitor," Line search: lambdas = [%g, %g, %g], obj = [%g, %g, %g]\n",
126: (double)lambda, (double)lambda_mid, (double)lambda_old, (double)fnrm, (double)fnrm_mid, (double)fnrm_old);
128: }
129: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
130: }
132: /* compute the search direction -- always go downhill */
133: if (del2Fnrm > 0.) lambda_update = lambda - delFnrm / del2Fnrm;
134: else lambda_update = lambda + delFnrm / del2Fnrm;
136: if (lambda_update < steptol) lambda_update = 0.5*(lambda + lambda_old);
138: if (PetscIsInfOrNanReal(lambda_update)) break;
140: if (lambda_update > maxstep) break;
142: /* compute the new state of the line search */
143: lambda_old = lambda;
144: lambda = lambda_update;
145: fnrm_old = fnrm;
146: lambda_mid = 0.5*(lambda + lambda_old);
147: }
148: /* construct the solution */
149: VecCopy(X, W);
150: VecAXPY(W, -lambda, Y);
151: if (linesearch->ops->viproject) {
152: (*linesearch->ops->viproject)(snes, W);
153: }
155: /* postcheck */
156: SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
157: if (changed_y) {
158: VecAXPY(X, -lambda, Y);
159: if (linesearch->ops->viproject) {
160: (*linesearch->ops->viproject)(snes, X);
161: }
162: } else {
163: VecCopy(W, X);
164: }
165: (*linesearch->ops->snesfunc)(snes,X,F);
167: SNESLineSearchSetLambda(linesearch, lambda);
168: SNESLineSearchComputeNorms(linesearch);
169: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
171: if (monitor) {
172: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
173: PetscViewerASCIIPrintf(monitor," Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm);
174: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
175: }
176: if (lambda <= steptol) {
177: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
178: }
179: return(0);
180: }
184: /*MC
185: SNESLINESEARCHL2 - Secant search in the L2 norm of the function or the objective function if it is provided with SNESSetObjective().
187: The function norm is evaluated at points in [0, damping] to construct
188: a polynomial fitting. This fitting is used to construct a new lambda
189: based upon secant descent. The process is repeated on the new
190: interval, [lambda, lambda_old], max_it - 1 times.
192: Options Database Keys:
193: + -snes_linesearch_max_it <maxit> - maximum number of iterations, default is 1
194: . -snes_linesearch_maxstep <length> - the algorithm insures that a step length is never longer than this value
195: . -snes_linesearch_damping <damping> - initial step is scaled back by this factor, default is 1.0
196: - -snes_linesearch_minlambda <minlambda> - minimum allowable lambda
198: Level: advanced
200: .keywords: SNES, nonlinear, line search, norm, secant
202: .seealso: SNESLineSearchBT, SNESLineSearchCP, SNESLineSearch
203: M*/
204: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_L2(SNESLineSearch linesearch)
205: {
207: linesearch->ops->apply = SNESLineSearchApply_L2;
208: linesearch->ops->destroy = NULL;
209: linesearch->ops->setfromoptions = NULL;
210: linesearch->ops->reset = NULL;
211: linesearch->ops->view = NULL;
212: linesearch->ops->setup = NULL;
214: linesearch->max_its = 1;
215: return(0);
216: }