1 #include <petsc/private/linesearchimpl.h> 2 #include <petsc/private/snesimpl.h> 3 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; 15 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)); 23 24 /* pre-check */ 25 PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y)); 26 27 /* compute ynorm to normalize search direction */ 28 PetscCall(VecNorm(Y, NORM_2, &ynorm)); 29 30 /* initialize interval for bisection */ 31 lambda_left = 0.0; 32 lambda_right = lambda; 33 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; 41 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; 69 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 } 86 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 } 96 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 } 106 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 } 117 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 } 127 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 } 138 139 /* bisect interval */ 140 lambda_old = lambda; 141 lambda = 0.5 * (lambda_left + lambda_right); 142 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 } 158 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 } 165 166 /* count up iteration */ 167 it++; 168 } 169 } 170 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 } 178 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 } 187 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. 192 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) 199 200 Level: intermediate 201 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()`. 207 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; 219 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 } 227