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 23 /* pre-check */ 24 PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y)); 25 26 /* compute ynorm to normalize search direction */ 27 PetscCall(VecNorm(Y, NORM_2, &ynorm)); 28 29 /* initialize interval for bisection */ 30 lambda_left = 0.0; 31 lambda_right = lambda; 32 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; 40 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; 68 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 } 85 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 } 95 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 } 105 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 } 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 186 /* finalization */ 187 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 188 PetscFunctionReturn(PETSC_SUCCESS); 189 } 190 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. 195 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 202 203 Level: intermediate 204 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()`. 210 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; 222 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 } 230