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