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 PetscReal 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 PetscCall(VecDotRealPart(F, Y, &fty_left)); 35 fty_initial = fty_left; 36 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)); 48 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; 61 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 } 74 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 } 84 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 } 94 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 } 106 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 } 116 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 } 127 128 /* bisect interval */ 129 lambda_old = lambda; 130 lambda = 0.5 * (lambda_left + lambda_right); 131 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)); 143 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 } 150 151 /* count up iteration */ 152 it++; 153 } 154 } 155 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 } 163 164 /* update solution*/ 165 PetscCall(VecCopy(W, X)); 166 PetscCall((*linesearch->ops->snesfunc)(snes, X, F)); 167 PetscCall(SNESLineSearchComputeNorms(linesearch)); 168 169 /* finalization */ 170 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 171 PetscFunctionReturn(PETSC_SUCCESS); 172 } 173 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. 178 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 185 186 Level: intermediate 187 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]. 191 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; 203 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 } 211