diff --git a/CMakeLists.txt b/CMakeLists.txt index 43b627b0..960ed3b3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -73,7 +73,6 @@ set(STA_SOURCE dcalc/DelayCalcBase.cc dcalc/DmpCeff.cc dcalc/DmpDelayCalc.cc - dcalc/FindRoot.cc dcalc/GraphDelayCalc.cc dcalc/LumpedCapDelayCalc.cc dcalc/NetCaps.cc diff --git a/dcalc/DmpCeff.cc b/dcalc/DmpCeff.cc index c00f85aa..fc8d4d95 100644 --- a/dcalc/DmpCeff.cc +++ b/dcalc/DmpCeff.cc @@ -92,13 +92,14 @@ gateModelRd(const LibertyCell *cell, double c1, const Pvt *pvt, bool pocv_enabled); +template static void newtonRaphson(const int max_iter, double x[], const int n, const double x_tol, // eval(state) is called to fill fvec and fjac. - function eval, + const Eval &eval, // Temporaries supplied by caller. double *fvec, double **fjac, @@ -494,9 +495,9 @@ DmpAlg::findVoCrossing(double vth, double t_lower, double t_upper) { - FindRootFunc vo_func = [&] (double t, - double &y, - double &dy) { + const auto vo_func = [&] (double t, + double &y, + double &dy) { double vo, vo_dt; Vo(t, vo, vo_dt); y = vo - vth; @@ -612,9 +613,9 @@ DmpAlg::findVlCrossing(double vth, double t_lower, double t_upper) { - FindRootFunc vl_func = [&] (double t, - double &y, - double &dy) { + const auto vl_func = [&] (double t, + double &y, + double &dy) { double vl, vl_dt; Vl(t, vl, vl_dt); y = vl - vth; @@ -1278,12 +1279,14 @@ DmpZeroC2::voCrossingUpperBound() // x_tol is percentage that all changes in x must be less than (1.0 = 100%). // Eval(state) is called to fill fvec and fjac (returns false if fails). // Return error msg on failure. +// Eval should be callable like void(void) +template static void newtonRaphson(const int max_iter, double x[], const int size, const double x_tol, - function eval, + const Eval &eval, // Temporaries supplied by caller. double *fvec, double **fjac, @@ -1291,6 +1294,7 @@ newtonRaphson(const int max_iter, double *p, double *scale) { + static_assert(std::is_invocable_v, "Must be a callable with no arguments"); for (int k = 0; k < max_iter; k++) { eval(); for (int i = 0; i < size; i++) diff --git a/dcalc/FindRoot.cc b/dcalc/FindRoot.cc deleted file mode 100644 index 78826d62..00000000 --- a/dcalc/FindRoot.cc +++ /dev/null @@ -1,106 +0,0 @@ -// OpenSTA, Static Timing Analyzer -// Copyright (c) 2023, Parallax Software, Inc. -// -// This program is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with this program. If not, see . - -#include "FindRoot.hh" - -#include // abs - -namespace sta { - -using std::abs; - -double -findRoot(FindRootFunc func, - double x1, - double x2, - double x_tol, - int max_iter, - // Return value. - bool &fail) -{ - double y1, y2, dy1; - func(x1, y1, dy1); - func(x2, y2, dy1); - return findRoot(func, x1, y1, x2, y2, x_tol, max_iter, fail); -} - -double -findRoot(FindRootFunc func, - double x1, - double y1, - double x2, - double y2, - double x_tol, - int max_iter, - // Return value. - bool &fail) -{ - if ((y1 > 0.0 && y2 > 0.0) || (y1 < 0.0 && y2 < 0.0)) { - // Initial bounds do not surround a root. - fail = true; - return 0.0; - } - - if (y1 == 0.0) { - fail = false; - return x1; - } - - if (y2 == 0.0) { - fail = false; - return x2; - } - - if (y1 > 0.0) - // Swap x1/x2 so func(x1) < 0. - std::swap(x1, x2); - double root = (x1 + x2) * 0.5; - double dx_prev = abs(x2 - x1); - double dx = dx_prev; - double y, dy; - func(root, y, dy); - for (int iter = 0; iter < max_iter; iter++) { - // Newton/raphson out of range. - if ((((root - x2) * dy - y) * ((root - x1) * dy - y) > 0.0) - // Not decreasing fast enough. - || (abs(2.0 * y) > abs(dx_prev * dy))) { - // Bisect x1/x2 interval. - dx_prev = dx; - dx = (x2 - x1) * 0.5; - root = x1 + dx; - } - else { - dx_prev = dx; - dx = y / dy; - root -= dx; - } - if (abs(dx) <= x_tol * abs(root)) { - // Converged. - fail = false; - return root; - } - - func(root, y, dy); - if (y < 0.0) - x1 = root; - else - x2 = root; - } - fail = true; - return root; -} - -} // namespace diff --git a/dcalc/FindRoot.hh b/dcalc/FindRoot.hh index 1a62df8b..85a828e1 100644 --- a/dcalc/FindRoot.hh +++ b/dcalc/FindRoot.hh @@ -16,26 +16,33 @@ #pragma once -#include +#include // abs namespace sta { -typedef const std::function FindRootFunc; +using std::abs; +// FindRootFunc should be callable like void(double, double&, double&) +template double -findRoot(FindRootFunc func, +findRoot(const FindRootFunc &func, double x1, double x2, double x_tol, int max_iter, // Return value. - bool &fail); + bool &fail) +{ + double y1, y2, dy1; + func(x1, y1, dy1); + func(x2, y2, dy1); + return findRoot(func, x1, y1, x2, y2, x_tol, max_iter, fail); +} +// FindRootFunc should be callable like void(double, double&, double&) +template double -findRoot(FindRootFunc func, +findRoot(const FindRootFunc &func, double x1, double y1, double x2, @@ -43,6 +50,61 @@ findRoot(FindRootFunc func, double x_tol, int max_iter, // Return value. - bool &fail); + bool &fail) +{ + if ((y1 > 0.0 && y2 > 0.0) || (y1 < 0.0 && y2 < 0.0)) { + // Initial bounds do not surround a root. + fail = true; + return 0.0; + } + + if (y1 == 0.0) { + fail = false; + return x1; + } + + if (y2 == 0.0) { + fail = false; + return x2; + } + + if (y1 > 0.0) + // Swap x1/x2 so func(x1) < 0. + std::swap(x1, x2); + double root = (x1 + x2) * 0.5; + double dx_prev = abs(x2 - x1); + double dx = dx_prev; + double y, dy; + func(root, y, dy); + for (int iter = 0; iter < max_iter; iter++) { + // Newton/raphson out of range. + if ((((root - x2) * dy - y) * ((root - x1) * dy - y) > 0.0) + // Not decreasing fast enough. + || (abs(2.0 * y) > abs(dx_prev * dy))) { + // Bisect x1/x2 interval. + dx_prev = dx; + dx = (x2 - x1) * 0.5; + root = x1 + dx; + } + else { + dx_prev = dx; + dx = y / dy; + root -= dx; + } + if (abs(dx) <= x_tol * abs(root)) { + // Converged. + fail = false; + return root; + } + + func(root, y, dy); + if (y < 0.0) + x1 = root; + else + x2 = root; + } + fail = true; + return root; +} } // namespace