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: PetscFunctionBegin;
16: PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL));
17: PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm));
18: PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
19: PetscCall(SNESLineSearchGetLambda(linesearch, &lambda));
20: PetscCall(SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, <ol, &max_its));
21: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
22: PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor));
24: /* precheck */
25: PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y));
26: lambda_old = 0.0;
28: if (linesearch->ops->vidirderiv) {
29: PetscCall((*linesearch->ops->vidirderiv)(snes, F, X, Y, &fty_old));
30: } else {
31: PetscCall(VecDot(F, Y, &fty_old));
32: }
33: if (PetscAbsScalar(fty_old) < atol * ynorm) {
34: if (monitor) {
35: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
36: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search terminated at initial point because dot(F,Y) = %g < atol*||y|| = %g\n", (double)PetscAbsScalar(fty_old), (double)(atol * ynorm)));
37: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
38: }
39: PetscCall(SNESSetConvergedReason(linesearch->snes, SNES_CONVERGED_FNORM_ABS));
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: fty_init = fty_old;
45: for (i = 0; i < max_its; i++) {
46: /* compute the norm at lambda */
47: PetscCall(VecWAXPY(W, -lambda, Y, X));
48: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
49: PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
50: if (linesearch->ops->vidirderiv) {
51: PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty));
52: } else {
53: PetscCall(VecDot(F, Y, &fty));
54: }
56: delLambda = lambda - lambda_old;
58: /* check for convergence */
59: if (PetscAbsReal(delLambda) < steptol * lambda) break;
60: if (PetscAbsScalar(fty) / PetscAbsScalar(fty_init) < rtol) break;
61: if (PetscAbsScalar(fty) < atol * ynorm && i > 0) break;
62: if (monitor) {
63: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
64: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: lambdas = [%g, %g], ftys = [%g, %g]\n", (double)lambda, (double)lambda_old, (double)PetscRealPart(fty), (double)PetscRealPart(fty_old)));
65: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
66: }
68: /* compute the search direction */
69: if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) {
70: s = (fty - fty_old) / delLambda;
71: } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
72: PetscCall(VecWAXPY(W, -0.5 * (lambda + lambda_old), Y, X));
73: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
74: PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
75: if (linesearch->ops->vidirderiv) {
76: PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty_mid1));
77: } else {
78: PetscCall(VecDot(F, Y, &fty_mid1));
79: }
80: s = (3. * fty - 4. * fty_mid1 + fty_old) / delLambda;
81: } else {
82: PetscCall(VecWAXPY(W, -0.5 * (lambda + lambda_old), Y, X));
83: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
84: PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
85: if (linesearch->ops->vidirderiv) {
86: PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty_mid1));
87: } else {
88: PetscCall(VecDot(F, Y, &fty_mid1));
89: }
90: PetscCall(VecWAXPY(W, -(lambda + 0.5 * (lambda - lambda_old)), Y, X));
91: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
92: PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
93: if (linesearch->ops->vidirderiv) {
94: PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty_mid2));
95: } else {
96: PetscCall(VecDot(F, Y, &fty_mid2));
97: }
98: s = (2. * fty_mid2 + 3. * fty - 6. * fty_mid1 + fty_old) / (3. * delLambda);
99: }
100: /* if the solve is going in the wrong direction, fix it */
101: if (PetscRealPart(s) > 0.) s = -s;
102: if (s == 0.0) break;
103: lambda_update = lambda - PetscRealPart(fty / s);
105: /* switch directions if we stepped out of bounds */
106: if (lambda_update < steptol) lambda_update = lambda + PetscRealPart(fty / s);
108: if (PetscIsInfOrNanReal(lambda_update)) break;
109: if (lambda_update > maxstep) break;
111: /* compute the new state of the line search */
112: lambda_old = lambda;
113: lambda = lambda_update;
114: fty_old = fty;
115: }
116: /* construct the solution */
117: PetscCall(VecWAXPY(W, -lambda, Y, X));
118: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
119: /* postcheck */
120: PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
121: PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
122: if (changed_y) {
123: if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X));
124: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
125: }
126: PetscCall(VecCopy(W, X));
127: PetscCall((*linesearch->ops->snesfunc)(snes, X, F));
129: PetscCall(SNESLineSearchComputeNorms(linesearch));
130: PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm));
132: if (monitor) {
133: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
134: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm));
135: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
136: }
137: if (lambda <= steptol) PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: /*MC
142: SNESLINESEARCHCP - Critical point line search. This line search assumes that there exists some
143: artificial $G(x)$ for which the `SNESFunctionFn` $ F(x) = grad G(x)$. Therefore, this line search seeks
144: to find roots of $ F^T 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 on entry to the line search, default is 1.0
150: - -snes_linesearch_max_it <max_it> - the maximum number of secant steps performed.
152: Level: advanced
154: Notes:
155: This method does NOT use the objective function if it is provided with `SNESSetObjective()`.
157: This method is the preferred line search for `SNESQN` and `SNESNCG`.
159: .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`, `SNESLINESEARCHBISECTION`
160: M*/
161: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_CP(SNESLineSearch linesearch)
162: {
163: PetscFunctionBegin;
164: linesearch->ops->apply = SNESLineSearchApply_CP;
165: linesearch->ops->destroy = NULL;
166: linesearch->ops->setfromoptions = NULL;
167: linesearch->ops->reset = NULL;
168: linesearch->ops->view = NULL;
169: linesearch->ops->setup = NULL;
170: linesearch->order = SNES_LINESEARCH_ORDER_LINEAR;
172: linesearch->max_its = 1;
173: PetscFunctionReturn(PETSC_SUCCESS);
174: }