Actual source code: linesearchcp.c
petsc-3.4.5 2014-06-29
1: #include <petsc-private/linesearchimpl.h>
2: #include <petscsnes.h>
6: static PetscErrorCode SNESLineSearchApply_CP(SNESLineSearch linesearch)
7: {
8: PetscBool changed_y, changed_w, domainerror;
10: Vec X, Y, F, W;
11: SNES snes;
12: PetscReal xnorm, ynorm, gnorm, steptol, atol, rtol, ltol, maxstep;
14: PetscReal lambda, lambda_old, lambda_update, delLambda;
15: PetscScalar fty, fty_init, fty_old, fty_mid1, fty_mid2, s;
16: PetscInt i, max_its;
18: PetscViewer monitor;
21: SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL);
22: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
23: SNESLineSearchGetSNES(linesearch, &snes);
24: SNESLineSearchGetLambda(linesearch, &lambda);
25: SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, <ol, &max_its);
26: SNESLineSearchSetSuccess(linesearch, PETSC_TRUE);
27: SNESLineSearchGetMonitor(linesearch, &monitor);
29: /* precheck */
30: SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);
31: lambda_old = 0.0;
32: VecDot(F, Y, &fty_old);
33: fty_init = fty_old;
35: for (i = 0; i < max_its; i++) {
37: /* compute the norm at lambda */
38: VecCopy(X, W);
39: VecAXPY(W, -lambda, Y);
40: if (linesearch->ops->viproject) {
41: (*linesearch->ops->viproject)(snes, W);
42: }
43: SNESComputeFunction(snes, W, F);
45: VecDot(F, Y, &fty);
47: delLambda = lambda - lambda_old;
49: /* check for convergence */
50: if (PetscAbsReal(delLambda) < steptol*lambda) break;
51: if (PetscAbsScalar(fty) / PetscAbsScalar(fty_init) < rtol) break;
52: if (PetscAbsScalar(fty) < atol && i > 0) break;
53: if (monitor) {
54: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
55: PetscViewerASCIIPrintf(monitor," Line search: lambdas = [%g, %g], ftys = [%g, %g]\n",(double)lambda, (double)lambda_old, (double)PetscRealPart(fty), (double)PetscRealPart(fty_old));
56: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
57: }
59: /* compute the search direction */
60: if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) {
61: s = (fty - fty_old) / delLambda;
62: } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
63: VecCopy(X, W);
64: VecAXPY(W, -0.5*(lambda + lambda_old), Y);
65: if (linesearch->ops->viproject) {
66: (*linesearch->ops->viproject)(snes, W);
67: }
68: SNESComputeFunction(snes, W, F);
69: VecDot(F, Y, &fty_mid1);
70: s = (3.*fty - 4.*fty_mid1 + fty_old) / delLambda;
71: } else {
72: VecCopy(X, W);
73: VecAXPY(W, -0.5*(lambda + lambda_old), Y);
74: if (linesearch->ops->viproject) {
75: (*linesearch->ops->viproject)(snes, W);
76: }
77: SNESComputeFunction(snes, W, F);
78: VecDot(F, Y, &fty_mid1);
79: VecCopy(X, W);
80: VecAXPY(W, -(lambda + 0.5*(lambda - lambda_old)), Y);
81: if (linesearch->ops->viproject) {
82: (*linesearch->ops->viproject)(snes, W);
83: }
84: SNESComputeFunction(snes, W, F);
85: VecDot(F, Y, &fty_mid2);
86: s = (2.*fty_mid2 + 3.*fty - 6.*fty_mid1 + fty_old) / (3.*delLambda);
87: }
88: /* if the solve is going in the wrong direction, fix it */
89: if (PetscRealPart(s) > 0.) s = -s;
90: lambda_update = lambda - PetscRealPart(fty / s);
92: /* switch directions if we stepped out of bounds */
93: if (lambda_update < steptol) lambda_update = lambda + PetscRealPart(fty / s);
95: if (PetscIsInfOrNanScalar(lambda_update)) break;
96: if (lambda_update > maxstep) break;
98: /* compute the new state of the line search */
99: lambda_old = lambda;
100: lambda = lambda_update;
101: fty_old = fty;
102: }
103: /* construct the solution */
104: VecCopy(X, W);
105: VecAXPY(W, -lambda, Y);
106: if (linesearch->ops->viproject) {
107: (*linesearch->ops->viproject)(snes, W);
108: }
109: /* postcheck */
110: SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
111: if (changed_y) {
112: VecAXPY(X, -lambda, Y);
113: if (linesearch->ops->viproject) {
114: (*linesearch->ops->viproject)(snes, X);
115: }
116: } else {
117: VecCopy(W, X);
118: }
119: SNESComputeFunction(snes,X,F);
120: SNESGetFunctionDomainError(snes, &domainerror);
121: if (domainerror) {
122: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
123: return(0);
124: }
126: SNESLineSearchComputeNorms(linesearch);
127: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
129: SNESLineSearchSetLambda(linesearch, lambda);
131: if (monitor) {
132: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
133: PetscViewerASCIIPrintf(monitor," Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm);
134: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
135: }
136: if (lambda <= steptol) {
137: SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);
138: }
139: return(0);
140: }
144: /*MC
145: SNESLINESEARCHCP - Critical point line search. This line search assumes that there exists some
146: artificial G(x) for which the SNESFunction F(x) = grad G(x). Therefore, this line search seeks
147: to find roots of dot(F, Y) via a secant method.
149: Options Database Keys:
150: + -snes_linesearch_minlambda - the minimum acceptable lambda
151: . -snes_linesearch_damping - initial trial step length
152: - -snes_linesearch_max_it - the maximum number of secant steps performed.
154: Notes:
155: This method is the preferred line search for SNESQN and SNESNCG.
157: Level: advanced
159: .keywords: SNES, SNESLineSearch, damping
161: .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
162: M*/
163: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_CP(SNESLineSearch linesearch)
164: {
166: linesearch->ops->apply = SNESLineSearchApply_CP;
167: linesearch->ops->destroy = NULL;
168: linesearch->ops->setfromoptions = NULL;
169: linesearch->ops->reset = NULL;
170: linesearch->ops->view = NULL;
171: linesearch->ops->setup = NULL;
172: linesearch->order = SNES_LINESEARCH_ORDER_LINEAR;
174: linesearch->max_its = 1;
175: return(0);
176: }