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