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", 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(VecCopy(X, G)); 213 PetscCall(VecAXPY(G, -1.0, Y)); 214 PetscCall(VecAXPY(G, -1.0, W)); 215 break; 216 } 217 218 /* Deuflhard suggests to add the following: 219 else if (lambdadash >= 4.0 * lambda) { 220 lambda = lambdadash; 221 } 222 to continue through the loop, i.e. go back to regularity test. 223 I deliberately exclude this, as I have practical experience of this 224 getting stuck in infinite loops (on e.g. an Allen--Cahn problem). */ 225 226 else { 227 /* accept iterate without adding on, i.e. don't use bar_delta_x; 228 again, I need to keep W for the next linesearch */ 229 PetscCall(VecWAXPY(G, -lambda, Y, X)); 230 break; 231 } 232 } 233 } 234 235 if (linesearch->ops->viproject) { 236 PetscCall((*linesearch->ops->viproject)(snes, G)); 237 } 238 239 /* W currently contains -bar_delta_u. Scale it so that it contains bar_delta_u. */ 240 PetscCall(VecScale(W, -1.0)); 241 242 /* postcheck */ 243 PetscCall(SNESLineSearchPostCheck(linesearch,X,Y,G,&changed_y,&changed_w)); 244 if (changed_y || changed_w) { 245 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_USER)); 246 PetscCall(PetscInfo(snes,"Changing the search direction here doesn't make sense.\n")); 247 PetscFunctionReturn(0); 248 } 249 250 /* copy the solution and information from this iteration over */ 251 nleqerr->norm_delta_x_prev = ynorm; 252 nleqerr->norm_bar_delta_x_prev = wnorm; 253 nleqerr->lambda_prev = lambda; 254 255 PetscCall(VecCopy(G, X)); 256 PetscCall(SNESComputeFunction(snes, X, F)); 257 PetscCall(VecNorm(X, NORM_2, &xnorm)); 258 PetscCall(VecNorm(F, NORM_2, &fnorm)); 259 PetscCall(SNESLineSearchSetLambda(linesearch, lambda)); 260 PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, (ynorm < 0 ? PETSC_INFINITY : ynorm))); 261 PetscFunctionReturn(0); 262 } 263 264 PetscErrorCode SNESLineSearchView_NLEQERR(SNESLineSearch linesearch, PetscViewer viewer) 265 { 266 PetscBool iascii; 267 SNESLineSearch_NLEQERR *nleqerr; 268 269 PetscFunctionBegin; 270 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 271 nleqerr = (SNESLineSearch_NLEQERR*)linesearch->data; 272 if (iascii) { 273 PetscCall(PetscViewerASCIIPrintf(viewer, " NLEQ-ERR affine-covariant linesearch")); 274 PetscCall(PetscViewerASCIIPrintf(viewer, " current local Lipschitz estimate omega=%e\n", (double)nleqerr->mu_curr)); 275 } 276 PetscFunctionReturn(0); 277 } 278 279 static PetscErrorCode SNESLineSearchDestroy_NLEQERR(SNESLineSearch linesearch) 280 { 281 PetscFunctionBegin; 282 PetscCall(PetscFree(linesearch->data)); 283 PetscFunctionReturn(0); 284 } 285 286 /*MC 287 SNESLINESEARCHNLEQERR - Error-oriented affine-covariant globalised Newton algorithm of Deuflhard (2011). 288 289 This linesearch is intended for Newton-type methods which are affine covariant. Affine covariance 290 means that Newton's method will give the same iterations for F(x) = 0 and AF(x) = 0 for a nonsingular 291 matrix A. This is a fundamental property; the philosophy of this linesearch is that globalisations 292 of Newton's method should carefully preserve it. 293 294 For a discussion of the theory behind this algorithm, see 295 296 @book{deuflhard2011, 297 title={Newton Methods for Nonlinear Problems}, 298 author={Deuflhard, P.}, 299 volume={35}, 300 year={2011}, 301 publisher={Springer-Verlag}, 302 address={Berlin, Heidelberg} 303 } 304 305 Pseudocode is given on page 148. 306 307 Options Database Keys: 308 + -snes_linesearch_damping<1.0> - initial step length 309 - -snes_linesearch_minlambda<1e-12> - minimum step length allowed 310 311 Contributed by Patrick Farrell <patrick.farrell@maths.ox.ac.uk> 312 313 Level: advanced 314 315 .seealso: SNESLineSearchCreate(), SNESLineSearchSetType() 316 M*/ 317 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_NLEQERR(SNESLineSearch linesearch) 318 { 319 SNESLineSearch_NLEQERR *nleqerr; 320 321 PetscFunctionBegin; 322 linesearch->ops->apply = SNESLineSearchApply_NLEQERR; 323 linesearch->ops->destroy = SNESLineSearchDestroy_NLEQERR; 324 linesearch->ops->setfromoptions = NULL; 325 linesearch->ops->reset = SNESLineSearchReset_NLEQERR; 326 linesearch->ops->view = SNESLineSearchView_NLEQERR; 327 linesearch->ops->setup = NULL; 328 329 PetscCall(PetscNewLog(linesearch,&nleqerr)); 330 331 linesearch->data = (void*)nleqerr; 332 linesearch->max_its = 40; 333 SNESLineSearchReset_NLEQERR(linesearch); 334 PetscFunctionReturn(0); 335 } 336