1 #include <private/linesearchimpl.h> 2 #include <private/snesimpl.h> 3 4 typedef enum {LINESEARCH_BT_QUADRATIC, LINESEARCH_BT_CUBIC} LineSearchBTOrder; 5 6 typedef struct { 7 PetscReal alpha; /* sufficient decrease parameter */ 8 LineSearchBTOrder order; 9 } LineSearch_BT; 10 11 /*MC 12 LineSearchBT - 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, LineSearch, damping 20 21 .seealso: LineSearchCreate(), LineSearchSetType() 22 M*/ 23 24 #undef __FUNCT__ 25 #define __FUNCT__ "LineSearchApply_BT" 26 27 PetscErrorCode LineSearchApply_BT(LineSearch 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 LineSearch_BT *bt; 43 Mat jac; 44 45 46 PetscFunctionBegin; 47 48 ierr = LineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G);CHKERRQ(ierr); 49 ierr = LineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr); 50 ierr = LineSearchGetLambda(linesearch, &lambda);CHKERRQ(ierr); 51 ierr = LineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr); 52 ierr = LineSearchGetMonitor(linesearch, &monitor);CHKERRQ(ierr); 53 ierr = LineSearchGetTolerances(linesearch, &steptol, &maxstep, PETSC_NULL, PETSC_NULL, PETSC_NULL, &max_its); 54 bt = (LineSearch_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, "LineSearchBT requires a Jacobian matrix"); 61 } 62 63 ierr = PetscLogEventBegin(LineSearch_Apply,snes,X,F,G);CHKERRQ(ierr); 64 65 /* precheck */ 66 ierr = LineSearchPreCheck(linesearch, &changed_y);CHKERRQ(ierr); 67 ierr = LineSearchSetSuccess(linesearch, PETSC_TRUE);CHKERRQ(ierr); 68 69 ierr = VecNorm(Y, NORM_2, &ynorm);CHKERRQ(ierr); 70 if (ynorm == 0.0) { 71 if (monitor) { 72 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 73 ierr = PetscViewerASCIIPrintf(monitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr); 74 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 75 } 76 ierr = VecCopy(X,W);CHKERRQ(ierr); 77 ierr = VecCopy(F,G);CHKERRQ(ierr); 78 ierr = LineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 79 PetscFunctionReturn(0); 80 } 81 if (ynorm > maxstep) { /* Step too big, so scale back */ 82 if (monitor) { 83 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 84 ierr = PetscViewerASCIIPrintf(monitor," Line search: Scaling step by %14.12e old ynorm %14.12e\n", (maxstep/ynorm),ynorm);CHKERRQ(ierr); 85 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 86 } 87 ierr = VecScale(Y,maxstep/(ynorm));CHKERRQ(ierr); 88 ynorm = maxstep; 89 } 90 ierr = VecMaxPointwiseDivide(Y,X,&rellength);CHKERRQ(ierr); 91 minlambda = steptol/rellength; 92 ierr = MatMult(jac,Y,W);CHKERRQ(ierr); 93 #if defined(PETSC_USE_COMPLEX) 94 ierr = VecDot(F,W,&cinitslope);CHKERRQ(ierr); 95 initslope = PetscRealPart(cinitslope); 96 #else 97 ierr = VecDot(F,W,&initslope);CHKERRQ(ierr); 98 #endif 99 if (initslope > 0.0) initslope = -initslope; 100 if (initslope == 0.0) initslope = -1.0; 101 102 ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr); 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 = LineSearchSetSuccess(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 = LineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 113 ierr = PetscLogEventEnd(SNES_LineSearch,snes,X,F,G);CHKERRQ(ierr); 114 PetscFunctionReturn(0); 115 } 116 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); 117 if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 118 ierr = PetscInfo2(snes,"Initial fnorm %14.12e gnorm %14.12e\n", fnorm, gnorm);CHKERRQ(ierr); 119 if (.5*gnorm*gnorm <= .5*fnorm*fnorm + lambda*alpha*initslope) { /* Sufficient reduction */ 120 if (monitor) { 121 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 122 ierr = PetscViewerASCIIPrintf(monitor," Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", fnorm, gnorm);CHKERRQ(ierr); 123 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 124 } 125 } else { 126 /* Fit points with quadratic */ 127 lambdatemp = -initslope/(gnorm*gnorm - fnorm*fnorm - 2.0*lambda*initslope); 128 lambdaprev = lambda; 129 gnormprev = gnorm; 130 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 131 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 132 else lambda = lambdatemp; 133 134 ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr); 135 if (snes->nfuncs >= snes->max_funcs) { 136 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 137 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 138 ierr = LineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 139 PetscFunctionReturn(0); 140 } 141 ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr); 142 if (snes->domainerror) { 143 ierr = PetscLogEventEnd(SNES_LineSearch,snes,X,F,G);CHKERRQ(ierr); 144 PetscFunctionReturn(0); 145 } 146 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); 147 if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 148 if (monitor) { 149 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 150 ierr = PetscViewerASCIIPrintf(monitor," Line search: gnorm after quadratic fit %14.12e\n",gnorm);CHKERRQ(ierr); 151 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 152 } 153 if (.5*gnorm*gnorm < .5*fnorm*fnorm + lambda*alpha*initslope) { /* sufficient reduction */ 154 if (monitor) { 155 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 156 ierr = PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda);CHKERRQ(ierr); 157 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 158 } 159 } else { 160 /* Fit points with cubic */ 161 for (count = 0; count < max_its; count++) { 162 if (lambda <= minlambda) { 163 if (monitor) { 164 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 165 ierr = PetscViewerASCIIPrintf(monitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 166 ierr = PetscViewerASCIIPrintf(monitor, 167 " Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", 168 fnorm, gnorm, ynorm, minlambda, lambda, initslope);CHKERRQ(ierr); 169 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 170 } 171 ierr = LineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 172 PetscFunctionReturn(0); 173 } 174 if (bt->order == LINESEARCH_BT_CUBIC) { 175 t1 = .5*(gnorm*gnorm - fnorm*fnorm) - lambda*initslope; 176 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 177 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 178 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 179 d = b*b - 3*a*initslope; 180 if (d < 0.0) d = 0.0; 181 if (a == 0.0) { 182 lambdatemp = -initslope/(2.0*b); 183 } else { 184 lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a); 185 } 186 } else if (bt->order == LINESEARCH_BT_QUADRATIC) { 187 lambdatemp = -initslope/(gnorm*gnorm - fnorm*fnorm - 2.0*initslope); 188 } 189 lambdaprev = lambda; 190 gnormprev = gnorm; 191 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 192 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 193 else lambda = lambdatemp; 194 ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr); 195 if (snes->nfuncs >= snes->max_funcs) { 196 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 197 ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", 198 fnorm,gnorm,ynorm,lambda,initslope);CHKERRQ(ierr); 199 ierr = LineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 200 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 201 PetscFunctionReturn(0); 202 } 203 ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr); 204 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 205 if (domainerror) { 206 ierr = PetscLogEventEnd(SNES_LineSearch,snes,X,F,G);CHKERRQ(ierr); 207 PetscFunctionReturn(0); 208 } 209 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); 210 if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 211 if (.5*gnorm*gnorm < .5*fnorm*fnorm + lambda*alpha*initslope) { /* is reduction enough? */ 212 if (monitor) { 213 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 214 if (bt->order == LINESEARCH_BT_CUBIC) { 215 ierr = PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr); 216 } else { 217 ierr = PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr); 218 } 219 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 220 } 221 break; 222 } else { 223 if (monitor) { 224 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 225 if (bt->order == LINESEARCH_BT_CUBIC) { 226 ierr = PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr); 227 } else { 228 ierr = PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr); 229 } 230 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 231 } 232 } 233 } 234 } 235 } 236 237 /* postcheck */ 238 ierr = LineSearchPostCheck(linesearch, &changed_y, &changed_w);CHKERRQ(ierr); 239 if (changed_y) { 240 ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr); 241 } 242 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 243 ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr); 244 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 245 if (domainerror) { 246 ierr = LineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 247 ierr = PetscLogEventEnd(SNES_LineSearch,linesearch,X,F,G);CHKERRQ(ierr); 248 PetscFunctionReturn(0); 249 } 250 ierr = VecNormBegin(G,NORM_2,&gnorm);CHKERRQ(ierr); 251 if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 252 ierr = VecNormBegin(Y,NORM_2,&ynorm);CHKERRQ(ierr); 253 ierr = VecNormEnd(G,NORM_2,&gnorm);CHKERRQ(ierr); 254 ierr = VecNormEnd(Y,NORM_2,&ynorm);CHKERRQ(ierr); 255 } 256 257 /* copy the solution over */ 258 ierr = VecCopy(W, X);CHKERRQ(ierr); 259 ierr = VecCopy(G, F);CHKERRQ(ierr); 260 ierr = VecNorm(X, NORM_2, &xnorm);CHKERRQ(ierr); 261 ierr = LineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr); 262 ierr = LineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);CHKERRQ(ierr); 263 PetscFunctionReturn(0); 264 } 265 266 #undef __FUNCT__ 267 #define __FUNCT__ "LineSearchDestroy_BT" 268 PetscErrorCode LineSearchDestroy_BT(LineSearch linesearch) { 269 270 PetscFunctionBegin; 271 PetscFunctionReturn(0); 272 273 } 274 275 276 #undef __FUNCT__ 277 #define __FUNCT__ "LineSearchSetFromOptions_BT" 278 static PetscErrorCode LineSearchSetFromOptions_BT(LineSearch linesearch) 279 { 280 281 PetscErrorCode ierr; 282 LineSearch_BT *bt; 283 const char *orders[] = {"quadratic", "cubic"}; 284 PetscInt indx = 0; 285 PetscBool flg; 286 PetscFunctionBegin; 287 288 bt = (LineSearch_BT*)linesearch->data; 289 290 ierr = PetscOptionsHead("LineSearch BT options");CHKERRQ(ierr); 291 ierr = PetscOptionsReal("-linesearch_bt_alpha", "Descent tolerance", "LineSearchBT", bt->alpha, &bt->alpha, PETSC_NULL);CHKERRQ(ierr); 292 ierr = PetscOptionsEList("-linesearch_bt_order", "Order of approximation", "LineSearchBT", orders,2,"quadratic",&indx,&flg);CHKERRQ(ierr); 293 if (flg) { 294 switch (indx) { 295 case 0: bt->order = LINESEARCH_BT_QUADRATIC; 296 break; 297 case 1: bt->order = LINESEARCH_BT_CUBIC; 298 break; 299 } 300 } 301 302 ierr = PetscOptionsTail();CHKERRQ(ierr); 303 PetscFunctionReturn(0); 304 } 305 306 307 EXTERN_C_BEGIN 308 #undef __FUNCT__ 309 #define __FUNCT__ "LineSearchCreate_BT" 310 PetscErrorCode LineSearchCreate_BT(LineSearch linesearch) 311 { 312 313 LineSearch_BT *bt; 314 PetscErrorCode ierr; 315 316 PetscFunctionBegin; 317 linesearch->ops->apply = LineSearchApply_BT; 318 linesearch->ops->destroy = LineSearchDestroy_BT; 319 linesearch->ops->setfromoptions = LineSearchSetFromOptions_BT; 320 linesearch->ops->reset = PETSC_NULL; 321 linesearch->ops->view = PETSC_NULL; 322 linesearch->ops->setup = PETSC_NULL; 323 324 ierr = PetscNewLog(linesearch, LineSearch_BT, &bt);CHKERRQ(ierr); 325 linesearch->data = (void *)bt; 326 linesearch->max_its = 40; 327 bt->order = LINESEARCH_BT_CUBIC; 328 bt->alpha = 1e-4; 329 330 PetscFunctionReturn(0); 331 } 332 EXTERN_C_END 333