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