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,-1.0,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 if (snes->domainerror) { 111 ierr = PetscLogEventEnd(SNES_LineSearch,snes,X,F,G);CHKERRQ(ierr); 112 PetscFunctionReturn(0); 113 } 114 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); 115 if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 116 ierr = PetscInfo2(snes,"Initial fnorm %14.12e gnorm %14.12e\n", fnorm, gnorm);CHKERRQ(ierr); 117 if (.5*gnorm*gnorm <= .5*fnorm*fnorm + lambda*alpha*initslope) { /* Sufficient reduction */ 118 if (monitor) { 119 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 120 ierr = PetscViewerASCIIPrintf(monitor," Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", fnorm, gnorm);CHKERRQ(ierr); 121 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 122 } 123 } else { 124 /* Fit points with quadratic */ 125 lambdatemp = -initslope/(gnorm*gnorm - fnorm*fnorm - 2.0*lambda*initslope); 126 lambdaprev = lambda; 127 gnormprev = gnorm; 128 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 129 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 130 else lambda = lambdatemp; 131 132 ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr); 133 if (snes->nfuncs >= snes->max_funcs) { 134 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 135 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 136 ierr = LineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 137 PetscFunctionReturn(0); 138 } 139 ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr); 140 if (snes->domainerror) { 141 ierr = PetscLogEventEnd(SNES_LineSearch,snes,X,F,G);CHKERRQ(ierr); 142 PetscFunctionReturn(0); 143 } 144 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); 145 if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 146 if (monitor) { 147 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 148 ierr = PetscViewerASCIIPrintf(monitor," Line search: gnorm after quadratic fit %14.12e\n",gnorm);CHKERRQ(ierr); 149 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 150 } 151 if (.5*gnorm*gnorm < .5*fnorm*fnorm + lambda*alpha*initslope) { /* sufficient reduction */ 152 if (monitor) { 153 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 154 ierr = PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda);CHKERRQ(ierr); 155 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 156 } 157 } else { 158 /* Fit points with cubic */ 159 for (count = 0; count < max_its; count++) { 160 if (lambda <= minlambda) { 161 if (monitor) { 162 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 163 ierr = PetscViewerASCIIPrintf(monitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 164 ierr = PetscViewerASCIIPrintf(monitor, 165 " Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", 166 fnorm, gnorm, ynorm, minlambda, lambda, initslope);CHKERRQ(ierr); 167 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 168 } 169 ierr = LineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 170 PetscFunctionReturn(0); 171 } 172 if (bt->order == LINESEARCH_BT_CUBIC) { 173 t1 = .5*(gnorm*gnorm - fnorm*fnorm) - lambda*initslope; 174 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 175 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 176 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 177 d = b*b - 3*a*initslope; 178 if (d < 0.0) d = 0.0; 179 if (a == 0.0) { 180 lambdatemp = -initslope/(2.0*b); 181 } else { 182 lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a); 183 } 184 } else if (bt->order == LINESEARCH_BT_QUADRATIC) { 185 lambdatemp = -initslope/(gnorm*gnorm - fnorm*fnorm - 2.0*initslope); 186 } 187 lambdaprev = lambda; 188 gnormprev = gnorm; 189 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 190 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 191 else lambda = lambdatemp; 192 ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr); 193 if (snes->nfuncs >= snes->max_funcs) { 194 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 195 ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", 196 fnorm,gnorm,ynorm,lambda,initslope);CHKERRQ(ierr); 197 ierr = LineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 198 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 199 PetscFunctionReturn(0); 200 } 201 ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr); 202 if (snes->domainerror) { 203 ierr = PetscLogEventEnd(SNES_LineSearch,snes,X,F,G);CHKERRQ(ierr); 204 PetscFunctionReturn(0); 205 } 206 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); 207 if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 208 if (.5*gnorm*gnorm < .5*fnorm*fnorm + lambda*alpha*initslope) { /* is reduction enough? */ 209 if (monitor) { 210 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 211 if (bt->order == LINESEARCH_BT_CUBIC) { 212 ierr = PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr); 213 } else { 214 ierr = PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr); 215 } 216 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 217 } 218 break; 219 } else { 220 if (monitor) { 221 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 222 if (bt->order == LINESEARCH_BT_CUBIC) { 223 ierr = PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr); 224 } else { 225 ierr = PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",gnorm,lambda);CHKERRQ(ierr); 226 } 227 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 228 } 229 } 230 } 231 } 232 } 233 234 /* postcheck */ 235 ierr = LineSearchPostCheck(linesearch, &changed_y, &changed_w);CHKERRQ(ierr); 236 if (changed_y) { 237 ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr); 238 } 239 ierr = SNESComputeFunction(snes,W,F);CHKERRQ(ierr); 240 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 241 if (domainerror) { 242 ierr = LineSearchSetSuccess(linesearch, PETSC_FALSE); 243 PetscFunctionReturn(0); 244 } 245 246 /* copy the solution over */ 247 ierr = VecCopy(W, X);CHKERRQ(ierr); 248 ierr = LineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr); 249 ierr = LineSearchComputeNorms(linesearch);CHKERRQ(ierr); 250 PetscFunctionReturn(0); 251 } 252 253 #undef __FUNCT__ 254 #define __FUNCT__ "LineSearchDestroy_BT" 255 PetscErrorCode LineSearchDestroy_BT(LineSearch linesearch) { 256 257 PetscFunctionBegin; 258 PetscFunctionReturn(0); 259 260 } 261 262 263 #undef __FUNCT__ 264 #define __FUNCT__ "LineSearchSetFromOptions_BT" 265 static PetscErrorCode LineSearchSetFromOptions_BT(LineSearch linesearch) 266 { 267 268 PetscErrorCode ierr; 269 LineSearch_BT *bt; 270 const char *orders[] = {"quadratic", "cubic"}; 271 PetscInt indx = 0; 272 PetscBool flg; 273 PetscFunctionBegin; 274 275 bt = (LineSearch_BT*)linesearch->data; 276 277 ierr = PetscOptionsHead("LineSearch BT options");CHKERRQ(ierr); 278 ierr = PetscOptionsReal("-linesearch_bt_alpha", "Descent tolerance", "LineSearchBT", bt->alpha, &bt->alpha, PETSC_NULL);CHKERRQ(ierr); 279 ierr = PetscOptionsEList("-linesearch_bt_order", "Order of approximation", "LineSearchBT", orders,2,"quadratic",&indx,&flg);CHKERRQ(ierr); 280 if (flg) { 281 switch (indx) { 282 case 0: bt->order = LINESEARCH_BT_QUADRATIC; 283 break; 284 case 1: bt->order = LINESEARCH_BT_CUBIC; 285 break; 286 } 287 } 288 289 ierr = PetscOptionsTail();CHKERRQ(ierr); 290 PetscFunctionReturn(0); 291 } 292 293 294 EXTERN_C_BEGIN 295 #undef __FUNCT__ 296 #define __FUNCT__ "LineSearchCreate_BT" 297 PetscErrorCode LineSearchCreate_BT(LineSearch linesearch) 298 { 299 300 LineSearch_BT *bt; 301 PetscErrorCode ierr; 302 303 PetscFunctionBegin; 304 linesearch->ops->apply = LineSearchApply_BT; 305 linesearch->ops->destroy = LineSearchDestroy_BT; 306 linesearch->ops->setfromoptions = LineSearchSetFromOptions_BT; 307 linesearch->ops->reset = PETSC_NULL; 308 linesearch->ops->view = PETSC_NULL; 309 linesearch->ops->setup = PETSC_NULL; 310 311 ierr = PetscNewLog(linesearch, LineSearch_BT, &bt);CHKERRQ(ierr); 312 linesearch->data = (void *)bt; 313 linesearch->max_its = 40; 314 bt->order = LINESEARCH_BT_CUBIC; 315 bt->alpha = 1e-4; 316 317 PetscFunctionReturn(0); 318 } 319 EXTERN_C_END 320