1af0996ceSBarry Smith #include <petsc/private/linesearchimpl.h> /*I "petscsnes.h" I*/ 2af0996ceSBarry Smith #include <petsc/private/snesimpl.h> 371b4ebd2SPeter Brune 471b4ebd2SPeter Brune typedef struct { 571b4ebd2SPeter Brune PetscReal alpha; /* sufficient decrease parameter */ 6f1c6b773SPeter Brune } SNESLineSearch_BT; 771b4ebd2SPeter Brune 82f4102e2SPeter Brune /*@ 92f4102e2SPeter Brune SNESLineSearchBTSetAlpha - Sets the descent parameter, alpha, in the BT linesearch variant. 102f4102e2SPeter Brune 112f4102e2SPeter Brune Input Parameters: 122f4102e2SPeter Brune + linesearch - linesearch context 132f4102e2SPeter Brune - alpha - The descent parameter 142f4102e2SPeter Brune 152f4102e2SPeter Brune Level: intermediate 162f4102e2SPeter Brune 17db781477SPatrick Sanan .seealso: `SNESLineSearchSetLambda()`, `SNESLineSearchGetTolerances()` `SNESLINESEARCHBT` 182f4102e2SPeter Brune @*/ 192f4102e2SPeter Brune PetscErrorCode SNESLineSearchBTSetAlpha(SNESLineSearch linesearch, PetscReal alpha) 202f4102e2SPeter Brune { 21d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data; 226e111a19SKarl Rupp 232f4102e2SPeter Brune PetscFunctionBegin; 242f4102e2SPeter Brune PetscValidHeaderSpecific(linesearch,SNESLINESEARCH_CLASSID,1); 252f4102e2SPeter Brune bt->alpha = alpha; 262f4102e2SPeter Brune PetscFunctionReturn(0); 272f4102e2SPeter Brune } 282f4102e2SPeter Brune 292f4102e2SPeter Brune /*@ 302f4102e2SPeter Brune SNESLineSearchBTGetAlpha - Gets the descent parameter, alpha, in the BT linesearch variant. 312f4102e2SPeter Brune 322f4102e2SPeter Brune Input Parameters: 332f4102e2SPeter Brune . linesearch - linesearch context 348e557f58SPeter Brune 358e557f58SPeter Brune Output Parameters: 362f4102e2SPeter Brune . alpha - The descent parameter 372f4102e2SPeter Brune 382f4102e2SPeter Brune Level: intermediate 392f4102e2SPeter Brune 40db781477SPatrick Sanan .seealso: `SNESLineSearchGetLambda()`, `SNESLineSearchGetTolerances()` `SNESLINESEARCHBT` 412f4102e2SPeter Brune @*/ 422f4102e2SPeter Brune PetscErrorCode SNESLineSearchBTGetAlpha(SNESLineSearch linesearch, PetscReal *alpha) 432f4102e2SPeter Brune { 44d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data; 456e111a19SKarl Rupp 462f4102e2SPeter Brune PetscFunctionBegin; 472f4102e2SPeter Brune PetscValidHeaderSpecific(linesearch,SNESLINESEARCH_CLASSID,1); 482f4102e2SPeter Brune *alpha = bt->alpha; 492f4102e2SPeter Brune PetscFunctionReturn(0); 502f4102e2SPeter Brune } 512f4102e2SPeter Brune 52f1c6b773SPeter Brune static PetscErrorCode SNESLineSearchApply_BT(SNESLineSearch linesearch) 5371b4ebd2SPeter Brune { 5471b4ebd2SPeter Brune PetscBool changed_y,changed_w; 5571b4ebd2SPeter Brune Vec X,F,Y,W,G; 5671b4ebd2SPeter Brune SNES snes; 578a430afbSPeter Brune PetscReal fnorm, xnorm, ynorm, gnorm; 583b288129SPeter Brune PetscReal lambda,lambdatemp,lambdaprev,minlambda,maxstep,initslope,alpha,stol; 59dfcd3828SPeter Brune PetscReal t1,t2,a,b,d; 608a430afbSPeter Brune PetscReal f; 618a430afbSPeter Brune PetscReal g,gprev; 6271b4ebd2SPeter Brune PetscViewer monitor; 6371b4ebd2SPeter Brune PetscInt max_its,count; 64d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data; 6571b4ebd2SPeter Brune Mat jac; 66bd4a8a71SBarry Smith PetscErrorCode (*objective)(SNES,Vec,PetscReal*,void*); 6771b4ebd2SPeter Brune 6871b4ebd2SPeter Brune PetscFunctionBegin; 699566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G)); 709566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm)); 719566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetLambda(linesearch, &lambda)); 729566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetSNES(linesearch, &snes)); 739566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor)); 749566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetTolerances(linesearch,&minlambda,&maxstep,NULL,NULL,NULL,&max_its)); 759566063dSJacob Faibussowitsch PetscCall(SNESGetTolerances(snes,NULL,NULL,&stol,NULL,NULL)); 769566063dSJacob Faibussowitsch PetscCall(SNESGetObjective(snes,&objective,NULL)); 7771b4ebd2SPeter Brune alpha = bt->alpha; 7871b4ebd2SPeter Brune 799566063dSJacob Faibussowitsch PetscCall(SNESGetJacobian(snes, &jac, NULL, NULL, NULL)); 8008401ef6SPierre Jolivet PetscCheck(jac || objective,PetscObjectComm((PetscObject)linesearch), PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix"); 81c69d1a72SBarry Smith 829566063dSJacob Faibussowitsch PetscCall(SNESLineSearchPreCheck(linesearch,X,Y,&changed_y)); 839566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 8471b4ebd2SPeter Brune 859566063dSJacob Faibussowitsch PetscCall(VecNormBegin(Y, NORM_2, &ynorm)); 869566063dSJacob Faibussowitsch PetscCall(VecNormBegin(X, NORM_2, &xnorm)); 879566063dSJacob Faibussowitsch PetscCall(VecNormEnd(Y, NORM_2, &ynorm)); 889566063dSJacob Faibussowitsch PetscCall(VecNormEnd(X, NORM_2, &xnorm)); 893b288129SPeter Brune 9071b4ebd2SPeter Brune if (ynorm == 0.0) { 9171b4ebd2SPeter Brune if (monitor) { 929566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 939566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Initial direction and size is 0\n")); 949566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 9571b4ebd2SPeter Brune } 969566063dSJacob Faibussowitsch PetscCall(VecCopy(X,W)); 979566063dSJacob Faibussowitsch PetscCall(VecCopy(F,G)); 989566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm)); 999566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 10071b4ebd2SPeter Brune PetscFunctionReturn(0); 10171b4ebd2SPeter Brune } 10271b4ebd2SPeter Brune if (ynorm > maxstep) { /* Step too big, so scale back */ 10371b4ebd2SPeter Brune if (monitor) { 1049566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 1059566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Scaling step by %14.12e old ynorm %14.12e\n", (double)(maxstep/ynorm),(double)ynorm)); 1069566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 10771b4ebd2SPeter Brune } 1089566063dSJacob Faibussowitsch PetscCall(VecScale(Y,maxstep/(ynorm))); 10971b4ebd2SPeter Brune ynorm = maxstep; 11071b4ebd2SPeter Brune } 1118a430afbSPeter Brune 1128a430afbSPeter Brune /* if the SNES has an objective set, use that instead of the function value */ 113bd4a8a71SBarry Smith if (objective) { 1149566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes,X,&f)); 1158a430afbSPeter Brune } else { 1168a430afbSPeter Brune f = fnorm*fnorm; 1178a430afbSPeter Brune } 1188a430afbSPeter Brune 1198a430afbSPeter Brune /* compute the initial slope */ 120bd4a8a71SBarry Smith if (objective) { 1218a430afbSPeter Brune /* slope comes from the function (assumed to be the gradient of the objective */ 1229566063dSJacob Faibussowitsch PetscCall(VecDotRealPart(Y,F,&initslope)); 1238a430afbSPeter Brune } else { 1248a430afbSPeter Brune /* slope comes from the normal equations */ 1259566063dSJacob Faibussowitsch PetscCall(MatMult(jac,Y,W)); 1269566063dSJacob Faibussowitsch PetscCall(VecDotRealPart(F,W,&initslope)); 12771b4ebd2SPeter Brune if (initslope > 0.0) initslope = -initslope; 12871b4ebd2SPeter Brune if (initslope == 0.0) initslope = -1.0; 1298a430afbSPeter Brune } 13071b4ebd2SPeter Brune 131df019d78SBarry Smith while (PETSC_TRUE) { 1329566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W,-lambda,Y,X)); 133*1baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 134e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 1359566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n")); 13671b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1379566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 13871b4ebd2SPeter Brune PetscFunctionReturn(0); 13971b4ebd2SPeter Brune } 1408a430afbSPeter Brune 141bd4a8a71SBarry Smith if (objective) { 1429566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes,W,&g)); 1438a430afbSPeter Brune } else { 1449566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes,W,G)); 1459bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 1469bd66eb0SPeter Brune gnorm = fnorm; 1479566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 1489bd66eb0SPeter Brune } else { 1499566063dSJacob Faibussowitsch PetscCall(VecNorm(G,NORM_2,&gnorm)); 1508a430afbSPeter Brune } 15179e7e81bSJed Brown g = PetscSqr(gnorm); 1529bd66eb0SPeter Brune } 1539566063dSJacob Faibussowitsch PetscCall(SNESLineSearchMonitor(linesearch)); 1549bd66eb0SPeter Brune 155df019d78SBarry Smith if (!PetscIsInfOrNanReal(g)) break; 156df019d78SBarry Smith if (monitor) { 1579566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 1589566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n",(double)lambda)); 1599566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 160df019d78SBarry Smith } 161df019d78SBarry Smith if (lambda <= minlambda) { 1627e49c775SBarry Smith SNESCheckFunctionNorm(snes,g); 163c98378a5SBarry Smith } 164df019d78SBarry Smith lambda = .5*lambda; 165df019d78SBarry Smith } 166df019d78SBarry Smith 167bd4a8a71SBarry Smith if (!objective) { 1689566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm)); 1698a430afbSPeter Brune } 1708a430afbSPeter Brune if (.5*g <= .5*f + lambda*alpha*initslope) { /* Sufficient reduction or step tolerance convergence */ 17171b4ebd2SPeter Brune if (monitor) { 1729566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 173bd4a8a71SBarry Smith if (!objective) { 1749566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm)); 1758a430afbSPeter Brune } else { 1769566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Using full step: obj %14.12e obj %14.12e\n", (double)f, (double)g)); 1778a430afbSPeter Brune } 1789566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 17971b4ebd2SPeter Brune } 18071b4ebd2SPeter Brune } else { 181c21ba15cSPeter Brune /* Since the full step didn't work and the step is tiny, quit */ 182c21ba15cSPeter Brune if (stol*xnorm > ynorm) { 1839566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm)); 1849566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 185c61ba1e2SPeter Brune if (monitor) { 1869566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 18763a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n",(double)ynorm,(double)(stol*xnorm))); 1889566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 189c61ba1e2SPeter Brune } 190c21ba15cSPeter Brune PetscFunctionReturn(0); 191c21ba15cSPeter Brune } 19271b4ebd2SPeter Brune /* Fit points with quadratic */ 1938a430afbSPeter Brune lambdatemp = -initslope/(g - f - 2.0*lambda*initslope); 19471b4ebd2SPeter Brune lambdaprev = lambda; 1958a430afbSPeter Brune gprev = g; 19671b4ebd2SPeter Brune if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 19771b4ebd2SPeter Brune if (lambdatemp <= .1*lambda) lambda = .1*lambda; 19871b4ebd2SPeter Brune else lambda = lambdatemp; 19971b4ebd2SPeter Brune 2009566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W,-lambda,Y,X)); 201*1baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 202e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 20363a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %" PetscInt_FMT " \n",snes->nfuncs)); 20471b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2059566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 20671b4ebd2SPeter Brune PetscFunctionReturn(0); 20771b4ebd2SPeter Brune } 208bd4a8a71SBarry Smith if (objective) { 2099566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes,W,&g)); 2108a430afbSPeter Brune } else { 2119566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes,W,G)); 2129bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 2139bd66eb0SPeter Brune gnorm = fnorm; 2149566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 2159bd66eb0SPeter Brune } else { 2169566063dSJacob Faibussowitsch PetscCall(VecNorm(G,NORM_2,&gnorm)); 2179bd66eb0SPeter Brune } 218f1ab69e3SPeter Brune g = PetscSqr(gnorm); 2198a430afbSPeter Brune } 220c98378a5SBarry Smith if (PetscIsInfOrNanReal(g)) { 2219566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 2229566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n")); 223c98378a5SBarry Smith PetscFunctionReturn(0); 224c98378a5SBarry Smith } 22571b4ebd2SPeter Brune if (monitor) { 2269566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 227bd4a8a71SBarry Smith if (!objective) { 2289566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: gnorm after quadratic fit %14.12e\n",(double)gnorm)); 2298a430afbSPeter Brune } else { 2309566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: obj after quadratic fit %14.12e\n",(double)g)); 2318a430afbSPeter Brune } 2329566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 23371b4ebd2SPeter Brune } 2348a430afbSPeter Brune if (.5*g < .5*f + lambda*alpha*initslope) { /* sufficient reduction */ 23571b4ebd2SPeter Brune if (monitor) { 2369566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 2379566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda)); 2389566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 23971b4ebd2SPeter Brune } 24071b4ebd2SPeter Brune } else { 24171b4ebd2SPeter Brune /* Fit points with cubic */ 24271b4ebd2SPeter Brune for (count = 0; count < max_its; count++) { 24371b4ebd2SPeter Brune if (lambda <= minlambda) { 24471b4ebd2SPeter Brune if (monitor) { 2459566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 24663a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: unable to find good step length! After %" PetscInt_FMT " tries \n",count)); 247bd4a8a71SBarry Smith if (!objective) { 248d0609cedSBarry Smith PetscCall(PetscViewerASCIIPrintf(monitor," Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", 249d0609cedSBarry Smith (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope)); 2508a430afbSPeter Brune } else { 251d0609cedSBarry Smith PetscCall(PetscViewerASCIIPrintf(monitor," Line search: obj(0)=%18.16e, obj=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", 252d0609cedSBarry Smith (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope)); 2538a430afbSPeter Brune } 2549566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 25571b4ebd2SPeter Brune } 2569566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 25771b4ebd2SPeter Brune PetscFunctionReturn(0); 25871b4ebd2SPeter Brune } 259b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 2608a430afbSPeter Brune t1 = .5*(g - f) - lambda*initslope; 2618a430afbSPeter Brune t2 = .5*(gprev - f) - lambdaprev*initslope; 26271b4ebd2SPeter Brune a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 26371b4ebd2SPeter Brune b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 26471b4ebd2SPeter Brune d = b*b - 3*a*initslope; 26571b4ebd2SPeter Brune if (d < 0.0) d = 0.0; 266f5af7f23SKarl Rupp if (a == 0.0) lambdatemp = -initslope/(2.0*b); 267f5af7f23SKarl Rupp else lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a); 268f5af7f23SKarl Rupp 269b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 2708a430afbSPeter Brune lambdatemp = -initslope/(g - f - 2.0*initslope); 271ce94432eSBarry Smith } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt"); 27271b4ebd2SPeter Brune lambdaprev = lambda; 2738a430afbSPeter Brune gprev = g; 27471b4ebd2SPeter Brune if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 27571b4ebd2SPeter Brune if (lambdatemp <= .1*lambda) lambda = .1*lambda; 27671b4ebd2SPeter Brune else lambda = lambdatemp; 2779566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W,-lambda,Y,X)); 278*1baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes,W)); 279e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 28063a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Exceeded maximum function evaluations, while looking for good step length! %" PetscInt_FMT " \n",count)); 281bd4a8a71SBarry Smith if (!objective) { 282d0609cedSBarry Smith PetscCall(PetscInfo(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",(double)fnorm,(double)gnorm,(double)ynorm,(double)lambda,(double)initslope)); 2838a430afbSPeter Brune } 2849566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 28571b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 28671b4ebd2SPeter Brune PetscFunctionReturn(0); 28771b4ebd2SPeter Brune } 288bd4a8a71SBarry Smith if (objective) { 2899566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes,W,&g)); 2908a430afbSPeter Brune } else { 2919566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes,W,G)); 2929bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 2939bd66eb0SPeter Brune gnorm = fnorm; 2949566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 2959bd66eb0SPeter Brune } else { 2969566063dSJacob Faibussowitsch PetscCall(VecNorm(G,NORM_2,&gnorm)); 2978a430afbSPeter Brune } 298f1ab69e3SPeter Brune g = PetscSqr(gnorm); 2999bd66eb0SPeter Brune } 3004a2f8832SBarry Smith if (PetscIsInfOrNanReal(g)) { 3019566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 3029566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n")); 303c98378a5SBarry Smith PetscFunctionReturn(0); 304c98378a5SBarry Smith } 3058a430afbSPeter Brune if (.5*g < .5*f + lambda*alpha*initslope) { /* is reduction enough? */ 30671b4ebd2SPeter Brune if (monitor) { 3079566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 308bd4a8a71SBarry Smith if (!objective) { 309b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3109566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda)); 31171b4ebd2SPeter Brune } else { 3129566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda)); 31371b4ebd2SPeter Brune } 3149566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 3158a430afbSPeter Brune } else { 3168a430afbSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3179566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda)); 3188a430afbSPeter Brune } else { 3199566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda)); 3208a430afbSPeter Brune } 3219566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 3228a430afbSPeter Brune } 32371b4ebd2SPeter Brune } 32471b4ebd2SPeter Brune break; 325f5af7f23SKarl Rupp } else if (monitor) { 3269566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 327bd4a8a71SBarry Smith if (!objective) { 328b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3299566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda)); 33071b4ebd2SPeter Brune } else { 3319566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda)); 33271b4ebd2SPeter Brune } 3339566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 3348a430afbSPeter Brune } else { 3358a430afbSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3369566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda)); 3378a430afbSPeter Brune } else { 3389566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda)); 3398a430afbSPeter Brune } 3409566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 3418a430afbSPeter Brune } 34271b4ebd2SPeter Brune } 34371b4ebd2SPeter Brune } 34471b4ebd2SPeter Brune } 34571b4ebd2SPeter Brune } 34671b4ebd2SPeter Brune 34771b4ebd2SPeter Brune /* postcheck */ 3480c1ac483SGlenn Hammond /* update Y to lambda*Y so that W is consistent with X - lambda*Y */ 3499566063dSJacob Faibussowitsch PetscCall(VecScale(Y,lambda)); 3509566063dSJacob Faibussowitsch PetscCall(SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w)); 35171b4ebd2SPeter Brune if (changed_y) { 3529566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W,-1.0,Y,X)); 353*1baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 35471b4ebd2SPeter Brune } 355bd4a8a71SBarry Smith if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */ 3569566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes,W,G)); 3579bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 3589bd66eb0SPeter Brune gnorm = fnorm; 3599566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 3609bd66eb0SPeter Brune } else { 3619566063dSJacob Faibussowitsch PetscCall(VecNorm(G,NORM_2,&gnorm)); 3629bd66eb0SPeter Brune } 3639566063dSJacob Faibussowitsch PetscCall(VecNorm(Y,NORM_2,&ynorm)); 364c98378a5SBarry Smith if (PetscIsInfOrNanReal(gnorm)) { 3659566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_FAILED_NANORINF)); 3669566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n")); 367c98378a5SBarry Smith PetscFunctionReturn(0); 368c98378a5SBarry Smith } 369c427506bSPeter Brune } 37071b4ebd2SPeter Brune 37171b4ebd2SPeter Brune /* copy the solution over */ 3729566063dSJacob Faibussowitsch PetscCall(VecCopy(W, X)); 3739566063dSJacob Faibussowitsch PetscCall(VecCopy(G, F)); 3749566063dSJacob Faibussowitsch PetscCall(VecNorm(X, NORM_2, &xnorm)); 3759566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetLambda(linesearch, lambda)); 3769566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm)); 37771b4ebd2SPeter Brune PetscFunctionReturn(0); 37871b4ebd2SPeter Brune } 37971b4ebd2SPeter Brune 3807f1410a3SPeter Brune PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer) 3817f1410a3SPeter Brune { 3827f1410a3SPeter Brune PetscBool iascii; 383d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data; 384075cc632SBarry Smith 3857f1410a3SPeter Brune PetscFunctionBegin; 3869566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 3877f1410a3SPeter Brune if (iascii) { 388b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3899566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: cubic\n")); 390b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 3919566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: quadratic\n")); 3927f1410a3SPeter Brune } 3939566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " alpha=%e\n", (double)bt->alpha)); 3947f1410a3SPeter Brune } 3957f1410a3SPeter Brune PetscFunctionReturn(0); 3967f1410a3SPeter Brune } 3977f1410a3SPeter Brune 398f1c6b773SPeter Brune static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch) 399731c067cSJed Brown { 40071b4ebd2SPeter Brune PetscFunctionBegin; 4019566063dSJacob Faibussowitsch PetscCall(PetscFree(linesearch->data)); 40271b4ebd2SPeter Brune PetscFunctionReturn(0); 40371b4ebd2SPeter Brune } 40471b4ebd2SPeter Brune 4054416b707SBarry Smith static PetscErrorCode SNESLineSearchSetFromOptions_BT(PetscOptionItems *PetscOptionsObject,SNESLineSearch linesearch) 40671b4ebd2SPeter Brune { 407eb23ba39SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data; 40871b4ebd2SPeter Brune 4096e111a19SKarl Rupp PetscFunctionBegin; 410d0609cedSBarry Smith PetscOptionsHeadBegin(PetscOptionsObject,"SNESLineSearch BT options"); 4119566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL)); 412d0609cedSBarry Smith PetscOptionsHeadEnd(); 41371b4ebd2SPeter Brune PetscFunctionReturn(0); 41471b4ebd2SPeter Brune } 41571b4ebd2SPeter Brune 416954494b2SJed Brown /*MC 417db609ea7SPeter Brune SNESLINESEARCHBT - Backtracking line search. 418954494b2SJed Brown 419db609ea7SPeter Brune This line search finds the minimum of a polynomial fitting of the L2 norm of the 420eb23ba39SBarry Smith function or the objective function if it is provided with SNESSetObjective(). If this fit does not satisfy the conditions for progress, the interval shrinks 421db609ea7SPeter Brune and the fit is reattempted at most max_it times or until lambda is below minlambda. 422954494b2SJed Brown 423cd7522eaSPeter Brune Options Database Keys: 4243eec2f6aSPierre Jolivet + -snes_linesearch_alpha <1e\-4> - slope descent parameter 425db609ea7SPeter Brune . -snes_linesearch_damping <1.0> - initial step length 426eb23ba39SBarry Smith . -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the 427e1bc860dSBarry Smith step is scaled back to be of this length at the beginning of the line search 428db609ea7SPeter Brune . -snes_linesearch_max_it <40> - maximum number of shrinking step 4293eec2f6aSPierre Jolivet . -snes_linesearch_minlambda <1e\-12> - minimum step length allowed 430cd7522eaSPeter Brune - -snes_linesearch_order <cubic,quadratic> - order of the approximation 431cd7522eaSPeter Brune 432954494b2SJed Brown Level: advanced 433954494b2SJed Brown 434db609ea7SPeter Brune Notes: 435db609ea7SPeter Brune This line search is taken from "Numerical Methods for Unconstrained 436db609ea7SPeter Brune Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325. 437e55accfeSBarryFSmith 438e55accfeSBarryFSmith This line search will always produce a step that is less than or equal to, in length, the full step size. 439db609ea7SPeter Brune 440db781477SPatrick Sanan .seealso: `SNESLineSearchCreate()`, `SNESLineSearchSetType()` 441954494b2SJed Brown M*/ 4428cc058d9SJed Brown PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch) 44371b4ebd2SPeter Brune { 44471b4ebd2SPeter Brune 445f1c6b773SPeter Brune SNESLineSearch_BT *bt; 44671b4ebd2SPeter Brune 44771b4ebd2SPeter Brune PetscFunctionBegin; 448f1c6b773SPeter Brune linesearch->ops->apply = SNESLineSearchApply_BT; 449f1c6b773SPeter Brune linesearch->ops->destroy = SNESLineSearchDestroy_BT; 450f1c6b773SPeter Brune linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT; 4510298fd71SBarry Smith linesearch->ops->reset = NULL; 4527f1410a3SPeter Brune linesearch->ops->view = SNESLineSearchView_BT; 4530298fd71SBarry Smith linesearch->ops->setup = NULL; 45471b4ebd2SPeter Brune 4559566063dSJacob Faibussowitsch PetscCall(PetscNewLog(linesearch,&bt)); 456f5af7f23SKarl Rupp 45771b4ebd2SPeter Brune linesearch->data = (void*)bt; 45871b4ebd2SPeter Brune linesearch->max_its = 40; 459b000cd8dSPeter Brune linesearch->order = SNES_LINESEARCH_ORDER_CUBIC; 46071b4ebd2SPeter Brune bt->alpha = 1e-4; 46171b4ebd2SPeter Brune PetscFunctionReturn(0); 46271b4ebd2SPeter Brune } 463