1 #include <private/linesearchimpl.h> 2 #include <private/snesimpl.h> 3 4 typedef enum {PETSCLINESEARCH_BT_QUADRATIC, PETSCLINESEARCH_BT_CUBIC} PetscLineSearchBTOrder; 5 6 typedef struct { 7 PetscReal alpha; /* sufficient decrease parameter */ 8 PetscLineSearchBTOrder order; 9 } PetscLineSearch_BT; 10 11 /*MC 12 PetscLineSearchBT - Backtracking line searches. 13 14 These linesearches try a polynomial fit for the L2 norm of the error 15 using the gradient. Failing that, they step back and try again. 16 17 Level: advanced 18 19 .keywords: SNES, PetscLineSearch, damping 20 21 .seealso: PetscLineSearchCreate(), PetscLineSearchSetType() 22 M*/ 23 24 #undef __FUNCT__ 25 #define __FUNCT__ "PetscLineSearchApply_BT" 26 27 PetscErrorCode PetscLineSearchApply_BT(PetscLineSearch linesearch) 28 { 29 PetscBool changed_y, changed_w; 30 PetscErrorCode ierr; 31 Vec X, F, Y, W, G; 32 SNES snes; 33 PetscReal fnorm, xnorm, ynorm, gnorm, gnormprev; 34 PetscReal lambda, lambdatemp, lambdaprev, minlambda, maxstep, rellength, initslope, alpha; 35 PetscReal t1, t2, a, b, d, steptol; 36 #if defined(PETSC_USE_COMPLEX) 37 PetscScalar cinitslope; 38 #endif 39 PetscBool domainerror; 40 PetscViewer monitor; 41 PetscInt max_its, count; 42 PetscLineSearch_BT *bt; 43 Mat jac; 44 45 46 PetscFunctionBegin; 47 48 ierr = PetscLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G);CHKERRQ(ierr); 49 ierr = PetscLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr); 50 ierr = PetscLineSearchGetLambda(linesearch, &lambda);CHKERRQ(ierr); 51 ierr = PetscLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr); 52 ierr = PetscLineSearchGetMonitor(linesearch, &monitor);CHKERRQ(ierr); 53 ierr = PetscLineSearchGetTolerances(linesearch, &steptol, &maxstep, PETSC_NULL, PETSC_NULL, PETSC_NULL, &max_its); 54 bt = (PetscLineSearch_BT *)linesearch->data; 55 56 alpha = bt->alpha; 57 58 ierr = SNESGetJacobian(snes, &jac, PETSC_NULL, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 59 if (!jac) { 60 SETERRQ(((PetscObject)linesearch)->comm, PETSC_ERR_USER, "PetscLineSearchBT requires a Jacobian matrix"); 61 } 62 /* precheck */ 63 ierr = PetscLineSearchPreCheck(linesearch, &changed_y);CHKERRQ(ierr); 64 ierr = PetscLineSearchSetSuccess(linesearch, PETSC_TRUE);CHKERRQ(ierr); 65 66 ierr = VecNorm(Y, NORM_2, &ynorm);CHKERRQ(ierr); 67 if (ynorm == 0.0) { 68 if (monitor) { 69 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 70 ierr = PetscViewerASCIIPrintf(monitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr); 71 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 72 } 73 ierr = VecCopy(X,W);CHKERRQ(ierr); 74 ierr = VecCopy(F,G);CHKERRQ(ierr); 75 ierr = PetscLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 76 PetscFunctionReturn(0); 77 } 78 if (ynorm > maxstep) { /* Step too big, so scale back */ 79 if (monitor) { 80 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 81 ierr = PetscViewerASCIIPrintf(monitor," Line search: Scaling step by %14.12e old ynorm %14.12e\n", (maxstep/ynorm),ynorm);CHKERRQ(ierr); 82 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 83 } 84 ierr = VecScale(Y,maxstep/(ynorm));CHKERRQ(ierr); 85 ynorm = maxstep; 86 } 87 ierr = VecMaxPointwiseDivide(Y,X,&rellength);CHKERRQ(ierr); 88 minlambda = steptol/rellength; 89 ierr = MatMult(jac,Y,W);CHKERRQ(ierr); 90 #if defined(PETSC_USE_COMPLEX) 91 ierr = VecDot(F,W,&cinitslope);CHKERRQ(ierr); 92 initslope = PetscRealPart(cinitslope); 93 #else 94 ierr = VecDot(F,W,&initslope);CHKERRQ(ierr); 95 #endif 96 if (initslope > 0.0) initslope = -initslope; 97 if (initslope == 0.0) initslope = -1.0; 98 99 ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr); 100 if (linesearch->ops->viproject) { 101 ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 102 } 103 if (snes->nfuncs >= snes->max_funcs) { 104 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 105 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 106 ierr = PetscLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 107 PetscFunctionReturn(0); 108 } 109 ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr); 110 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 111 if (domainerror) { 112 ierr = PetscLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 113 PetscFunctionReturn(0); 114 } 115 if (linesearch->ops->vinorm) { 116 gnorm = fnorm; 117 ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr); 118 } else { 119 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); 120 } 121 122 if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 123 ierr = PetscInfo2(snes,"Initial fnorm %14.12e gnorm %14.12e\n", fnorm, gnorm);CHKERRQ(ierr); 124 if (.5*gnorm*gnorm <= .5*fnorm*fnorm + lambda*alpha*initslope) { /* Sufficient reduction */ 125 if (monitor) { 126 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 127 ierr = PetscViewerASCIIPrintf(monitor," Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", fnorm, gnorm);CHKERRQ(ierr); 128 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 129 } 130 } else { 131 /* Fit points with quadratic */ 132 lambdatemp = -initslope/(gnorm*gnorm - fnorm*fnorm - 2.0*lambda*initslope); 133 lambdaprev = lambda; 134 gnormprev = gnorm; 135 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 136 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 137 else lambda = lambdatemp; 138 139 ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr); 140 if (linesearch->ops->viproject) { 141 ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 142 } 143 if (snes->nfuncs >= snes->max_funcs) { 144 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 145 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 146 ierr = PetscLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 147 PetscFunctionReturn(0); 148 } 149 ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr); 150 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 151 if (domainerror) { 152 PetscFunctionReturn(0); 153 } 154 if (linesearch->ops->vinorm) { 155 gnorm = fnorm; 156 ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr); 157 } else { 158 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); 159 } 160 if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 161 if (monitor) { 162 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 163 ierr = PetscViewerASCIIPrintf(monitor," Line search: gnorm after quadratic fit %14.12e\n",gnorm);CHKERRQ(ierr); 164 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 165 } 166 if (.5*gnorm*gnorm < .5*fnorm*fnorm + lambda*alpha*initslope) { /* sufficient reduction */ 167 if (monitor) { 168 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 169 ierr = PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda);CHKERRQ(ierr); 170 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 171 } 172 } else { 173 /* Fit points with cubic */ 174 for (count = 0; count < max_its; count++) { 175 if (lambda <= minlambda) { 176 if (monitor) { 177 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 178 ierr = PetscViewerASCIIPrintf(monitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 179 ierr = PetscViewerASCIIPrintf(monitor, 180 " Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", 181 fnorm, gnorm, ynorm, minlambda, lambda, initslope);CHKERRQ(ierr); 182 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 183 } 184 ierr = PetscLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 185 PetscFunctionReturn(0); 186 } 187 if (bt->order == PETSCLINESEARCH_BT_CUBIC) { 188 t1 = .5*(gnorm*gnorm - fnorm*fnorm) - lambda*initslope; 189 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 190 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 191 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 192 d = b*b - 3*a*initslope; 193 if (d < 0.0) d = 0.0; 194 if (a == 0.0) { 195 lambdatemp = -initslope/(2.0*b); 196 } else { 197 lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a); 198 } 199 } else if (bt->order == PETSCLINESEARCH_BT_QUADRATIC) { 200 lambdatemp = -initslope/(gnorm*gnorm - fnorm*fnorm - 2.0*initslope); 201 } 202 lambdaprev = lambda; 203 gnormprev = gnorm; 204 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 205 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 206 else lambda = lambdatemp; 207 ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr); 208 if (snes->nfuncs >= snes->max_funcs) { 209 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 210 ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", 211 fnorm,gnorm,ynorm,lambda,initslope);CHKERRQ(ierr); 212 ierr = PetscLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 213 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 214 PetscFunctionReturn(0); 215 } 216 ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr); 217 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 218 if (domainerror) { 219 ierr = PetscLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 220 PetscFunctionReturn(0); 221 } 222 if (linesearch->ops->vinorm) { 223 gnorm = fnorm; 224 ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr); 225 } else { 226 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); 227 } 228 if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 229 if (.5*gnorm*gnorm < .5*fnorm*fnorm + lambda*alpha*initslope) { /* is reduction enough? */ 230 if (monitor) { 231 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 232 if (bt->order == PETSCLINESEARCH_BT_CUBIC) { 233 ierr = PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr); 234 } else { 235 ierr = PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr); 236 } 237 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 238 } 239 break; 240 } else { 241 if (monitor) { 242 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 243 if (bt->order == PETSCLINESEARCH_BT_CUBIC) { 244 ierr = PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr); 245 } else { 246 ierr = PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr); 247 } 248 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 249 } 250 } 251 } 252 } 253 } 254 255 /* postcheck */ 256 ierr = PetscLineSearchPostCheck(linesearch, &changed_y, &changed_w);CHKERRQ(ierr); 257 if (changed_y) { 258 ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr); 259 if (linesearch->ops->viproject) { 260 ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 261 } 262 } 263 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 264 ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr); 265 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 266 if (domainerror) { 267 ierr = PetscLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 268 PetscFunctionReturn(0); 269 } 270 if (linesearch->ops->vinorm) { 271 gnorm = fnorm; 272 ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr); 273 } else { 274 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); 275 } 276 ierr = VecNorm(Y,NORM_2,&ynorm);CHKERRQ(ierr); 277 if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 278 279 } 280 281 /* copy the solution over */ 282 ierr = VecCopy(W, X);CHKERRQ(ierr); 283 ierr = VecCopy(G, F);CHKERRQ(ierr); 284 ierr = VecNorm(X, NORM_2, &xnorm);CHKERRQ(ierr); 285 ierr = PetscLineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr); 286 ierr = PetscLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);CHKERRQ(ierr); 287 PetscFunctionReturn(0); 288 } 289 290 #undef __FUNCT__ 291 #define __FUNCT__ "PetscLineSearchDestroy_BT" 292 PetscErrorCode PetscLineSearchDestroy_BT(PetscLineSearch linesearch) 293 { 294 PetscErrorCode ierr; 295 296 PetscFunctionBegin; 297 ierr = PetscFree(linesearch->data);CHKERRQ(ierr); 298 PetscFunctionReturn(0); 299 } 300 301 302 #undef __FUNCT__ 303 #define __FUNCT__ "PetscLineSearchSetFromOptions_BT" 304 static PetscErrorCode PetscLineSearchSetFromOptions_BT(PetscLineSearch linesearch) 305 { 306 307 PetscErrorCode ierr; 308 PetscLineSearch_BT *bt; 309 const char *orders[] = {"quadratic", "cubic"}; 310 PetscInt indx = 0; 311 PetscBool flg; 312 PetscFunctionBegin; 313 314 bt = (PetscLineSearch_BT*)linesearch->data; 315 316 ierr = PetscOptionsHead("PetscLineSearch BT options");CHKERRQ(ierr); 317 ierr = PetscOptionsReal("-linesearch_bt_alpha", "Descent tolerance", "PetscLineSearchBT", bt->alpha, &bt->alpha, PETSC_NULL);CHKERRQ(ierr); 318 ierr = PetscOptionsEList("-linesearch_bt_order", "Order of approximation", "PetscLineSearchBT", orders,2,"quadratic",&indx,&flg);CHKERRQ(ierr); 319 if (flg) { 320 switch (indx) { 321 case 0: bt->order = PETSCLINESEARCH_BT_QUADRATIC; 322 break; 323 case 1: bt->order = PETSCLINESEARCH_BT_CUBIC; 324 break; 325 } 326 } 327 328 ierr = PetscOptionsTail();CHKERRQ(ierr); 329 PetscFunctionReturn(0); 330 } 331 332 333 EXTERN_C_BEGIN 334 #undef __FUNCT__ 335 #define __FUNCT__ "PetscLineSearchCreate_BT" 336 PetscErrorCode PetscLineSearchCreate_BT(PetscLineSearch linesearch) 337 { 338 339 PetscLineSearch_BT *bt; 340 PetscErrorCode ierr; 341 342 PetscFunctionBegin; 343 linesearch->ops->apply = PetscLineSearchApply_BT; 344 linesearch->ops->destroy = PetscLineSearchDestroy_BT; 345 linesearch->ops->setfromoptions = PetscLineSearchSetFromOptions_BT; 346 linesearch->ops->reset = PETSC_NULL; 347 linesearch->ops->view = PETSC_NULL; 348 linesearch->ops->setup = PETSC_NULL; 349 350 ierr = PetscNewLog(linesearch, PetscLineSearch_BT, &bt);CHKERRQ(ierr); 351 linesearch->data = (void *)bt; 352 linesearch->max_its = 40; 353 bt->order = PETSCLINESEARCH_BT_CUBIC; 354 bt->alpha = 1e-4; 355 356 PetscFunctionReturn(0); 357 } 358 EXTERN_C_END 359