1 #include <petsc/private/linesearchimpl.h> /*I "petscsnes.h" I*/ 2 #include <petsc/private/snesimpl.h> 3 4 typedef struct { 5 PetscReal norm_delta_x_prev; /* norm of previous update */ 6 PetscReal norm_bar_delta_x_prev; /* norm of previous bar update */ 7 PetscReal mu_curr; /* current local Lipschitz estimate */ 8 PetscReal lambda_prev; /* previous step length: for some reason SNESLineSearchGetLambda returns 1 instead of the previous step length */ 9 } SNESLineSearch_NLEQERR; 10 11 static PetscBool NLEQERR_cited = PETSC_FALSE; 12 static const char NLEQERR_citation[] = "@book{deuflhard2011,\n" 13 " title = {Newton Methods for Nonlinear Problems},\n" 14 " author = {Peter Deuflhard},\n" 15 " volume = 35,\n" 16 " year = 2011,\n" 17 " isbn = {978-3-642-23898-7},\n" 18 " doi = {10.1007/978-3-642-23899-4},\n" 19 " publisher = {Springer-Verlag},\n" 20 " address = {Berlin, Heidelberg}\n}\n"; 21 22 static PetscErrorCode SNESLineSearchReset_NLEQERR(SNESLineSearch linesearch) 23 { 24 SNESLineSearch_NLEQERR *nleqerr = (SNESLineSearch_NLEQERR *)linesearch->data; 25 26 PetscFunctionBegin; 27 nleqerr->mu_curr = 0.0; 28 nleqerr->norm_delta_x_prev = -1.0; 29 nleqerr->norm_bar_delta_x_prev = -1.0; 30 PetscFunctionReturn(PETSC_SUCCESS); 31 } 32 33 static PetscErrorCode SNESLineSearchApply_NLEQERR(SNESLineSearch linesearch) 34 { 35 PetscBool changed_y, changed_w; 36 Vec X, F, Y, W, G; 37 SNES snes; 38 PetscReal fnorm, xnorm, ynorm, gnorm, wnorm; 39 PetscReal lambda, minlambda, stol; 40 PetscViewer monitor; 41 PetscInt max_its, count, snes_iteration; 42 PetscReal theta, mudash, lambdadash; 43 SNESLineSearch_NLEQERR *nleqerr = (SNESLineSearch_NLEQERR *)linesearch->data; 44 KSPConvergedReason kspreason; 45 46 PetscFunctionBegin; 47 PetscCall(PetscCitationsRegister(NLEQERR_citation, &NLEQERR_cited)); 48 49 PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G)); 50 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm)); 51 PetscCall(SNESLineSearchGetLambda(linesearch, &lambda)); 52 PetscCall(SNESLineSearchGetSNES(linesearch, &snes)); 53 PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor)); 54 PetscCall(SNESLineSearchGetTolerances(linesearch, &minlambda, NULL, NULL, NULL, NULL, &max_its)); 55 PetscCall(SNESGetTolerances(snes, NULL, NULL, &stol, NULL, NULL)); 56 57 /* reset the state of the Lipschitz estimates */ 58 PetscCall(SNESGetIterationNumber(snes, &snes_iteration)); 59 if (!snes_iteration) PetscCall(SNESLineSearchReset_NLEQERR(linesearch)); 60 61 /* precheck */ 62 PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y)); 63 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 64 65 PetscCall(VecNormBegin(Y, NORM_2, &ynorm)); 66 PetscCall(VecNormBegin(X, NORM_2, &xnorm)); 67 PetscCall(VecNormEnd(Y, NORM_2, &ynorm)); 68 PetscCall(VecNormEnd(X, NORM_2, &xnorm)); 69 70 /* Note: Y is *minus* the Newton step. For whatever reason PETSc doesn't solve with the minus on the RHS. */ 71 72 if (ynorm == 0.0) { 73 if (monitor) { 74 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 75 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Initial direction and size is 0\n")); 76 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 77 } 78 PetscCall(VecCopy(X, W)); 79 PetscCall(VecCopy(F, G)); 80 PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm)); 81 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 82 PetscFunctionReturn(PETSC_SUCCESS); 83 } 84 85 /* At this point, we've solved the Newton system for delta_x, and we assume that 86 its norm is greater than the solution tolerance (otherwise we wouldn't be in 87 here). So let's go ahead and estimate the Lipschitz constant. 88 89 W contains bar_delta_x_prev at this point. */ 90 91 if (monitor) { 92 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 93 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: norm of Newton step: %14.12e\n", (double)ynorm)); 94 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 95 } 96 97 /* this needs information from a previous iteration, so can't do it on the first one */ 98 if (nleqerr->norm_delta_x_prev > 0 && nleqerr->norm_bar_delta_x_prev > 0) { 99 PetscCall(VecWAXPY(G, +1.0, Y, W)); /* bar_delta_x - delta_x; +1 because Y is -delta_x */ 100 PetscCall(VecNormBegin(G, NORM_2, &gnorm)); 101 PetscCall(VecNormEnd(G, NORM_2, &gnorm)); 102 103 nleqerr->mu_curr = nleqerr->lambda_prev * (nleqerr->norm_delta_x_prev * nleqerr->norm_bar_delta_x_prev) / (gnorm * ynorm); 104 lambda = PetscMin(1.0, nleqerr->mu_curr); 105 106 if (monitor) { 107 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 108 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Lipschitz estimate: %14.12e; lambda: %14.12e\n", (double)nleqerr->mu_curr, (double)lambda)); 109 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 110 } 111 } else { 112 lambda = linesearch->damping; 113 } 114 115 /* The main while loop of the algorithm. 116 At the end of this while loop, G should have the accepted new X in it. */ 117 118 count = 0; 119 while (PETSC_TRUE) { 120 if (monitor) { 121 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 122 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: entering iteration with lambda: %14.12e\n", (double)lambda)); 123 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 124 } 125 126 /* Check that we haven't performed too many iterations */ 127 count += 1; 128 if (count >= max_its) { 129 if (monitor) { 130 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 131 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: maximum iterations reached\n")); 132 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 133 } 134 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 135 PetscFunctionReturn(PETSC_SUCCESS); 136 } 137 138 /* Now comes the Regularity Test. */ 139 if (lambda <= minlambda) { 140 /* This isn't what is suggested by Deuflhard, but it works better in my experience */ 141 if (monitor) { 142 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 143 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: lambda has reached lambdamin, taking full Newton step\n")); 144 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 145 } 146 lambda = 1.0; 147 PetscCall(VecWAXPY(G, -lambda, Y, X)); 148 149 /* and clean up the state for next time */ 150 PetscCall(SNESLineSearchReset_NLEQERR(linesearch)); 151 /* 152 The clang static analyzer detected a problem here; once the loop is broken the values 153 nleqerr->norm_delta_x_prev = ynorm; 154 nleqerr->norm_bar_delta_x_prev = wnorm; 155 are set, but wnorm has not even been computed. 156 I don't know if this is the correct fix but by setting ynorm and wnorm to -1.0 at 157 least the linesearch object is kept in the state set by the SNESLineSearchReset_NLEQERR() call above 158 */ 159 ynorm = wnorm = -1.0; 160 break; 161 } 162 163 /* Compute new trial iterate */ 164 PetscCall(VecWAXPY(W, -lambda, Y, X)); 165 PetscCall(SNESComputeFunction(snes, W, G)); 166 167 /* Solve linear system for bar_delta_x_curr: old Jacobian, new RHS. Note absence of minus sign, compared to Deuflhard, in keeping with PETSc convention */ 168 PetscCall(KSPSolve(snes->ksp, G, W)); 169 PetscCall(KSPGetConvergedReason(snes->ksp, &kspreason)); 170 if (kspreason < 0) PetscCall(PetscInfo(snes, "Solution for \\bar{delta x}^{k+1} failed.\n")); 171 172 /* W now contains -bar_delta_x_curr. */ 173 174 PetscCall(VecNorm(W, NORM_2, &wnorm)); 175 if (monitor) { 176 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 177 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: norm of simplified Newton update: %14.12e\n", (double)wnorm)); 178 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 179 } 180 181 /* compute the monitoring quantities theta and mudash. */ 182 183 theta = wnorm / ynorm; 184 185 PetscCall(VecWAXPY(G, -(1.0 - lambda), Y, W)); 186 PetscCall(VecNorm(G, NORM_2, &gnorm)); 187 188 mudash = (0.5 * ynorm * lambda * lambda) / gnorm; 189 190 /* Check for termination of the linesearch */ 191 if (theta >= 1.0) { 192 /* need to go around again with smaller lambda */ 193 if (monitor) { 194 PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 195 PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: monotonicity check failed, ratio: %14.12e\n", (double)theta)); 196 PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 197 } 198 lambda = PetscMin(mudash, 0.5 * lambda); 199 lambda = PetscMax(lambda, minlambda); 200 /* continue through the loop, i.e. go back to regularity test */ 201 } else { 202 /* linesearch terminated */ 203 lambdadash = PetscMin(1.0, mudash); 204 205 if (lambdadash == 1.0 && lambda == 1.0 && wnorm <= stol) { 206 /* store the updated state, X - Y - W, in G: 207 I need to keep W for the next linesearch */ 208 PetscCall(VecWAXPY(G, -1.0, Y, X)); 209 PetscCall(VecAXPY(G, -1.0, W)); 210 break; 211 } 212 213 /* Deuflhard suggests to add the following: 214 else if (lambdadash >= 4.0 * lambda) { 215 lambda = lambdadash; 216 } 217 to continue through the loop, i.e. go back to regularity test. 218 I deliberately exclude this, as I have practical experience of this 219 getting stuck in infinite loops (on e.g. an Allen--Cahn problem). */ 220 221 else { 222 /* accept iterate without adding on, i.e. don't use bar_delta_x; 223 again, I need to keep W for the next linesearch */ 224 PetscCall(VecWAXPY(G, -lambda, Y, X)); 225 break; 226 } 227 } 228 } 229 230 if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, G)); 231 232 /* W currently contains -bar_delta_u. Scale it so that it contains bar_delta_u. */ 233 PetscCall(VecScale(W, -1.0)); 234 235 /* postcheck */ 236 PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, G, &changed_y, &changed_w)); 237 if (changed_y || changed_w) { 238 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_USER)); 239 PetscCall(PetscInfo(snes, "Changing the search direction here doesn't make sense.\n")); 240 PetscFunctionReturn(PETSC_SUCCESS); 241 } 242 243 /* copy the solution and information from this iteration over */ 244 nleqerr->norm_delta_x_prev = ynorm; 245 nleqerr->norm_bar_delta_x_prev = wnorm; 246 nleqerr->lambda_prev = lambda; 247 248 PetscCall(VecCopy(G, X)); 249 PetscCall(SNESComputeFunction(snes, X, F)); 250 PetscCall(VecNorm(X, NORM_2, &xnorm)); 251 PetscCall(VecNorm(F, NORM_2, &fnorm)); 252 PetscCall(SNESLineSearchSetLambda(linesearch, lambda)); 253 PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, (ynorm < 0 ? PETSC_INFINITY : ynorm))); 254 PetscFunctionReturn(PETSC_SUCCESS); 255 } 256 257 PetscErrorCode SNESLineSearchView_NLEQERR(SNESLineSearch linesearch, PetscViewer viewer) 258 { 259 PetscBool iascii; 260 SNESLineSearch_NLEQERR *nleqerr; 261 262 PetscFunctionBegin; 263 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 264 nleqerr = (SNESLineSearch_NLEQERR *)linesearch->data; 265 if (iascii) { 266 PetscCall(PetscViewerASCIIPrintf(viewer, " NLEQ-ERR affine-covariant linesearch")); 267 PetscCall(PetscViewerASCIIPrintf(viewer, " current local Lipschitz estimate omega=%e\n", (double)nleqerr->mu_curr)); 268 } 269 PetscFunctionReturn(PETSC_SUCCESS); 270 } 271 272 static PetscErrorCode SNESLineSearchDestroy_NLEQERR(SNESLineSearch linesearch) 273 { 274 PetscFunctionBegin; 275 PetscCall(PetscFree(linesearch->data)); 276 PetscFunctionReturn(PETSC_SUCCESS); 277 } 278 279 /*MC 280 SNESLINESEARCHNLEQERR - Error-oriented affine-covariant globalised Newton algorithm of Deuflhard (2011). 281 282 This linesearch is intended for Newton-type methods which are affine covariant. Affine covariance 283 means that Newton's method will give the same iterations for F(x) = 0 and AF(x) = 0 for a nonsingular 284 matrix A. This is a fundamental property; the philosophy of this linesearch is that globalisations 285 of Newton's method should carefully preserve it. 286 287 Options Database Keys: 288 + -snes_linesearch_damping<1.0> - initial step length 289 - -snes_linesearch_minlambda<1e-12> - minimum step length allowed 290 291 Level: advanced 292 293 Note: 294 Contributed by Patrick Farrell <patrick.farrell@maths.ox.ac.uk> 295 296 Reference: 297 . - * - Newton Methods for Nonlinear Problems, Deuflhard, P. 2011, Springer-Verlag, page 148 298 299 .seealso: `SNESLineSearch`, `SNES`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()` 300 M*/ 301 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_NLEQERR(SNESLineSearch linesearch) 302 { 303 SNESLineSearch_NLEQERR *nleqerr; 304 305 PetscFunctionBegin; 306 linesearch->ops->apply = SNESLineSearchApply_NLEQERR; 307 linesearch->ops->destroy = SNESLineSearchDestroy_NLEQERR; 308 linesearch->ops->setfromoptions = NULL; 309 linesearch->ops->reset = SNESLineSearchReset_NLEQERR; 310 linesearch->ops->view = SNESLineSearchView_NLEQERR; 311 linesearch->ops->setup = NULL; 312 313 PetscCall(PetscNew(&nleqerr)); 314 315 linesearch->data = (void *)nleqerr; 316 linesearch->max_its = 40; 317 PetscCall(SNESLineSearchReset_NLEQERR(linesearch)); 318 PetscFunctionReturn(PETSC_SUCCESS); 319 } 320