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 /*@ 9420bcc1bSBarry Smith SNESLineSearchBTSetAlpha - Sets the descent parameter, `alpha`, in the `SNESLINESEARCHBT` `SNESLineSearch` 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 17420bcc1bSBarry Smith .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchSetLambda()`, `SNESLineSearchGetTolerances()`, `SNESLINESEARCHBT`, `SNESLineSearchBTGetAlpha()` 182f4102e2SPeter Brune @*/ 19d71ae5a4SJacob Faibussowitsch PetscErrorCode SNESLineSearchBTSetAlpha(SNESLineSearch linesearch, PetscReal alpha) 20d71ae5a4SJacob Faibussowitsch { 21d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 22d1b38757SStefano Zampini PetscBool isbt; 236e111a19SKarl Rupp 242f4102e2SPeter Brune PetscFunctionBegin; 252f4102e2SPeter Brune PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 1); 26d1b38757SStefano Zampini PetscCall(PetscObjectTypeCompare((PetscObject)linesearch, SNESLINESEARCHBT, &isbt)); 27d1b38757SStefano Zampini if (isbt) bt->alpha = alpha; 283ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 292f4102e2SPeter Brune } 302f4102e2SPeter Brune 312f4102e2SPeter Brune /*@ 32420bcc1bSBarry Smith SNESLineSearchBTGetAlpha - Gets the descent parameter, `alpha`, in the `SNESLINESEARCHBT` variant that was set with `SNESLineSearchBTSetAlpha()` 332f4102e2SPeter Brune 34f6dfbefdSBarry Smith Input Parameter: 352f4102e2SPeter Brune . linesearch - linesearch context 368e557f58SPeter Brune 37f6dfbefdSBarry Smith Output Parameter: 382f4102e2SPeter Brune . alpha - The descent parameter 392f4102e2SPeter Brune 402f4102e2SPeter Brune Level: intermediate 412f4102e2SPeter Brune 420241b274SPierre Jolivet .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchGetLambda()`, `SNESLineSearchGetTolerances()`, `SNESLINESEARCHBT`, `SNESLineSearchBTSetAlpha()` 432f4102e2SPeter Brune @*/ 44d71ae5a4SJacob Faibussowitsch PetscErrorCode SNESLineSearchBTGetAlpha(SNESLineSearch linesearch, PetscReal *alpha) 45d71ae5a4SJacob Faibussowitsch { 46d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 47d1b38757SStefano Zampini PetscBool isbt; 486e111a19SKarl Rupp 492f4102e2SPeter Brune PetscFunctionBegin; 502f4102e2SPeter Brune PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 1); 51d1b38757SStefano Zampini PetscCall(PetscObjectTypeCompare((PetscObject)linesearch, SNESLINESEARCHBT, &isbt)); 52d1b38757SStefano Zampini PetscCheck(isbt, PetscObjectComm((PetscObject)linesearch), PETSC_ERR_USER, "Not for type %s", ((PetscObject)linesearch)->type_name); 532f4102e2SPeter Brune *alpha = bt->alpha; 543ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 552f4102e2SPeter Brune } 562f4102e2SPeter Brune 57d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchApply_BT(SNESLineSearch linesearch) 58d71ae5a4SJacob Faibussowitsch { 596b72add0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 6071b4ebd2SPeter Brune PetscBool changed_y, changed_w; 6171b4ebd2SPeter Brune Vec X, F, Y, W, G; 6271b4ebd2SPeter Brune SNES snes; 638a430afbSPeter Brune PetscReal fnorm, xnorm, ynorm, gnorm; 64a99ef635SJonas Heinzmann PetscReal lambda, lambdatemp, lambdaprev, minlambda, initslope, alpha, stol; 65dfcd3828SPeter Brune PetscReal t1, t2, a, b, d; 668a430afbSPeter Brune PetscReal f; 678a430afbSPeter Brune PetscReal g, gprev; 6871b4ebd2SPeter Brune PetscViewer monitor; 69a99ef635SJonas Heinzmann PetscInt max_it, count; 7071b4ebd2SPeter Brune Mat jac; 716b72add0SBarry Smith SNESObjectiveFn *objective; 7232661dccSStefano Zampini const char *const ordStr[] = {"Linear", "Quadratic", "Cubic"}; 7371b4ebd2SPeter Brune 7471b4ebd2SPeter Brune PetscFunctionBegin; 759566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G)); 76a68bbae5SBarry Smith PetscCall(SNESLineSearchGetNorms(linesearch, NULL, &fnorm, NULL)); 779566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetLambda(linesearch, &lambda)); 789566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetSNES(linesearch, &snes)); 799566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor)); 80a99ef635SJonas Heinzmann PetscCall(SNESLineSearchGetTolerances(linesearch, &minlambda, NULL, NULL, NULL, NULL, &max_it)); 819566063dSJacob Faibussowitsch PetscCall(SNESGetTolerances(snes, NULL, NULL, &stol, NULL, NULL)); 829566063dSJacob Faibussowitsch PetscCall(SNESGetObjective(snes, &objective, NULL)); 8371b4ebd2SPeter Brune alpha = bt->alpha; 8471b4ebd2SPeter Brune 859566063dSJacob Faibussowitsch PetscCall(SNESGetJacobian(snes, &jac, NULL, NULL, NULL)); 8608401ef6SPierre Jolivet PetscCheck(jac || objective, PetscObjectComm((PetscObject)linesearch), PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix"); 87c69d1a72SBarry Smith 889566063dSJacob Faibussowitsch PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y)); 899566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 9071b4ebd2SPeter Brune 919566063dSJacob Faibussowitsch PetscCall(VecNormBegin(Y, NORM_2, &ynorm)); 929566063dSJacob Faibussowitsch PetscCall(VecNormBegin(X, NORM_2, &xnorm)); 939566063dSJacob Faibussowitsch PetscCall(VecNormEnd(Y, NORM_2, &ynorm)); 949566063dSJacob Faibussowitsch PetscCall(VecNormEnd(X, NORM_2, &xnorm)); 953b288129SPeter Brune 9671b4ebd2SPeter Brune if (ynorm == 0.0) { 9771b4ebd2SPeter Brune if (monitor) { 989566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 999566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Initial direction and size is 0\n")); 1009566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 10171b4ebd2SPeter Brune } 1029566063dSJacob Faibussowitsch PetscCall(VecCopy(X, W)); 1039566063dSJacob Faibussowitsch PetscCall(VecCopy(F, G)); 1049566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm)); 1059566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 1063ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 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)); 112*76c63389SBarry Smith SNESLineSearchCheckObjectiveDomainError(snes, f); 1138a430afbSPeter Brune } else { 114e3ce9940SStefano Zampini f = 0.5 * PetscSqr(fnorm); 1158a430afbSPeter Brune } 1168a430afbSPeter Brune 1178a430afbSPeter Brune /* compute the initial slope */ 118bd4a8a71SBarry Smith if (objective) { 11932661dccSStefano Zampini /* slope comes from the function (assumed to be the gradient of the objective) */ 1209566063dSJacob Faibussowitsch PetscCall(VecDotRealPart(Y, F, &initslope)); 1218a430afbSPeter Brune } else { 1228a430afbSPeter Brune /* slope comes from the normal equations */ 1239566063dSJacob Faibussowitsch PetscCall(MatMult(jac, Y, W)); 1249566063dSJacob Faibussowitsch PetscCall(VecDotRealPart(F, W, &initslope)); 12571b4ebd2SPeter Brune if (initslope > 0.0) initslope = -initslope; 12671b4ebd2SPeter Brune if (initslope == 0.0) initslope = -1.0; 1278a430afbSPeter Brune } 12871b4ebd2SPeter Brune 129*76c63389SBarry Smith /* repeatedly cut lambda until the norm or objective function is not infinity or NaN or lambda is too small */ 130df019d78SBarry Smith while (PETSC_TRUE) { 1319566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -lambda, Y, X)); 1321baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 133e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 1349566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while checking full step length!\n")); 13571b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1369566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 1373ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 13871b4ebd2SPeter Brune } 1398a430afbSPeter Brune 140bd4a8a71SBarry Smith if (objective) { 1419566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes, W, &g)); 1428a430afbSPeter Brune } else { 1439566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 1449bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 1459bd66eb0SPeter Brune gnorm = fnorm; 1469566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 1479bd66eb0SPeter Brune } else { 1489566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 1498a430afbSPeter Brune } 150e3ce9940SStefano Zampini g = 0.5 * PetscSqr(gnorm); 1519bd66eb0SPeter Brune } 1529566063dSJacob Faibussowitsch PetscCall(SNESLineSearchMonitor(linesearch)); 1539bd66eb0SPeter Brune 154df019d78SBarry Smith if (!PetscIsInfOrNanReal(g)) break; 155df019d78SBarry Smith if (monitor) { 1569566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 157*76c63389SBarry Smith PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: objective function at lambdas = %g is infinity or NaN, cutting lambda\n", (double)lambda)); 1589566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 159df019d78SBarry Smith } 160*76c63389SBarry Smith if (lambda <= minlambda) SNESCheckFunctionDomainError(snes, g); 16132661dccSStefano Zampini lambda *= .5; 162df019d78SBarry Smith } 163df019d78SBarry Smith 16448a46eb9SPierre Jolivet if (!objective) PetscCall(PetscInfo(snes, "Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm)); 165e3ce9940SStefano Zampini if (g <= f + lambda * alpha * initslope) { /* Sufficient reduction or step tolerance convergence */ 16671b4ebd2SPeter Brune if (monitor) { 1679566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 168bd4a8a71SBarry Smith if (!objective) { 1699566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm)); 1708a430afbSPeter Brune } else { 171e3ce9940SStefano Zampini PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Using full step: old obj %14.12e new obj %14.12e\n", (double)f, (double)g)); 1728a430afbSPeter Brune } 1739566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 17471b4ebd2SPeter Brune } 17571b4ebd2SPeter Brune } else { 176c21ba15cSPeter Brune if (stol * xnorm > ynorm) { 177a68bbae5SBarry Smith /* Since the full step didn't give sufficient decrease and the step is tiny, exit */ 1789566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm)); 1799566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 180c61ba1e2SPeter Brune if (monitor) { 1819566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 18263a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n", (double)ynorm, (double)(stol * xnorm))); 1839566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 184c61ba1e2SPeter Brune } 1853ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 186c21ba15cSPeter Brune } 18732661dccSStefano Zampini /* Here to avoid -Wmaybe-uninitiliazed warnings */ 18871b4ebd2SPeter Brune lambdaprev = lambda; 1898a430afbSPeter Brune gprev = g; 19032661dccSStefano Zampini if (linesearch->order != SNES_LINESEARCH_ORDER_LINEAR) { 19132661dccSStefano Zampini /* Fit points with quadratic */ 192e3ce9940SStefano Zampini lambdatemp = -initslope * PetscSqr(lambda) / (2.0 * (g - f - lambda * initslope)); 193e3ce9940SStefano Zampini lambda = PetscClipInterval(lambdatemp, .1 * lambda, .5 * lambda); 19471b4ebd2SPeter Brune 1959566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -lambda, Y, X)); 1961baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 197e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 19863a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while attempting quadratic backtracking! %" PetscInt_FMT " \n", snes->nfuncs)); 19971b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2009566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 2013ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 20271b4ebd2SPeter Brune } 203bd4a8a71SBarry Smith if (objective) { 2049566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes, W, &g)); 205*76c63389SBarry Smith SNESLineSearchCheckObjectiveDomainError(snes, g); 2068a430afbSPeter Brune } else { 2079566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 2089bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 2099bd66eb0SPeter Brune gnorm = fnorm; 2109566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 2119bd66eb0SPeter Brune } else { 2129566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 2139bd66eb0SPeter Brune } 214*76c63389SBarry Smith SNESLineSearchCheckFunctionDomainError(snes, linesearch, gnorm); 215e3ce9940SStefano Zampini g = 0.5 * PetscSqr(gnorm); 2168a430afbSPeter Brune } 217c98378a5SBarry Smith if (PetscIsInfOrNanReal(g)) { 2189566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 219*76c63389SBarry Smith PetscCall(PetscInfo(snes, "Aborted due to infinity or NaN in function evaluation\n")); 2203ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 221c98378a5SBarry Smith } 22271b4ebd2SPeter Brune if (monitor) { 2239566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 224bd4a8a71SBarry Smith if (!objective) { 2259566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: gnorm after quadratic fit %14.12e\n", (double)gnorm)); 2268a430afbSPeter Brune } else { 2279566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: obj after quadratic fit %14.12e\n", (double)g)); 2288a430afbSPeter Brune } 2299566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 23071b4ebd2SPeter Brune } 23132661dccSStefano Zampini } 232e3ce9940SStefano Zampini if (linesearch->order != SNES_LINESEARCH_ORDER_LINEAR && g <= f + lambda * alpha * initslope) { /* sufficient reduction */ 23371b4ebd2SPeter Brune if (monitor) { 2349566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 2359566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Quadratically determined step, lambda=%18.16e\n", (double)lambda)); 2369566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 23771b4ebd2SPeter Brune } 23871b4ebd2SPeter Brune } else { 239a99ef635SJonas Heinzmann for (count = 0; count < max_it; count++) { 24071b4ebd2SPeter Brune if (lambda <= minlambda) { 24171b4ebd2SPeter Brune if (monitor) { 2429566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 24363a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: unable to find good step length! After %" PetscInt_FMT " tries \n", count)); 244bd4a8a71SBarry Smith if (!objective) { 2459371c9d4SSatish 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)); 2468a430afbSPeter Brune } else { 247e3ce9940SStefano Zampini 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)); 2488a430afbSPeter Brune } 2499566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 25071b4ebd2SPeter Brune } 2519566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 2523ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 25371b4ebd2SPeter Brune } 254b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 255e3ce9940SStefano Zampini /* Fit points with cubic */ 256e3ce9940SStefano Zampini t1 = g - f - lambda * initslope; 257e3ce9940SStefano Zampini t2 = gprev - f - lambdaprev * initslope; 25871b4ebd2SPeter Brune a = (t1 / (lambda * lambda) - t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev); 25971b4ebd2SPeter Brune b = (-lambdaprev * t1 / (lambda * lambda) + lambda * t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev); 26071b4ebd2SPeter Brune d = b * b - 3 * a * initslope; 26171b4ebd2SPeter Brune if (d < 0.0) d = 0.0; 262f5af7f23SKarl Rupp if (a == 0.0) lambdatemp = -initslope / (2.0 * b); 263f5af7f23SKarl Rupp else lambdatemp = (-b + PetscSqrtReal(d)) / (3.0 * a); 264b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 265e3ce9940SStefano Zampini lambdatemp = -initslope * PetscSqr(lambda) / (2.0 * (g - f - lambda * initslope)); 26632661dccSStefano Zampini } else if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) { /* Just backtrack */ 26732661dccSStefano Zampini lambdatemp = .5 * lambda; 26832661dccSStefano Zampini } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "Line search order %" PetscInt_FMT " for type bt", linesearch->order); 26971b4ebd2SPeter Brune lambdaprev = lambda; 2708a430afbSPeter Brune gprev = g; 271e3ce9940SStefano Zampini 272e3ce9940SStefano Zampini lambda = PetscClipInterval(lambdatemp, .1 * lambda, .5 * lambda); 2739566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -lambda, Y, X)); 2741baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 275e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 27663a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while looking for good step length! %" PetscInt_FMT " \n", count)); 27748a46eb9SPierre 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)); 2789566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 27971b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2803ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 28171b4ebd2SPeter Brune } 282bd4a8a71SBarry Smith if (objective) { 2839566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes, W, &g)); 284*76c63389SBarry Smith SNESLineSearchCheckObjectiveDomainError(snes, g); 2858a430afbSPeter Brune } else { 2869566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 2879bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 2889bd66eb0SPeter Brune gnorm = fnorm; 2899566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 2909bd66eb0SPeter Brune } else { 2919566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 2928a430afbSPeter Brune } 293*76c63389SBarry Smith SNESLineSearchCheckFunctionDomainError(snes, linesearch, gnorm); 294e3ce9940SStefano Zampini g = 0.5 * PetscSqr(gnorm); 2959bd66eb0SPeter Brune } 296e3ce9940SStefano Zampini if (g <= f + lambda * alpha * initslope) { /* is reduction enough? */ 29771b4ebd2SPeter Brune if (monitor) { 2989566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 299bd4a8a71SBarry Smith if (!objective) { 30032661dccSStefano Zampini PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: %s step, current gnorm %14.12e lambda=%18.16e\n", ordStr[linesearch->order - 1], (double)gnorm, (double)lambda)); 3019566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3028a430afbSPeter Brune } else { 30332661dccSStefano Zampini PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: %s step, obj %14.12e lambda=%18.16e\n", ordStr[linesearch->order - 1], (double)g, (double)lambda)); 3049566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3058a430afbSPeter Brune } 30671b4ebd2SPeter Brune } 30771b4ebd2SPeter Brune break; 308f5af7f23SKarl Rupp } else if (monitor) { 3099566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 310bd4a8a71SBarry Smith if (!objective) { 31132661dccSStefano Zampini PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: %s step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n", ordStr[linesearch->order - 1], (double)gnorm, (double)lambda)); 3129566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3138a430afbSPeter Brune } else { 31432661dccSStefano Zampini PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: %s step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n", ordStr[linesearch->order - 1], (double)g, (double)lambda)); 3159566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3168a430afbSPeter Brune } 31771b4ebd2SPeter Brune } 31871b4ebd2SPeter Brune } 31971b4ebd2SPeter Brune } 32071b4ebd2SPeter Brune } 32171b4ebd2SPeter Brune 32271b4ebd2SPeter Brune /* postcheck */ 3236b095a85SStefano Zampini PetscCall(SNESLineSearchSetLambda(linesearch, lambda)); 3249566063dSJacob Faibussowitsch PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w)); 32571b4ebd2SPeter Brune if (changed_y) { 3266b095a85SStefano Zampini if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X)); 3271baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 32871b4ebd2SPeter Brune } 329bd4a8a71SBarry Smith if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */ 3309566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 3319bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 3329bd66eb0SPeter Brune gnorm = fnorm; 3339566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 3349bd66eb0SPeter Brune } else { 3359566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 3369bd66eb0SPeter Brune } 337*76c63389SBarry Smith SNESLineSearchCheckFunctionDomainError(snes, linesearch, gnorm); 3389566063dSJacob Faibussowitsch PetscCall(VecNorm(Y, NORM_2, &ynorm)); 339c427506bSPeter Brune } 34071b4ebd2SPeter Brune 34171b4ebd2SPeter Brune /* copy the solution over */ 3429566063dSJacob Faibussowitsch PetscCall(VecCopy(W, X)); 3439566063dSJacob Faibussowitsch PetscCall(VecCopy(G, F)); 3449566063dSJacob Faibussowitsch PetscCall(VecNorm(X, NORM_2, &xnorm)); 3459566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm)); 3463ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 34771b4ebd2SPeter Brune } 34871b4ebd2SPeter Brune 34966976f2fSJacob Faibussowitsch static PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer) 350d71ae5a4SJacob Faibussowitsch { 3519f196a02SMartin Diehl PetscBool isascii; 352d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 353075cc632SBarry Smith 3547f1410a3SPeter Brune PetscFunctionBegin; 3559f196a02SMartin Diehl PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 3569f196a02SMartin Diehl if (isascii) { 357b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3589566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: cubic\n")); 359b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 3609566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: quadratic\n")); 3617f1410a3SPeter Brune } 3629566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " alpha=%e\n", (double)bt->alpha)); 3637f1410a3SPeter Brune } 3643ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 3657f1410a3SPeter Brune } 3667f1410a3SPeter Brune 367d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch) 368d71ae5a4SJacob Faibussowitsch { 36971b4ebd2SPeter Brune PetscFunctionBegin; 3709566063dSJacob Faibussowitsch PetscCall(PetscFree(linesearch->data)); 3713ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 37271b4ebd2SPeter Brune } 37371b4ebd2SPeter Brune 374ce78bad3SBarry Smith static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch, PetscOptionItems PetscOptionsObject) 375d71ae5a4SJacob Faibussowitsch { 376eb23ba39SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 37771b4ebd2SPeter Brune 3786e111a19SKarl Rupp PetscFunctionBegin; 379d0609cedSBarry Smith PetscOptionsHeadBegin(PetscOptionsObject, "SNESLineSearch BT options"); 3809566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL)); 381d0609cedSBarry Smith PetscOptionsHeadEnd(); 3823ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 38371b4ebd2SPeter Brune } 38471b4ebd2SPeter Brune 385954494b2SJed Brown /*MC 3861d27aa22SBarry Smith SNESLINESEARCHBT - Backtracking line search {cite}`dennis:83`. 387954494b2SJed Brown 388a99ef635SJonas Heinzmann This line search finds the minimum of a polynomial fitting either $1/2 ||F(x_k + \lambda Y_k)||_2^2$, 389a99ef635SJonas Heinzmann or the objective function $G(x_k + \lambda Y_k)$ if it is provided with `SNESSetObjective()`. 390a99ef635SJonas Heinzmann If this fit does not satisfy the sufficient decrease conditions, the interval shrinks 39189dfbbd5SBarry Smith and the fit is reattempted at most `max_it` times or until $\lambda$ is below `minlambda`. 392954494b2SJed Brown 393cd7522eaSPeter Brune Options Database Keys: 3943eec2f6aSPierre Jolivet + -snes_linesearch_alpha <1e\-4> - slope descent parameter 395a99ef635SJonas Heinzmann . -snes_linesearch_damping <1.0> - initial `lambda` on entry to the line search 396a99ef635SJonas Heinzmann . -snes_linesearch_max_it <40> - maximum number of shrinking iterations in the line search 397a99ef635SJonas Heinzmann . -snes_linesearch_minlambda <1e\-12> - minimum `lambda` (scaling of solution update) allowed 398a99ef635SJonas Heinzmann - -snes_linesearch_order <3> - order of the polynomial fit, must be 1, 2, or 3. With order 1, it performs a simple backtracking without any curve fitting 399cd7522eaSPeter Brune 400954494b2SJed Brown Level: advanced 401954494b2SJed Brown 402f6dfbefdSBarry Smith Note: 403e55accfeSBarryFSmith This line search will always produce a step that is less than or equal to, in length, the full step size. 404db609ea7SPeter Brune 405420bcc1bSBarry Smith .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()` 406954494b2SJed Brown M*/ 407d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch) 408d71ae5a4SJacob Faibussowitsch { 409f1c6b773SPeter Brune SNESLineSearch_BT *bt; 41071b4ebd2SPeter Brune 41171b4ebd2SPeter Brune PetscFunctionBegin; 412f1c6b773SPeter Brune linesearch->ops->apply = SNESLineSearchApply_BT; 413f1c6b773SPeter Brune linesearch->ops->destroy = SNESLineSearchDestroy_BT; 414f1c6b773SPeter Brune linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT; 4150298fd71SBarry Smith linesearch->ops->reset = NULL; 4167f1410a3SPeter Brune linesearch->ops->view = SNESLineSearchView_BT; 4170298fd71SBarry Smith linesearch->ops->setup = NULL; 41871b4ebd2SPeter Brune 4194dfa11a4SJacob Faibussowitsch PetscCall(PetscNew(&bt)); 420f5af7f23SKarl Rupp 42171b4ebd2SPeter Brune linesearch->data = (void *)bt; 422a99ef635SJonas Heinzmann linesearch->max_it = 40; 423b000cd8dSPeter Brune linesearch->order = SNES_LINESEARCH_ORDER_CUBIC; 42471b4ebd2SPeter Brune bt->alpha = 1e-4; 4253ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 42671b4ebd2SPeter Brune } 427