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: PetscScalar 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: 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 NaN or Inf */
71: if (PetscIsInfOrNanScalar(fty)) {
72: if (monitor) {
73: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
74: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search fty is NaN or Inf!\n"));
75: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
76: }
77: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
78: PetscFunctionReturn(PETSC_SUCCESS);
79: break;
80: }
82: /* check absolute tolerance */
83: if (PetscAbsScalar(fty) <= atol * ynorm) {
84: if (monitor) {
85: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
86: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(fty)/||y|| = %g <= atol = %g\n", (double)(PetscAbsScalar(fty) / ynorm), (double)atol));
87: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
88: }
89: break;
90: }
92: /* check relative tolerance */
93: if (PetscAbsScalar(fty) / PetscAbsScalar(fty_initial) <= rtol) {
94: if (monitor) {
95: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
96: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(fty/fty_initial) = %g <= rtol = %g\n", (double)(PetscAbsScalar(fty) / PetscAbsScalar(fty_initial)), (double)rtol));
97: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
98: }
99: break;
100: }
102: /* check maximum number of iterations */
103: if (it > max_its) {
104: if (monitor) {
105: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
106: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: maximum iterations reached\n"));
107: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
108: }
109: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
110: PetscFunctionReturn(PETSC_SUCCESS);
111: break;
112: }
114: /* check change of lambda tolerance */
115: if (PetscAbsReal(lambda - lambda_old) < ltol) {
116: if (monitor) {
117: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
118: PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(dlambda) = %g < ltol = %g\n", (double)PetscAbsReal(lambda - lambda_old), (double)ltol));
119: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
120: }
121: break;
122: }
124: /* determine direction of bisection (not necessary for 0th iteration) */
125: if (it > 0) {
126: if (PetscRealPart(fty * fty_left) <= 0.0) {
127: lambda_right = lambda;
128: } else {
129: lambda_left = lambda;
130: /* also update fty_left for direction check in next iteration */
131: fty_left = fty;
132: }
133: }
135: /* bisect interval */
136: lambda_old = lambda;
137: lambda = 0.5 * (lambda_left + lambda_right);
139: /* compute fty at new lambda */
140: PetscCall(VecWAXPY(W, -lambda, Y, X));
141: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
142: PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
143: if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
144: PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations during line search!\n"));
145: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
146: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION));
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
149: if (linesearch->ops->vidirderiv) {
150: PetscCall((*linesearch->ops->vidirderiv)(snes, G, W, Y, &fty));
151: } else {
152: PetscCall(VecDot(G, Y, &fty));
153: }
155: /* print iteration information */
156: if (monitor) {
157: PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
158: PetscCall(PetscViewerASCIIPrintf(monitor, " %3" PetscInt_FMT " Line search: fty/||y|| = %g, lambda = %g\n", it, (double)(PetscRealPart(fty) / ynorm), (double)lambda));
159: PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
160: }
162: /* count up iteration */
163: it++;
164: }
165: }
167: /* post-check */
168: PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
169: PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
170: if (changed_y) {
171: if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X));
172: if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
173: }
175: /* update solution*/
176: PetscCall(VecCopy(W, X));
177: PetscCall((*linesearch->ops->snesfunc)(snes, X, F));
178: PetscCall(SNESLineSearchComputeNorms(linesearch));
180: /* finalization */
181: PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
185: /*MC
186: SNESLINESEARCHBISECTION - Bisection line search.
187: 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)$.
188: This line search seeks to find the root of the directional derivative along the search direction $F^T Y$ through bisection.
190: Options Database Keys:
191: + -snes_linesearch_max_it <50> - maximum number of iterations for the line search
192: . -snes_linesearch_damping <1.0> - initial trial step length on entry to the line search
193: . -snes_linesearch_rtol <1e\-8> - relative tolerance for the directional derivative
194: . -snes_linesearch_atol <1e\-6> - absolute tolerance for the directional derivative
195: - -snes_linesearch_ltol <1e\-6> - minimum absolute change in lambda allowed
197: Level: intermediate
199: Note:
200: This method does NOT use the objective function if it is provided with `SNESSetObjective()`.
201: This line search will always give a step size in the interval [0, damping].
203: .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`, `SNESLINESEARCHCP`
204: M*/
205: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_Bisection(SNESLineSearch linesearch)
206: {
207: PetscFunctionBegin;
208: linesearch->ops->apply = SNESLineSearchApply_Bisection;
209: linesearch->ops->destroy = NULL;
210: linesearch->ops->setfromoptions = NULL;
211: linesearch->ops->reset = NULL;
212: linesearch->ops->view = NULL;
213: linesearch->ops->setup = NULL;
215: /* set default option values */
216: linesearch->max_its = 50;
217: linesearch->rtol = 1e-8;
218: linesearch->atol = 1e-6;
219: linesearch->ltol = 1e-6;
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }