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