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