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