1 #include <petsc/private/linesearchimpl.h> 2 #include <petscsnes.h> 3 4 static PetscErrorCode SNESLineSearchApply_CP(SNESLineSearch linesearch) 5 { 6 PetscBool changed_y, changed_w; 7 Vec X, Y, F, W; 8 SNES snes; 9 PetscReal xnorm, ynorm, gnorm, minlambda, maxlambda, rtol, atol, ltol; 10 PetscReal lambda, lambda_old, lambda_update, delLambda; 11 PetscScalar fty, fty_init, fty_old, fty_mid1, fty_mid2, s; 12 PetscInt i, max_it; 13 PetscViewer monitor; 14 15 PetscFunctionBegin; 16 PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL)); 17 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm)); 18 PetscCall(SNESLineSearchGetSNES(linesearch, &snes)); 19 PetscCall(SNESLineSearchGetLambda(linesearch, &lambda)); 20 PetscCall(SNESLineSearchGetTolerances(linesearch, &minlambda, &maxlambda, &rtol, &atol, <ol, &max_it)); 21 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 22 PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor)); 23 24 /* precheck */ 25 PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y)); 26 lambda_old = 0.0; 27 28 /* evaluate initial point */ 29 if (linesearch->ops->vidirderiv) { 30 PetscCall((*linesearch->ops->vidirderiv)(snes, F, X, Y, &fty_old)); 31 } else { 32 PetscCall(VecDot(F, Y, &fty_old)); 33 } 34 if (PetscAbsScalar(fty_old) < atol * ynorm) { 35 if (monitor) { 36 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 37 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search terminated at initial point because dot(F,Y) = %g < atol*||y|| = %g\n", (double)PetscAbsScalar(fty_old), (double)(atol * ynorm))); 38 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 39 } 40 PetscCall(SNESSetConvergedReason(linesearch->snes, SNES_CONVERGED_FNORM_ABS)); 41 PetscFunctionReturn(PETSC_SUCCESS); 42 } 43 fty_init = fty_old; 44 45 for (i = 0; i < max_it; i++) { 46 /* compute the norm at lambda */ 47 PetscCall(VecWAXPY(W, -lambda, Y, X)); 48 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 49 PetscCall((*linesearch->ops->snesfunc)(snes, W, F)); 50 if (linesearch->ops->vidirderiv) { 51 PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty)); 52 } else { 53 PetscCall(VecDot(F, Y, &fty)); 54 } 55 56 /* compute change of lambda */ 57 delLambda = lambda - lambda_old; 58 59 /* check change of lambda tolerance */ 60 if (PetscAbsReal(delLambda) < ltol) { 61 if (monitor) { 62 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 63 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(dlambda) = %g < ltol = %g\n", (double)PetscAbsReal(delLambda), (double)ltol)); 64 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 65 } 66 break; 67 } 68 69 /* check relative tolerance */ 70 if (PetscAbsScalar(fty) / PetscAbsScalar(fty_init) < rtol) { 71 if (monitor) { 72 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 73 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(fty/fty_init) = %g <= rtol = %g\n", (double)(PetscAbsScalar(fty) / PetscAbsScalar(fty_init)), (double)rtol)); 74 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 75 } 76 break; 77 } 78 79 /* check absolute tolerance */ 80 if (PetscAbsScalar(fty) < atol * ynorm && i > 0) { 81 if (monitor) { 82 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 83 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: abs(fty)/||y|| = %g <= atol = %g\n", (double)(PetscAbsScalar(fty) / ynorm), (double)atol)); 84 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 85 } 86 break; 87 } 88 89 /* print iteration information */ 90 if (monitor) { 91 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 92 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: lambdas = [%g, %g], ftys = [%g, %g]\n", (double)lambda, (double)lambda_old, (double)PetscRealPart(fty), (double)PetscRealPart(fty_old))); 93 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 94 } 95 96 /* approximate the second derivative */ 97 if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) { 98 /* first order finite difference approximation */ 99 s = (fty - fty_old) / delLambda; 100 } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 101 /* evaluate function at midpoint 0.5 * (lambda + lambda_old) */ 102 PetscCall(VecWAXPY(W, -0.5 * (lambda + lambda_old), Y, X)); 103 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 104 PetscCall((*linesearch->ops->snesfunc)(snes, W, F)); 105 if (linesearch->ops->vidirderiv) { 106 PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty_mid1)); 107 } else { 108 PetscCall(VecDot(F, Y, &fty_mid1)); 109 } 110 111 /* second order finite difference approximation */ 112 s = (3. * fty - 4. * fty_mid1 + fty_old) / delLambda; 113 114 } else { 115 /* evaluate function at midpoint 0.5 * (lambda + lambda_old) */ 116 PetscCall(VecWAXPY(W, -0.5 * (lambda + lambda_old), Y, X)); 117 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 118 PetscCall((*linesearch->ops->snesfunc)(snes, W, F)); 119 if (linesearch->ops->vidirderiv) { 120 PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty_mid1)); 121 } else { 122 PetscCall(VecDot(F, Y, &fty_mid1)); 123 } 124 125 /* evaluate function at lambda + 0.5 * (lambda - lambda_old) */ 126 PetscCall(VecWAXPY(W, -(lambda + 0.5 * (lambda - lambda_old)), Y, X)); 127 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 128 PetscCall((*linesearch->ops->snesfunc)(snes, W, F)); 129 if (linesearch->ops->vidirderiv) { 130 PetscCall((*linesearch->ops->vidirderiv)(snes, F, W, Y, &fty_mid2)); 131 } else { 132 PetscCall(VecDot(F, Y, &fty_mid2)); 133 } 134 135 /* third order finite difference approximation */ 136 s = (2. * fty_mid2 + 3. * fty - 6. * fty_mid1 + fty_old) / (3. * delLambda); 137 } 138 139 /* compute secant update (modifying the search direction if necessary) */ 140 if (PetscRealPart(s) > 0.) s = -s; 141 if (s == 0.0) { 142 if (monitor) { 143 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 144 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: del2Fnrm = 0\n")); 145 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 146 } 147 break; 148 } 149 lambda_update = lambda - PetscRealPart(fty / s); 150 151 /* if secant method would step out of bounds, exit with the respective bound */ 152 if (lambda_update < minlambda) { 153 lambda_update = minlambda; 154 break; 155 } 156 if (lambda_update > maxlambda) { 157 lambda_update = maxlambda; 158 break; 159 } 160 161 /* if lambda is NaN or Inf, do not accept update but exit with previous lambda */ 162 if (PetscIsInfOrNanReal(lambda_update)) { 163 if (monitor) { 164 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 165 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: lambda_update is NaN or Inf\n")); 166 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 167 } 168 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 169 break; 170 } 171 172 /* compute the new state of the line search */ 173 lambda_old = lambda; 174 lambda = lambda_update; 175 fty_old = fty; 176 } 177 178 /* construct the solution */ 179 PetscCall(VecWAXPY(W, -lambda, Y, X)); 180 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 181 /* postcheck */ 182 PetscCall(SNESLineSearchSetLambda(linesearch, lambda)); 183 PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w)); 184 if (changed_y) { 185 if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X)); 186 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 187 } 188 PetscCall(VecCopy(W, X)); 189 PetscCall((*linesearch->ops->snesfunc)(snes, X, F)); 190 191 PetscCall(SNESLineSearchComputeNorms(linesearch)); 192 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm)); 193 194 if (monitor) { 195 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 196 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm)); 197 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 198 } 199 PetscFunctionReturn(PETSC_SUCCESS); 200 } 201 202 /*MC 203 SNESLINESEARCHCP - Critical point line search. This line search assumes that there exists some 204 artificial $G(x)$ for which the `SNESFunctionFn` $F(x) = grad G(x)$. Therefore, this line search seeks 205 to find roots of the directional derivative via a secant method, that is $F(x_k - \lambda Y_k) \cdot Y_k / ||Y_k|| = 0$. 206 207 Options Database Keys: 208 + -snes_linesearch_minlambda <1e\-12> - the minimum acceptable `lambda` (scaling of solution update) 209 . -snes_linesearch_maxlambda <1.0> - the algorithm ensures that `lambda` is never larger than this value 210 . -snes_linesearch_damping <1.0> - initial `lambda` on entry to the line search 211 . -snes_linesearch_order <1> - order of the approximation in the secant method, must be 1, 2, or 3 212 . -snes_linesearch_max_it <1> - the maximum number of secant iterations performed 213 . -snes_linesearch_rtol <1e\-8> - relative tolerance for the directional derivative 214 . -snes_linesearch_atol <1e\-15> - absolute tolerance for the directional derivative 215 - -snes_linesearch_ltol <1e\-8> - minimum absolute change in `lambda` allowed 216 217 Level: advanced 218 219 Notes: 220 This method does NOT use the objective function if it is provided with `SNESSetObjective()`. 221 222 This method is the preferred line search for `SNESQN` and `SNESNCG`. 223 224 .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`, `SNESLINESEARCHBISECTION` 225 M*/ 226 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_CP(SNESLineSearch linesearch) 227 { 228 PetscFunctionBegin; 229 linesearch->ops->apply = SNESLineSearchApply_CP; 230 linesearch->ops->destroy = NULL; 231 linesearch->ops->setfromoptions = NULL; 232 linesearch->ops->reset = NULL; 233 linesearch->ops->view = NULL; 234 linesearch->ops->setup = NULL; 235 linesearch->order = SNES_LINESEARCH_ORDER_LINEAR; 236 237 linesearch->max_it = 1; 238 PetscFunctionReturn(PETSC_SUCCESS); 239 } 240