Actual source code: linesearchbisection.c
1: #include <petsc/private/linesearchimpl.h>
2: #include <petsc/private/snesimpl.h>
4: static PetscErrorCode SNESLineSearchApply_Bisection(SNESLineSearch linesearch)
5: {
6: PetscBool changed_y, changed_w;
7: Vec X, F, Y, W, G;
8: SNES snes;
9: PetscReal ynorm;
10: PetscReal lambda_left, lambda, lambda_right, lambda_old, fnorm;
11: PetscScalar fty_left, fty, fty_initial;
12: PetscViewer monitor;
13: PetscReal rtol, atol, ltol;
14: PetscInt it, max_it;
16: PetscFunctionBegin;
17: PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G));
18: PetscCall(SNESLineSearchGetLambda(linesearch, &lambda));
19: PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
20: PetscCall(SNESLineSearchGetTolerances(linesearch, NULL, NULL, &rtol, &atol, <ol, &max_it));
21: PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor));
22: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
24: /* pre-check */
25: PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y));
27: /* compute ynorm to normalize search direction */
28: PetscCall(VecNorm(Y, NORM_2, &ynorm));
30: /* initialize interval for bisection */
31: lambda_left = 0.0;
32: lambda_right = lambda;
34: /* compute fty at left end of interval */
35: if (linesearch->ops->vidirderiv) {
36: PetscCall((*linesearch->ops->vidirderiv)(snes, F, X, Y, &fty_left));
37: } else {
38: PetscCall(VecDot(F, Y, &fty_left));
39: }
40: fty_initial = fty_left;
42: /* compute fty at right end of interval (initial lambda) */
43: PetscCall(VecWAXPY(W, -lambda, Y, X));
44: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
45: PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
46: if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
47: PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations during line search!\n"));
48: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
49: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
52: if (linesearch->ops->vidirderiv) {
53: PetscCall((*linesearch->ops->vidirderiv)(snes, G, W, Y, &fty));
54: } else {
55: PetscCall(VecDot(G, Y, &fty));
56: }
57: /* check whether sign changes in interval */
58: if (!PetscIsInfOrNanScalar(fty) && (PetscRealPart(fty_left * fty) > 0.0)) {
59: /* no change of sign: accept full step */
60: if (monitor) {
61: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
62: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: sign of fty does not change in step interval, accepting full step\n"));
63: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
64: }
65: } else {
66: /* change of sign: iteratively bisect interval */
67: lambda_old = 0.0;
68: it = 0;
70: while (PETSC_TRUE) {
71: /* check for infinity or NaN */
72: if (PetscIsInfOrNanScalar(fty)) {
73: if (monitor) {
74: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
75: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search fty is infinity or NaN!\n"));
76: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
77: }
78: PetscCheck(!snes->errorifnotconverged, PetscObjectComm((PetscObject)snes), PETSC_ERR_CONV_FAILED, "infinity or NaN in function evaluation");
79: if (snes->functiondomainerror) {
80: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION_DOMAIN));
81: snes->functiondomainerror = PETSC_FALSE;
82: } else PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: break;
85: }
87: /* check absolute tolerance */
88: if (PetscAbsScalar(fty) <= atol * ynorm) {
89: if (monitor) {
90: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
91: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(fty)/||y|| = %g <= atol = %g\n", (double)(PetscAbsScalar(fty) / ynorm), (double)atol));
92: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
93: }
94: break;
95: }
97: /* check relative tolerance */
98: if (PetscAbsScalar(fty) / PetscAbsScalar(fty_initial) <= rtol) {
99: if (monitor) {
100: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
101: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(fty/fty_initial) = %g <= rtol = %g\n", (double)(PetscAbsScalar(fty) / PetscAbsScalar(fty_initial)), (double)rtol));
102: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
103: }
104: break;
105: }
107: /* check maximum number of iterations */
108: if (it > max_it) {
109: if (monitor) {
110: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
111: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: maximum iterations reached\n"));
112: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
113: }
114: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
115: break;
116: }
118: /* check change of lambda tolerance */
119: if (PetscAbsReal(lambda - lambda_old) < ltol) {
120: if (monitor) {
121: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
122: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(dlambda) = %g < ltol = %g\n", (double)PetscAbsReal(lambda - lambda_old), (double)ltol));
123: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
124: }
125: break;
126: }
128: /* determine direction of bisection (not necessary for 0th iteration) */
129: if (it > 0) {
130: if (PetscRealPart(fty * fty_left) <= 0.0) {
131: lambda_right = lambda;
132: } else {
133: lambda_left = lambda;
134: /* also update fty_left for direction check in next iteration */
135: fty_left = fty;
136: }
137: }
139: /* bisect interval */
140: lambda_old = lambda;
141: lambda = 0.5 * (lambda_left + lambda_right);
143: /* compute fty at new lambda */
144: PetscCall(VecWAXPY(W, -lambda, Y, X));
145: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
146: PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
147: if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
148: PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations during line search!\n"));
149: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
150: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
153: if (linesearch->ops->vidirderiv) {
154: PetscCall((*linesearch->ops->vidirderiv)(snes, G, W, Y, &fty));
155: } else {
156: PetscCall(VecDot(G, Y, &fty));
157: }
159: /* print iteration information */
160: if (monitor) {
161: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
162: PetscCall(PetscViewerASCIIPrintf(monitor, " %3" PetscInt_FMT " Line search: fty/||y|| = %g, lambda = %g\n", it, (double)(PetscRealPart(fty) / ynorm), (double)lambda));
163: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
164: }
166: /* count up iteration */
167: it++;
168: }
169: }
171: /* post-check */
172: PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
173: PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
174: if (changed_y) {
175: if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X));
176: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
177: }
179: /* update solution*/
180: PetscCall(VecCopy(W, X));
181: PetscCall((*linesearch->ops->snesfunc)(snes, X, F));
182: PetscCall(SNESLineSearchComputeNorms(linesearch));
183: PetscCall(SNESLineSearchGetNorms(linesearch, NULL, &fnorm, NULL));
184: SNESLineSearchCheckFunctionDomainError(snes, linesearch, fnorm);
185: PetscFunctionReturn(PETSC_SUCCESS);
186: }
188: /*MC
189: SNESLINESEARCHBISECTION - Bisection line search.
190: Similar to the critical point line search, `SNESLINESEARCHCP`, the bisection line search assumes that there exists some $G(x)$ for which the `SNESFunctionFn` $F(x) = grad G(x)$.
191: This line search seeks to find the root of the directional derivative, that is $F(x_k - \lambda Y_k) \cdot Y_k / ||Y_k|| = 0$, along the search direction $Y_k$ through bisection.
193: Options Database Keys:
194: + -snes_linesearch_max_it <50> - maximum number of bisection iterations for the line search
195: . -snes_linesearch_damping <1.0> - initial `lambda` on entry to the line search
196: . -snes_linesearch_rtol <1e\-8> - relative tolerance for the directional derivative
197: . -snes_linesearch_atol <1e\-6> - absolute tolerance for the directional derivative
198: - -snes_linesearch_ltol <1e\-6> - minimum absolute change in `lambda` allowed (this is an alternative to setting a maximum number of iterations)
200: Level: intermediate
202: Notes:
203: `lambda` is the scaling of the search direction (vector) that is computed by this algorithm.
204: If there is no change of sign in the directional derivative from $\lambda=0$ to the initial `lambda` (the damping), then the initial `lambda` will be used.
205: Hence, this line search will always give a `lambda` in the interval $[0, damping]$.
206: This method does NOT use the objective function if it is provided with `SNESSetObjective()`.
208: .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`, `SNESLINESEARCHCP`
209: M*/
210: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_Bisection(SNESLineSearch linesearch)
211: {
212: PetscFunctionBegin;
213: linesearch->ops->apply = SNESLineSearchApply_Bisection;
214: linesearch->ops->destroy = NULL;
215: linesearch->ops->setfromoptions = NULL;
216: linesearch->ops->reset = NULL;
217: linesearch->ops->view = NULL;
218: linesearch->ops->setup = NULL;
220: /* set default option values */
221: linesearch->max_it = 50;
222: linesearch->rtol = 1e-8;
223: linesearch->atol = 1e-6;
224: linesearch->ltol = 1e-6;
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }