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 } 1096b095a85SStefano Zampini 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 */ 3266b095a85SStefano Zampini PetscCall(SNESLineSearchSetLambda(linesearch, lambda)); 3279566063dSJacob Faibussowitsch PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w)); 32871b4ebd2SPeter Brune if (changed_y) { 3296b095a85SStefano Zampini if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X)); 3301baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 33171b4ebd2SPeter Brune } 332bd4a8a71SBarry Smith if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */ 3339566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 3349bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 3359bd66eb0SPeter Brune gnorm = fnorm; 3369566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 3379bd66eb0SPeter Brune } else { 3389566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 3399bd66eb0SPeter Brune } 3409566063dSJacob Faibussowitsch PetscCall(VecNorm(Y, NORM_2, &ynorm)); 341c98378a5SBarry Smith if (PetscIsInfOrNanReal(gnorm)) { 3429566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 3439566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n")); 3443ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 345c98378a5SBarry Smith } 346c427506bSPeter Brune } 34771b4ebd2SPeter Brune 34871b4ebd2SPeter Brune /* copy the solution over */ 3499566063dSJacob Faibussowitsch PetscCall(VecCopy(W, X)); 3509566063dSJacob Faibussowitsch PetscCall(VecCopy(G, F)); 3519566063dSJacob Faibussowitsch PetscCall(VecNorm(X, NORM_2, &xnorm)); 3529566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm)); 3533ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 35471b4ebd2SPeter Brune } 35571b4ebd2SPeter Brune 35666976f2fSJacob Faibussowitsch static PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer) 357d71ae5a4SJacob Faibussowitsch { 3587f1410a3SPeter Brune PetscBool iascii; 359d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 360075cc632SBarry Smith 3617f1410a3SPeter Brune PetscFunctionBegin; 3629566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 3637f1410a3SPeter Brune if (iascii) { 364b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3659566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: cubic\n")); 366b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 3679566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: quadratic\n")); 3687f1410a3SPeter Brune } 3699566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " alpha=%e\n", (double)bt->alpha)); 3707f1410a3SPeter Brune } 3713ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 3727f1410a3SPeter Brune } 3737f1410a3SPeter Brune 374d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch) 375d71ae5a4SJacob Faibussowitsch { 37671b4ebd2SPeter Brune PetscFunctionBegin; 3779566063dSJacob Faibussowitsch PetscCall(PetscFree(linesearch->data)); 3783ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 37971b4ebd2SPeter Brune } 38071b4ebd2SPeter Brune 381d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch, PetscOptionItems *PetscOptionsObject) 382d71ae5a4SJacob Faibussowitsch { 383eb23ba39SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 38471b4ebd2SPeter Brune 3856e111a19SKarl Rupp PetscFunctionBegin; 386d0609cedSBarry Smith PetscOptionsHeadBegin(PetscOptionsObject, "SNESLineSearch BT options"); 3879566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL)); 388d0609cedSBarry Smith PetscOptionsHeadEnd(); 3893ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 39071b4ebd2SPeter Brune } 39171b4ebd2SPeter Brune 392954494b2SJed Brown /*MC 3931d27aa22SBarry Smith SNESLINESEARCHBT - Backtracking line search {cite}`dennis:83`. 394954494b2SJed Brown 395db609ea7SPeter Brune This line search finds the minimum of a polynomial fitting of the L2 norm of the 39632661dccSStefano Zampini function or the objective function if it is provided with `SNESSetObjective()`. 39732661dccSStefano Zampini If this fit does not satisfy the conditions for progress, the interval shrinks 398*89dfbbd5SBarry Smith and the fit is reattempted at most `max_it` times or until $\lambda$ is below `minlambda`. 399954494b2SJed Brown 400cd7522eaSPeter Brune Options Database Keys: 4013eec2f6aSPierre Jolivet + -snes_linesearch_alpha <1e\-4> - slope descent parameter 402*89dfbbd5SBarry Smith . -snes_linesearch_damping <1.0> - scaling of initial step length on entry to the line search 403eb23ba39SBarry Smith . -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the 404e1bc860dSBarry Smith step is scaled back to be of this length at the beginning of the line search 405*89dfbbd5SBarry Smith . -snes_linesearch_max_it <40> - maximum number of shrinking steps 4063eec2f6aSPierre Jolivet . -snes_linesearch_minlambda <1e\-12> - minimum step length allowed 40732661dccSStefano Zampini - -snes_linesearch_order <1,2,3> - order of the approximation. With order 1, it performs a simple backtracking without any curve fitting 408cd7522eaSPeter Brune 409954494b2SJed Brown Level: advanced 410954494b2SJed Brown 411f6dfbefdSBarry Smith Note: 412e55accfeSBarryFSmith This line search will always produce a step that is less than or equal to, in length, the full step size. 413db609ea7SPeter Brune 414420bcc1bSBarry Smith .seealso: [](ch_snes), `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()` 415954494b2SJed Brown M*/ 416d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch) 417d71ae5a4SJacob Faibussowitsch { 418f1c6b773SPeter Brune SNESLineSearch_BT *bt; 41971b4ebd2SPeter Brune 42071b4ebd2SPeter Brune PetscFunctionBegin; 421f1c6b773SPeter Brune linesearch->ops->apply = SNESLineSearchApply_BT; 422f1c6b773SPeter Brune linesearch->ops->destroy = SNESLineSearchDestroy_BT; 423f1c6b773SPeter Brune linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT; 4240298fd71SBarry Smith linesearch->ops->reset = NULL; 4257f1410a3SPeter Brune linesearch->ops->view = SNESLineSearchView_BT; 4260298fd71SBarry Smith linesearch->ops->setup = NULL; 42771b4ebd2SPeter Brune 4284dfa11a4SJacob Faibussowitsch PetscCall(PetscNew(&bt)); 429f5af7f23SKarl Rupp 43071b4ebd2SPeter Brune linesearch->data = (void *)bt; 43171b4ebd2SPeter Brune linesearch->max_its = 40; 432b000cd8dSPeter Brune linesearch->order = SNES_LINESEARCH_ORDER_CUBIC; 43371b4ebd2SPeter Brune bt->alpha = 1e-4; 4343ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 43571b4ebd2SPeter Brune } 436