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; 22*d1b38757SStefano Zampini PetscBool isbt; 236e111a19SKarl Rupp 242f4102e2SPeter Brune PetscFunctionBegin; 252f4102e2SPeter Brune PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 1); 26*d1b38757SStefano Zampini PetscCall(PetscObjectTypeCompare((PetscObject)linesearch, SNESLINESEARCHBT, &isbt)); 27*d1b38757SStefano 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; 47*d1b38757SStefano Zampini PetscBool isbt; 486e111a19SKarl Rupp 492f4102e2SPeter Brune PetscFunctionBegin; 502f4102e2SPeter Brune PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 1); 51*d1b38757SStefano Zampini PetscCall(PetscObjectTypeCompare((PetscObject)linesearch, SNESLINESEARCHBT, &isbt)); 52*d1b38757SStefano 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)); 1128a430afbSPeter Brune } else { 113e3ce9940SStefano Zampini f = 0.5 * PetscSqr(fnorm); 1148a430afbSPeter Brune } 1158a430afbSPeter Brune 1168a430afbSPeter Brune /* compute the initial slope */ 117bd4a8a71SBarry Smith if (objective) { 11832661dccSStefano Zampini /* 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)); 1353ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 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 } 148e3ce9940SStefano Zampini g = 0.5 * 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 } 158ad540459SPierre Jolivet if (lambda <= minlambda) SNESCheckFunctionNorm(snes, g); 15932661dccSStefano Zampini lambda *= .5; 160df019d78SBarry Smith } 161df019d78SBarry Smith 16248a46eb9SPierre Jolivet if (!objective) PetscCall(PetscInfo(snes, "Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm)); 163e3ce9940SStefano Zampini if (g <= 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 { 169e3ce9940SStefano Zampini PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Using full step: old obj %14.12e new 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 if (stol * xnorm > ynorm) { 175a68bbae5SBarry Smith /* Since the full step didn't give sufficient decrease and the step is tiny, exit */ 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 } 1833ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 184c21ba15cSPeter Brune } 18532661dccSStefano Zampini /* Here to avoid -Wmaybe-uninitiliazed warnings */ 18671b4ebd2SPeter Brune lambdaprev = lambda; 1878a430afbSPeter Brune gprev = g; 18832661dccSStefano Zampini if (linesearch->order != SNES_LINESEARCH_ORDER_LINEAR) { 18932661dccSStefano Zampini /* Fit points with quadratic */ 190e3ce9940SStefano Zampini lambdatemp = -initslope * PetscSqr(lambda) / (2.0 * (g - f - lambda * initslope)); 191e3ce9940SStefano Zampini lambda = PetscClipInterval(lambdatemp, .1 * lambda, .5 * lambda); 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)); 1993ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 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 } 211e3ce9940SStefano Zampini g = 0.5 * 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")); 2163ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 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 } 22732661dccSStefano Zampini } 228e3ce9940SStefano Zampini if (linesearch->order != SNES_LINESEARCH_ORDER_LINEAR && g <= f + lambda * alpha * initslope) { /* sufficient reduction */ 22971b4ebd2SPeter Brune if (monitor) { 2309566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 2319566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Quadratically determined step, lambda=%18.16e\n", (double)lambda)); 2329566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 23371b4ebd2SPeter Brune } 23471b4ebd2SPeter Brune } else { 235a99ef635SJonas Heinzmann for (count = 0; count < max_it; 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 { 243e3ce9940SStefano 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)); 2448a430afbSPeter Brune } 2459566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 24671b4ebd2SPeter Brune } 2479566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 2483ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 24971b4ebd2SPeter Brune } 250b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 251e3ce9940SStefano Zampini /* Fit points with cubic */ 252e3ce9940SStefano Zampini t1 = g - f - lambda * initslope; 253e3ce9940SStefano Zampini t2 = gprev - f - lambdaprev * initslope; 25471b4ebd2SPeter Brune a = (t1 / (lambda * lambda) - t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev); 25571b4ebd2SPeter Brune b = (-lambdaprev * t1 / (lambda * lambda) + lambda * t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev); 25671b4ebd2SPeter Brune d = b * b - 3 * a * initslope; 25771b4ebd2SPeter Brune if (d < 0.0) d = 0.0; 258f5af7f23SKarl Rupp if (a == 0.0) lambdatemp = -initslope / (2.0 * b); 259f5af7f23SKarl Rupp else lambdatemp = (-b + PetscSqrtReal(d)) / (3.0 * a); 260b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 261e3ce9940SStefano Zampini lambdatemp = -initslope * PetscSqr(lambda) / (2.0 * (g - f - lambda * initslope)); 26232661dccSStefano Zampini } else if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) { /* Just backtrack */ 26332661dccSStefano Zampini lambdatemp = .5 * lambda; 26432661dccSStefano Zampini } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "Line search order %" PetscInt_FMT " for type bt", linesearch->order); 26571b4ebd2SPeter Brune lambdaprev = lambda; 2668a430afbSPeter Brune gprev = g; 267e3ce9940SStefano Zampini 268e3ce9940SStefano Zampini lambda = PetscClipInterval(lambdatemp, .1 * lambda, .5 * lambda); 2699566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -lambda, Y, X)); 2701baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 271e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 27263a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while looking for good step length! %" PetscInt_FMT " \n", count)); 27348a46eb9SPierre 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)); 2749566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 27571b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2763ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 27771b4ebd2SPeter Brune } 278bd4a8a71SBarry Smith if (objective) { 2799566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes, W, &g)); 2808a430afbSPeter Brune } else { 2819566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 2829bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 2839bd66eb0SPeter Brune gnorm = fnorm; 2849566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 2859bd66eb0SPeter Brune } else { 2869566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 2878a430afbSPeter Brune } 288e3ce9940SStefano Zampini g = 0.5 * PetscSqr(gnorm); 2899bd66eb0SPeter Brune } 2904a2f8832SBarry Smith if (PetscIsInfOrNanReal(g)) { 2919566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 2929566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n")); 2933ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 294c98378a5SBarry Smith } 295e3ce9940SStefano Zampini if (g <= f + lambda * alpha * initslope) { /* is reduction enough? */ 29671b4ebd2SPeter Brune if (monitor) { 2979566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 298bd4a8a71SBarry Smith if (!objective) { 29932661dccSStefano Zampini PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: %s step, current gnorm %14.12e lambda=%18.16e\n", ordStr[linesearch->order - 1], (double)gnorm, (double)lambda)); 3009566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3018a430afbSPeter Brune } else { 30232661dccSStefano Zampini PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: %s step, obj %14.12e lambda=%18.16e\n", ordStr[linesearch->order - 1], (double)g, (double)lambda)); 3039566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3048a430afbSPeter Brune } 30571b4ebd2SPeter Brune } 30671b4ebd2SPeter Brune break; 307f5af7f23SKarl Rupp } else if (monitor) { 3089566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 309bd4a8a71SBarry Smith if (!objective) { 31032661dccSStefano 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)); 3119566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3128a430afbSPeter Brune } else { 31332661dccSStefano 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)); 3149566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3158a430afbSPeter Brune } 31671b4ebd2SPeter Brune } 31771b4ebd2SPeter Brune } 31871b4ebd2SPeter Brune } 31971b4ebd2SPeter Brune } 32071b4ebd2SPeter Brune 32171b4ebd2SPeter Brune /* postcheck */ 3226b095a85SStefano Zampini PetscCall(SNESLineSearchSetLambda(linesearch, lambda)); 3239566063dSJacob Faibussowitsch PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w)); 32471b4ebd2SPeter Brune if (changed_y) { 3256b095a85SStefano Zampini if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X)); 3261baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 32771b4ebd2SPeter Brune } 328bd4a8a71SBarry Smith if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */ 3299566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 3309bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 3319bd66eb0SPeter Brune gnorm = fnorm; 3329566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 3339bd66eb0SPeter Brune } else { 3349566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 3359bd66eb0SPeter Brune } 3369566063dSJacob Faibussowitsch PetscCall(VecNorm(Y, NORM_2, &ynorm)); 337c98378a5SBarry Smith if (PetscIsInfOrNanReal(gnorm)) { 3389566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 3399566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n")); 3403ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 341c98378a5SBarry Smith } 342c427506bSPeter Brune } 34371b4ebd2SPeter Brune 34471b4ebd2SPeter Brune /* copy the solution over */ 3459566063dSJacob Faibussowitsch PetscCall(VecCopy(W, X)); 3469566063dSJacob Faibussowitsch PetscCall(VecCopy(G, F)); 3479566063dSJacob Faibussowitsch PetscCall(VecNorm(X, NORM_2, &xnorm)); 3489566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm)); 3493ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 35071b4ebd2SPeter Brune } 35171b4ebd2SPeter Brune 35266976f2fSJacob Faibussowitsch static PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer) 353d71ae5a4SJacob Faibussowitsch { 3549f196a02SMartin Diehl PetscBool isascii; 355d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 356075cc632SBarry Smith 3577f1410a3SPeter Brune PetscFunctionBegin; 3589f196a02SMartin Diehl PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 3599f196a02SMartin Diehl if (isascii) { 360b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3619566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: cubic\n")); 362b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 3639566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: quadratic\n")); 3647f1410a3SPeter Brune } 3659566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " alpha=%e\n", (double)bt->alpha)); 3667f1410a3SPeter Brune } 3673ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 3687f1410a3SPeter Brune } 3697f1410a3SPeter Brune 370d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch) 371d71ae5a4SJacob Faibussowitsch { 37271b4ebd2SPeter Brune PetscFunctionBegin; 3739566063dSJacob Faibussowitsch PetscCall(PetscFree(linesearch->data)); 3743ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 37571b4ebd2SPeter Brune } 37671b4ebd2SPeter Brune 377ce78bad3SBarry Smith static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch, PetscOptionItems PetscOptionsObject) 378d71ae5a4SJacob Faibussowitsch { 379eb23ba39SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 38071b4ebd2SPeter Brune 3816e111a19SKarl Rupp PetscFunctionBegin; 382d0609cedSBarry Smith PetscOptionsHeadBegin(PetscOptionsObject, "SNESLineSearch BT options"); 3839566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL)); 384d0609cedSBarry Smith PetscOptionsHeadEnd(); 3853ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 38671b4ebd2SPeter Brune } 38771b4ebd2SPeter Brune 388954494b2SJed Brown /*MC 3891d27aa22SBarry Smith SNESLINESEARCHBT - Backtracking line search {cite}`dennis:83`. 390954494b2SJed Brown 391a99ef635SJonas Heinzmann This line search finds the minimum of a polynomial fitting either $1/2 ||F(x_k + \lambda Y_k)||_2^2$, 392a99ef635SJonas Heinzmann or the objective function $G(x_k + \lambda Y_k)$ if it is provided with `SNESSetObjective()`. 393a99ef635SJonas Heinzmann If this fit does not satisfy the sufficient decrease conditions, the interval shrinks 39489dfbbd5SBarry Smith and the fit is reattempted at most `max_it` times or until $\lambda$ is below `minlambda`. 395954494b2SJed Brown 396cd7522eaSPeter Brune Options Database Keys: 3973eec2f6aSPierre Jolivet + -snes_linesearch_alpha <1e\-4> - slope descent parameter 398a99ef635SJonas Heinzmann . -snes_linesearch_damping <1.0> - initial `lambda` on entry to the line search 399a99ef635SJonas Heinzmann . -snes_linesearch_max_it <40> - maximum number of shrinking iterations in the line search 400a99ef635SJonas Heinzmann . -snes_linesearch_minlambda <1e\-12> - minimum `lambda` (scaling of solution update) allowed 401a99ef635SJonas 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 402cd7522eaSPeter Brune 403954494b2SJed Brown Level: advanced 404954494b2SJed Brown 405f6dfbefdSBarry Smith Note: 406e55accfeSBarryFSmith This line search will always produce a step that is less than or equal to, in length, the full step size. 407db609ea7SPeter Brune 408420bcc1bSBarry Smith .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()` 409954494b2SJed Brown M*/ 410d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch) 411d71ae5a4SJacob Faibussowitsch { 412f1c6b773SPeter Brune SNESLineSearch_BT *bt; 41371b4ebd2SPeter Brune 41471b4ebd2SPeter Brune PetscFunctionBegin; 415f1c6b773SPeter Brune linesearch->ops->apply = SNESLineSearchApply_BT; 416f1c6b773SPeter Brune linesearch->ops->destroy = SNESLineSearchDestroy_BT; 417f1c6b773SPeter Brune linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT; 4180298fd71SBarry Smith linesearch->ops->reset = NULL; 4197f1410a3SPeter Brune linesearch->ops->view = SNESLineSearchView_BT; 4200298fd71SBarry Smith linesearch->ops->setup = NULL; 42171b4ebd2SPeter Brune 4224dfa11a4SJacob Faibussowitsch PetscCall(PetscNew(&bt)); 423f5af7f23SKarl Rupp 42471b4ebd2SPeter Brune linesearch->data = (void *)bt; 425a99ef635SJonas Heinzmann linesearch->max_it = 40; 426b000cd8dSPeter Brune linesearch->order = SNES_LINESEARCH_ORDER_CUBIC; 42771b4ebd2SPeter Brune bt->alpha = 1e-4; 4283ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 42971b4ebd2SPeter Brune } 430