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 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)); 1339bd66eb0SPeter Brune if (linesearch->ops->viproject) { 1349566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->viproject)(snes, W)); 1359bd66eb0SPeter Brune } 136e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 1379566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n")); 13871b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1399566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 14071b4ebd2SPeter Brune PetscFunctionReturn(0); 14171b4ebd2SPeter Brune } 1428a430afbSPeter Brune 143bd4a8a71SBarry Smith if (objective) { 1449566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes,W,&g)); 1458a430afbSPeter Brune } else { 1469566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes,W,G)); 1479bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 1489bd66eb0SPeter Brune gnorm = fnorm; 1499566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 1509bd66eb0SPeter Brune } else { 1519566063dSJacob Faibussowitsch PetscCall(VecNorm(G,NORM_2,&gnorm)); 1528a430afbSPeter Brune } 15379e7e81bSJed Brown g = PetscSqr(gnorm); 1549bd66eb0SPeter Brune } 1559566063dSJacob Faibussowitsch PetscCall(SNESLineSearchMonitor(linesearch)); 1569bd66eb0SPeter Brune 157df019d78SBarry Smith if (!PetscIsInfOrNanReal(g)) break; 158df019d78SBarry Smith if (monitor) { 1599566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 1609566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n",(double)lambda)); 1619566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 162df019d78SBarry Smith } 163df019d78SBarry Smith if (lambda <= minlambda) { 1647e49c775SBarry Smith SNESCheckFunctionNorm(snes,g); 165c98378a5SBarry Smith } 166df019d78SBarry Smith lambda = .5*lambda; 167df019d78SBarry Smith } 168df019d78SBarry Smith 169bd4a8a71SBarry Smith if (!objective) { 1709566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm)); 1718a430afbSPeter Brune } 1728a430afbSPeter Brune if (.5*g <= .5*f + lambda*alpha*initslope) { /* Sufficient reduction or step tolerance convergence */ 17371b4ebd2SPeter Brune if (monitor) { 1749566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 175bd4a8a71SBarry Smith if (!objective) { 1769566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm)); 1778a430afbSPeter Brune } else { 1789566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Using full step: obj %14.12e obj %14.12e\n", (double)f, (double)g)); 1798a430afbSPeter Brune } 1809566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 18171b4ebd2SPeter Brune } 18271b4ebd2SPeter Brune } else { 183c21ba15cSPeter Brune /* Since the full step didn't work and the step is tiny, quit */ 184c21ba15cSPeter Brune if (stol*xnorm > ynorm) { 1859566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch,xnorm,fnorm,ynorm)); 1869566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 187c61ba1e2SPeter Brune if (monitor) { 1889566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 189*63a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n",(double)ynorm,(double)(stol*xnorm))); 1909566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 191c61ba1e2SPeter Brune } 192c21ba15cSPeter Brune PetscFunctionReturn(0); 193c21ba15cSPeter Brune } 19471b4ebd2SPeter Brune /* Fit points with quadratic */ 1958a430afbSPeter Brune lambdatemp = -initslope/(g - f - 2.0*lambda*initslope); 19671b4ebd2SPeter Brune lambdaprev = lambda; 1978a430afbSPeter Brune gprev = g; 19871b4ebd2SPeter Brune if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 19971b4ebd2SPeter Brune if (lambdatemp <= .1*lambda) lambda = .1*lambda; 20071b4ebd2SPeter Brune else lambda = lambdatemp; 20171b4ebd2SPeter Brune 2029566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W,-lambda,Y,X)); 2039bd66eb0SPeter Brune if (linesearch->ops->viproject) { 2049566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->viproject)(snes, W)); 2059bd66eb0SPeter Brune } 206e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 207*63a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %" PetscInt_FMT " \n",snes->nfuncs)); 20871b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2099566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 21071b4ebd2SPeter Brune PetscFunctionReturn(0); 21171b4ebd2SPeter Brune } 212bd4a8a71SBarry Smith if (objective) { 2139566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes,W,&g)); 2148a430afbSPeter Brune } else { 2159566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes,W,G)); 2169bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 2179bd66eb0SPeter Brune gnorm = fnorm; 2189566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 2199bd66eb0SPeter Brune } else { 2209566063dSJacob Faibussowitsch PetscCall(VecNorm(G,NORM_2,&gnorm)); 2219bd66eb0SPeter Brune } 222f1ab69e3SPeter Brune g = PetscSqr(gnorm); 2238a430afbSPeter Brune } 224c98378a5SBarry Smith if (PetscIsInfOrNanReal(g)) { 2259566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 2269566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n")); 227c98378a5SBarry Smith PetscFunctionReturn(0); 228c98378a5SBarry Smith } 22971b4ebd2SPeter Brune if (monitor) { 2309566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 231bd4a8a71SBarry Smith if (!objective) { 2329566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: gnorm after quadratic fit %14.12e\n",(double)gnorm)); 2338a430afbSPeter Brune } else { 2349566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: obj after quadratic fit %14.12e\n",(double)g)); 2358a430afbSPeter Brune } 2369566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 23771b4ebd2SPeter Brune } 2388a430afbSPeter Brune if (.5*g < .5*f + lambda*alpha*initslope) { /* sufficient reduction */ 23971b4ebd2SPeter Brune if (monitor) { 2409566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 2419566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda)); 2429566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 24371b4ebd2SPeter Brune } 24471b4ebd2SPeter Brune } else { 24571b4ebd2SPeter Brune /* Fit points with cubic */ 24671b4ebd2SPeter Brune for (count = 0; count < max_its; count++) { 24771b4ebd2SPeter Brune if (lambda <= minlambda) { 24871b4ebd2SPeter Brune if (monitor) { 2499566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 250*63a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: unable to find good step length! After %" PetscInt_FMT " tries \n",count)); 251bd4a8a71SBarry Smith if (!objective) { 252d0609cedSBarry 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", 253d0609cedSBarry Smith (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope)); 2548a430afbSPeter Brune } else { 255d0609cedSBarry 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", 256d0609cedSBarry Smith (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope)); 2578a430afbSPeter Brune } 2589566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 25971b4ebd2SPeter Brune } 2609566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 26171b4ebd2SPeter Brune PetscFunctionReturn(0); 26271b4ebd2SPeter Brune } 263b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 2648a430afbSPeter Brune t1 = .5*(g - f) - lambda*initslope; 2658a430afbSPeter Brune t2 = .5*(gprev - f) - lambdaprev*initslope; 26671b4ebd2SPeter Brune a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 26771b4ebd2SPeter Brune b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 26871b4ebd2SPeter Brune d = b*b - 3*a*initslope; 26971b4ebd2SPeter Brune if (d < 0.0) d = 0.0; 270f5af7f23SKarl Rupp if (a == 0.0) lambdatemp = -initslope/(2.0*b); 271f5af7f23SKarl Rupp else lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a); 272f5af7f23SKarl Rupp 273b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 2748a430afbSPeter Brune lambdatemp = -initslope/(g - f - 2.0*initslope); 275ce94432eSBarry Smith } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt"); 27671b4ebd2SPeter Brune lambdaprev = lambda; 2778a430afbSPeter Brune gprev = g; 27871b4ebd2SPeter Brune if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 27971b4ebd2SPeter Brune if (lambdatemp <= .1*lambda) lambda = .1*lambda; 28071b4ebd2SPeter Brune else lambda = lambdatemp; 2819566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W,-lambda,Y,X)); 282b7b2e573SPeter Brune if (linesearch->ops->viproject) { 2839566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->viproject)(snes,W)); 284b7b2e573SPeter Brune } 285e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 286*63a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Exceeded maximum function evaluations, while looking for good step length! %" PetscInt_FMT " \n",count)); 287bd4a8a71SBarry Smith if (!objective) { 288d0609cedSBarry 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)); 2898a430afbSPeter Brune } 2909566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 29171b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 29271b4ebd2SPeter Brune PetscFunctionReturn(0); 29371b4ebd2SPeter Brune } 294bd4a8a71SBarry Smith if (objective) { 2959566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes,W,&g)); 2968a430afbSPeter Brune } else { 2979566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes,W,G)); 2989bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 2999bd66eb0SPeter Brune gnorm = fnorm; 3009566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 3019bd66eb0SPeter Brune } else { 3029566063dSJacob Faibussowitsch PetscCall(VecNorm(G,NORM_2,&gnorm)); 3038a430afbSPeter Brune } 304f1ab69e3SPeter Brune g = PetscSqr(gnorm); 3059bd66eb0SPeter Brune } 3064a2f8832SBarry Smith if (PetscIsInfOrNanReal(g)) { 3079566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 3089566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n")); 309c98378a5SBarry Smith PetscFunctionReturn(0); 310c98378a5SBarry Smith } 3118a430afbSPeter Brune if (.5*g < .5*f + lambda*alpha*initslope) { /* is reduction enough? */ 31271b4ebd2SPeter Brune if (monitor) { 3139566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 314bd4a8a71SBarry Smith if (!objective) { 315b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3169566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda)); 31771b4ebd2SPeter Brune } else { 3189566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda)); 31971b4ebd2SPeter Brune } 3209566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 3218a430afbSPeter Brune } else { 3228a430afbSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3239566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda)); 3248a430afbSPeter Brune } else { 3259566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n",(double)g,(double)lambda)); 3268a430afbSPeter Brune } 3279566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 3288a430afbSPeter Brune } 32971b4ebd2SPeter Brune } 33071b4ebd2SPeter Brune break; 331f5af7f23SKarl Rupp } else if (monitor) { 3329566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 333bd4a8a71SBarry Smith if (!objective) { 334b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3359566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda)); 33671b4ebd2SPeter Brune } else { 3379566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda)); 33871b4ebd2SPeter Brune } 3399566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 3408a430afbSPeter Brune } else { 3418a430afbSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3429566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda)); 3438a430afbSPeter Brune } else { 3449566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n",(double)g,(double)lambda)); 3458a430afbSPeter Brune } 3469566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 3478a430afbSPeter Brune } 34871b4ebd2SPeter Brune } 34971b4ebd2SPeter Brune } 35071b4ebd2SPeter Brune } 35171b4ebd2SPeter Brune } 35271b4ebd2SPeter Brune 35371b4ebd2SPeter Brune /* postcheck */ 3540c1ac483SGlenn Hammond /* update Y to lambda*Y so that W is consistent with X - lambda*Y */ 3559566063dSJacob Faibussowitsch PetscCall(VecScale(Y,lambda)); 3569566063dSJacob Faibussowitsch PetscCall(SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w)); 35771b4ebd2SPeter Brune if (changed_y) { 3589566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W,-1.0,Y,X)); 3599bd66eb0SPeter Brune if (linesearch->ops->viproject) { 3609566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->viproject)(snes, W)); 3619bd66eb0SPeter Brune } 36271b4ebd2SPeter Brune } 363bd4a8a71SBarry Smith if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */ 3649566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes,W,G)); 3659bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 3669bd66eb0SPeter Brune gnorm = fnorm; 3679566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 3689bd66eb0SPeter Brune } else { 3699566063dSJacob Faibussowitsch PetscCall(VecNorm(G,NORM_2,&gnorm)); 3709bd66eb0SPeter Brune } 3719566063dSJacob Faibussowitsch PetscCall(VecNorm(Y,NORM_2,&ynorm)); 372c98378a5SBarry Smith if (PetscIsInfOrNanReal(gnorm)) { 3739566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_FAILED_NANORINF)); 3749566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes,"Aborted due to Nan or Inf in function evaluation\n")); 375c98378a5SBarry Smith PetscFunctionReturn(0); 376c98378a5SBarry Smith } 377c427506bSPeter Brune } 37871b4ebd2SPeter Brune 37971b4ebd2SPeter Brune /* copy the solution over */ 3809566063dSJacob Faibussowitsch PetscCall(VecCopy(W, X)); 3819566063dSJacob Faibussowitsch PetscCall(VecCopy(G, F)); 3829566063dSJacob Faibussowitsch PetscCall(VecNorm(X, NORM_2, &xnorm)); 3839566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetLambda(linesearch, lambda)); 3849566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm)); 38571b4ebd2SPeter Brune PetscFunctionReturn(0); 38671b4ebd2SPeter Brune } 38771b4ebd2SPeter Brune 3887f1410a3SPeter Brune PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer) 3897f1410a3SPeter Brune { 3907f1410a3SPeter Brune PetscBool iascii; 391d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data; 392075cc632SBarry Smith 3937f1410a3SPeter Brune PetscFunctionBegin; 3949566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 3957f1410a3SPeter Brune if (iascii) { 396b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3979566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: cubic\n")); 398b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 3999566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: quadratic\n")); 4007f1410a3SPeter Brune } 4019566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " alpha=%e\n", (double)bt->alpha)); 4027f1410a3SPeter Brune } 4037f1410a3SPeter Brune PetscFunctionReturn(0); 4047f1410a3SPeter Brune } 4057f1410a3SPeter Brune 406f1c6b773SPeter Brune static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch) 407731c067cSJed Brown { 40871b4ebd2SPeter Brune PetscFunctionBegin; 4099566063dSJacob Faibussowitsch PetscCall(PetscFree(linesearch->data)); 41071b4ebd2SPeter Brune PetscFunctionReturn(0); 41171b4ebd2SPeter Brune } 41271b4ebd2SPeter Brune 4134416b707SBarry Smith static PetscErrorCode SNESLineSearchSetFromOptions_BT(PetscOptionItems *PetscOptionsObject,SNESLineSearch linesearch) 41471b4ebd2SPeter Brune { 415eb23ba39SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT*)linesearch->data; 41671b4ebd2SPeter Brune 4176e111a19SKarl Rupp PetscFunctionBegin; 418d0609cedSBarry Smith PetscOptionsHeadBegin(PetscOptionsObject,"SNESLineSearch BT options"); 4199566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL)); 420d0609cedSBarry Smith PetscOptionsHeadEnd(); 42171b4ebd2SPeter Brune PetscFunctionReturn(0); 42271b4ebd2SPeter Brune } 42371b4ebd2SPeter Brune 424954494b2SJed Brown /*MC 425db609ea7SPeter Brune SNESLINESEARCHBT - Backtracking line search. 426954494b2SJed Brown 427db609ea7SPeter Brune This line search finds the minimum of a polynomial fitting of the L2 norm of the 428eb23ba39SBarry 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 429db609ea7SPeter Brune and the fit is reattempted at most max_it times or until lambda is below minlambda. 430954494b2SJed Brown 431cd7522eaSPeter Brune Options Database Keys: 4323eec2f6aSPierre Jolivet + -snes_linesearch_alpha <1e\-4> - slope descent parameter 433db609ea7SPeter Brune . -snes_linesearch_damping <1.0> - initial step length 434eb23ba39SBarry Smith . -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the 435e1bc860dSBarry Smith step is scaled back to be of this length at the beginning of the line search 436db609ea7SPeter Brune . -snes_linesearch_max_it <40> - maximum number of shrinking step 4373eec2f6aSPierre Jolivet . -snes_linesearch_minlambda <1e\-12> - minimum step length allowed 438cd7522eaSPeter Brune - -snes_linesearch_order <cubic,quadratic> - order of the approximation 439cd7522eaSPeter Brune 440954494b2SJed Brown Level: advanced 441954494b2SJed Brown 442db609ea7SPeter Brune Notes: 443db609ea7SPeter Brune This line search is taken from "Numerical Methods for Unconstrained 444db609ea7SPeter Brune Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325. 445e55accfeSBarryFSmith 446e55accfeSBarryFSmith This line search will always produce a step that is less than or equal to, in length, the full step size. 447db609ea7SPeter Brune 448f1c6b773SPeter Brune .seealso: SNESLineSearchCreate(), SNESLineSearchSetType() 449954494b2SJed Brown M*/ 4508cc058d9SJed Brown PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch) 45171b4ebd2SPeter Brune { 45271b4ebd2SPeter Brune 453f1c6b773SPeter Brune SNESLineSearch_BT *bt; 45471b4ebd2SPeter Brune 45571b4ebd2SPeter Brune PetscFunctionBegin; 456f1c6b773SPeter Brune linesearch->ops->apply = SNESLineSearchApply_BT; 457f1c6b773SPeter Brune linesearch->ops->destroy = SNESLineSearchDestroy_BT; 458f1c6b773SPeter Brune linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT; 4590298fd71SBarry Smith linesearch->ops->reset = NULL; 4607f1410a3SPeter Brune linesearch->ops->view = SNESLineSearchView_BT; 4610298fd71SBarry Smith linesearch->ops->setup = NULL; 46271b4ebd2SPeter Brune 4639566063dSJacob Faibussowitsch PetscCall(PetscNewLog(linesearch,&bt)); 464f5af7f23SKarl Rupp 46571b4ebd2SPeter Brune linesearch->data = (void*)bt; 46671b4ebd2SPeter Brune linesearch->max_its = 40; 467b000cd8dSPeter Brune linesearch->order = SNES_LINESEARCH_ORDER_CUBIC; 46871b4ebd2SPeter Brune bt->alpha = 1e-4; 46971b4ebd2SPeter Brune PetscFunctionReturn(0); 47071b4ebd2SPeter Brune } 471