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; 11 PetscScalar fty_left, fty, fty_initial; 12 PetscViewer monitor; 13 PetscReal rtol, atol, ltol; 14 PetscInt it, max_its; 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_its)); 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 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 } 81 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 } 91 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 } 101 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 } 113 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 } 123 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 } 134 135 /* bisect interval */ 136 lambda_old = lambda; 137 lambda = 0.5 * (lambda_left + lambda_right); 138 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 } 154 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 } 161 162 /* count up iteration */ 163 it++; 164 } 165 } 166 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 } 174 175 /* update solution*/ 176 PetscCall(VecCopy(W, X)); 177 PetscCall((*linesearch->ops->snesfunc)(snes, X, F)); 178 PetscCall(SNESLineSearchComputeNorms(linesearch)); 179 180 /* finalization */ 181 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 182 PetscFunctionReturn(PETSC_SUCCESS); 183 } 184 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. 189 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 196 197 Level: intermediate 198 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]. 202 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; 214 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 } 222