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));
23: /* pre-check */
24: PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y));
26: /* compute ynorm to normalize search direction */
27: PetscCall(VecNorm(Y, NORM_2, &ynorm));
29: /* initialize interval for bisection */
30: lambda_left = 0.0;
31: lambda_right = lambda;
33: /* compute fty at left end of interval */
34: if (linesearch->ops->vidirderiv) {
35: PetscCall((*linesearch->ops->vidirderiv)(snes, F, X, Y, &fty_left));
36: } else {
37: PetscCall(VecDot(F, Y, &fty_left));
38: }
39: fty_initial = fty_left;
41: /* compute fty at right end of interval (initial lambda) */
42: PetscCall(VecWAXPY(W, -lambda, Y, X));
43: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
44: PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
45: if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
46: PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations during line search!\n"));
47: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
48: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
49: PetscFunctionReturn(PETSC_SUCCESS);
50: }
51: if (linesearch->ops->vidirderiv) {
52: PetscCall((*linesearch->ops->vidirderiv)(snes, G, W, Y, &fty));
53: } else {
54: PetscCall(VecDot(G, Y, &fty));
55: }
56: /* check whether sign changes in interval */
57: if (!PetscIsInfOrNanScalar(fty) && (PetscRealPart(fty_left * fty) > 0.0)) {
58: /* no change of sign: accept full step */
59: if (monitor) {
60: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
61: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: sign of fty does not change in step interval, accepting full step\n"));
62: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
63: }
64: } else {
65: /* change of sign: iteratively bisect interval */
66: lambda_old = 0.0;
67: it = 0;
69: while (PETSC_TRUE) {
70: /* check for infinity or NaN */
71: if (PetscIsInfOrNanScalar(fty)) {
72: if (monitor) {
73: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
74: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search fty is infinity or NaN!\n"));
75: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
76: }
77: PetscCheck(!snes->errorifnotconverged, PetscObjectComm((PetscObject)snes), PETSC_ERR_CONV_FAILED, "infinity or NaN in function evaluation");
78: if (snes->functiondomainerror) {
79: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION_DOMAIN));
80: snes->functiondomainerror = PETSC_FALSE;
81: } else PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
82: PetscFunctionReturn(PETSC_SUCCESS);
83: break;
84: }
86: /* check absolute tolerance */
87: if (PetscAbsScalar(fty) <= atol * ynorm) {
88: if (monitor) {
89: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
90: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(fty)/||y|| = %g <= atol = %g\n", (double)(PetscAbsScalar(fty) / ynorm), (double)atol));
91: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
92: }
93: break;
94: }
96: /* check relative tolerance */
97: if (PetscAbsScalar(fty) / PetscAbsScalar(fty_initial) <= rtol) {
98: if (monitor) {
99: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
100: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(fty/fty_initial) = %g <= rtol = %g\n", (double)(PetscAbsScalar(fty) / PetscAbsScalar(fty_initial)), (double)rtol));
101: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
102: }
103: break;
104: }
106: /* check maximum number of iterations */
107: if (it > max_it) {
108: if (monitor) {
109: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
110: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: maximum iterations reached\n"));
111: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
112: }
113: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
114: PetscFunctionReturn(PETSC_SUCCESS);
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);
186: /* finalization */
187: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: /*MC
192: SNESLINESEARCHBISECTION - Bisection line search.
193: 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)$.
194: 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.
196: Options Database Keys:
197: + -snes_linesearch_max_it <50> - maximum number of bisection iterations for the line search
198: . -snes_linesearch_damping <1.0> - initial `lambda` on entry to the line search
199: . -snes_linesearch_rtol <1e\-8> - relative tolerance for the directional derivative
200: . -snes_linesearch_atol <1e\-6> - absolute tolerance for the directional derivative
201: - -snes_linesearch_ltol <1e\-6> - minimum absolute change in `lambda` allowed
203: Level: intermediate
205: Notes:
206: `lambda` is the scaling of the search direction (vector) that is computed by this algorithm.
207: 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.
208: Hence, this line search will always give a `lambda` in the interval $[0, damping]$.
209: This method does NOT use the objective function if it is provided with `SNESSetObjective()`.
211: .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`, `SNESLINESEARCHCP`
212: M*/
213: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_Bisection(SNESLineSearch linesearch)
214: {
215: PetscFunctionBegin;
216: linesearch->ops->apply = SNESLineSearchApply_Bisection;
217: linesearch->ops->destroy = NULL;
218: linesearch->ops->setfromoptions = NULL;
219: linesearch->ops->reset = NULL;
220: linesearch->ops->view = NULL;
221: linesearch->ops->setup = NULL;
223: /* set default option values */
224: linesearch->max_it = 50;
225: linesearch->rtol = 1e-8;
226: linesearch->atol = 1e-6;
227: linesearch->ltol = 1e-6;
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }