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 @*/ 199371c9d4SSatish Balay PetscErrorCode SNESLineSearchBTSetAlpha(SNESLineSearch linesearch, PetscReal alpha) { 20d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 216e111a19SKarl Rupp 222f4102e2SPeter Brune PetscFunctionBegin; 232f4102e2SPeter Brune PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 1); 242f4102e2SPeter Brune bt->alpha = alpha; 252f4102e2SPeter Brune PetscFunctionReturn(0); 262f4102e2SPeter Brune } 272f4102e2SPeter Brune 282f4102e2SPeter Brune /*@ 292f4102e2SPeter Brune SNESLineSearchBTGetAlpha - Gets the descent parameter, alpha, in the BT linesearch variant. 302f4102e2SPeter Brune 312f4102e2SPeter Brune Input Parameters: 322f4102e2SPeter Brune . linesearch - linesearch context 338e557f58SPeter Brune 348e557f58SPeter Brune Output Parameters: 352f4102e2SPeter Brune . alpha - The descent parameter 362f4102e2SPeter Brune 372f4102e2SPeter Brune Level: intermediate 382f4102e2SPeter Brune 39db781477SPatrick Sanan .seealso: `SNESLineSearchGetLambda()`, `SNESLineSearchGetTolerances()` `SNESLINESEARCHBT` 402f4102e2SPeter Brune @*/ 419371c9d4SSatish Balay PetscErrorCode SNESLineSearchBTGetAlpha(SNESLineSearch linesearch, PetscReal *alpha) { 42d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 436e111a19SKarl Rupp 442f4102e2SPeter Brune PetscFunctionBegin; 452f4102e2SPeter Brune PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 1); 462f4102e2SPeter Brune *alpha = bt->alpha; 472f4102e2SPeter Brune PetscFunctionReturn(0); 482f4102e2SPeter Brune } 492f4102e2SPeter Brune 509371c9d4SSatish Balay static PetscErrorCode SNESLineSearchApply_BT(SNESLineSearch linesearch) { 5171b4ebd2SPeter Brune PetscBool changed_y, changed_w; 5271b4ebd2SPeter Brune Vec X, F, Y, W, G; 5371b4ebd2SPeter Brune SNES snes; 548a430afbSPeter Brune PetscReal fnorm, xnorm, ynorm, gnorm; 553b288129SPeter Brune PetscReal lambda, lambdatemp, lambdaprev, minlambda, maxstep, initslope, alpha, stol; 56dfcd3828SPeter Brune PetscReal t1, t2, a, b, d; 578a430afbSPeter Brune PetscReal f; 588a430afbSPeter Brune PetscReal g, gprev; 5971b4ebd2SPeter Brune PetscViewer monitor; 6071b4ebd2SPeter Brune PetscInt max_its, count; 61d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 6271b4ebd2SPeter Brune Mat jac; 63bd4a8a71SBarry Smith PetscErrorCode (*objective)(SNES, Vec, PetscReal *, void *); 6471b4ebd2SPeter Brune 6571b4ebd2SPeter Brune PetscFunctionBegin; 669566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G)); 679566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm)); 689566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetLambda(linesearch, &lambda)); 699566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetSNES(linesearch, &snes)); 709566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor)); 719566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetTolerances(linesearch, &minlambda, &maxstep, NULL, NULL, NULL, &max_its)); 729566063dSJacob Faibussowitsch PetscCall(SNESGetTolerances(snes, NULL, NULL, &stol, NULL, NULL)); 739566063dSJacob Faibussowitsch PetscCall(SNESGetObjective(snes, &objective, NULL)); 7471b4ebd2SPeter Brune alpha = bt->alpha; 7571b4ebd2SPeter Brune 769566063dSJacob Faibussowitsch PetscCall(SNESGetJacobian(snes, &jac, NULL, NULL, NULL)); 7708401ef6SPierre Jolivet PetscCheck(jac || objective, PetscObjectComm((PetscObject)linesearch), PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix"); 78c69d1a72SBarry Smith 799566063dSJacob Faibussowitsch PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y)); 809566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 8171b4ebd2SPeter Brune 829566063dSJacob Faibussowitsch PetscCall(VecNormBegin(Y, NORM_2, &ynorm)); 839566063dSJacob Faibussowitsch PetscCall(VecNormBegin(X, NORM_2, &xnorm)); 849566063dSJacob Faibussowitsch PetscCall(VecNormEnd(Y, NORM_2, &ynorm)); 859566063dSJacob Faibussowitsch PetscCall(VecNormEnd(X, NORM_2, &xnorm)); 863b288129SPeter Brune 8771b4ebd2SPeter Brune if (ynorm == 0.0) { 8871b4ebd2SPeter Brune if (monitor) { 899566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 909566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Initial direction and size is 0\n")); 919566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 9271b4ebd2SPeter Brune } 939566063dSJacob Faibussowitsch PetscCall(VecCopy(X, W)); 949566063dSJacob Faibussowitsch PetscCall(VecCopy(F, G)); 959566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm)); 969566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 9771b4ebd2SPeter Brune PetscFunctionReturn(0); 9871b4ebd2SPeter Brune } 9971b4ebd2SPeter Brune if (ynorm > maxstep) { /* Step too big, so scale back */ 10071b4ebd2SPeter Brune if (monitor) { 1019566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 1029566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Scaling step by %14.12e old ynorm %14.12e\n", (double)(maxstep / ynorm), (double)ynorm)); 1039566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 10471b4ebd2SPeter Brune } 1059566063dSJacob Faibussowitsch PetscCall(VecScale(Y, maxstep / (ynorm))); 10671b4ebd2SPeter Brune ynorm = maxstep; 10771b4ebd2SPeter Brune } 1088a430afbSPeter Brune 1098a430afbSPeter Brune /* if the SNES has an objective set, use that instead of the function value */ 110bd4a8a71SBarry Smith if (objective) { 1119566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes, X, &f)); 1128a430afbSPeter Brune } else { 1138a430afbSPeter Brune f = fnorm * fnorm; 1148a430afbSPeter Brune } 1158a430afbSPeter Brune 1168a430afbSPeter Brune /* compute the initial slope */ 117bd4a8a71SBarry Smith if (objective) { 1188a430afbSPeter Brune /* slope comes from the function (assumed to be the gradient of the objective */ 1199566063dSJacob Faibussowitsch PetscCall(VecDotRealPart(Y, F, &initslope)); 1208a430afbSPeter Brune } else { 1218a430afbSPeter Brune /* slope comes from the normal equations */ 1229566063dSJacob Faibussowitsch PetscCall(MatMult(jac, Y, W)); 1239566063dSJacob Faibussowitsch PetscCall(VecDotRealPart(F, W, &initslope)); 12471b4ebd2SPeter Brune if (initslope > 0.0) initslope = -initslope; 12571b4ebd2SPeter Brune if (initslope == 0.0) initslope = -1.0; 1268a430afbSPeter Brune } 12771b4ebd2SPeter Brune 128df019d78SBarry Smith while (PETSC_TRUE) { 1299566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -lambda, Y, X)); 1301baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 131e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 1329566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while checking full step length!\n")); 13371b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1349566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 13571b4ebd2SPeter Brune PetscFunctionReturn(0); 13671b4ebd2SPeter Brune } 1378a430afbSPeter Brune 138bd4a8a71SBarry Smith if (objective) { 1399566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes, W, &g)); 1408a430afbSPeter Brune } else { 1419566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 1429bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 1439bd66eb0SPeter Brune gnorm = fnorm; 1449566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 1459bd66eb0SPeter Brune } else { 1469566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 1478a430afbSPeter Brune } 14879e7e81bSJed Brown g = PetscSqr(gnorm); 1499bd66eb0SPeter Brune } 1509566063dSJacob Faibussowitsch PetscCall(SNESLineSearchMonitor(linesearch)); 1519bd66eb0SPeter Brune 152df019d78SBarry Smith if (!PetscIsInfOrNanReal(g)) break; 153df019d78SBarry Smith if (monitor) { 1549566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 1559566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n", (double)lambda)); 1569566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 157df019d78SBarry Smith } 1589371c9d4SSatish Balay if (lambda <= minlambda) { SNESCheckFunctionNorm(snes, g); } 159df019d78SBarry Smith lambda = .5 * lambda; 160df019d78SBarry Smith } 161df019d78SBarry Smith 162*48a46eb9SPierre Jolivet if (!objective) PetscCall(PetscInfo(snes, "Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm)); 1638a430afbSPeter Brune if (.5 * g <= .5 * f + lambda * alpha * initslope) { /* Sufficient reduction or step tolerance convergence */ 16471b4ebd2SPeter Brune if (monitor) { 1659566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 166bd4a8a71SBarry Smith if (!objective) { 1679566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm)); 1688a430afbSPeter Brune } else { 1699566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Using full step: obj %14.12e obj %14.12e\n", (double)f, (double)g)); 1708a430afbSPeter Brune } 1719566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 17271b4ebd2SPeter Brune } 17371b4ebd2SPeter Brune } else { 174c21ba15cSPeter Brune /* Since the full step didn't work and the step is tiny, quit */ 175c21ba15cSPeter Brune if (stol * xnorm > ynorm) { 1769566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm)); 1779566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 178c61ba1e2SPeter Brune if (monitor) { 1799566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 18063a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n", (double)ynorm, (double)(stol * xnorm))); 1819566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 182c61ba1e2SPeter Brune } 183c21ba15cSPeter Brune PetscFunctionReturn(0); 184c21ba15cSPeter Brune } 18571b4ebd2SPeter Brune /* Fit points with quadratic */ 1868a430afbSPeter Brune lambdatemp = -initslope / (g - f - 2.0 * lambda * initslope); 18771b4ebd2SPeter Brune lambdaprev = lambda; 1888a430afbSPeter Brune gprev = g; 18971b4ebd2SPeter Brune if (lambdatemp > .5 * lambda) lambdatemp = .5 * lambda; 19071b4ebd2SPeter Brune if (lambdatemp <= .1 * lambda) lambda = .1 * lambda; 19171b4ebd2SPeter Brune else lambda = lambdatemp; 19271b4ebd2SPeter Brune 1939566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -lambda, Y, X)); 1941baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 195e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 19663a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while attempting quadratic backtracking! %" PetscInt_FMT " \n", snes->nfuncs)); 19771b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1989566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 19971b4ebd2SPeter Brune PetscFunctionReturn(0); 20071b4ebd2SPeter Brune } 201bd4a8a71SBarry Smith if (objective) { 2029566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes, W, &g)); 2038a430afbSPeter Brune } else { 2049566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 2059bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 2069bd66eb0SPeter Brune gnorm = fnorm; 2079566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 2089bd66eb0SPeter Brune } else { 2099566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 2109bd66eb0SPeter Brune } 211f1ab69e3SPeter Brune g = PetscSqr(gnorm); 2128a430afbSPeter Brune } 213c98378a5SBarry Smith if (PetscIsInfOrNanReal(g)) { 2149566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 2159566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n")); 216c98378a5SBarry Smith PetscFunctionReturn(0); 217c98378a5SBarry Smith } 21871b4ebd2SPeter Brune if (monitor) { 2199566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 220bd4a8a71SBarry Smith if (!objective) { 2219566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: gnorm after quadratic fit %14.12e\n", (double)gnorm)); 2228a430afbSPeter Brune } else { 2239566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: obj after quadratic fit %14.12e\n", (double)g)); 2248a430afbSPeter Brune } 2259566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 22671b4ebd2SPeter Brune } 2278a430afbSPeter Brune if (.5 * g < .5 * f + lambda * alpha * initslope) { /* sufficient reduction */ 22871b4ebd2SPeter Brune if (monitor) { 2299566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 2309566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Quadratically determined step, lambda=%18.16e\n", (double)lambda)); 2319566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 23271b4ebd2SPeter Brune } 23371b4ebd2SPeter Brune } else { 23471b4ebd2SPeter Brune /* Fit points with cubic */ 23571b4ebd2SPeter Brune for (count = 0; count < max_its; count++) { 23671b4ebd2SPeter Brune if (lambda <= minlambda) { 23771b4ebd2SPeter Brune if (monitor) { 2389566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 23963a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: unable to find good step length! After %" PetscInt_FMT " tries \n", count)); 240bd4a8a71SBarry Smith if (!objective) { 2419371c9d4SSatish Balay 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", (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope)); 2428a430afbSPeter Brune } else { 2439371c9d4SSatish Balay 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", (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope)); 2448a430afbSPeter Brune } 2459566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 24671b4ebd2SPeter Brune } 2479566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 24871b4ebd2SPeter Brune PetscFunctionReturn(0); 24971b4ebd2SPeter Brune } 250b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 2518a430afbSPeter Brune t1 = .5 * (g - f) - lambda * initslope; 2528a430afbSPeter Brune t2 = .5 * (gprev - f) - lambdaprev * initslope; 25371b4ebd2SPeter Brune a = (t1 / (lambda * lambda) - t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev); 25471b4ebd2SPeter Brune b = (-lambdaprev * t1 / (lambda * lambda) + lambda * t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev); 25571b4ebd2SPeter Brune d = b * b - 3 * a * initslope; 25671b4ebd2SPeter Brune if (d < 0.0) d = 0.0; 257f5af7f23SKarl Rupp if (a == 0.0) lambdatemp = -initslope / (2.0 * b); 258f5af7f23SKarl Rupp else lambdatemp = (-b + PetscSqrtReal(d)) / (3.0 * a); 259f5af7f23SKarl Rupp 260b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 2618a430afbSPeter Brune lambdatemp = -initslope / (g - f - 2.0 * initslope); 262ce94432eSBarry Smith } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt"); 26371b4ebd2SPeter Brune lambdaprev = lambda; 2648a430afbSPeter Brune gprev = g; 26571b4ebd2SPeter Brune if (lambdatemp > .5 * lambda) lambdatemp = .5 * lambda; 26671b4ebd2SPeter Brune if (lambdatemp <= .1 * lambda) lambda = .1 * lambda; 26771b4ebd2SPeter Brune else lambda = lambdatemp; 2689566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -lambda, Y, X)); 2691baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 270e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 27163a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while looking for good step length! %" PetscInt_FMT " \n", count)); 272*48a46eb9SPierre Jolivet if (!objective) 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)); 2739566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 27471b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 27571b4ebd2SPeter Brune PetscFunctionReturn(0); 27671b4ebd2SPeter Brune } 277bd4a8a71SBarry Smith if (objective) { 2789566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes, W, &g)); 2798a430afbSPeter Brune } else { 2809566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 2819bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 2829bd66eb0SPeter Brune gnorm = fnorm; 2839566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 2849bd66eb0SPeter Brune } else { 2859566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 2868a430afbSPeter Brune } 287f1ab69e3SPeter Brune g = PetscSqr(gnorm); 2889bd66eb0SPeter Brune } 2894a2f8832SBarry Smith if (PetscIsInfOrNanReal(g)) { 2909566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 2919566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n")); 292c98378a5SBarry Smith PetscFunctionReturn(0); 293c98378a5SBarry Smith } 2948a430afbSPeter Brune if (.5 * g < .5 * f + lambda * alpha * initslope) { /* is reduction enough? */ 29571b4ebd2SPeter Brune if (monitor) { 2969566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 297bd4a8a71SBarry Smith if (!objective) { 298b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 2999566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n", (double)gnorm, (double)lambda)); 30071b4ebd2SPeter Brune } else { 3019566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n", (double)gnorm, (double)lambda)); 30271b4ebd2SPeter Brune } 3039566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3048a430afbSPeter Brune } else { 3058a430afbSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3069566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n", (double)g, (double)lambda)); 3078a430afbSPeter Brune } else { 3089566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n", (double)g, (double)lambda)); 3098a430afbSPeter Brune } 3109566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3118a430afbSPeter Brune } 31271b4ebd2SPeter Brune } 31371b4ebd2SPeter Brune break; 314f5af7f23SKarl Rupp } else 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: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n", (double)gnorm, (double)lambda)); 31971b4ebd2SPeter Brune } else { 3209566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Quadratic step no good, shrinking lambda, current gnorm %12.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: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n", (double)g, (double)lambda)); 3268a430afbSPeter Brune } else { 3279566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Quadratic step no good, shrinking lambda, obj %12.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 } 33371b4ebd2SPeter Brune } 33471b4ebd2SPeter Brune } 33571b4ebd2SPeter Brune 33671b4ebd2SPeter Brune /* postcheck */ 3370c1ac483SGlenn Hammond /* update Y to lambda*Y so that W is consistent with X - lambda*Y */ 3389566063dSJacob Faibussowitsch PetscCall(VecScale(Y, lambda)); 3399566063dSJacob Faibussowitsch PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w)); 34071b4ebd2SPeter Brune if (changed_y) { 3419566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -1.0, Y, X)); 3421baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 34371b4ebd2SPeter Brune } 344bd4a8a71SBarry Smith if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */ 3459566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 3469bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 3479bd66eb0SPeter Brune gnorm = fnorm; 3489566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 3499bd66eb0SPeter Brune } else { 3509566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 3519bd66eb0SPeter Brune } 3529566063dSJacob Faibussowitsch PetscCall(VecNorm(Y, NORM_2, &ynorm)); 353c98378a5SBarry Smith if (PetscIsInfOrNanReal(gnorm)) { 3549566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 3559566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n")); 356c98378a5SBarry Smith PetscFunctionReturn(0); 357c98378a5SBarry Smith } 358c427506bSPeter Brune } 35971b4ebd2SPeter Brune 36071b4ebd2SPeter Brune /* copy the solution over */ 3619566063dSJacob Faibussowitsch PetscCall(VecCopy(W, X)); 3629566063dSJacob Faibussowitsch PetscCall(VecCopy(G, F)); 3639566063dSJacob Faibussowitsch PetscCall(VecNorm(X, NORM_2, &xnorm)); 3649566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetLambda(linesearch, lambda)); 3659566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm)); 36671b4ebd2SPeter Brune PetscFunctionReturn(0); 36771b4ebd2SPeter Brune } 36871b4ebd2SPeter Brune 3699371c9d4SSatish Balay PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer) { 3707f1410a3SPeter Brune PetscBool iascii; 371d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 372075cc632SBarry Smith 3737f1410a3SPeter Brune PetscFunctionBegin; 3749566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 3757f1410a3SPeter Brune if (iascii) { 376b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3779566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: cubic\n")); 378b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 3799566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: quadratic\n")); 3807f1410a3SPeter Brune } 3819566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " alpha=%e\n", (double)bt->alpha)); 3827f1410a3SPeter Brune } 3837f1410a3SPeter Brune PetscFunctionReturn(0); 3847f1410a3SPeter Brune } 3857f1410a3SPeter Brune 3869371c9d4SSatish Balay static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch) { 38771b4ebd2SPeter Brune PetscFunctionBegin; 3889566063dSJacob Faibussowitsch PetscCall(PetscFree(linesearch->data)); 38971b4ebd2SPeter Brune PetscFunctionReturn(0); 39071b4ebd2SPeter Brune } 39171b4ebd2SPeter Brune 3929371c9d4SSatish Balay static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch, PetscOptionItems *PetscOptionsObject) { 393eb23ba39SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 39471b4ebd2SPeter Brune 3956e111a19SKarl Rupp PetscFunctionBegin; 396d0609cedSBarry Smith PetscOptionsHeadBegin(PetscOptionsObject, "SNESLineSearch BT options"); 3979566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL)); 398d0609cedSBarry Smith PetscOptionsHeadEnd(); 39971b4ebd2SPeter Brune PetscFunctionReturn(0); 40071b4ebd2SPeter Brune } 40171b4ebd2SPeter Brune 402954494b2SJed Brown /*MC 403db609ea7SPeter Brune SNESLINESEARCHBT - Backtracking line search. 404954494b2SJed Brown 405db609ea7SPeter Brune This line search finds the minimum of a polynomial fitting of the L2 norm of the 406eb23ba39SBarry 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 407db609ea7SPeter Brune and the fit is reattempted at most max_it times or until lambda is below minlambda. 408954494b2SJed Brown 409cd7522eaSPeter Brune Options Database Keys: 4103eec2f6aSPierre Jolivet + -snes_linesearch_alpha <1e\-4> - slope descent parameter 411db609ea7SPeter Brune . -snes_linesearch_damping <1.0> - initial step length 412eb23ba39SBarry Smith . -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the 413e1bc860dSBarry Smith step is scaled back to be of this length at the beginning of the line search 414db609ea7SPeter Brune . -snes_linesearch_max_it <40> - maximum number of shrinking step 4153eec2f6aSPierre Jolivet . -snes_linesearch_minlambda <1e\-12> - minimum step length allowed 416cd7522eaSPeter Brune - -snes_linesearch_order <cubic,quadratic> - order of the approximation 417cd7522eaSPeter Brune 418954494b2SJed Brown Level: advanced 419954494b2SJed Brown 420db609ea7SPeter Brune Notes: 421db609ea7SPeter Brune This line search is taken from "Numerical Methods for Unconstrained 422db609ea7SPeter Brune Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325. 423e55accfeSBarryFSmith 424e55accfeSBarryFSmith This line search will always produce a step that is less than or equal to, in length, the full step size. 425db609ea7SPeter Brune 426db781477SPatrick Sanan .seealso: `SNESLineSearchCreate()`, `SNESLineSearchSetType()` 427954494b2SJed Brown M*/ 4289371c9d4SSatish Balay PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch) { 429f1c6b773SPeter Brune SNESLineSearch_BT *bt; 43071b4ebd2SPeter Brune 43171b4ebd2SPeter Brune PetscFunctionBegin; 432f1c6b773SPeter Brune linesearch->ops->apply = SNESLineSearchApply_BT; 433f1c6b773SPeter Brune linesearch->ops->destroy = SNESLineSearchDestroy_BT; 434f1c6b773SPeter Brune linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT; 4350298fd71SBarry Smith linesearch->ops->reset = NULL; 4367f1410a3SPeter Brune linesearch->ops->view = SNESLineSearchView_BT; 4370298fd71SBarry Smith linesearch->ops->setup = NULL; 43871b4ebd2SPeter Brune 4399566063dSJacob Faibussowitsch PetscCall(PetscNewLog(linesearch, &bt)); 440f5af7f23SKarl Rupp 44171b4ebd2SPeter Brune linesearch->data = (void *)bt; 44271b4ebd2SPeter Brune linesearch->max_its = 40; 443b000cd8dSPeter Brune linesearch->order = SNES_LINESEARCH_ORDER_CUBIC; 44471b4ebd2SPeter Brune bt->alpha = 1e-4; 44571b4ebd2SPeter Brune PetscFunctionReturn(0); 44671b4ebd2SPeter Brune } 447