1 #include <petsc/private/linesearchimpl.h> 2 #include <petsc/private/snesimpl.h> 3 #include <petscsnes.h> 4 5 static PetscErrorCode SNESLineSearchApply_Secant(SNESLineSearch linesearch) 6 { 7 PetscBool changed_y, changed_w; 8 Vec X; 9 Vec F; 10 Vec Y; 11 Vec W; 12 SNES snes; 13 PetscReal gnorm; 14 PetscReal ynorm; 15 PetscReal xnorm; 16 PetscReal minlambda, maxlambda, atol, ltol; 17 PetscViewer monitor; 18 PetscReal lambda, lambda_old, lambda_mid, lambda_update, delLambda; 19 PetscReal fnrm, fnrm_old, fnrm_mid; 20 PetscReal delFnrm, delFnrm_old, del2Fnrm; 21 PetscInt i, max_it; 22 SNESObjectiveFn *objective; 23 24 PetscFunctionBegin; 25 PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL)); 26 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm)); 27 PetscCall(SNESLineSearchGetLambda(linesearch, &lambda)); 28 PetscCall(SNESLineSearchGetSNES(linesearch, &snes)); 29 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 30 PetscCall(SNESLineSearchGetTolerances(linesearch, &minlambda, &maxlambda, NULL, &atol, <ol, &max_it)); 31 PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor)); 32 33 PetscCall(SNESGetObjective(snes, &objective, NULL)); 34 35 /* precheck */ 36 PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y)); 37 lambda_old = 0.0; 38 if (!objective) { 39 fnrm_old = gnorm * gnorm; 40 } else { 41 PetscCall(SNESComputeObjective(snes, X, &fnrm_old)); 42 SNESLineSearchCheckObjectiveDomainError(snes, fnrm_old); 43 } 44 lambda_mid = 0.5 * (lambda + lambda_old); 45 46 for (i = 0; i < max_it; i++) { 47 /* repeatedly cut lambda until the norm or objective function is not infinity or NaN or lambda is too small */ 48 while (PETSC_TRUE) { 49 PetscCall(VecWAXPY(W, -lambda_mid, Y, X)); 50 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 51 if (!objective) { 52 /* compute the norm at the midpoint */ 53 PetscCall((*linesearch->ops->snesfunc)(snes, W, F)); 54 if (linesearch->ops->vinorm) { 55 fnrm_mid = gnorm; 56 PetscCall((*linesearch->ops->vinorm)(snes, F, W, &fnrm_mid)); 57 } else { 58 PetscCall(VecNorm(F, NORM_2, &fnrm_mid)); 59 } 60 61 /* compute the norm at the new endpoint */ 62 PetscCall(VecWAXPY(W, -lambda, Y, X)); 63 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 64 PetscCall((*linesearch->ops->snesfunc)(snes, W, F)); 65 if (linesearch->ops->vinorm) { 66 fnrm = gnorm; 67 PetscCall((*linesearch->ops->vinorm)(snes, F, W, &fnrm)); 68 } else { 69 PetscCall(VecNorm(F, NORM_2, &fnrm)); 70 } 71 fnrm_mid = fnrm_mid * fnrm_mid; 72 fnrm = fnrm * fnrm; 73 } else { 74 /* compute the objective at the midpoint */ 75 PetscCall(SNESComputeObjective(snes, W, &fnrm_mid)); 76 77 /* compute the objective at the new endpoint */ 78 PetscCall(VecWAXPY(W, -lambda, Y, X)); 79 PetscCall(SNESComputeObjective(snes, W, &fnrm)); 80 } 81 82 /* if new endpoint is viable, exit */ 83 if (!PetscIsInfOrNanReal(fnrm)) break; 84 if (monitor) { 85 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 86 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: objective function at lambda = %g is infinity or NaN, cutting lambda\n", (double)lambda)); 87 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 88 } 89 90 /* if smallest allowable lambda gives NaN or Inf, exit line search */ 91 if (lambda <= minlambda) { 92 if (monitor) { 93 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 94 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: objective function at lambda = %g <= lambda_min = %g is infinity or NaN, can not further cut lambda\n", (double)lambda, (double)lambda)); 95 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 96 } 97 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 98 PetscFunctionReturn(PETSC_SUCCESS); 99 } 100 101 /* forbid the search from ever going back to the "failed" length that generates infinity or NaN */ 102 maxlambda = .95 * lambda; 103 104 /* shrink lambda towards the previous one which was viable */ 105 lambda = .5 * (lambda + lambda_old); 106 lambda_mid = .5 * (lambda + lambda_old); 107 } 108 109 /* monitoring output */ 110 if (monitor) { 111 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 112 if (!objective) { 113 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: lambdas = [%g, %g, %g], fnorms = [%g, %g, %g]\n", (double)lambda, (double)lambda_mid, (double)lambda_old, (double)PetscSqrtReal(fnrm), (double)PetscSqrtReal(fnrm_mid), (double)PetscSqrtReal(fnrm_old))); 114 } else { 115 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: lambdas = [%g, %g, %g], obj = [%g, %g, %g]\n", (double)lambda, (double)lambda_mid, (double)lambda_old, (double)fnrm, (double)fnrm_mid, (double)fnrm_old)); 116 } 117 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 118 } 119 120 /* determine change of lambda */ 121 delLambda = lambda - lambda_old; 122 123 /* check change of lambda tolerance */ 124 if (PetscAbsReal(delLambda) < ltol) { 125 if (monitor) { 126 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 127 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(delLambda) = %g < ltol = %g\n", (double)PetscAbsReal(delLambda), (double)ltol)); 128 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 129 } 130 break; 131 } 132 133 /* compute f'() at the end points using second order one sided differencing */ 134 delFnrm = (3. * fnrm - 4. * fnrm_mid + 1. * fnrm_old) / delLambda; 135 delFnrm_old = (-3. * fnrm_old + 4. * fnrm_mid - 1. * fnrm) / delLambda; 136 137 /* compute f''() at the midpoint using centered differencing */ 138 del2Fnrm = (delFnrm - delFnrm_old) / delLambda; 139 140 /* check absolute tolerance */ 141 if (PetscAbsReal(delFnrm) <= atol) { 142 if (monitor) { 143 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 144 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(delFnrm) = %g <= atol = %g\n", (double)PetscAbsReal(delFnrm), (double)atol)); 145 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 146 } 147 break; 148 } 149 150 /* compute the secant (Newton) update -- always go downhill */ 151 if (del2Fnrm > 0.) lambda_update = lambda - delFnrm / del2Fnrm; 152 else if (del2Fnrm < 0.) lambda_update = lambda + delFnrm / del2Fnrm; 153 else { 154 if (monitor) { 155 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 156 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: del2Fnrm = 0\n")); 157 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 158 } 159 break; 160 } 161 162 /* prevent secant method from stepping out of bounds */ 163 if (lambda_update < minlambda) lambda_update = 0.5 * (lambda + lambda_old); 164 if (lambda_update > maxlambda) { 165 lambda_update = maxlambda; 166 break; 167 } 168 169 /* if lambda is NaN or Inf, do not accept update but exit with previous lambda */ 170 if (PetscIsInfOrNanReal(lambda_update)) { 171 if (monitor) { 172 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 173 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: lambda_update is NaN or Inf\n")); 174 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 175 } 176 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 177 break; 178 } 179 180 /* update the endpoints and the midpoint of the bracketed secant region */ 181 lambda_old = lambda; 182 lambda = lambda_update; 183 fnrm_old = fnrm; 184 lambda_mid = 0.5 * (lambda + lambda_old); 185 186 if ((i == max_it - 1) && monitor) { 187 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 188 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: reached maximum number of iterations!\n")); 189 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 190 } 191 } 192 193 /* construct the solution */ 194 PetscCall(VecWAXPY(W, -lambda, Y, X)); 195 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 196 197 /* postcheck */ 198 PetscCall(SNESLineSearchSetLambda(linesearch, lambda)); 199 PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w)); 200 if (changed_y) { 201 if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X)); 202 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 203 } 204 PetscCall(VecCopy(W, X)); 205 PetscCall((*linesearch->ops->snesfunc)(snes, X, F)); 206 207 PetscCall(SNESLineSearchComputeNorms(linesearch)); 208 PetscCall(SNESLineSearchGetNorms(linesearch, NULL, &gnorm, NULL)); 209 SNESLineSearchCheckFunctionDomainError(snes, linesearch, gnorm); 210 if (monitor) { 211 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 212 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search terminated: lambda = %g, fnorm = %g\n", (double)lambda, (double)gnorm)); 213 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 214 } 215 PetscFunctionReturn(PETSC_SUCCESS); 216 } 217 218 /*MC 219 SNESLINESEARCHSECANT - Secant search in the L2 norm of the function or the objective function. 220 221 Attempts to solve $ \min_{\lambda} f(x_k + \lambda Y_k) $ using the secant method with the initial bracketing of $ \lambda $ between [0,damping]. 222 $f(x_k + \lambda Y_k)$ is either the squared L2-norm of the function $||F(x_k + \lambda Y_k)||_2^2$, 223 or the objective function $G(x_k + \lambda Y_k)$ if it is provided with `SNESSetObjective()` 224 Differences of $f()$ are used to approximate the first and second derivative of $f()$ with respect to 225 $\lambda$, $f'()$ and $f''()$. 226 227 When an objective function is provided $f(w)$ is the objective function otherwise $f(w) = ||F(w)||^2$. 228 $x$ is the current step and $y$ is the search direction. 229 230 Options Database Keys: 231 + -snes_linesearch_max_it <1> - maximum number of iterations within the line search 232 . -snes_linesearch_damping <1.0> - initial `lambda` on entry to the line search 233 . -snes_linesearch_minlambda <1e\-12> - minimum allowable `lambda` 234 . -snes_linesearch_maxlambda <1.0> - maximum `lambda` (scaling of solution update) allowed 235 . -snes_linesearch_atol <1e\-15> - absolute tolerance for the secant method $ f'() < atol $ 236 - -snes_linesearch_ltol <1e\-8> - minimum absolute change in `lambda` allowed 237 238 Level: advanced 239 240 .seealso: [](ch_snes), `SNESLINESEARCHBT`, `SNESLINESEARCHCP`, `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()` 241 M*/ 242 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_Secant(SNESLineSearch linesearch) 243 { 244 PetscFunctionBegin; 245 linesearch->ops->apply = SNESLineSearchApply_Secant; 246 linesearch->ops->destroy = NULL; 247 linesearch->ops->setfromoptions = NULL; 248 linesearch->ops->reset = NULL; 249 linesearch->ops->view = NULL; 250 linesearch->ops->setup = NULL; 251 252 linesearch->max_it = 1; 253 PetscFunctionReturn(PETSC_SUCCESS); 254 } 255