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;
11: PetscReal fty_left, fty, fty_initial;
12: PetscViewer monitor;
13: PetscReal rtol, atol, ltol;
14: PetscInt it, max_its;
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_its));
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: PetscCall(VecDotRealPart(F, Y, &fty_left));
35: fty_initial = fty_left;
37: /* compute fty at right end of interval (initial lambda) */
38: PetscCall(VecWAXPY(W, -lambda, Y, X));
39: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
40: PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
41: if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
42: PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations during line search!\n"));
43: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
44: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
45: PetscFunctionReturn(PETSC_SUCCESS);
46: }
47: PetscCall(VecDotRealPart(G, Y, &fty));
49: /* check whether sign changes in interval */
50: if (!PetscIsInfOrNanReal(fty) && (fty_left * fty) > 0.0) {
51: /* no change of sign: accept full step */
52: if (monitor) {
53: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
54: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: sign of fty does not change in step intervall, accepting full step\n"));
55: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
56: }
57: } else {
58: /* change of sign: iteratively bisect interval */
59: lambda_old = 0.0;
60: it = 0;
62: while (PETSC_TRUE) {
63: /* check for NaN or Inf */
64: if (PetscIsInfOrNanReal(fty)) {
65: if (monitor) {
66: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
67: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search fty is NaN or Inf!\n"));
68: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
69: }
70: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
71: PetscFunctionReturn(PETSC_SUCCESS);
72: break;
73: }
75: /* check absolute tolerance */
76: if (PetscAbsReal(fty) <= atol * ynorm) {
77: if (monitor) {
78: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
79: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(fty)/||y|| = %g <= atol = %g\n", (double)PetscAbsReal(fty / ynorm), (double)(atol)));
80: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
81: }
82: break;
83: }
85: /* check relative tolerance */
86: if (PetscAbsReal(fty / fty_initial) <= rtol) {
87: if (monitor) {
88: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
89: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(fty/fty_initial) = %g <= rtol = %g\n", (double)PetscAbsReal(fty / fty_initial), (double)rtol));
90: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
91: }
92: break;
93: }
95: /* check maximum number of iterations */
96: if (it > max_its) {
97: if (monitor) {
98: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
99: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: maximum iterations reached\n"));
100: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
101: }
102: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
103: PetscFunctionReturn(PETSC_SUCCESS);
104: break;
105: }
107: /* check change of lambda tolerance */
108: if (PetscAbsReal(lambda - lambda_old) < ltol) {
109: if (monitor) {
110: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
111: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(dlambda) = %g < ltol = %g\n", (double)PetscAbsReal(lambda - lambda_old), (double)ltol));
112: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
113: }
114: break;
115: }
117: /* determine direction of bisection (not necessary for 0th iteration) */
118: if (it > 0) {
119: if (fty * fty_left <= 0.0) {
120: lambda_right = lambda;
121: } else {
122: lambda_left = lambda;
123: /* also update fty_left for direction check in next iteration */
124: fty_left = fty;
125: }
126: }
128: /* bisect interval */
129: lambda_old = lambda;
130: lambda = 0.5 * (lambda_left + lambda_right);
132: /* compute fty at new lambda */
133: PetscCall(VecWAXPY(W, -lambda, Y, X));
134: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
135: PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
136: if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
137: PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations during line search!\n"));
138: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
139: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
142: PetscCall(VecDotRealPart(G, Y, &fty));
144: /* print iteration information */
145: if (monitor) {
146: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
147: PetscCall(PetscViewerASCIIPrintf(monitor, " %3" PetscInt_FMT " Line search: fty/||y|| = %g, lambda = %g\n", it, (double)(fty / ynorm), (double)lambda));
148: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
149: }
151: /* count up iteration */
152: it++;
153: }
154: }
156: /* post-check */
157: PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
158: PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
159: if (changed_y) {
160: if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X));
161: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
162: }
164: /* update solution*/
165: PetscCall(VecCopy(W, X));
166: PetscCall((*linesearch->ops->snesfunc)(snes, X, F));
167: PetscCall(SNESLineSearchComputeNorms(linesearch));
169: /* finalization */
170: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: /*MC
175: SNESLINESEARCHBISECTION - Bisection line search.
176: Similar to the critical point line search, `SNESLINESEARCHCP`, the bisection line search assumes that there exists some $G(x)$ for which the `SNESFunction` $F(x) = grad G(x)$.
177: This line search seeks to find the root of the directional derivative along the search direction $F^T Y$ through bisection.
179: Options Database Keys:
180: + -snes_linesearch_max_it <50> - maximum number of iterations for the line search
181: . -snes_linesearch_damping <1.0> - initial trial step length on entry to the line search
182: - -snes_linesearch_rtol <1e\-8> - relative tolerance for the directional derivative
183: - -snes_linesearch_atol <1e\-6> - absolute tolerance for the directional derivative
184: . -snes_linesearch_ltol <1e\-6> - minimum absolute change in lambda allowed
186: Level: intermediate
188: Note:
189: This method does NOT use the objective function if it is provided with `SNESSetObjective()`.
190: This line search will always give a step size in the interval [0, damping].
192: .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`, `SNESLINESEARCHCP`
193: M*/
194: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_Bisection(SNESLineSearch linesearch)
195: {
196: PetscFunctionBegin;
197: linesearch->ops->apply = SNESLineSearchApply_Bisection;
198: linesearch->ops->destroy = NULL;
199: linesearch->ops->setfromoptions = NULL;
200: linesearch->ops->reset = NULL;
201: linesearch->ops->view = NULL;
202: linesearch->ops->setup = NULL;
204: /* set default option values */
205: linesearch->max_its = 50;
206: linesearch->rtol = 1e-8;
207: linesearch->atol = 1e-6;
208: linesearch->ltol = 1e-6;
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }