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, minlambda, maxlambda, rtol, atol, ltol;
10: PetscReal lambda, lambda_old, lambda_update, delLambda;
11: PetscScalar fty, fty_init, fty_old, fty_mid1, fty_mid2, s;
12: PetscInt i, max_it;
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, &minlambda, &maxlambda, &rtol, &atol, <ol, &max_it));
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: /* evaluate initial point */
29: if (linesearch->ops->vidirderiv) {
30: PetscCall((*linesearch->ops->vidirderiv)(snes, F, X, Y, &fty_old));
31: } else {
32: PetscCall(VecDot(F, Y, &fty_old));
33: }
34: if (PetscAbsScalar(fty_old) < atol * ynorm) {
35: if (monitor) {
36: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
37: 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)));
38: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
39: }
40: PetscCall(SNESSetConvergedReason(linesearch->snes, SNES_CONVERGED_FNORM_ABS));
41: PetscFunctionReturn(PETSC_SUCCESS);
42: }
43: fty_init = fty_old;
45: for (i = 0; i < max_it; 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: /* compute change of lambda */
57: delLambda = lambda - lambda_old;
59: /* check change of lambda tolerance */
60: if (PetscAbsReal(delLambda) < ltol) {
61: if (monitor) {
62: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
63: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(dlambda) = %g < ltol = %g\n", (double)PetscAbsReal(delLambda), (double)ltol));
64: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
65: }
66: break;
67: }
69: /* check relative tolerance */
70: if (PetscAbsScalar(fty) / PetscAbsScalar(fty_init) < rtol) {
71: if (monitor) {
72: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
73: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(fty/fty_init) = %g <= rtol = %g\n", (double)(PetscAbsScalar(fty) / PetscAbsScalar(fty_init)), (double)rtol));
74: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
75: }
76: break;
77: }
79: /* check absolute tolerance */
80: if (PetscAbsScalar(fty) < atol * ynorm && i > 0) {
81: if (monitor) {
82: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
83: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(fty)/||y|| = %g <= atol = %g\n", (double)(PetscAbsScalar(fty) / ynorm), (double)atol));
84: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
85: }
86: break;
87: }
89: /* print iteration information */
90: if (monitor) {
91: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
92: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: lambdas = [%g, %g], ftys = [%g, %g]\n", (double)lambda, (double)lambda_old, (double)PetscRealPart(fty), (double)PetscRealPart(fty_old)));
93: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
94: }
96: /* approximate the second derivative */
97: if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) {
98: /* first order finite difference approximation */
99: s = (fty - fty_old) / delLambda;
100: } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
101: /* evaluate function at midpoint 0.5 * (lambda + lambda_old) */
102: PetscCall(VecWAXPY(W, -0.5 * (lambda + lambda_old), Y, X));
103: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
104: PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
105: if (linesearch->ops->vidirderiv) {
106: PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty_mid1));
107: } else {
108: PetscCall(VecDot(F, Y, &fty_mid1));
109: }
111: /* second order finite difference approximation */
112: s = (3. * fty - 4. * fty_mid1 + fty_old) / delLambda;
114: } else {
115: /* evaluate function at midpoint 0.5 * (lambda + lambda_old) */
116: PetscCall(VecWAXPY(W, -0.5 * (lambda + lambda_old), Y, X));
117: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
118: PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
119: if (linesearch->ops->vidirderiv) {
120: PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty_mid1));
121: } else {
122: PetscCall(VecDot(F, Y, &fty_mid1));
123: }
125: /* evaluate function at lambda + 0.5 * (lambda - lambda_old) */
126: PetscCall(VecWAXPY(W, -(lambda + 0.5 * (lambda - lambda_old)), Y, X));
127: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
128: PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
129: if (linesearch->ops->vidirderiv) {
130: PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty_mid2));
131: } else {
132: PetscCall(VecDot(F, Y, &fty_mid2));
133: }
135: /* third order finite difference approximation */
136: s = (2. * fty_mid2 + 3. * fty - 6. * fty_mid1 + fty_old) / (3. * delLambda);
137: }
139: /* compute secant update (modifying the search direction if necessary) */
140: if (PetscRealPart(s) > 0.) s = -s;
141: if (s == 0.0) {
142: if (monitor) {
143: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
144: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: del2Fnrm = 0\n"));
145: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
146: }
147: break;
148: }
149: lambda_update = lambda - PetscRealPart(fty / s);
151: /* if secant method would step out of bounds, exit with the respective bound */
152: if (lambda_update < minlambda) {
153: lambda_update = minlambda;
154: break;
155: }
156: if (lambda_update > maxlambda) {
157: lambda_update = maxlambda;
158: break;
159: }
161: /* if lambda is NaN or Inf, do not accept update but exit with previous lambda */
162: if (PetscIsInfOrNanReal(lambda_update)) {
163: if (monitor) {
164: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
165: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: lambda_update is NaN or Inf\n"));
166: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
167: }
168: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
169: break;
170: }
172: /* compute the new state of the line search */
173: lambda_old = lambda;
174: lambda = lambda_update;
175: fty_old = fty;
176: }
178: /* construct the solution */
179: PetscCall(VecWAXPY(W, -lambda, Y, X));
180: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
181: /* postcheck */
182: PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
183: PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
184: if (changed_y) {
185: if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X));
186: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
187: }
188: PetscCall(VecCopy(W, X));
189: PetscCall((*linesearch->ops->snesfunc)(snes, X, F));
191: PetscCall(SNESLineSearchComputeNorms(linesearch));
192: PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm));
194: if (monitor) {
195: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
196: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm));
197: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
198: }
199: PetscFunctionReturn(PETSC_SUCCESS);
200: }
202: /*MC
203: SNESLINESEARCHCP - Critical point line search. This line search assumes that there exists some
204: artificial $G(x)$ for which the `SNESFunctionFn` $F(x) = grad G(x)$. Therefore, this line search seeks
205: to find roots of the directional derivative via a secant method, that is $F(x_k - \lambda Y_k) \cdot Y_k / ||Y_k|| = 0$.
207: Options Database Keys:
208: + -snes_linesearch_minlambda <1e\-12> - the minimum acceptable `lambda` (scaling of solution update)
209: . -snes_linesearch_maxlambda <1.0> - the algorithm ensures that `lambda` is never larger than this value
210: . -snes_linesearch_damping <1.0> - initial `lambda` on entry to the line search
211: . -snes_linesearch_order <1> - order of the approximation in the secant method, must be 1, 2, or 3
212: . -snes_linesearch_max_it <1> - the maximum number of secant iterations performed
213: . -snes_linesearch_rtol <1e\-8> - relative tolerance for the directional derivative
214: . -snes_linesearch_atol <1e\-15> - absolute tolerance for the directional derivative
215: - -snes_linesearch_ltol <1e\-8> - minimum absolute change in `lambda` allowed
217: Level: advanced
219: Notes:
220: This method does NOT use the objective function if it is provided with `SNESSetObjective()`.
222: This method is the preferred line search for `SNESQN` and `SNESNCG`.
224: .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`, `SNESLINESEARCHBISECTION`
225: M*/
226: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_CP(SNESLineSearch linesearch)
227: {
228: PetscFunctionBegin;
229: linesearch->ops->apply = SNESLineSearchApply_CP;
230: linesearch->ops->destroy = NULL;
231: linesearch->ops->setfromoptions = NULL;
232: linesearch->ops->reset = NULL;
233: linesearch->ops->view = NULL;
234: linesearch->ops->setup = NULL;
235: linesearch->order = SNES_LINESEARCH_ORDER_LINEAR;
237: linesearch->max_it = 1;
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }