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; 263ba16761SJacob 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; 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)); 71*a68bbae5SBarry 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 { 11732661dccSStefano Zampini f = 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)); 12432661dccSStefano Zampini initslope *= 2.0; /* all the fitting assumes f = ||grad||^2, not f = 0.5 ||grad||^2 */ 1258a430afbSPeter Brune } else { 1268a430afbSPeter Brune /* slope comes from the normal equations */ 1279566063dSJacob Faibussowitsch PetscCall(MatMult(jac, Y, W)); 1289566063dSJacob Faibussowitsch PetscCall(VecDotRealPart(F, W, &initslope)); 12971b4ebd2SPeter Brune if (initslope > 0.0) initslope = -initslope; 13071b4ebd2SPeter Brune if (initslope == 0.0) initslope = -1.0; 1318a430afbSPeter Brune } 13271b4ebd2SPeter Brune 133df019d78SBarry Smith while (PETSC_TRUE) { 1349566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -lambda, Y, X)); 1351baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 136e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 1379566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while checking full step length!\n")); 13871b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1399566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 1403ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 14171b4ebd2SPeter Brune } 1428a430afbSPeter Brune 143bd4a8a71SBarry Smith if (objective) { 1449566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes, W, &g)); 1458a430afbSPeter Brune } else { 1469566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 1479bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 1489bd66eb0SPeter Brune gnorm = fnorm; 1499566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 1509bd66eb0SPeter Brune } else { 1519566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 1528a430afbSPeter Brune } 15379e7e81bSJed Brown g = PetscSqr(gnorm); 1549bd66eb0SPeter Brune } 1559566063dSJacob Faibussowitsch PetscCall(SNESLineSearchMonitor(linesearch)); 1569bd66eb0SPeter Brune 157df019d78SBarry Smith if (!PetscIsInfOrNanReal(g)) break; 158df019d78SBarry Smith if (monitor) { 1599566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 1609566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: objective function at lambdas = %g is Inf or Nan, cutting lambda\n", (double)lambda)); 1619566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 162df019d78SBarry Smith } 163ad540459SPierre Jolivet if (lambda <= minlambda) SNESCheckFunctionNorm(snes, g); 16432661dccSStefano Zampini lambda *= .5; 165df019d78SBarry Smith } 166df019d78SBarry Smith 16748a46eb9SPierre Jolivet if (!objective) PetscCall(PetscInfo(snes, "Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm)); 1688a430afbSPeter Brune if (.5 * g <= .5 * f + lambda * alpha * initslope) { /* Sufficient reduction or step tolerance convergence */ 16971b4ebd2SPeter Brune if (monitor) { 1709566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 171bd4a8a71SBarry Smith if (!objective) { 1729566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm)); 1738a430afbSPeter Brune } else { 1749566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Using full step: obj %14.12e obj %14.12e\n", (double)f, (double)g)); 1758a430afbSPeter Brune } 1769566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 17771b4ebd2SPeter Brune } 17871b4ebd2SPeter Brune } else { 179c21ba15cSPeter Brune if (stol * xnorm > ynorm) { 180*a68bbae5SBarry Smith /* Since the full step didn't give sufficient decrease and the step is tiny, exit */ 1819566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, fnorm, ynorm)); 1829566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 183c61ba1e2SPeter Brune if (monitor) { 1849566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 18563a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Ended due to ynorm < stol*xnorm (%14.12e < %14.12e).\n", (double)ynorm, (double)(stol * xnorm))); 1869566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 187c61ba1e2SPeter Brune } 1883ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 189c21ba15cSPeter Brune } 19032661dccSStefano Zampini /* Here to avoid -Wmaybe-uninitiliazed warnings */ 19171b4ebd2SPeter Brune lambdaprev = lambda; 1928a430afbSPeter Brune gprev = g; 19332661dccSStefano Zampini if (linesearch->order != SNES_LINESEARCH_ORDER_LINEAR) { 19432661dccSStefano Zampini /* Fit points with quadratic */ 19532661dccSStefano Zampini lambdatemp = -initslope / (g - f - 2.0 * lambda * initslope); 19671b4ebd2SPeter Brune if (lambdatemp > .5 * lambda) lambdatemp = .5 * lambda; 19771b4ebd2SPeter Brune if (lambdatemp <= .1 * lambda) lambda = .1 * lambda; 19871b4ebd2SPeter Brune else lambda = lambdatemp; 19971b4ebd2SPeter Brune 2009566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -lambda, Y, X)); 2011baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 202e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 20363a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while attempting quadratic backtracking! %" PetscInt_FMT " \n", snes->nfuncs)); 20471b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2059566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 2063ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 20771b4ebd2SPeter Brune } 208bd4a8a71SBarry Smith if (objective) { 2099566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes, W, &g)); 2108a430afbSPeter Brune } else { 2119566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 2129bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 2139bd66eb0SPeter Brune gnorm = fnorm; 2149566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 2159bd66eb0SPeter Brune } else { 2169566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 2179bd66eb0SPeter Brune } 218f1ab69e3SPeter Brune g = PetscSqr(gnorm); 2198a430afbSPeter Brune } 220c98378a5SBarry Smith if (PetscIsInfOrNanReal(g)) { 2219566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 2229566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n")); 2233ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 224c98378a5SBarry Smith } 22571b4ebd2SPeter Brune if (monitor) { 2269566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 227bd4a8a71SBarry Smith if (!objective) { 2289566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: gnorm after quadratic fit %14.12e\n", (double)gnorm)); 2298a430afbSPeter Brune } else { 2309566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: obj after quadratic fit %14.12e\n", (double)g)); 2318a430afbSPeter Brune } 2329566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 23371b4ebd2SPeter Brune } 23432661dccSStefano Zampini } 23532661dccSStefano Zampini if (linesearch->order != SNES_LINESEARCH_ORDER_LINEAR && .5 * g < .5 * f + lambda * alpha * initslope) { /* sufficient reduction */ 23671b4ebd2SPeter Brune if (monitor) { 2379566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 2389566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: Quadratically determined step, lambda=%18.16e\n", (double)lambda)); 2399566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 24071b4ebd2SPeter Brune } 24171b4ebd2SPeter Brune } else { 24271b4ebd2SPeter Brune /* Fit points with cubic */ 24371b4ebd2SPeter Brune for (count = 0; count < max_its; count++) { 24471b4ebd2SPeter Brune if (lambda <= minlambda) { 24571b4ebd2SPeter Brune if (monitor) { 2469566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 24763a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: unable to find good step length! After %" PetscInt_FMT " tries \n", count)); 248bd4a8a71SBarry Smith if (!objective) { 2499371c9d4SSatish 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)); 2508a430afbSPeter Brune } else { 25132661dccSStefano 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 / 2.0)); 2528a430afbSPeter Brune } 2539566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 25471b4ebd2SPeter Brune } 2559566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 2563ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 25771b4ebd2SPeter Brune } 258b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 2598a430afbSPeter Brune t1 = .5 * (g - f) - lambda * initslope; 2608a430afbSPeter Brune t2 = .5 * (gprev - f) - lambdaprev * initslope; 26171b4ebd2SPeter Brune a = (t1 / (lambda * lambda) - t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev); 26271b4ebd2SPeter Brune b = (-lambdaprev * t1 / (lambda * lambda) + lambda * t2 / (lambdaprev * lambdaprev)) / (lambda - lambdaprev); 26371b4ebd2SPeter Brune d = b * b - 3 * a * initslope; 26471b4ebd2SPeter Brune if (d < 0.0) d = 0.0; 265f5af7f23SKarl Rupp if (a == 0.0) lambdatemp = -initslope / (2.0 * b); 266f5af7f23SKarl Rupp else lambdatemp = (-b + PetscSqrtReal(d)) / (3.0 * a); 267b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 2688a430afbSPeter Brune lambdatemp = -initslope / (g - f - 2.0 * initslope); 26932661dccSStefano Zampini } else if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) { /* Just backtrack */ 27032661dccSStefano Zampini lambdatemp = .5 * lambda; 27132661dccSStefano Zampini } else SETERRQ(PetscObjectComm((PetscObject)linesearch), PETSC_ERR_SUP, "Line search order %" PetscInt_FMT " for type bt", linesearch->order); 27271b4ebd2SPeter Brune lambdaprev = lambda; 2738a430afbSPeter Brune gprev = g; 27471b4ebd2SPeter Brune if (lambdatemp > .5 * lambda) lambdatemp = .5 * lambda; 27571b4ebd2SPeter Brune if (lambdatemp <= .1 * lambda) lambda = .1 * lambda; 27671b4ebd2SPeter Brune else lambda = lambdatemp; 2779566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -lambda, Y, X)); 2781baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 279e71169deSBarry Smith if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 28063a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Exceeded maximum function evaluations, while looking for good step length! %" PetscInt_FMT " \n", count)); 28148a46eb9SPierre 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)); 2829566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_FUNCTION)); 28371b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2843ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 28571b4ebd2SPeter Brune } 286bd4a8a71SBarry Smith if (objective) { 2879566063dSJacob Faibussowitsch PetscCall(SNESComputeObjective(snes, W, &g)); 2888a430afbSPeter Brune } else { 2899566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 2909bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 2919bd66eb0SPeter Brune gnorm = fnorm; 2929566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 2939bd66eb0SPeter Brune } else { 2949566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 2958a430afbSPeter Brune } 296f1ab69e3SPeter Brune g = PetscSqr(gnorm); 2979bd66eb0SPeter Brune } 2984a2f8832SBarry Smith if (PetscIsInfOrNanReal(g)) { 2999566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 3009566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n")); 3013ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 302c98378a5SBarry Smith } 3038a430afbSPeter Brune if (.5 * g < .5 * f + lambda * alpha * initslope) { /* is reduction enough? */ 30471b4ebd2SPeter Brune if (monitor) { 3059566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 306bd4a8a71SBarry Smith if (!objective) { 30732661dccSStefano Zampini PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: %s step, current gnorm %14.12e lambda=%18.16e\n", ordStr[linesearch->order - 1], (double)gnorm, (double)lambda)); 3089566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3098a430afbSPeter Brune } else { 31032661dccSStefano Zampini PetscCall(PetscViewerASCIIPrintf(monitor, " Line search: %s step, obj %14.12e lambda=%18.16e\n", ordStr[linesearch->order - 1], (double)g, (double)lambda)); 3119566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3128a430afbSPeter Brune } 31371b4ebd2SPeter Brune } 31471b4ebd2SPeter Brune break; 315f5af7f23SKarl Rupp } else if (monitor) { 3169566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel)); 317bd4a8a71SBarry Smith if (!objective) { 31832661dccSStefano 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)); 3199566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3208a430afbSPeter Brune } else { 32132661dccSStefano 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)); 3229566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel)); 3238a430afbSPeter Brune } 32471b4ebd2SPeter Brune } 32571b4ebd2SPeter Brune } 32671b4ebd2SPeter Brune } 32771b4ebd2SPeter Brune } 32871b4ebd2SPeter Brune 32971b4ebd2SPeter Brune /* postcheck */ 3300c1ac483SGlenn Hammond /* update Y to lambda*Y so that W is consistent with X - lambda*Y */ 3319566063dSJacob Faibussowitsch PetscCall(VecScale(Y, lambda)); 3329566063dSJacob Faibussowitsch PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w)); 33371b4ebd2SPeter Brune if (changed_y) { 3349566063dSJacob Faibussowitsch PetscCall(VecWAXPY(W, -1.0, Y, X)); 3351baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 33671b4ebd2SPeter Brune } 337bd4a8a71SBarry Smith if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */ 3389566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 3399bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 3409bd66eb0SPeter Brune gnorm = fnorm; 3419566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 3429bd66eb0SPeter Brune } else { 3439566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 3449bd66eb0SPeter Brune } 3459566063dSJacob Faibussowitsch PetscCall(VecNorm(Y, NORM_2, &ynorm)); 346c98378a5SBarry Smith if (PetscIsInfOrNanReal(gnorm)) { 3479566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 3489566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n")); 3493ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 350c98378a5SBarry Smith } 351c427506bSPeter Brune } 35271b4ebd2SPeter Brune 35371b4ebd2SPeter Brune /* copy the solution over */ 3549566063dSJacob Faibussowitsch PetscCall(VecCopy(W, X)); 3559566063dSJacob Faibussowitsch PetscCall(VecCopy(G, F)); 3569566063dSJacob Faibussowitsch PetscCall(VecNorm(X, NORM_2, &xnorm)); 3579566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetLambda(linesearch, lambda)); 3589566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm)); 3593ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 36071b4ebd2SPeter Brune } 36171b4ebd2SPeter Brune 362d71ae5a4SJacob Faibussowitsch PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer) 363d71ae5a4SJacob Faibussowitsch { 3647f1410a3SPeter Brune PetscBool iascii; 365d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 366075cc632SBarry Smith 3677f1410a3SPeter Brune PetscFunctionBegin; 3689566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 3697f1410a3SPeter Brune if (iascii) { 370b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3719566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: cubic\n")); 372b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 3739566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: quadratic\n")); 3747f1410a3SPeter Brune } 3759566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " alpha=%e\n", (double)bt->alpha)); 3767f1410a3SPeter Brune } 3773ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 3787f1410a3SPeter Brune } 3797f1410a3SPeter Brune 380d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch) 381d71ae5a4SJacob Faibussowitsch { 38271b4ebd2SPeter Brune PetscFunctionBegin; 3839566063dSJacob Faibussowitsch PetscCall(PetscFree(linesearch->data)); 3843ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 38571b4ebd2SPeter Brune } 38671b4ebd2SPeter Brune 387d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch, PetscOptionItems *PetscOptionsObject) 388d71ae5a4SJacob Faibussowitsch { 389eb23ba39SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 39071b4ebd2SPeter Brune 3916e111a19SKarl Rupp PetscFunctionBegin; 392d0609cedSBarry Smith PetscOptionsHeadBegin(PetscOptionsObject, "SNESLineSearch BT options"); 3939566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL)); 394d0609cedSBarry Smith PetscOptionsHeadEnd(); 3953ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 39671b4ebd2SPeter Brune } 39771b4ebd2SPeter Brune 398954494b2SJed Brown /*MC 399db609ea7SPeter Brune SNESLINESEARCHBT - Backtracking line search. 400954494b2SJed Brown 401db609ea7SPeter Brune This line search finds the minimum of a polynomial fitting of the L2 norm of the 40232661dccSStefano Zampini function or the objective function if it is provided with `SNESSetObjective()`. 40332661dccSStefano Zampini If this fit does not satisfy the conditions for progress, the interval shrinks 404db609ea7SPeter Brune and the fit is reattempted at most max_it times or until lambda is below minlambda. 405954494b2SJed Brown 406cd7522eaSPeter Brune Options Database Keys: 4073eec2f6aSPierre Jolivet + -snes_linesearch_alpha <1e\-4> - slope descent parameter 408db609ea7SPeter Brune . -snes_linesearch_damping <1.0> - initial step length 409eb23ba39SBarry Smith . -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the 410e1bc860dSBarry Smith step is scaled back to be of this length at the beginning of the line search 411db609ea7SPeter Brune . -snes_linesearch_max_it <40> - maximum number of shrinking step 4123eec2f6aSPierre Jolivet . -snes_linesearch_minlambda <1e\-12> - minimum step length allowed 41332661dccSStefano Zampini - -snes_linesearch_order <1,2,3> - order of the approximation. With order 1, it performs a simple backtracking without any curve fitting 414cd7522eaSPeter Brune 415954494b2SJed Brown Level: advanced 416954494b2SJed Brown 417f6dfbefdSBarry Smith Note: 418e55accfeSBarryFSmith This line search will always produce a step that is less than or equal to, in length, the full step size. 419db609ea7SPeter Brune 420f6dfbefdSBarry Smith Reference: 421f6dfbefdSBarry Smith . - * - "Numerical Methods for Unconstrained Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325. 422f6dfbefdSBarry Smith 423f6dfbefdSBarry Smith .seealso: `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()` 424954494b2SJed Brown M*/ 425d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch) 426d71ae5a4SJacob Faibussowitsch { 427f1c6b773SPeter Brune SNESLineSearch_BT *bt; 42871b4ebd2SPeter Brune 42971b4ebd2SPeter Brune PetscFunctionBegin; 430f1c6b773SPeter Brune linesearch->ops->apply = SNESLineSearchApply_BT; 431f1c6b773SPeter Brune linesearch->ops->destroy = SNESLineSearchDestroy_BT; 432f1c6b773SPeter Brune linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT; 4330298fd71SBarry Smith linesearch->ops->reset = NULL; 4347f1410a3SPeter Brune linesearch->ops->view = SNESLineSearchView_BT; 4350298fd71SBarry Smith linesearch->ops->setup = NULL; 43671b4ebd2SPeter Brune 4374dfa11a4SJacob Faibussowitsch PetscCall(PetscNew(&bt)); 438f5af7f23SKarl Rupp 43971b4ebd2SPeter Brune linesearch->data = (void *)bt; 44071b4ebd2SPeter Brune linesearch->max_its = 40; 441b000cd8dSPeter Brune linesearch->order = SNES_LINESEARCH_ORDER_CUBIC; 44271b4ebd2SPeter Brune bt->alpha = 1e-4; 4433ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 44471b4ebd2SPeter Brune } 445