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