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 /*@ 9f6dfbefdSBarry Smith SNESLineSearchBTSetAlpha - Sets the descent parameter, alpha, in the `SNESLINESEARCHBT` 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 17f6dfbefdSBarry Smith .seealso: `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; 26*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 272f4102e2SPeter Brune } 282f4102e2SPeter Brune 292f4102e2SPeter Brune /*@ 30f6dfbefdSBarry Smith SNESLineSearchBTGetAlpha - Gets the descent parameter, alpha, in the `SNESLINESEARCHBT` variant. 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 40f6dfbefdSBarry Smith .seealso: `SNESLineSearchGetLambda()`, `SNESLineSearchGetTolerances()` `SNESLINESEARCHBT`, `SNESLineSearchBTGetAlpha()` 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; 49*3ba16761SJacob 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; 66bd4a8a71SBarry Smith PetscErrorCode (*objective)(SNES, Vec, PetscReal *, void *); 6771b4ebd2SPeter Brune 6871b4ebd2SPeter Brune PetscFunctionBegin; 699566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G)); 709566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm)); 719566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetLambda(linesearch, &lambda)); 729566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetSNES(linesearch, &snes)); 739566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor)); 749566063dSJacob Faibussowitsch PetscCall(SNESLineSearchGetTolerances(linesearch, &minlambda, &maxstep, NULL, NULL, NULL, &max_its)); 759566063dSJacob Faibussowitsch PetscCall(SNESGetTolerances(snes, NULL, NULL, &stol, NULL, NULL)); 769566063dSJacob Faibussowitsch PetscCall(SNESGetObjective(snes, &objective, NULL)); 7771b4ebd2SPeter Brune alpha = bt->alpha; 7871b4ebd2SPeter Brune 799566063dSJacob Faibussowitsch PetscCall(SNESGetJacobian(snes, &jac, NULL, NULL, NULL)); 8008401ef6SPierre Jolivet PetscCheck(jac || objective, PetscObjectComm((PetscObject)linesearch), PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix"); 81c69d1a72SBarry Smith 829566063dSJacob Faibussowitsch PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y)); 839566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 8471b4ebd2SPeter Brune 859566063dSJacob Faibussowitsch PetscCall(VecNormBegin(Y, NORM_2, &ynorm)); 869566063dSJacob Faibussowitsch PetscCall(VecNormBegin(X, NORM_2, &xnorm)); 879566063dSJacob Faibussowitsch PetscCall(VecNormEnd(Y, NORM_2, &ynorm)); 889566063dSJacob Faibussowitsch PetscCall(VecNormEnd(X, NORM_2, &xnorm)); 893b288129SPeter Brune 9071b4ebd2SPeter Brune if (ynorm == 0.0) { 9171b4ebd2SPeter Brune if (monitor) { 929566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 939566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Initial direction and size is 0\n")); 949566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 9571b4ebd2SPeter Brune } 969566063dSJacob Faibussowitsch PetscCall(VecCopy(X, W)); 979566063dSJacob Faibussowitsch PetscCall(VecCopy(F, G)); 989566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm)); 999566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 100*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 10171b4ebd2SPeter Brune } 10271b4ebd2SPeter Brune if (ynorm > maxstep) { /* Step too big, so scale back */ 10371b4ebd2SPeter Brune if (monitor) { 1049566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 1059566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Scaling step by %14.12e old ynorm %14.12e\n", (double)(maxstep / ynorm), (double)ynorm)); 1069566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 10771b4ebd2SPeter Brune } 1089566063dSJacob Faibussowitsch PetscCall(VecScale(Y, maxstep / (ynorm))); 10971b4ebd2SPeter Brune ynorm = maxstep; 11071b4ebd2SPeter Brune } 1118a430afbSPeter Brune 1128a430afbSPeter Brune /* if the SNES has an objective set, use that instead of the function value */ 113bd4a8a71SBarry Smith if (objective) { 1149566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes, X, &f)); 1158a430afbSPeter Brune } else { 1168a430afbSPeter Brune f = fnorm * fnorm; 1178a430afbSPeter Brune } 1188a430afbSPeter Brune 1198a430afbSPeter Brune /* compute the initial slope */ 120bd4a8a71SBarry Smith if (objective) { 1218a430afbSPeter Brune /* slope comes from the function (assumed to be the gradient of the objective */ 1229566063dSJacob Faibussowitsch PetscCall(VecDotRealPart(Y, F, &initslope)); 1238a430afbSPeter Brune } else { 1248a430afbSPeter Brune /* slope comes from the normal equations */ 1259566063dSJacob Faibussowitsch PetscCall(MatMult(jac, Y, W)); 1269566063dSJacob Faibussowitsch PetscCall(VecDotRealPart(F, W, &initslope)); 12771b4ebd2SPeter Brune if (initslope > 0.0) initslope = -initslope; 12871b4ebd2SPeter Brune if (initslope == 0.0) initslope = -1.0; 1298a430afbSPeter Brune } 13071b4ebd2SPeter Brune 131df019d78SBarry Smith while (PETSC_TRUE) { 1329566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -lambda, Y, X)); 1331baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 134e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 1359566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while checking full step length!\n")); 13671b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1379566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 138*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 13971b4ebd2SPeter Brune } 1408a430afbSPeter Brune 141bd4a8a71SBarry Smith if (objective) { 1429566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes, W, &g)); 1438a430afbSPeter Brune } else { 1449566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 1459bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 1469bd66eb0SPeter Brune gnorm = fnorm; 1479566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 1489bd66eb0SPeter Brune } else { 1499566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 1508a430afbSPeter Brune } 15179e7e81bSJed Brown g = PetscSqr(gnorm); 1529bd66eb0SPeter Brune } 1539566063dSJacob Faibussowitsch PetscCall(SNESLineSearchMonitor(linesearch)); 1549bd66eb0SPeter Brune 155df019d78SBarry Smith if (!PetscIsInfOrNanReal(g)) break; 156df019d78SBarry Smith if (monitor) { 1579566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 1589566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n", (double)lambda)); 1599566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 160df019d78SBarry Smith } 161ad540459SPierre Jolivet if (lambda <= minlambda) SNESCheckFunctionNorm(snes, g); 162df019d78SBarry Smith lambda = .5 * lambda; 163df019d78SBarry Smith } 164df019d78SBarry Smith 16548a46eb9SPierre Jolivet if (!objective) PetscCall(PetscInfo(snes, "Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm)); 1668a430afbSPeter Brune if (.5 * g <= .5 * f + lambda * alpha * initslope) { /* Sufficient reduction or step tolerance convergence */ 16771b4ebd2SPeter Brune if (monitor) { 1689566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 169bd4a8a71SBarry Smith if (!objective) { 1709566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm)); 1718a430afbSPeter Brune } else { 1729566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Using full step: obj %14.12e obj %14.12e\n", (double)f, (double)g)); 1738a430afbSPeter Brune } 1749566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 17571b4ebd2SPeter Brune } 17671b4ebd2SPeter Brune } else { 177c21ba15cSPeter Brune /* Since the full step didn't work and the step is tiny, quit */ 178c21ba15cSPeter Brune if (stol * xnorm > ynorm) { 1799566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm)); 1809566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 181c61ba1e2SPeter Brune if (monitor) { 1829566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 18363a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n", (double)ynorm, (double)(stol * xnorm))); 1849566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 185c61ba1e2SPeter Brune } 186*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 187c21ba15cSPeter Brune } 18871b4ebd2SPeter Brune /* Fit points with quadratic */ 1898a430afbSPeter Brune lambdatemp = -initslope / (g - f - 2.0 * lambda * initslope); 19071b4ebd2SPeter Brune lambdaprev = lambda; 1918a430afbSPeter Brune gprev = g; 19271b4ebd2SPeter Brune if (lambdatemp > .5 * lambda) lambdatemp = .5 * lambda; 19371b4ebd2SPeter Brune if (lambdatemp <= .1 * lambda) lambda = .1 * lambda; 19471b4ebd2SPeter Brune else lambda = lambdatemp; 19571b4ebd2SPeter Brune 1969566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -lambda, Y, X)); 1971baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 198e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 19963a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while attempting quadratic backtracking! %" PetscInt_FMT " \n", snes->nfuncs)); 20071b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2019566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 202*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 20371b4ebd2SPeter Brune } 204bd4a8a71SBarry Smith if (objective) { 2059566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes, W, &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 } 214f1ab69e3SPeter Brune g = PetscSqr(gnorm); 2158a430afbSPeter Brune } 216c98378a5SBarry Smith if (PetscIsInfOrNanReal(g)) { 2179566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 2189566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n")); 219*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 220c98378a5SBarry Smith } 22171b4ebd2SPeter Brune if (monitor) { 2229566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 223bd4a8a71SBarry Smith if (!objective) { 2249566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: gnorm after quadratic fit %14.12e\n", (double)gnorm)); 2258a430afbSPeter Brune } else { 2269566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: obj after quadratic fit %14.12e\n", (double)g)); 2278a430afbSPeter Brune } 2289566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 22971b4ebd2SPeter Brune } 2308a430afbSPeter Brune if (.5 * g < .5 * f + lambda * alpha * initslope) { /* sufficient reduction */ 23171b4ebd2SPeter Brune if (monitor) { 2329566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 2339566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Quadratically determined step, lambda=%18.16e\n", (double)lambda)); 2349566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 23571b4ebd2SPeter Brune } 23671b4ebd2SPeter Brune } else { 23771b4ebd2SPeter Brune /* Fit points with cubic */ 23871b4ebd2SPeter Brune for (count = 0; count < max_its; count++) { 23971b4ebd2SPeter Brune if (lambda <= minlambda) { 24071b4ebd2SPeter Brune if (monitor) { 2419566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 24263a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: unable to find good step length! After %" PetscInt_FMT " tries \n", count)); 243bd4a8a71SBarry Smith if (!objective) { 2449371c9d4SSatish 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)); 2458a430afbSPeter Brune } else { 2469371c9d4SSatish Balay PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: obj(0)=%18.16e, obj=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", (double)f, (double)g, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope)); 2478a430afbSPeter Brune } 2489566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 24971b4ebd2SPeter Brune } 2509566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 251*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 25271b4ebd2SPeter Brune } 253b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 2548a430afbSPeter Brune t1 = .5 * (g - f) - lambda * initslope; 2558a430afbSPeter Brune t2 = .5 * (gprev - f) - lambdaprev * initslope; 25671b4ebd2SPeter Brune a = (t1 / (lambda * lambda) - t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev); 25771b4ebd2SPeter Brune b = (-lambdaprev * t1 / (lambda * lambda) + lambda * t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev); 25871b4ebd2SPeter Brune d = b * b - 3 * a * initslope; 25971b4ebd2SPeter Brune if (d < 0.0) d = 0.0; 260f5af7f23SKarl Rupp if (a == 0.0) lambdatemp = -initslope / (2.0 * b); 261f5af7f23SKarl Rupp else lambdatemp = (-b + PetscSqrtReal(d)) / (3.0 * a); 262f5af7f23SKarl Rupp 263b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 2648a430afbSPeter Brune lambdatemp = -initslope / (g - f - 2.0 * initslope); 265ce94432eSBarry Smith } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "unsupported line search order for type bt"); 26671b4ebd2SPeter Brune lambdaprev = lambda; 2678a430afbSPeter Brune gprev = g; 26871b4ebd2SPeter Brune if (lambdatemp > .5 * lambda) lambdatemp = .5 * lambda; 26971b4ebd2SPeter Brune if (lambdatemp <= .1 * lambda) lambda = .1 * lambda; 27071b4ebd2SPeter Brune else lambda = lambdatemp; 2719566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -lambda, Y, X)); 2721baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 273e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 27463a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while looking for good step length! %" PetscInt_FMT " \n", count)); 27548a46eb9SPierre 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)); 2769566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 27771b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 278*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 27971b4ebd2SPeter Brune } 280bd4a8a71SBarry Smith if (objective) { 2819566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes, W, &g)); 2828a430afbSPeter Brune } else { 2839566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 2849bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 2859bd66eb0SPeter Brune gnorm = fnorm; 2869566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 2879bd66eb0SPeter Brune } else { 2889566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 2898a430afbSPeter Brune } 290f1ab69e3SPeter Brune g = PetscSqr(gnorm); 2919bd66eb0SPeter Brune } 2924a2f8832SBarry Smith if (PetscIsInfOrNanReal(g)) { 2939566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 2949566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n")); 295*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 296c98378a5SBarry Smith } 2978a430afbSPeter Brune if (.5 * g < .5 * f + lambda * alpha * initslope) { /* is reduction enough? */ 29871b4ebd2SPeter Brune if (monitor) { 2999566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 300bd4a8a71SBarry Smith if (!objective) { 301b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3029566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n", (double)gnorm, (double)lambda)); 30371b4ebd2SPeter Brune } else { 3049566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n", (double)gnorm, (double)lambda)); 30571b4ebd2SPeter Brune } 3069566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3078a430afbSPeter Brune } else { 3088a430afbSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3099566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Cubically determined step, obj %14.12e lambda=%18.16e\n", (double)g, (double)lambda)); 3108a430afbSPeter Brune } else { 3119566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Quadratically determined step, obj %14.12e lambda=%18.16e\n", (double)g, (double)lambda)); 3128a430afbSPeter Brune } 3139566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3148a430afbSPeter Brune } 31571b4ebd2SPeter Brune } 31671b4ebd2SPeter Brune break; 317f5af7f23SKarl Rupp } else if (monitor) { 3189566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 319bd4a8a71SBarry Smith if (!objective) { 320b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3219566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n", (double)gnorm, (double)lambda)); 32271b4ebd2SPeter Brune } else { 3239566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n", (double)gnorm, (double)lambda)); 32471b4ebd2SPeter Brune } 3259566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3268a430afbSPeter Brune } else { 3278a430afbSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3289566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Cubic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n", (double)g, (double)lambda)); 3298a430afbSPeter Brune } else { 3309566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Quadratic step no good, shrinking lambda, obj %12.12e lambda=%18.16e\n", (double)g, (double)lambda)); 3318a430afbSPeter Brune } 3329566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3338a430afbSPeter Brune } 33471b4ebd2SPeter Brune } 33571b4ebd2SPeter Brune } 33671b4ebd2SPeter Brune } 33771b4ebd2SPeter Brune } 33871b4ebd2SPeter Brune 33971b4ebd2SPeter Brune /* postcheck */ 3400c1ac483SGlenn Hammond /* update Y to lambda*Y so that W is consistent with X - lambda*Y */ 3419566063dSJacob Faibussowitsch PetscCall(VecScale(Y, lambda)); 3429566063dSJacob Faibussowitsch PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w)); 34371b4ebd2SPeter Brune if (changed_y) { 3449566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -1.0, Y, X)); 3451baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 34671b4ebd2SPeter Brune } 347bd4a8a71SBarry Smith if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */ 3489566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 3499bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 3509bd66eb0SPeter Brune gnorm = fnorm; 3519566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 3529bd66eb0SPeter Brune } else { 3539566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 3549bd66eb0SPeter Brune } 3559566063dSJacob Faibussowitsch PetscCall(VecNorm(Y, NORM_2, &ynorm)); 356c98378a5SBarry Smith if (PetscIsInfOrNanReal(gnorm)) { 3579566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 3589566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n")); 359*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 360c98378a5SBarry Smith } 361c427506bSPeter Brune } 36271b4ebd2SPeter Brune 36371b4ebd2SPeter Brune /* copy the solution over */ 3649566063dSJacob Faibussowitsch PetscCall(VecCopy(W, X)); 3659566063dSJacob Faibussowitsch PetscCall(VecCopy(G, F)); 3669566063dSJacob Faibussowitsch PetscCall(VecNorm(X, NORM_2, &xnorm)); 3679566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetLambda(linesearch, lambda)); 3689566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm)); 369*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 37071b4ebd2SPeter Brune } 37171b4ebd2SPeter Brune 372d71ae5a4SJacob Faibussowitsch PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer) 373d71ae5a4SJacob Faibussowitsch { 3747f1410a3SPeter Brune PetscBool iascii; 375d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 376075cc632SBarry Smith 3777f1410a3SPeter Brune PetscFunctionBegin; 3789566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 3797f1410a3SPeter Brune if (iascii) { 380b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3819566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: cubic\n")); 382b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 3839566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: quadratic\n")); 3847f1410a3SPeter Brune } 3859566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " alpha=%e\n", (double)bt->alpha)); 3867f1410a3SPeter Brune } 387*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 3887f1410a3SPeter Brune } 3897f1410a3SPeter Brune 390d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch) 391d71ae5a4SJacob Faibussowitsch { 39271b4ebd2SPeter Brune PetscFunctionBegin; 3939566063dSJacob Faibussowitsch PetscCall(PetscFree(linesearch->data)); 394*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 39571b4ebd2SPeter Brune } 39671b4ebd2SPeter Brune 397d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch, PetscOptionItems *PetscOptionsObject) 398d71ae5a4SJacob Faibussowitsch { 399eb23ba39SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 40071b4ebd2SPeter Brune 4016e111a19SKarl Rupp PetscFunctionBegin; 402d0609cedSBarry Smith PetscOptionsHeadBegin(PetscOptionsObject, "SNESLineSearch BT options"); 4039566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL)); 404d0609cedSBarry Smith PetscOptionsHeadEnd(); 405*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 40671b4ebd2SPeter Brune } 40771b4ebd2SPeter Brune 408954494b2SJed Brown /*MC 409db609ea7SPeter Brune SNESLINESEARCHBT - Backtracking line search. 410954494b2SJed Brown 411db609ea7SPeter Brune This line search finds the minimum of a polynomial fitting of the L2 norm of the 412f6dfbefdSBarry Smith function or the objective function if it is provided with `SNESSetObjective()`. If this fit does not satisfy the conditions for progress, the interval shrinks 413db609ea7SPeter Brune and the fit is reattempted at most max_it times or until lambda is below minlambda. 414954494b2SJed Brown 415cd7522eaSPeter Brune Options Database Keys: 4163eec2f6aSPierre Jolivet + -snes_linesearch_alpha <1e\-4> - slope descent parameter 417db609ea7SPeter Brune . -snes_linesearch_damping <1.0> - initial step length 418eb23ba39SBarry Smith . -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the 419e1bc860dSBarry Smith step is scaled back to be of this length at the beginning of the line search 420db609ea7SPeter Brune . -snes_linesearch_max_it <40> - maximum number of shrinking step 4213eec2f6aSPierre Jolivet . -snes_linesearch_minlambda <1e\-12> - minimum step length allowed 422cd7522eaSPeter Brune - -snes_linesearch_order <cubic,quadratic> - order of the approximation 423cd7522eaSPeter Brune 424954494b2SJed Brown Level: advanced 425954494b2SJed Brown 426f6dfbefdSBarry Smith Note: 427e55accfeSBarryFSmith This line search will always produce a step that is less than or equal to, in length, the full step size. 428db609ea7SPeter Brune 429f6dfbefdSBarry Smith Reference: 430f6dfbefdSBarry Smith . - * - "Numerical Methods for Unconstrained Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325. 431f6dfbefdSBarry Smith 432f6dfbefdSBarry Smith .seealso: `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()` 433954494b2SJed Brown M*/ 434d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch) 435d71ae5a4SJacob Faibussowitsch { 436f1c6b773SPeter Brune SNESLineSearch_BT *bt; 43771b4ebd2SPeter Brune 43871b4ebd2SPeter Brune PetscFunctionBegin; 439f1c6b773SPeter Brune linesearch->ops->apply = SNESLineSearchApply_BT; 440f1c6b773SPeter Brune linesearch->ops->destroy = SNESLineSearchDestroy_BT; 441f1c6b773SPeter Brune linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT; 4420298fd71SBarry Smith linesearch->ops->reset = NULL; 4437f1410a3SPeter Brune linesearch->ops->view = SNESLineSearchView_BT; 4440298fd71SBarry Smith linesearch->ops->setup = NULL; 44571b4ebd2SPeter Brune 4464dfa11a4SJacob Faibussowitsch PetscCall(PetscNew(&bt)); 447f5af7f23SKarl Rupp 44871b4ebd2SPeter Brune linesearch->data = (void *)bt; 44971b4ebd2SPeter Brune linesearch->max_its = 40; 450b000cd8dSPeter Brune linesearch->order = SNES_LINESEARCH_ORDER_CUBIC; 45171b4ebd2SPeter Brune bt->alpha = 1e-4; 452*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 45371b4ebd2SPeter Brune } 454