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 4042747ad1SJacob Faibussowitsch .seealso: `SNESLineSearchGetLambda()`, `SNESLineSearchGetTolerances()` `SNESLINESEARCHBT` 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 } 109*6b095a85SStefano 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 { 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) { 180a68bbae5SBarry 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 */ 330*6b095a85SStefano Zampini PetscCall(SNESLineSearchSetLambda(linesearch, lambda)); 3319566063dSJacob Faibussowitsch PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w)); 33271b4ebd2SPeter Brune if (changed_y) { 333*6b095a85SStefano Zampini if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X)); 3341baa6e33SBarry Smith if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W)); 33571b4ebd2SPeter Brune } 336bd4a8a71SBarry Smith if (changed_y || changed_w || objective) { /* recompute the function norm if the step has changed or the objective isn't the norm */ 3379566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->snesfunc)(snes, W, G)); 3389bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 3399bd66eb0SPeter Brune gnorm = fnorm; 3409566063dSJacob Faibussowitsch PetscCall((*linesearch->ops->vinorm)(snes, G, W, &gnorm)); 3419bd66eb0SPeter Brune } else { 3429566063dSJacob Faibussowitsch PetscCall(VecNorm(G, NORM_2, &gnorm)); 3439bd66eb0SPeter Brune } 3449566063dSJacob Faibussowitsch PetscCall(VecNorm(Y, NORM_2, &ynorm)); 345c98378a5SBarry Smith if (PetscIsInfOrNanReal(gnorm)) { 3469566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF)); 3479566063dSJacob Faibussowitsch PetscCall(PetscInfo(snes, "Aborted due to Nan or Inf in function evaluation\n")); 3483ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 349c98378a5SBarry Smith } 350c427506bSPeter Brune } 35171b4ebd2SPeter Brune 35271b4ebd2SPeter Brune /* copy the solution over */ 3539566063dSJacob Faibussowitsch PetscCall(VecCopy(W, X)); 3549566063dSJacob Faibussowitsch PetscCall(VecCopy(G, F)); 3559566063dSJacob Faibussowitsch PetscCall(VecNorm(X, NORM_2, &xnorm)); 3569566063dSJacob Faibussowitsch PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm)); 3573ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 35871b4ebd2SPeter Brune } 35971b4ebd2SPeter Brune 36066976f2fSJacob Faibussowitsch static PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer) 361d71ae5a4SJacob Faibussowitsch { 3627f1410a3SPeter Brune PetscBool iascii; 363d7ca06c0SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 364075cc632SBarry Smith 3657f1410a3SPeter Brune PetscFunctionBegin; 3669566063dSJacob Faibussowitsch PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 3677f1410a3SPeter Brune if (iascii) { 368b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3699566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: cubic\n")); 370b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 3719566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation: quadratic\n")); 3727f1410a3SPeter Brune } 3739566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(viewer, " alpha=%e\n", (double)bt->alpha)); 3747f1410a3SPeter Brune } 3753ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 3767f1410a3SPeter Brune } 3777f1410a3SPeter Brune 378d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch) 379d71ae5a4SJacob Faibussowitsch { 38071b4ebd2SPeter Brune PetscFunctionBegin; 3819566063dSJacob Faibussowitsch PetscCall(PetscFree(linesearch->data)); 3823ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 38371b4ebd2SPeter Brune } 38471b4ebd2SPeter Brune 385d71ae5a4SJacob Faibussowitsch static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch, PetscOptionItems *PetscOptionsObject) 386d71ae5a4SJacob Faibussowitsch { 387eb23ba39SBarry Smith SNESLineSearch_BT *bt = (SNESLineSearch_BT *)linesearch->data; 38871b4ebd2SPeter Brune 3896e111a19SKarl Rupp PetscFunctionBegin; 390d0609cedSBarry Smith PetscOptionsHeadBegin(PetscOptionsObject, "SNESLineSearch BT options"); 3919566063dSJacob Faibussowitsch PetscCall(PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, NULL)); 392d0609cedSBarry Smith PetscOptionsHeadEnd(); 3933ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 39471b4ebd2SPeter Brune } 39571b4ebd2SPeter Brune 396954494b2SJed Brown /*MC 397db609ea7SPeter Brune SNESLINESEARCHBT - Backtracking line search. 398954494b2SJed Brown 399db609ea7SPeter Brune This line search finds the minimum of a polynomial fitting of the L2 norm of the 40032661dccSStefano Zampini function or the objective function if it is provided with `SNESSetObjective()`. 40132661dccSStefano Zampini If this fit does not satisfy the conditions for progress, the interval shrinks 402db609ea7SPeter Brune and the fit is reattempted at most max_it times or until lambda is below minlambda. 403954494b2SJed Brown 404cd7522eaSPeter Brune Options Database Keys: 4053eec2f6aSPierre Jolivet + -snes_linesearch_alpha <1e\-4> - slope descent parameter 406db609ea7SPeter Brune . -snes_linesearch_damping <1.0> - initial step length 407eb23ba39SBarry Smith . -snes_linesearch_maxstep <length> - if the length the initial step is larger than this then the 408e1bc860dSBarry Smith step is scaled back to be of this length at the beginning of the line search 409db609ea7SPeter Brune . -snes_linesearch_max_it <40> - maximum number of shrinking step 4103eec2f6aSPierre Jolivet . -snes_linesearch_minlambda <1e\-12> - minimum step length allowed 41132661dccSStefano Zampini - -snes_linesearch_order <1,2,3> - order of the approximation. With order 1, it performs a simple backtracking without any curve fitting 412cd7522eaSPeter Brune 413954494b2SJed Brown Level: advanced 414954494b2SJed Brown 415f6dfbefdSBarry Smith Note: 416e55accfeSBarryFSmith This line search will always produce a step that is less than or equal to, in length, the full step size. 417db609ea7SPeter Brune 418f6dfbefdSBarry Smith Reference: 419f6dfbefdSBarry Smith . - * - "Numerical Methods for Unconstrained Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325. 420f6dfbefdSBarry Smith 421f6dfbefdSBarry Smith .seealso: `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()` 422954494b2SJed Brown M*/ 423d71ae5a4SJacob Faibussowitsch PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch) 424d71ae5a4SJacob Faibussowitsch { 425f1c6b773SPeter Brune SNESLineSearch_BT *bt; 42671b4ebd2SPeter Brune 42771b4ebd2SPeter Brune PetscFunctionBegin; 428f1c6b773SPeter Brune linesearch->ops->apply = SNESLineSearchApply_BT; 429f1c6b773SPeter Brune linesearch->ops->destroy = SNESLineSearchDestroy_BT; 430f1c6b773SPeter Brune linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT; 4310298fd71SBarry Smith linesearch->ops->reset = NULL; 4327f1410a3SPeter Brune linesearch->ops->view = SNESLineSearchView_BT; 4330298fd71SBarry Smith linesearch->ops->setup = NULL; 43471b4ebd2SPeter Brune 4354dfa11a4SJacob Faibussowitsch PetscCall(PetscNew(&bt)); 436f5af7f23SKarl Rupp 43771b4ebd2SPeter Brune linesearch->data = (void *)bt; 43871b4ebd2SPeter Brune linesearch->max_its = 40; 439b000cd8dSPeter Brune linesearch->order = SNES_LINESEARCH_ORDER_CUBIC; 44071b4ebd2SPeter Brune bt->alpha = 1e-4; 4413ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 44271b4ebd2SPeter Brune } 443