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: VecWAXPY(W, -lambda, Y, X);
43: if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, W);
44: (*linesearch->ops->snesfunc)(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 * ynorm && 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: VecWAXPY(W, -0.5 * (lambda + lambda_old), Y, X);
64: if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, W);
65: (*linesearch->ops->snesfunc)(snes, W, F);
66: VecDot(F, Y, &fty_mid1);
67: s = (3. * fty - 4. * fty_mid1 + fty_old) / delLambda;
68: } else {
69: VecWAXPY(W, -0.5 * (lambda + lambda_old), Y, X);
70: if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, W);
71: (*linesearch->ops->snesfunc)(snes, W, F);
72: VecDot(F, Y, &fty_mid1);
73: VecWAXPY(W, -(lambda + 0.5 * (lambda - lambda_old)), Y, X);
74: if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, W);
75: (*linesearch->ops->snesfunc)(snes, W, F);
76: VecDot(F, Y, &fty_mid2);
77: s = (2. * fty_mid2 + 3. * fty - 6. * fty_mid1 + fty_old) / (3. * delLambda);
78: }
79: /* if the solve is going in the wrong direction, fix it */
80: if (PetscRealPart(s) > 0.) s = -s;
81: if (s == 0.0) break;
82: lambda_update = lambda - PetscRealPart(fty / s);
84: /* switch directions if we stepped out of bounds */
85: if (lambda_update < steptol) lambda_update = lambda + PetscRealPart(fty / s);
87: if (PetscIsInfOrNanReal(lambda_update)) break;
88: if (lambda_update > maxstep) break;
90: /* compute the new state of the line search */
91: lambda_old = lambda;
92: lambda = lambda_update;
93: fty_old = fty;
94: }
95: /* construct the solution */
96: VecWAXPY(W, -lambda, Y, X);
97: if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, W);
98: /* postcheck */
99: SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w);
100: if (changed_y) {
101: VecAXPY(X, -lambda, Y);
102: if (linesearch->ops->viproject) (*linesearch->ops->viproject)(snes, X);
103: } else {
104: VecCopy(W, X);
105: }
106: (*linesearch->ops->snesfunc)(snes, X, F);
108: SNESLineSearchComputeNorms(linesearch);
109: SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);
111: SNESLineSearchSetLambda(linesearch, lambda);
113: if (monitor) {
114: PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel);
115: PetscViewerASCIIPrintf(monitor, " Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm);
116: PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel);
117: }
118: if (lambda <= steptol) SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);
119: return 0;
120: }
122: /*MC
123: SNESLINESEARCHCP - Critical point line search. This line search assumes that there exists some
124: artificial G(x) for which the `SNESFunction` F(x) = grad G(x). Therefore, this line search seeks
125: to find roots of dot(F, Y) via a secant method.
127: Options Database Keys:
128: + -snes_linesearch_minlambda <minlambda> - the minimum acceptable lambda
129: . -snes_linesearch_maxstep <length> - the algorithm insures that a step length is never longer than this value
130: . -snes_linesearch_damping <damping> - initial trial step length is scaled by this factor, default is 1.0
131: - -snes_linesearch_max_it <max_it> - the maximum number of secant steps performed.
133: Notes:
134: This method does NOT use the objective function if it is provided with `SNESSetObjective()`.
136: This method is the preferred line search for `SNESQN` and `SNESNCG`.
138: Level: advanced
140: .seealso: `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
141: M*/
142: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_CP(SNESLineSearch linesearch)
143: {
144: linesearch->ops->apply = SNESLineSearchApply_CP;
145: linesearch->ops->destroy = NULL;
146: linesearch->ops->setfromoptions = NULL;
147: linesearch->ops->reset = NULL;
148: linesearch->ops->view = NULL;
149: linesearch->ops->setup = NULL;
150: linesearch->order = SNES_LINESEARCH_ORDER_LINEAR;
152: linesearch->max_its = 1;
153: return 0;
154: }