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