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; 226e111a19SKarl Rupp 232f4102e2SPeter Brune PetscFunctionBegin; 242f4102e2SPeter Brune PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 1); 252f4102e2SPeter Brune bt->alpha = alpha; 263ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 272f4102e2SPeter Brune } 282f4102e2SPeter Brune 292f4102e2SPeter Brune /*@ 30420bcc1bSBarry Smith SNESLineSearchBTGetAlpha - Gets the descent parameter, `alpha`, in the `SNESLINESEARCHBT` variant that was set with `SNESLineSearchBTSetAlpha()` 312f4102e2SPeter Brune 32f6dfbefdSBarry Smith Input Parameter: 332f4102e2SPeter Brune . linesearch - linesearch context 348e557f58SPeter Brune 35f6dfbefdSBarry Smith Output Parameter: 362f4102e2SPeter Brune . alpha - The descent parameter 372f4102e2SPeter Brune 382f4102e2SPeter Brune Level: intermediate 392f4102e2SPeter Brune 40420bcc1bSBarry Smith .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchGetLambda()`, `SNESLineSearchGetTolerances()` `SNESLINESEARCHBT`, `SNESLineSearchBTSetAlpha()` 412f4102e2SPeter Brune @*/ 42d71ae5a4SJacob Faibussowitsch PetscErrorCode SNESLineSearchBTGetAlpha(SNESLineSearch linesearch, PetscReal *alpha) 43d71ae5a4SJacob Faibussowitsch { 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; 493ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 502f4102e2SPeter Brune } 512f4102e2SPeter Brune 52d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchApply_BT(SNESLineSearch linesearch) 53d71ae5a4SJacob Faibussowitsch { 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; 6632661dccSStefano Zampini const char *const ordStr[] = {"Linear", "Quadratic", "Cubic"}; 67bd4a8a71SBarry Smith PetscErrorCode (*objective)(SNES, Vec, PetscReal *, void *); 6871b4ebd2SPeter Brune 6971b4ebd2SPeter Brune PetscFunctionBegin; 709566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G)); 71a68bbae5SBarry Smith PetscCall(SNESLineSearchGetNorms(linesearch, NULL, &fnorm, NULL)); 729566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetLambda(linesearch, &lambda)); 739566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetSNES(linesearch, &snes)); 749566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor)); 759566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetTolerances(linesearch, &minlambda, &maxstep, NULL, NULL, NULL, &max_its)); 769566063dSJacob Faibussowitsch PetscCall(SNESGetTolerances(snes, NULL, NULL, &stol, NULL, NULL)); 779566063dSJacob Faibussowitsch PetscCall(SNESGetObjective(snes, &objective, NULL)); 7871b4ebd2SPeter Brune alpha = bt->alpha; 7971b4ebd2SPeter Brune 809566063dSJacob Faibussowitsch PetscCall(SNESGetJacobian(snes, &jac, NULL, NULL, NULL)); 8108401ef6SPierre Jolivet PetscCheck(jac || objective, PetscObjectComm((PetscObject)linesearch), PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix"); 82c69d1a72SBarry Smith 839566063dSJacob Faibussowitsch PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y)); 849566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 8571b4ebd2SPeter Brune 869566063dSJacob Faibussowitsch PetscCall(VecNormBegin(Y, NORM_2, &ynorm)); 879566063dSJacob Faibussowitsch PetscCall(VecNormBegin(X, NORM_2, &xnorm)); 889566063dSJacob Faibussowitsch PetscCall(VecNormEnd(Y, NORM_2, &ynorm)); 899566063dSJacob Faibussowitsch PetscCall(VecNormEnd(X, NORM_2, &xnorm)); 903b288129SPeter Brune 9171b4ebd2SPeter Brune if (ynorm == 0.0) { 9271b4ebd2SPeter Brune if (monitor) { 939566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 949566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Initial direction and size is 0\n")); 959566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 9671b4ebd2SPeter Brune } 979566063dSJacob Faibussowitsch PetscCall(VecCopy(X, W)); 989566063dSJacob Faibussowitsch PetscCall(VecCopy(F, G)); 999566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm)); 1009566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 1013ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 10271b4ebd2SPeter Brune } 10371b4ebd2SPeter Brune if (ynorm > maxstep) { /* Step too big, so scale back */ 10471b4ebd2SPeter Brune if (monitor) { 1059566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 1069566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Scaling step by %14.12e old ynorm %14.12e\n", (double)(maxstep / ynorm), (double)ynorm)); 1079566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 10871b4ebd2SPeter Brune } 1099566063dSJacob Faibussowitsch PetscCall(VecScale(Y, maxstep / (ynorm))); 11071b4ebd2SPeter Brune ynorm = maxstep; 11171b4ebd2SPeter Brune } 1128a430afbSPeter Brune 1138a430afbSPeter Brune /* if the SNES has an objective set, use that instead of the function value */ 114bd4a8a71SBarry Smith if (objective) { 1159566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes, X, &f)); 1168a430afbSPeter Brune } else { 117e3ce9940SStefano Zampini f = 0.5 * PetscSqr(fnorm); 1188a430afbSPeter Brune } 1198a430afbSPeter Brune 1208a430afbSPeter Brune /* compute the initial slope */ 121bd4a8a71SBarry Smith if (objective) { 12232661dccSStefano Zampini /* slope comes from the function (assumed to be the gradient of the objective) */ 1239566063dSJacob Faibussowitsch PetscCall(VecDotRealPart(Y, F, &initslope)); 1248a430afbSPeter Brune } else { 1258a430afbSPeter Brune /* slope comes from the normal equations */ 1269566063dSJacob Faibussowitsch PetscCall(MatMult(jac, Y, W)); 1279566063dSJacob Faibussowitsch PetscCall(VecDotRealPart(F, W, &initslope)); 12871b4ebd2SPeter Brune if (initslope > 0.0) initslope = -initslope; 12971b4ebd2SPeter Brune if (initslope == 0.0) initslope = -1.0; 1308a430afbSPeter Brune } 13171b4ebd2SPeter Brune 132df019d78SBarry Smith while (PETSC_TRUE) { 1339566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -lambda, Y, X)); 1341baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 135e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 1369566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while checking full step length!\n")); 13771b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1389566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 1393ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 14071b4ebd2SPeter Brune } 1418a430afbSPeter Brune 142bd4a8a71SBarry Smith if (objective) { 1439566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes, W, &g)); 1448a430afbSPeter Brune } else { 1459566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 1469bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 1479bd66eb0SPeter Brune gnorm = fnorm; 1489566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 1499bd66eb0SPeter Brune } else { 1509566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 1518a430afbSPeter Brune } 152e3ce9940SStefano Zampini g = 0.5 * PetscSqr(gnorm); 1539bd66eb0SPeter Brune } 1549566063dSJacob Faibussowitsch PetscCall(SNESLineSearchMonitor(linesearch)); 1559bd66eb0SPeter Brune 156df019d78SBarry Smith if (!PetscIsInfOrNanReal(g)) break; 157df019d78SBarry Smith if (monitor) { 1589566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 1599566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n", (double)lambda)); 1609566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 161df019d78SBarry Smith } 162ad540459SPierre Jolivet if (lambda <= minlambda) SNESCheckFunctionNorm(snes, g); 16332661dccSStefano Zampini lambda *= .5; 164df019d78SBarry Smith } 165df019d78SBarry Smith 16648a46eb9SPierre Jolivet if (!objective) PetscCall(PetscInfo(snes, "Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm)); 167e3ce9940SStefano Zampini if (g <= f + lambda * alpha * initslope) { /* Sufficient reduction or step tolerance convergence */ 16871b4ebd2SPeter Brune if (monitor) { 1699566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 170bd4a8a71SBarry Smith if (!objective) { 1719566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm)); 1728a430afbSPeter Brune } else { 173e3ce9940SStefano Zampini PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Using full step: old obj %14.12e new obj %14.12e\n", (double)f, (double)g)); 1748a430afbSPeter Brune } 1759566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 17671b4ebd2SPeter Brune } 17771b4ebd2SPeter Brune } else { 178c21ba15cSPeter Brune if (stol * xnorm > ynorm) { 179a68bbae5SBarry Smith /* Since the full step didn't give sufficient decrease and the step is tiny, exit */ 1809566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm)); 1819566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 182c61ba1e2SPeter Brune if (monitor) { 1839566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 18463a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n", (double)ynorm, (double)(stol * xnorm))); 1859566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 186c61ba1e2SPeter Brune } 1873ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 188c21ba15cSPeter Brune } 18932661dccSStefano Zampini /* Here to avoid -Wmaybe-uninitiliazed warnings */ 19071b4ebd2SPeter Brune lambdaprev = lambda; 1918a430afbSPeter Brune gprev = g; 19232661dccSStefano Zampini if (linesearch->order != SNES_LINESEARCH_ORDER_LINEAR) { 19332661dccSStefano Zampini /* Fit points with quadratic */ 194e3ce9940SStefano Zampini lambdatemp = -initslope * PetscSqr(lambda) / (2.0 * (g - f - lambda * initslope)); 195e3ce9940SStefano Zampini lambda = PetscClipInterval(lambdatemp, .1 * lambda, .5 * lambda); 19671b4ebd2SPeter Brune 1979566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -lambda, Y, X)); 1981baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 199e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 20063a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while attempting quadratic backtracking! %" PetscInt_FMT " \n", snes->nfuncs)); 20171b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2029566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 2033ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 20471b4ebd2SPeter Brune } 205bd4a8a71SBarry Smith if (objective) { 2069566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes, W, &g)); 2078a430afbSPeter Brune } else { 2089566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 2099bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 2109bd66eb0SPeter Brune gnorm = fnorm; 2119566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 2129bd66eb0SPeter Brune } else { 2139566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 2149bd66eb0SPeter Brune } 215e3ce9940SStefano Zampini g = 0.5 * PetscSqr(gnorm); 2168a430afbSPeter Brune } 217c98378a5SBarry Smith if (PetscIsInfOrNanReal(g)) { 2189566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 2199566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf 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 { 23971b4ebd2SPeter Brune for (count = 0; count < max_its; 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)); 2848a430afbSPeter Brune } else { 2859566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 2869bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 2879bd66eb0SPeter Brune gnorm = fnorm; 2889566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 2899bd66eb0SPeter Brune } else { 2909566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 2918a430afbSPeter Brune } 292e3ce9940SStefano Zampini g = 0.5 * PetscSqr(gnorm); 2939bd66eb0SPeter Brune } 2944a2f8832SBarry Smith if (PetscIsInfOrNanReal(g)) { 2959566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 2969566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n")); 2973ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 298c98378a5SBarry Smith } 299e3ce9940SStefano Zampini if (g <= f + lambda * alpha * initslope) { /* is reduction enough? */ 30071b4ebd2SPeter Brune if (monitor) { 3019566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 302bd4a8a71SBarry Smith if (!objective) { 30332661dccSStefano Zampini PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: %s step, current gnorm %14.12e lambda=%18.16e\n", ordStr[linesearch->order - 1], (double)gnorm, (double)lambda)); 3049566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3058a430afbSPeter Brune } else { 30632661dccSStefano Zampini PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: %s step, obj %14.12e lambda=%18.16e\n", ordStr[linesearch->order - 1], (double)g, (double)lambda)); 3079566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3088a430afbSPeter Brune } 30971b4ebd2SPeter Brune } 31071b4ebd2SPeter Brune break; 311f5af7f23SKarl Rupp } else if (monitor) { 3129566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 313bd4a8a71SBarry Smith if (!objective) { 31432661dccSStefano 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)); 3159566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3168a430afbSPeter Brune } else { 31732661dccSStefano 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)); 3189566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3198a430afbSPeter Brune } 32071b4ebd2SPeter Brune } 32171b4ebd2SPeter Brune } 32271b4ebd2SPeter Brune } 32371b4ebd2SPeter Brune } 32471b4ebd2SPeter Brune 32571b4ebd2SPeter Brune /* postcheck */ 3260c1ac483SGlenn Hammond /* update Y to lambda*Y so that W is consistent with X - lambda*Y */ 3279566063dSJacob Faibussowitsch PetscCall(VecScale(Y, lambda)); 3289566063dSJacob Faibussowitsch PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w)); 32971b4ebd2SPeter Brune if (changed_y) { 3309566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -1.0, Y, X)); 3311baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 33271b4ebd2SPeter Brune } 333bd4a8a71SBarry Smith if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */ 3349566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 3359bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 3369bd66eb0SPeter Brune gnorm = fnorm; 3379566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 3389bd66eb0SPeter Brune } else { 3399566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 3409bd66eb0SPeter Brune } 3419566063dSJacob Faibussowitsch PetscCall(VecNorm(Y, NORM_2, &ynorm)); 342c98378a5SBarry Smith if (PetscIsInfOrNanReal(gnorm)) { 3439566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 3449566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n")); 3453ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 346c98378a5SBarry Smith } 347c427506bSPeter Brune } 34871b4ebd2SPeter Brune 34971b4ebd2SPeter Brune /* copy the solution over */ 3509566063dSJacob Faibussowitsch PetscCall(VecCopy(W, X)); 3519566063dSJacob Faibussowitsch PetscCall(VecCopy(G, F)); 3529566063dSJacob Faibussowitsch PetscCall(VecNorm(X, NORM_2, &xnorm)); 3539566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetLambda(linesearch, lambda)); 3549566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm)); 3553ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 35671b4ebd2SPeter Brune } 35771b4ebd2SPeter Brune 35866976f2fSJacob Faibussowitsch static PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer) 359d71ae5a4SJacob Faibussowitsch { 3607f1410a3SPeter Brune PetscBool iascii; 361d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 362075cc632SBarry Smith 3637f1410a3SPeter Brune PetscFunctionBegin; 3649566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 3657f1410a3SPeter Brune if (iascii) { 366b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3679566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: cubic\n")); 368b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 3699566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: quadratic\n")); 3707f1410a3SPeter Brune } 3719566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " alpha=%e\n", (double)bt->alpha)); 3727f1410a3SPeter Brune } 3733ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 3747f1410a3SPeter Brune } 3757f1410a3SPeter Brune 376d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch) 377d71ae5a4SJacob Faibussowitsch { 37871b4ebd2SPeter Brune PetscFunctionBegin; 3799566063dSJacob Faibussowitsch PetscCall(PetscFree(linesearch->data)); 3803ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 38171b4ebd2SPeter Brune } 38271b4ebd2SPeter Brune 383d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch, PetscOptionItems *PetscOptionsObject) 384d71ae5a4SJacob Faibussowitsch { 385eb23ba39SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 38671b4ebd2SPeter Brune 3876e111a19SKarl Rupp PetscFunctionBegin; 388d0609cedSBarry Smith PetscOptionsHeadBegin(PetscOptionsObject, "SNESLineSearch BT options"); 3899566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL)); 390d0609cedSBarry Smith PetscOptionsHeadEnd(); 3913ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 39271b4ebd2SPeter Brune } 39371b4ebd2SPeter Brune 394954494b2SJed Brown /*MC 395*1d27aa22SBarry Smith SNESLINESEARCHBT - Backtracking line search {cite}`dennis:83`. 396954494b2SJed Brown 397db609ea7SPeter Brune This line search finds the minimum of a polynomial fitting of the L2 norm of the 39832661dccSStefano Zampini function or the objective function if it is provided with `SNESSetObjective()`. 39932661dccSStefano Zampini If this fit does not satisfy the conditions for progress, the interval shrinks 400db609ea7SPeter Brune and the fit is reattempted at most max_it times or until lambda is below minlambda. 401954494b2SJed Brown 402cd7522eaSPeter Brune Options Database Keys: 4033eec2f6aSPierre Jolivet + -snes_linesearch_alpha <1e\-4> - slope descent parameter 404db609ea7SPeter Brune . -snes_linesearch_damping <1.0> - initial step length 405eb23ba39SBarry Smith . -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the 406e1bc860dSBarry Smith step is scaled back to be of this length at the beginning of the line search 407db609ea7SPeter Brune . -snes_linesearch_max_it <40> - maximum number of shrinking step 4083eec2f6aSPierre Jolivet . -snes_linesearch_minlambda <1e\-12> - minimum step length allowed 40932661dccSStefano Zampini - -snes_linesearch_order <1,2,3> - order of the approximation. With order 1, it performs a simple backtracking without any curve fitting 410cd7522eaSPeter Brune 411954494b2SJed Brown Level: advanced 412954494b2SJed Brown 413f6dfbefdSBarry Smith Note: 414e55accfeSBarryFSmith This line search will always produce a step that is less than or equal to, in length, the full step size. 415db609ea7SPeter Brune 416420bcc1bSBarry Smith .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()` 417954494b2SJed Brown M*/ 418d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch) 419d71ae5a4SJacob Faibussowitsch { 420f1c6b773SPeter Brune SNESLineSearch_BT *bt; 42171b4ebd2SPeter Brune 42271b4ebd2SPeter Brune PetscFunctionBegin; 423f1c6b773SPeter Brune linesearch->ops->apply = SNESLineSearchApply_BT; 424f1c6b773SPeter Brune linesearch->ops->destroy = SNESLineSearchDestroy_BT; 425f1c6b773SPeter Brune linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT; 4260298fd71SBarry Smith linesearch->ops->reset = NULL; 4277f1410a3SPeter Brune linesearch->ops->view = SNESLineSearchView_BT; 4280298fd71SBarry Smith linesearch->ops->setup = NULL; 42971b4ebd2SPeter Brune 4304dfa11a4SJacob Faibussowitsch PetscCall(PetscNew(&bt)); 431f5af7f23SKarl Rupp 43271b4ebd2SPeter Brune linesearch->data = (void *)bt; 43371b4ebd2SPeter Brune linesearch->max_its = 40; 434b000cd8dSPeter Brune linesearch->order = SNES_LINESEARCH_ORDER_CUBIC; 43571b4ebd2SPeter Brune bt->alpha = 1e-4; 4363ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 43771b4ebd2SPeter Brune } 438