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