1 #include <petsc-private/linesearchimpl.h> /*I "petscsnes.h" I*/ 2 #include <petsc-private/snesimpl.h> 3 4 typedef struct { 5 PetscReal alpha; /* sufficient decrease parameter */ 6 } SNESLineSearch_BT; 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "SNESLineSearchBTSetAlpha" 10 /*@ 11 SNESLineSearchBTSetAlpha - Sets the descent parameter, alpha, in the BT linesearch variant. 12 13 Input Parameters: 14 + linesearch - linesearch context 15 - alpha - The descent parameter 16 17 Level: intermediate 18 19 .seealso: SNESLineSearchSetLambda(), SNESLineSearchGetTolerances() SNESLINESEARCHBT 20 @*/ 21 PetscErrorCode SNESLineSearchBTSetAlpha(SNESLineSearch linesearch, PetscReal alpha) 22 { 23 SNESLineSearch_BT *bt; 24 PetscFunctionBegin; 25 PetscValidHeaderSpecific(linesearch,SNESLINESEARCH_CLASSID,1); 26 bt = (SNESLineSearch_BT *)linesearch->data; 27 bt->alpha = alpha; 28 PetscFunctionReturn(0); 29 } 30 31 32 #undef __FUNCT__ 33 #define __FUNCT__ "SNESLineSearchBTGetAlpha" 34 /*@ 35 SNESLineSearchBTGetAlpha - Gets the descent parameter, alpha, in the BT linesearch variant. 36 37 Input Parameters: 38 . linesearch - linesearch context 39 40 Output Parameters: 41 . alpha - The descent parameter 42 43 Level: intermediate 44 45 .seealso: SNESLineSearchGetLambda(), SNESLineSearchGetTolerances() SNESLINESEARCHBT 46 @*/ 47 PetscErrorCode SNESLineSearchBTGetAlpha(SNESLineSearch linesearch, PetscReal *alpha) 48 { 49 SNESLineSearch_BT *bt; 50 PetscFunctionBegin; 51 PetscValidHeaderSpecific(linesearch,SNESLINESEARCH_CLASSID,1); 52 bt = (SNESLineSearch_BT *)linesearch->data; 53 *alpha = bt->alpha; 54 PetscFunctionReturn(0); 55 } 56 57 #undef __FUNCT__ 58 #define __FUNCT__ "SNESLineSearchApply_BT" 59 static PetscErrorCode SNESLineSearchApply_BT(SNESLineSearch linesearch) 60 { 61 PetscBool changed_y, changed_w; 62 PetscErrorCode ierr; 63 Vec X, F, Y, W, G; 64 SNES snes; 65 PetscReal fnorm, xnorm, ynorm, gnorm, gnormprev; 66 PetscReal lambda, lambdatemp, lambdaprev, minlambda, maxstep, rellength, initslope, alpha; 67 PetscReal t1, t2, a, b, d, steptol; 68 #if defined(PETSC_USE_COMPLEX) 69 PetscScalar cinitslope; 70 #endif 71 PetscBool domainerror; 72 PetscViewer monitor; 73 PetscInt max_its, count; 74 SNESLineSearch_BT *bt; 75 Mat jac; 76 77 78 PetscFunctionBegin; 79 80 ierr = SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G);CHKERRQ(ierr); 81 ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr); 82 ierr = SNESLineSearchGetLambda(linesearch, &lambda);CHKERRQ(ierr); 83 ierr = SNESLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr); 84 ierr = SNESLineSearchGetMonitor(linesearch, &monitor);CHKERRQ(ierr); 85 ierr = SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, PETSC_NULL, PETSC_NULL, PETSC_NULL, &max_its); 86 bt = (SNESLineSearch_BT *)linesearch->data; 87 88 alpha = bt->alpha; 89 90 ierr = SNESGetJacobian(snes, &jac, PETSC_NULL, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 91 if (!jac) { 92 SETERRQ(((PetscObject)linesearch)->comm, PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix"); 93 } 94 /* precheck */ 95 ierr = SNESLineSearchPreCheck(linesearch, &changed_y);CHKERRQ(ierr); 96 ierr = SNESLineSearchSetSuccess(linesearch, PETSC_TRUE);CHKERRQ(ierr); 97 98 ierr = VecNorm(Y, NORM_2, &ynorm);CHKERRQ(ierr); 99 if (ynorm == 0.0) { 100 if (monitor) { 101 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 102 ierr = PetscViewerASCIIPrintf(monitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr); 103 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 104 } 105 ierr = VecCopy(X,W);CHKERRQ(ierr); 106 ierr = VecCopy(F,G);CHKERRQ(ierr); 107 ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 108 PetscFunctionReturn(0); 109 } 110 if (ynorm > maxstep) { /* Step too big, so scale back */ 111 if (monitor) { 112 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 113 ierr = PetscViewerASCIIPrintf(monitor," Line search: Scaling step by %14.12e old ynorm %14.12e\n", (maxstep/ynorm),ynorm);CHKERRQ(ierr); 114 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 115 } 116 ierr = VecScale(Y,maxstep/(ynorm));CHKERRQ(ierr); 117 ynorm = maxstep; 118 } 119 ierr = VecMaxPointwiseDivide(Y,X,&rellength);CHKERRQ(ierr); 120 minlambda = steptol/rellength; 121 ierr = MatMult(jac,Y,W);CHKERRQ(ierr); 122 #if defined(PETSC_USE_COMPLEX) 123 ierr = VecDot(F,W,&cinitslope);CHKERRQ(ierr); 124 initslope = PetscRealPart(cinitslope); 125 #else 126 ierr = VecDot(F,W,&initslope);CHKERRQ(ierr); 127 #endif 128 if (initslope > 0.0) initslope = -initslope; 129 if (initslope == 0.0) initslope = -1.0; 130 131 ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr); 132 if (linesearch->ops->viproject) { 133 ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 134 } 135 if (snes->nfuncs >= snes->max_funcs) { 136 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 137 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 138 ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 139 PetscFunctionReturn(0); 140 } 141 ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr); 142 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 143 if (domainerror) { 144 ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 145 PetscFunctionReturn(0); 146 } 147 if (linesearch->ops->vinorm) { 148 gnorm = fnorm; 149 ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr); 150 } else { 151 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); 152 } 153 154 if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 155 ierr = PetscInfo2(snes,"Initial fnorm %14.12e gnorm %14.12e\n", fnorm, gnorm);CHKERRQ(ierr); 156 if (.5*gnorm*gnorm <= .5*fnorm*fnorm + lambda*alpha*initslope) { /* Sufficient reduction */ 157 if (monitor) { 158 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 159 ierr = PetscViewerASCIIPrintf(monitor," Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", fnorm, gnorm);CHKERRQ(ierr); 160 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 161 } 162 } else { 163 /* Fit points with quadratic */ 164 lambdatemp = -initslope/(gnorm*gnorm - fnorm*fnorm - 2.0*lambda*initslope); 165 lambdaprev = lambda; 166 gnormprev = gnorm; 167 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 168 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 169 else lambda = lambdatemp; 170 171 ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr); 172 if (linesearch->ops->viproject) { 173 ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 174 } 175 if (snes->nfuncs >= snes->max_funcs) { 176 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 177 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 178 ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 179 PetscFunctionReturn(0); 180 } 181 ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr); 182 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 183 if (domainerror) { 184 PetscFunctionReturn(0); 185 } 186 if (linesearch->ops->vinorm) { 187 gnorm = fnorm; 188 ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr); 189 } else { 190 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); 191 } 192 if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 193 if (monitor) { 194 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 195 ierr = PetscViewerASCIIPrintf(monitor," Line search: gnorm after quadratic fit %14.12e\n",gnorm);CHKERRQ(ierr); 196 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 197 } 198 if (.5*gnorm*gnorm < .5*fnorm*fnorm + lambda*alpha*initslope) { /* sufficient reduction */ 199 if (monitor) { 200 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 201 ierr = PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda);CHKERRQ(ierr); 202 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 203 } 204 } else { 205 /* Fit points with cubic */ 206 for (count = 0; count < max_its; count++) { 207 if (lambda <= minlambda) { 208 if (monitor) { 209 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 210 ierr = PetscViewerASCIIPrintf(monitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 211 ierr = PetscViewerASCIIPrintf(monitor, 212 " Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", 213 fnorm, gnorm, ynorm, minlambda, lambda, initslope);CHKERRQ(ierr); 214 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 215 } 216 ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 217 PetscFunctionReturn(0); 218 } 219 if (linesearch->order == SNES_LINESEARCH_CUBIC) { 220 t1 = .5*(gnorm*gnorm - fnorm*fnorm) - lambda*initslope; 221 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 222 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 223 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 224 d = b*b - 3*a*initslope; 225 if (d < 0.0) d = 0.0; 226 if (a == 0.0) { 227 lambdatemp = -initslope/(2.0*b); 228 } else { 229 lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a); 230 } 231 } else if (linesearch->order == SNES_LINESEARCH_QUADRATIC) { 232 lambdatemp = -initslope/(gnorm*gnorm - fnorm*fnorm - 2.0*initslope); 233 } else { 234 SETERRQ(((PetscObject)linesearch)->comm, PETSC_ERR_SUP, "unsupported line search order for type bt"); 235 } 236 lambdaprev = lambda; 237 gnormprev = gnorm; 238 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 239 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 240 else lambda = lambdatemp; 241 ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr); 242 if (snes->nfuncs >= snes->max_funcs) { 243 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 244 ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", 245 fnorm,gnorm,ynorm,lambda,initslope);CHKERRQ(ierr); 246 ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 247 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 248 PetscFunctionReturn(0); 249 } 250 ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr); 251 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 252 if (domainerror) { 253 ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 254 PetscFunctionReturn(0); 255 } 256 if (linesearch->ops->vinorm) { 257 gnorm = fnorm; 258 ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr); 259 } else { 260 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); 261 } 262 if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 263 if (.5*gnorm*gnorm < .5*fnorm*fnorm + lambda*alpha*initslope) { /* is reduction enough? */ 264 if (monitor) { 265 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 266 if (linesearch->order == SNES_LINESEARCH_CUBIC) { 267 ierr = PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr); 268 } else { 269 ierr = PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr); 270 } 271 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 272 } 273 break; 274 } else { 275 if (monitor) { 276 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 277 if (linesearch->order == SNES_LINESEARCH_CUBIC) { 278 ierr = PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr); 279 } else { 280 ierr = PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr); 281 } 282 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 283 } 284 } 285 } 286 } 287 } 288 289 /* postcheck */ 290 ierr = SNESLineSearchPostCheck(linesearch, &changed_y, &changed_w);CHKERRQ(ierr); 291 if (changed_y) { 292 ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr); 293 if (linesearch->ops->viproject) { 294 ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 295 } 296 } 297 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 298 ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr); 299 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 300 if (domainerror) { 301 ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 302 PetscFunctionReturn(0); 303 } 304 if (linesearch->ops->vinorm) { 305 gnorm = fnorm; 306 ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr); 307 } else { 308 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); 309 } 310 ierr = VecNorm(Y,NORM_2,&ynorm);CHKERRQ(ierr); 311 if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 312 313 } 314 315 /* copy the solution over */ 316 ierr = VecCopy(W, X);CHKERRQ(ierr); 317 ierr = VecCopy(G, F);CHKERRQ(ierr); 318 ierr = VecNorm(X, NORM_2, &xnorm);CHKERRQ(ierr); 319 ierr = SNESLineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr); 320 ierr = SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);CHKERRQ(ierr); 321 PetscFunctionReturn(0); 322 } 323 324 #undef __FUNCT__ 325 #define __FUNCT__ "SNESLineSearchView_BT" 326 PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer) 327 { 328 PetscErrorCode ierr; 329 PetscBool iascii; 330 SNESLineSearch_BT *bt; 331 PetscFunctionBegin; 332 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 333 bt = (SNESLineSearch_BT*)linesearch->data; 334 if (iascii) { 335 if (linesearch->order == SNES_LINESEARCH_CUBIC) { 336 ierr = PetscViewerASCIIPrintf(viewer, " interpolation: cubic\n");CHKERRQ(ierr); 337 } else if (linesearch->order == SNES_LINESEARCH_QUADRATIC) { 338 ierr = PetscViewerASCIIPrintf(viewer, " interpolation: quadratic\n");CHKERRQ(ierr); 339 } 340 ierr = PetscViewerASCIIPrintf(viewer, " alpha=%G\n", bt->alpha);CHKERRQ(ierr); 341 } 342 PetscFunctionReturn(0); 343 } 344 345 346 #undef __FUNCT__ 347 #define __FUNCT__ "SNESLineSearchDestroy_BT" 348 static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch) 349 { 350 PetscErrorCode ierr; 351 352 PetscFunctionBegin; 353 ierr = PetscFree(linesearch->data);CHKERRQ(ierr); 354 PetscFunctionReturn(0); 355 } 356 357 358 #undef __FUNCT__ 359 #define __FUNCT__ "SNESLineSearchSetFromOptions_BT" 360 static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch) 361 { 362 363 PetscErrorCode ierr; 364 SNESLineSearch_BT *bt; 365 PetscFunctionBegin; 366 367 bt = (SNESLineSearch_BT*)linesearch->data; 368 369 ierr = PetscOptionsHead("SNESLineSearch BT options");CHKERRQ(ierr); 370 ierr = PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, PETSC_NULL);CHKERRQ(ierr); 371 372 ierr = PetscOptionsTail();CHKERRQ(ierr); 373 PetscFunctionReturn(0); 374 } 375 376 377 #undef __FUNCT__ 378 #define __FUNCT__ "SNESLineSearchCreate_BT" 379 /*MC 380 SNESLINESEARCHBT - Backtracking line searches. 381 382 These linesearches try a polynomial fit for the L2 norm of the error 383 using the gradient. Failing that, they step back and try again. 384 385 Options Database Keys: 386 + -snes_linesearch_alpha<1e-4> - slope descent parameter 387 . -snes_linesearch_damping<1.0> - full steplength 388 - -snes_linesearch_order<cubic, quadratic> - order of the approximation 389 390 Level: advanced 391 392 .keywords: SNES, SNESLineSearch, damping 393 394 .seealso: SNESLineSearchCreate(), SNESLineSearchSetType() 395 M*/ 396 PETSC_EXTERN_C PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch) 397 { 398 399 SNESLineSearch_BT *bt; 400 PetscErrorCode ierr; 401 402 PetscFunctionBegin; 403 linesearch->ops->apply = SNESLineSearchApply_BT; 404 linesearch->ops->destroy = SNESLineSearchDestroy_BT; 405 linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT; 406 linesearch->ops->reset = PETSC_NULL; 407 linesearch->ops->view = SNESLineSearchView_BT; 408 linesearch->ops->setup = PETSC_NULL; 409 410 ierr = PetscNewLog(linesearch, SNESLineSearch_BT, &bt);CHKERRQ(ierr); 411 linesearch->data = (void *)bt; 412 linesearch->max_its = 40; 413 linesearch->order = SNES_LINESEARCH_CUBIC; 414 bt->alpha = 1e-4; 415 416 PetscFunctionReturn(0); 417 } 418