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