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 171a4f838cSPeter Brune .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 401a4f838cSPeter Brune .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 PetscErrorCode ierr; 5671b4ebd2SPeter Brune Vec X,F,Y,W,G; 5771b4ebd2SPeter Brune SNES snes; 588a430afbSPeter Brune PetscReal fnorm, xnorm, ynorm, gnorm; 593b288129SPeter Brune PetscReal lambda,lambdatemp,lambdaprev,minlambda,maxstep,initslope,alpha,stol; 60dfcd3828SPeter Brune PetscReal t1,t2,a,b,d; 618a430afbSPeter Brune PetscReal f; 628a430afbSPeter Brune PetscReal g,gprev; 6371b4ebd2SPeter Brune PetscViewer monitor; 6471b4ebd2SPeter Brune PetscInt max_its,count; 65d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data; 6671b4ebd2SPeter Brune Mat jac; 67bd4a8a71SBarry Smith PetscErrorCode (*objective)(SNES,Vec,PetscReal*,void*); 6871b4ebd2SPeter Brune 6971b4ebd2SPeter Brune PetscFunctionBegin; 709566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G)); 719566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm)); 729566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetLambda(linesearch, &lambda)); 739566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetSNES(linesearch, &snes)); 749566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor)); 759566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetTolerances(linesearch,&minlambda,&maxstep,NULL,NULL,NULL,&max_its)); 769566063dSJacob Faibussowitsch PetscCall(SNESGetTolerances(snes,NULL,NULL,&stol,NULL,NULL)); 779566063dSJacob Faibussowitsch PetscCall(SNESGetObjective(snes,&objective,NULL)); 7871b4ebd2SPeter Brune alpha = bt->alpha; 7971b4ebd2SPeter Brune 809566063dSJacob Faibussowitsch PetscCall(SNESGetJacobian(snes, &jac, NULL, NULL, NULL)); 81*08401ef6SPierre Jolivet PetscCheck(jac || objective,PetscObjectComm((PetscObject)linesearch), PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix"); 82c69d1a72SBarry Smith 839566063dSJacob Faibussowitsch PetscCall(SNESLineSearchPreCheck(linesearch,X,Y,&changed_y)); 849566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 8571b4ebd2SPeter Brune 869566063dSJacob Faibussowitsch PetscCall(VecNormBegin(Y, NORM_2, &ynorm)); 879566063dSJacob Faibussowitsch PetscCall(VecNormBegin(X, NORM_2, &xnorm)); 889566063dSJacob Faibussowitsch PetscCall(VecNormEnd(Y, NORM_2, &ynorm)); 899566063dSJacob Faibussowitsch PetscCall(VecNormEnd(X, NORM_2, &xnorm)); 903b288129SPeter Brune 9171b4ebd2SPeter Brune if (ynorm == 0.0) { 9271b4ebd2SPeter Brune if (monitor) { 939566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 949566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Initial direction and size is 0\n")); 959566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 9671b4ebd2SPeter Brune } 979566063dSJacob Faibussowitsch PetscCall(VecCopy(X,W)); 989566063dSJacob Faibussowitsch PetscCall(VecCopy(F,G)); 999566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm)); 1009566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 10171b4ebd2SPeter Brune PetscFunctionReturn(0); 10271b4ebd2SPeter Brune } 10371b4ebd2SPeter Brune if (ynorm > maxstep) { /* Step too big, so scale back */ 10471b4ebd2SPeter Brune if (monitor) { 1059566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 1069566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Scaling step by %14.12e old ynorm %14.12e\n", (double)(maxstep/ynorm),(double)ynorm)); 1079566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 10871b4ebd2SPeter Brune } 1099566063dSJacob Faibussowitsch PetscCall(VecScale(Y,maxstep/(ynorm))); 11071b4ebd2SPeter Brune ynorm = maxstep; 11171b4ebd2SPeter Brune } 1128a430afbSPeter Brune 1138a430afbSPeter Brune /* if the SNES has an objective set, use that instead of the function value */ 114bd4a8a71SBarry Smith if (objective) { 1159566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes,X,&f)); 1168a430afbSPeter Brune } else { 1178a430afbSPeter Brune f = fnorm*fnorm; 1188a430afbSPeter Brune } 1198a430afbSPeter Brune 1208a430afbSPeter Brune /* compute the initial slope */ 121bd4a8a71SBarry Smith if (objective) { 1228a430afbSPeter Brune /* slope comes from the function (assumed to be the gradient of the objective */ 1239566063dSJacob Faibussowitsch PetscCall(VecDotRealPart(Y,F,&initslope)); 1248a430afbSPeter Brune } else { 1258a430afbSPeter Brune /* slope comes from the normal equations */ 1269566063dSJacob Faibussowitsch PetscCall(MatMult(jac,Y,W)); 1279566063dSJacob Faibussowitsch PetscCall(VecDotRealPart(F,W,&initslope)); 12871b4ebd2SPeter Brune if (initslope > 0.0) initslope = -initslope; 12971b4ebd2SPeter Brune if (initslope == 0.0) initslope = -1.0; 1308a430afbSPeter Brune } 13171b4ebd2SPeter Brune 132df019d78SBarry Smith while (PETSC_TRUE) { 1339566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W,-lambda,Y,X)); 1349bd66eb0SPeter Brune if (linesearch->ops->viproject) { 1359566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->viproject)(snes, W)); 1369bd66eb0SPeter Brune } 137e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 1389566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n")); 13971b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1409566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 14171b4ebd2SPeter Brune PetscFunctionReturn(0); 14271b4ebd2SPeter Brune } 1438a430afbSPeter Brune 144bd4a8a71SBarry Smith if (objective) { 1459566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes,W,&g)); 1468a430afbSPeter Brune } else { 1479566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes,W,G)); 1489bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 1499bd66eb0SPeter Brune gnorm = fnorm; 1509566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 1519bd66eb0SPeter Brune } else { 1529566063dSJacob Faibussowitsch PetscCall(VecNorm(G,NORM_2,&gnorm)); 1538a430afbSPeter Brune } 15479e7e81bSJed Brown g = PetscSqr(gnorm); 1559bd66eb0SPeter Brune } 1569566063dSJacob Faibussowitsch PetscCall(SNESLineSearchMonitor(linesearch)); 1579bd66eb0SPeter Brune 158df019d78SBarry Smith if (!PetscIsInfOrNanReal(g)) break; 159df019d78SBarry Smith if (monitor) { 1609566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 1619566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n",(double)lambda)); 1629566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 163df019d78SBarry Smith } 164df019d78SBarry Smith if (lambda <= minlambda) { 1657e49c775SBarry Smith SNESCheckFunctionNorm(snes,g); 166c98378a5SBarry Smith } 167df019d78SBarry Smith lambda = .5*lambda; 168df019d78SBarry Smith } 169df019d78SBarry Smith 170bd4a8a71SBarry Smith if (!objective) { 1719566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm)); 1728a430afbSPeter Brune } 1738a430afbSPeter Brune if (.5*g <= .5*f + lambda*alpha*initslope) { /* Sufficient reduction or step tolerance convergence */ 17471b4ebd2SPeter Brune if (monitor) { 1759566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 176bd4a8a71SBarry Smith if (!objective) { 1779566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm)); 1788a430afbSPeter Brune } else { 1799566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Using full step: obj %14.12e obj %14.12e\n", (double)f, (double)g)); 1808a430afbSPeter Brune } 1819566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 18271b4ebd2SPeter Brune } 18371b4ebd2SPeter Brune } else { 184c21ba15cSPeter Brune /* Since the full step didn't work and the step is tiny, quit */ 185c21ba15cSPeter Brune if (stol*xnorm > ynorm) { 1869566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm)); 1879566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 188c61ba1e2SPeter Brune if (monitor) { 1899566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 1909566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n",(double)ynorm,(double)stol*xnorm)); 1919566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 192c61ba1e2SPeter Brune } 193c21ba15cSPeter Brune PetscFunctionReturn(0); 194c21ba15cSPeter Brune } 19571b4ebd2SPeter Brune /* Fit points with quadratic */ 1968a430afbSPeter Brune lambdatemp = -initslope/(g - f - 2.0*lambda*initslope); 19771b4ebd2SPeter Brune lambdaprev = lambda; 1988a430afbSPeter Brune gprev = g; 19971b4ebd2SPeter Brune if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 20071b4ebd2SPeter Brune if (lambdatemp <= .1*lambda) lambda = .1*lambda; 20171b4ebd2SPeter Brune else lambda = lambdatemp; 20271b4ebd2SPeter Brune 2039566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W,-lambda,Y,X)); 2049bd66eb0SPeter Brune if (linesearch->ops->viproject) { 2059566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->viproject)(snes, W)); 2069bd66eb0SPeter Brune } 207e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 2089566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs)); 20971b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2109566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 21171b4ebd2SPeter Brune PetscFunctionReturn(0); 21271b4ebd2SPeter Brune } 213bd4a8a71SBarry Smith if (objective) { 2149566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes,W,&g)); 2158a430afbSPeter Brune } else { 2169566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes,W,G)); 2179bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 2189bd66eb0SPeter Brune gnorm = fnorm; 2199566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 2209bd66eb0SPeter Brune } else { 2219566063dSJacob Faibussowitsch PetscCall(VecNorm(G,NORM_2,&gnorm)); 2229bd66eb0SPeter Brune } 223f1ab69e3SPeter Brune g = PetscSqr(gnorm); 2248a430afbSPeter Brune } 225c98378a5SBarry Smith if (PetscIsInfOrNanReal(g)) { 2269566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 2279566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n")); 228c98378a5SBarry Smith PetscFunctionReturn(0); 229c98378a5SBarry Smith } 23071b4ebd2SPeter Brune if (monitor) { 2319566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 232bd4a8a71SBarry Smith if (!objective) { 2339566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: gnorm after quadratic fit %14.12e\n",(double)gnorm)); 2348a430afbSPeter Brune } else { 2359566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: obj after quadratic fit %14.12e\n",(double)g)); 2368a430afbSPeter Brune } 2379566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 23871b4ebd2SPeter Brune } 2398a430afbSPeter Brune if (.5*g < .5*f + lambda*alpha*initslope) { /* sufficient reduction */ 24071b4ebd2SPeter Brune if (monitor) { 2419566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 2429566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda)); 2439566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 24471b4ebd2SPeter Brune } 24571b4ebd2SPeter Brune } else { 24671b4ebd2SPeter Brune /* Fit points with cubic */ 24771b4ebd2SPeter Brune for (count = 0; count < max_its; count++) { 24871b4ebd2SPeter Brune if (lambda <= minlambda) { 24971b4ebd2SPeter Brune if (monitor) { 2509566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 2519566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: unable to find good step length! After %D tries \n",count)); 252bd4a8a71SBarry Smith if (!objective) { 2534a2f8832SBarry Smith ierr = PetscViewerASCIIPrintf(monitor," Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", 2549566063dSJacob Faibussowitsch (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);PetscCall(ierr); 2558a430afbSPeter Brune } else { 2564a2f8832SBarry Smith ierr = 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", 2579566063dSJacob Faibussowitsch (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);PetscCall(ierr); 2588a430afbSPeter Brune } 2599566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 26071b4ebd2SPeter Brune } 2619566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 26271b4ebd2SPeter Brune PetscFunctionReturn(0); 26371b4ebd2SPeter Brune } 264b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 2658a430afbSPeter Brune t1 = .5*(g - f) - lambda*initslope; 2668a430afbSPeter Brune t2 = .5*(gprev - f) - lambdaprev*initslope; 26771b4ebd2SPeter Brune a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 26871b4ebd2SPeter Brune b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 26971b4ebd2SPeter Brune d = b*b - 3*a*initslope; 27071b4ebd2SPeter Brune if (d < 0.0) d = 0.0; 271f5af7f23SKarl Rupp if (a == 0.0) lambdatemp = -initslope/(2.0*b); 272f5af7f23SKarl Rupp else lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a); 273f5af7f23SKarl Rupp 274b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 2758a430afbSPeter Brune lambdatemp = -initslope/(g - f - 2.0*initslope); 276ce94432eSBarry Smith } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt"); 27771b4ebd2SPeter Brune lambdaprev = lambda; 2788a430afbSPeter Brune gprev = g; 27971b4ebd2SPeter Brune if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 28071b4ebd2SPeter Brune if (lambdatemp <= .1*lambda) lambda = .1*lambda; 28171b4ebd2SPeter Brune else lambda = lambdatemp; 2829566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W,-lambda,Y,X)); 283b7b2e573SPeter Brune if (linesearch->ops->viproject) { 2849566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->viproject)(snes,W)); 285b7b2e573SPeter Brune } 286e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 2879566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count)); 288bd4a8a71SBarry Smith if (!objective) { 2897d3de750SJacob Faibussowitsch ierr = PetscInfo(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", 2909566063dSJacob Faibussowitsch (double)fnorm,(double)gnorm,(double)ynorm,(double)lambda,(double)initslope);PetscCall(ierr); 2918a430afbSPeter Brune } 2929566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 29371b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 29471b4ebd2SPeter Brune PetscFunctionReturn(0); 29571b4ebd2SPeter Brune } 296bd4a8a71SBarry Smith if (objective) { 2979566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes,W,&g)); 2988a430afbSPeter Brune } else { 2999566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes,W,G)); 3009bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 3019bd66eb0SPeter Brune gnorm = fnorm; 3029566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 3039bd66eb0SPeter Brune } else { 3049566063dSJacob Faibussowitsch PetscCall(VecNorm(G,NORM_2,&gnorm)); 3058a430afbSPeter Brune } 306f1ab69e3SPeter Brune g = PetscSqr(gnorm); 3079bd66eb0SPeter Brune } 3084a2f8832SBarry Smith if (PetscIsInfOrNanReal(g)) { 3099566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 3109566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n")); 311c98378a5SBarry Smith PetscFunctionReturn(0); 312c98378a5SBarry Smith } 3138a430afbSPeter Brune if (.5*g < .5*f + lambda*alpha*initslope) { /* is reduction enough? */ 31471b4ebd2SPeter Brune if (monitor) { 3159566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 316bd4a8a71SBarry Smith if (!objective) { 317b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3189566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda)); 31971b4ebd2SPeter Brune } else { 3209566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda)); 32171b4ebd2SPeter Brune } 3229566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 3238a430afbSPeter Brune } else { 3248a430afbSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3259566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda)); 3268a430afbSPeter Brune } else { 3279566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda)); 3288a430afbSPeter Brune } 3299566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 3308a430afbSPeter Brune } 33171b4ebd2SPeter Brune } 33271b4ebd2SPeter Brune break; 333f5af7f23SKarl Rupp } else if (monitor) { 3349566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 335bd4a8a71SBarry Smith if (!objective) { 336b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3379566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda)); 33871b4ebd2SPeter Brune } else { 3399566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda)); 34071b4ebd2SPeter Brune } 3419566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 3428a430afbSPeter Brune } else { 3438a430afbSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3449566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda)); 3458a430afbSPeter Brune } else { 3469566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda)); 3478a430afbSPeter Brune } 3489566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 3498a430afbSPeter Brune } 35071b4ebd2SPeter Brune } 35171b4ebd2SPeter Brune } 35271b4ebd2SPeter Brune } 35371b4ebd2SPeter Brune } 35471b4ebd2SPeter Brune 35571b4ebd2SPeter Brune /* postcheck */ 3560c1ac483SGlenn Hammond /* update Y to lambda*Y so that W is consistent with X - lambda*Y */ 3579566063dSJacob Faibussowitsch PetscCall(VecScale(Y,lambda)); 3589566063dSJacob Faibussowitsch PetscCall(SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w)); 35971b4ebd2SPeter Brune if (changed_y) { 3609566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W,-1.0,Y,X)); 3619bd66eb0SPeter Brune if (linesearch->ops->viproject) { 3629566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->viproject)(snes, W)); 3639bd66eb0SPeter Brune } 36471b4ebd2SPeter Brune } 365bd4a8a71SBarry Smith if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */ 3669566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes,W,G)); 3679bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 3689bd66eb0SPeter Brune gnorm = fnorm; 3699566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 3709bd66eb0SPeter Brune } else { 3719566063dSJacob Faibussowitsch PetscCall(VecNorm(G,NORM_2,&gnorm)); 3729bd66eb0SPeter Brune } 3739566063dSJacob Faibussowitsch PetscCall(VecNorm(Y,NORM_2,&ynorm)); 374c98378a5SBarry Smith if (PetscIsInfOrNanReal(gnorm)) { 3759566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_FAILED_NANORINF)); 3769566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n")); 377c98378a5SBarry Smith PetscFunctionReturn(0); 378c98378a5SBarry Smith } 379c427506bSPeter Brune } 38071b4ebd2SPeter Brune 38171b4ebd2SPeter Brune /* copy the solution over */ 3829566063dSJacob Faibussowitsch PetscCall(VecCopy(W, X)); 3839566063dSJacob Faibussowitsch PetscCall(VecCopy(G, F)); 3849566063dSJacob Faibussowitsch PetscCall(VecNorm(X, NORM_2, &xnorm)); 3859566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetLambda(linesearch, lambda)); 3869566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm)); 38771b4ebd2SPeter Brune PetscFunctionReturn(0); 38871b4ebd2SPeter Brune } 38971b4ebd2SPeter Brune 3907f1410a3SPeter Brune PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer) 3917f1410a3SPeter Brune { 3927f1410a3SPeter Brune PetscBool iascii; 393d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data; 394075cc632SBarry Smith 3957f1410a3SPeter Brune PetscFunctionBegin; 3969566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 3977f1410a3SPeter Brune if (iascii) { 398b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3999566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: cubic\n")); 400b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 4019566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: quadratic\n")); 4027f1410a3SPeter Brune } 4039566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " alpha=%e\n", (double)bt->alpha)); 4047f1410a3SPeter Brune } 4057f1410a3SPeter Brune PetscFunctionReturn(0); 4067f1410a3SPeter Brune } 4077f1410a3SPeter Brune 408f1c6b773SPeter Brune static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch) 409731c067cSJed Brown { 41071b4ebd2SPeter Brune PetscFunctionBegin; 4119566063dSJacob Faibussowitsch PetscCall(PetscFree(linesearch->data)); 41271b4ebd2SPeter Brune PetscFunctionReturn(0); 41371b4ebd2SPeter Brune } 41471b4ebd2SPeter Brune 4154416b707SBarry Smith static PetscErrorCode SNESLineSearchSetFromOptions_BT(PetscOptionItems *PetscOptionsObject,SNESLineSearch linesearch) 41671b4ebd2SPeter Brune { 417eb23ba39SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data; 41871b4ebd2SPeter Brune 4196e111a19SKarl Rupp PetscFunctionBegin; 4209566063dSJacob Faibussowitsch PetscCall(PetscOptionsHead(PetscOptionsObject,"SNESLineSearch BT options")); 4219566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL)); 4229566063dSJacob Faibussowitsch PetscCall(PetscOptionsTail()); 42371b4ebd2SPeter Brune PetscFunctionReturn(0); 42471b4ebd2SPeter Brune } 42571b4ebd2SPeter Brune 426954494b2SJed Brown /*MC 427db609ea7SPeter Brune SNESLINESEARCHBT - Backtracking line search. 428954494b2SJed Brown 429db609ea7SPeter Brune This line search finds the minimum of a polynomial fitting of the L2 norm of the 430eb23ba39SBarry 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 431db609ea7SPeter Brune and the fit is reattempted at most max_it times or until lambda is below minlambda. 432954494b2SJed Brown 433cd7522eaSPeter Brune Options Database Keys: 4343eec2f6aSPierre Jolivet + -snes_linesearch_alpha <1e\-4> - slope descent parameter 435db609ea7SPeter Brune . -snes_linesearch_damping <1.0> - initial step length 436eb23ba39SBarry Smith . -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the 437e1bc860dSBarry Smith step is scaled back to be of this length at the beginning of the line search 438db609ea7SPeter Brune . -snes_linesearch_max_it <40> - maximum number of shrinking step 4393eec2f6aSPierre Jolivet . -snes_linesearch_minlambda <1e\-12> - minimum step length allowed 440cd7522eaSPeter Brune - -snes_linesearch_order <cubic,quadratic> - order of the approximation 441cd7522eaSPeter Brune 442954494b2SJed Brown Level: advanced 443954494b2SJed Brown 444db609ea7SPeter Brune Notes: 445db609ea7SPeter Brune This line search is taken from "Numerical Methods for Unconstrained 446db609ea7SPeter Brune Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325. 447e55accfeSBarryFSmith 448e55accfeSBarryFSmith This line search will always produce a step that is less than or equal to, in length, the full step size. 449db609ea7SPeter Brune 450f1c6b773SPeter Brune .seealso: SNESLineSearchCreate(), SNESLineSearchSetType() 451954494b2SJed Brown M*/ 4528cc058d9SJed Brown PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch) 45371b4ebd2SPeter Brune { 45471b4ebd2SPeter Brune 455f1c6b773SPeter Brune SNESLineSearch_BT *bt; 45671b4ebd2SPeter Brune 45771b4ebd2SPeter Brune PetscFunctionBegin; 458f1c6b773SPeter Brune linesearch->ops->apply = SNESLineSearchApply_BT; 459f1c6b773SPeter Brune linesearch->ops->destroy = SNESLineSearchDestroy_BT; 460f1c6b773SPeter Brune linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT; 4610298fd71SBarry Smith linesearch->ops->reset = NULL; 4627f1410a3SPeter Brune linesearch->ops->view = SNESLineSearchView_BT; 4630298fd71SBarry Smith linesearch->ops->setup = NULL; 46471b4ebd2SPeter Brune 4659566063dSJacob Faibussowitsch PetscCall(PetscNewLog(linesearch,&bt)); 466f5af7f23SKarl Rupp 46771b4ebd2SPeter Brune linesearch->data = (void*)bt; 46871b4ebd2SPeter Brune linesearch->max_its = 40; 469b000cd8dSPeter Brune linesearch->order = SNES_LINESEARCH_ORDER_CUBIC; 47071b4ebd2SPeter Brune bt->alpha = 1e-4; 47171b4ebd2SPeter Brune PetscFunctionReturn(0); 47271b4ebd2SPeter Brune } 473