Actual source code: linesearchcp.c
1: #include <petsc/private/linesearchimpl.h>
2: #include <petscsnes.h>
4: static PetscErrorCode SNESLineSearchApply_CP(SNESLineSearch linesearch)
5: {
6: PetscBool changed_y, changed_w;
7: Vec X, Y, F, W;
8: SNES snes;
9: PetscReal xnorm, ynorm, gnorm, steptol, atol, rtol, ltol, maxstep;
10: PetscReal lambda, lambda_old, lambda_update, delLambda;
11: PetscScalar fty, fty_init, fty_old, fty_mid1, fty_mid2, s;
12: PetscInt i, max_its;
13: PetscViewer monitor;
15: SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL);
16: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
17: SNESLineSearchGetSNES(linesearch, &snes);
18: SNESLineSearchGetLambda(linesearch, &lambda);
19: SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, <ol, &max_its);
20: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);
21: SNESLineSearchGetDefaultMonitor(linesearch, &monitor);
23: /* precheck */
24: SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);
25: lambda_old = 0.0;
27: VecDot(F,Y,&fty_old);
28: if (PetscAbsScalar(fty_old) < atol * ynorm) {
29: if (monitor) {
30: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
31: PetscViewerASCIIPrintf(monitor," Line search terminated at initial point because dot(F,Y) = %g < atol*||y|| = %g\n",(double)PetscAbsScalar(fty_old), (double)atol*ynorm);
32: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
33: }
34: SNESSetConvergedReason(linesearch->snes,SNES_CONVERGED_FNORM_ABS);
35: return 0;
36: }
38: fty_init = fty_old;
40: for (i = 0; i < max_its; i++) {
41: /* compute the norm at lambda */
42: VecCopy(X, W);
43: VecAXPY(W, -lambda, Y);
44: if (linesearch->ops->viproject) {
45: (*linesearch->ops->viproject)(snes, W);
46: }
47: (*linesearch->ops->snesfunc)(snes,W,F);
48: VecDot(F,Y,&fty);
50: delLambda = lambda - lambda_old;
52: /* check for convergence */
53: if (PetscAbsReal(delLambda) < steptol*lambda) break;
54: if (PetscAbsScalar(fty) / PetscAbsScalar(fty_init) < rtol) break;
55: if (PetscAbsScalar(fty) < atol * ynorm && i > 0) break;
56: if (monitor) {
57: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
58: PetscViewerASCIIPrintf(monitor," Line search: lambdas = [%g, %g], ftys = [%g, %g]\n",(double)lambda, (double)lambda_old, (double)PetscRealPart(fty), (double)PetscRealPart(fty_old));
59: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
60: }
62: /* compute the search direction */
63: if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) {
64: s = (fty - fty_old) / delLambda;
65: } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
66: VecCopy(X, W);
67: VecAXPY(W, -0.5*(lambda + lambda_old), Y);
68: if (linesearch->ops->viproject) {
69: (*linesearch->ops->viproject)(snes, W);
70: }
71: (*linesearch->ops->snesfunc)(snes,W,F);
72: VecDot(F, Y, &fty_mid1);
73: s = (3.*fty - 4.*fty_mid1 + fty_old) / delLambda;
74: } else {
75: VecCopy(X, W);
76: VecAXPY(W, -0.5*(lambda + lambda_old), Y);
77: if (linesearch->ops->viproject) {
78: (*linesearch->ops->viproject)(snes, W);
79: }
80: (*linesearch->ops->snesfunc)(snes,W,F);
81: VecDot(F, Y, &fty_mid1);
82: VecCopy(X, W);
83: VecAXPY(W, -(lambda + 0.5*(lambda - lambda_old)), Y);
84: if (linesearch->ops->viproject) {
85: (*linesearch->ops->viproject)(snes, W);
86: }
87: (*linesearch->ops->snesfunc)(snes, W, F);
88: VecDot(F, Y, &fty_mid2);
89: s = (2.*fty_mid2 + 3.*fty - 6.*fty_mid1 + fty_old) / (3.*delLambda);
90: }
91: /* if the solve is going in the wrong direction, fix it */
92: if (PetscRealPart(s) > 0.) s = -s;
93: if (s == 0.0) break;
94: lambda_update = lambda - PetscRealPart(fty / s);
96: /* switch directions if we stepped out of bounds */
97: if (lambda_update < steptol) lambda_update = lambda + PetscRealPart(fty / s);
99: if (PetscIsInfOrNanReal(lambda_update)) break;
100: if (lambda_update > maxstep) break;
102: /* compute the new state of the line search */
103: lambda_old = lambda;
104: lambda = lambda_update;
105: fty_old = fty;
106: }
107: /* construct the solution */
108: VecCopy(X, W);
109: VecAXPY(W, -lambda, Y);
110: if (linesearch->ops->viproject) {
111: (*linesearch->ops->viproject)(snes, W);
112: }
113: /* postcheck */
114: SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);
115: if (changed_y) {
116: VecAXPY(X, -lambda, Y);
117: if (linesearch->ops->viproject) {
118: (*linesearch->ops->viproject)(snes, X);
119: }
120: } else {
121: VecCopy(W, X);
122: }
123: (*linesearch->ops->snesfunc)(snes,X,F);
125: SNESLineSearchComputeNorms(linesearch);
126: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
128: SNESLineSearchSetLambda(linesearch, lambda);
130: if (monitor) {
131: PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);
132: PetscViewerASCIIPrintf(monitor," Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm);
133: PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);
134: }
135: if (lambda <= steptol) {
136: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
137: }
138: return 0;
139: }
141: /*MC
142: SNESLINESEARCHCP - Critical point line search. This line search assumes that there exists some
143: artificial G(x) for which the SNESFunction F(x) = grad G(x). Therefore, this line search seeks
144: to find roots of dot(F, Y) via a secant method.
146: Options Database Keys:
147: + -snes_linesearch_minlambda <minlambda> - the minimum acceptable lambda
148: . -snes_linesearch_maxstep <length> - the algorithm insures that a step length is never longer than this value
149: . -snes_linesearch_damping <damping> - initial trial step length is scaled by this factor, default is 1.0
150: - -snes_linesearch_max_it <max_it> - the maximum number of secant steps performed.
152: Notes:
153: This method does NOT use the objective function if it is provided with SNESSetObjective().
155: This method is the preferred line search for SNESQN and SNESNCG.
157: Level: advanced
159: .seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
160: M*/
161: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_CP(SNESLineSearch linesearch)
162: {
163: linesearch->ops->apply = SNESLineSearchApply_CP;
164: linesearch->ops->destroy = NULL;
165: linesearch->ops->setfromoptions = NULL;
166: linesearch->ops->reset = NULL;
167: linesearch->ops->view = NULL;
168: linesearch->ops->setup = NULL;
169: linesearch->order = SNES_LINESEARCH_ORDER_LINEAR;
171: linesearch->max_its = 1;
172: return 0;
173: }