Refactor DMP solver to use Eigen and implement Determinant Guarded solver (#450)

This commit refactors the Dartu-Menezes-Pileggi (DMP) effective capacitance
algorithm solver to use the Eigen library, and implements a Determinant Guarded
solver to handle singular and near-singular Jacobians safely and efficiently.

Key changes:
1. Refactored the DMP solver to use Eigen stack-based data structures,
   replacing the custom Crout LU decomposition and solver.
2. Extracted the linear solver logic into a dedicated helper function
   DmpAlg::solveNewtonStep.
3. Implemented a "Determinant Guarded" solver:
   - Manually checks the determinant of the Jacobian (safety guard).
   - Throws a DmpError if the determinant is dangerously close to zero (|det| < 1e-12),
     safely triggering the lumped capacitance fallback.
   - Otherwise, solves using Eigen's highly optimized analytical inverse (fast path).
4. Added documentation in the comments on how to easily swap back to LU
   decomposition with partial pivoting (PartialPivLU) if any numerical
   precision issues arise in the future.

Benchmark Results (Optimized Release Build, 5,288 DMP calls):

| Test Case | Calls | 1. Original (Crout LU) | 2. Eigen (PartialPivLU) | 3. Eigen (Determinant Guarded) | Speedup (3 vs 1) | Speedup (3 vs 2) |
| :--- | :---: | :---: | :---: | :---: | :---: | :---: |
| power | 2,608 | 1.8822 ms | 1.6791 ms | 1.2702 ms | +32.5% | +24.4% |
| power_vcd | 2,608 | 1.8729 ms | 1.6704 ms | 1.2768 ms | +31.8% | +23.6% |
| spef_parasitics | 24 | 0.0193 ms | 0.0187 ms | 0.0140 ms | +27.5% | +25.1% |
| mcmm3 | 48 | 0.0728 ms | 0.0756 ms | 0.0817 ms | -12.2% | -8.1% |
| Total DMP Time | 5,288 | 3.8472 ms | 3.4438 ms | 2.6427 ms | +31.3% | +23.3% |
This commit is contained in:
Drew Lewis 2026-06-18 16:52:14 -04:00 committed by GitHub
parent 133580eeae
commit 832c62b579
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
2 changed files with 194 additions and 202 deletions

View File

@ -31,6 +31,7 @@
// slew voltage is matched instead of y20 in eqn 12.
#include "DmpCeff.hh"
#include <Eigen/Dense>
#include <algorithm>
#include <array>
@ -137,20 +138,23 @@ DmpAlg::init(const LibertyLibrary *drvr_library,
void
DmpAlg::findDriverParams(double ceff)
{
Eigen::Vector3d x = Eigen::Vector3d::Zero();
if (nr_order_ == 3)
x_[DmpParam::ceff] = ceff;
x[DmpParam::ceff] = ceff;
auto [t_vth, t_vl, slew] = gateDelays(ceff);
// Scale slew to 0-100%
double dt = slew / (vh_ - vl_);
double t0 = t_vth + std::log(1.0 - vth_) * rd_ * ceff - vth_ * dt;
x_[DmpParam::dt] = dt;
x_[DmpParam::t0] = t0;
newtonRaphson();
t0_ = x_[DmpParam::t0];
dt_ = x_[DmpParam::dt];
x[DmpParam::dt] = dt;
x[DmpParam::t0] = t0;
newtonRaphson(x);
t0_ = x[DmpParam::t0];
dt_ = x[DmpParam::dt];
if (nr_order_ == 3)
ceff_ = x[DmpParam::ceff];
debugPrint(debug_, "dmp_ceff", 3, " t0 = {} dt = {} ceff = {}",
units_->timeUnit()->asString(t0_), units_->timeUnit()->asString(dt_),
units_->capacitanceUnit()->asString(x_[DmpParam::ceff]));
units_->capacitanceUnit()->asString(ceff_));
if (debug_->check("dmp_ceff", 4))
showVo();
}
@ -241,21 +245,21 @@ DmpAlg::y0dcl(double t,
}
void
DmpAlg::showX()
DmpAlg::showX(const Eigen::Vector3d &x)
{
for (int i = 0; i < nr_order_; i++)
report_->report("{:4} {:12.3e}", dmp_param_index_strings[i], x_[i]);
report_->report("{:4} {:12.3e}", dmp_param_index_strings[i], x[i]);
}
void
DmpAlg::showFvec()
DmpAlg::showFvec(const Eigen::Vector3d &fvec)
{
for (int i = 0; i < nr_order_; i++)
report_->report("{:4} {:12.3e}", dmp_func_index_strings[i], fvec_[i]);
report_->report("{:4} {:12.3e}", dmp_func_index_strings[i], fvec[i]);
}
void
DmpAlg::showJacobian()
DmpAlg::showJacobian(const Eigen::Matrix3d &fjac)
{
std::string line = " ";
for (int j = 0; j < nr_order_; j++)
@ -265,7 +269,7 @@ DmpAlg::showJacobian()
line.clear();
line += sta::format("{:4} ", dmp_func_index_strings[i]);
for (int j = 0; j < nr_order_; j++)
line += sta::format("{:12.3e} ", fjac_[i][j]);
line += sta::format("{:12.3e} ", fjac(i, j));
report_->reportLine(line);
}
}
@ -505,7 +509,9 @@ DmpCap::loadDelaySlew(const Pin *,
}
void
DmpCap::evalDmpEqns()
DmpCap::evalDmpEqns(Eigen::Vector3d &,
Eigen::Vector3d &,
Eigen::Matrix3d &)
{
}
@ -584,7 +590,6 @@ DmpPi::gateDelaySlew()
double slew = 0.0;
try {
findDriverParamsPi();
ceff_ = x_[DmpParam::ceff];
auto [table_delay, table_slew] = gateCapDelaySlew(ceff_);
delay = table_delay;
// slew = table_slew;
@ -623,60 +628,85 @@ DmpPi::findDriverParamsPi()
// Given x_ as a vector of input parameters, fill fvec_ with the
// equations evaluated at x_ and fjac_ with the jacobian evaluated at x_.
void
DmpPi::evalDmpEqns()
DmpPi::evalDmpEqns(Eigen::Vector3d &x,
Eigen::Vector3d &fvec,
Eigen::Matrix3d &fjac)
{
double t0 = x_[DmpParam::t0];
double dt = x_[DmpParam::dt];
double ceff = x_[DmpParam::ceff];
const double t0 = x[DmpParam::t0];
const double dt = x[DmpParam::dt];
const double ceff = x[DmpParam::ceff];
if (ceff < 0.0)
// Validate bounds to prevent mathematical domain errors.
if (ceff < 0.0) {
throw DmpError("eqn eval failed: ceff < 0");
if (ceff > (c1_ + c2_))
}
if (ceff > (c1_ + c2_)) {
throw DmpError("eqn eval failed: ceff > c2 + c1");
}
if (dt <= 0.0) {
throw DmpError("eqn eval failed: dt < 0");
}
auto [t_vth, t_vl, slew] = gateDelays(ceff);
if (slew == 0.0)
if (slew == 0.0) {
throw DmpError("eqn eval failed: slew = 0");
}
double ceff_time = slew / (vh_ - vl_);
ceff_time = std::min(ceff_time, 1.4 * dt);
// ceff_time is bounded by 1.4 * dt.
const double ceff_time = std::min(slew / (vh_ - vl_), 1.4 * dt);
if (dt <= 0.0)
throw DmpError("eqn eval failed: dt < 0");
// Pre-calculate exponential terms to avoid redundant calls to
// transcendental functions.
const double exp_p1_dt = exp2(-p1_ * dt);
const double exp_p2_dt = exp2(-p2_ * dt);
const double exp_dt_rd_ceff = exp2(-dt / (rd_ * ceff));
double exp_p1_dt = exp2(-p1_ * dt);
double exp_p2_dt = exp2(-p2_ * dt);
double exp_dt_rd_ceff = exp2(-dt / (rd_ * ceff));
// Evaluate function values (residuals).
const double y50 = y(t_vth, t0, dt, ceff).first;
const double y20 = y(t_vl, t0, dt, ceff).first;
double y50 = y(t_vth, t0, dt, ceff).first;
// Match Vl.
double y20 = y(t_vl, t0, dt, ceff).first;
fvec_[DmpFunc::ipi] = ipiIceff(t0, dt, ceff_time, ceff);
fvec_[DmpFunc::y50] = y50 - vth_;
fvec_[DmpFunc::y20] = y20 - vl_;
fjac_[DmpFunc::ipi][DmpParam::t0] = 0.0;
fjac_[DmpFunc::ipi][DmpParam::dt] =
(-A_ * dt + B_ * dt * exp_p1_dt - (2 * B_ / p1_) * (1.0 - exp_p1_dt)
+ D_ * dt * exp_p2_dt - (2 * D_ / p2_) * (1.0 - exp_p2_dt)
+ rd_ * ceff
* (dt + dt * exp_dt_rd_ceff - 2 * rd_ * ceff * (1.0 - exp_dt_rd_ceff)))
/ (rd_ * dt * dt * dt);
fjac_[DmpFunc::ipi][DmpParam::ceff] =
(2 * rd_ * ceff - dt - (2 * rd_ * ceff + dt) * exp2(-dt / (rd_ * ceff)))
/ (dt * dt);
fvec[DmpFunc::ipi] = ipiIceff(t0, dt, ceff_time, ceff);
fvec[DmpFunc::y50] = y50 - vth_;
fvec[DmpFunc::y20] = y20 - vl_;
std::tie(fjac_[DmpFunc::y20][DmpParam::t0],
fjac_[DmpFunc::y20][DmpParam::dt],
fjac_[DmpFunc::y20][DmpParam::ceff]) = dy(t_vl, t0, dt, ceff);
// Pre-calculate common sub-expressions for the Jacobian derivatives.
const double b_div_p1 = B_ / p1_;
const double d_div_p2 = D_ / p2_;
const double rd_ceff = rd_ * ceff;
std::tie(fjac_[DmpFunc::y50][DmpParam::t0],
fjac_[DmpFunc::y50][DmpParam::dt],
fjac_[DmpFunc::y50][DmpParam::ceff]) = dy(t_vth, t0, dt, ceff);
// Row 1 (Ipi derivatives).
fjac(DmpFunc::ipi, DmpParam::t0) = 0.0;
// Derivative w.r.t dt (broken down into physical terms).
const double term_a = -A_ * dt;
const double term_b =
B_ * dt * exp_p1_dt - 2.0 * b_div_p1 * (1.0 - exp_p1_dt);
const double term_d =
D_ * dt * exp_p2_dt - 2.0 * d_div_p2 * (1.0 - exp_p2_dt);
const double term_rd = rd_ceff
* (dt + dt * exp_dt_rd_ceff - 2.0 * rd_ceff * (1.0 - exp_dt_rd_ceff));
fjac(DmpFunc::ipi, DmpParam::dt) =
(term_a + term_b + term_d + term_rd) / (rd_ * dt * dt * dt);
// Derivative w.r.t ceff (reusing exp_dt_rd_ceff).
const double two_rd_ceff = 2.0 * rd_ceff;
fjac(DmpFunc::ipi, DmpParam::ceff) =
(two_rd_ceff - dt - (two_rd_ceff + dt) * exp_dt_rd_ceff) / (dt * dt);
// Rows 2 & 3 (y20 and y50 derivatives).
std::tie(fjac(DmpFunc::y20, DmpParam::t0),
fjac(DmpFunc::y20, DmpParam::dt),
fjac(DmpFunc::y20, DmpParam::ceff)) = dy(t_vl, t0, dt, ceff);
std::tie(fjac(DmpFunc::y50, DmpParam::t0),
fjac(DmpFunc::y50, DmpParam::dt),
fjac(DmpFunc::y50, DmpParam::ceff)) = dy(t_vth, t0, dt, ceff);
if (debug_->check("dmp_ceff", 4)) {
showX();
showFvec();
showJacobian();
showX(x);
showFvec(fvec);
showJacobian(fjac);
report_->report(".................");
}
}
@ -741,35 +771,37 @@ DmpOnePole::DmpOnePole(StaState *sta) :
}
void
DmpOnePole::evalDmpEqns()
DmpOnePole::evalDmpEqns(Eigen::Vector3d &x,
Eigen::Vector3d &fvec,
Eigen::Matrix3d &fjac)
{
double t0 = x_[DmpParam::t0];
double dt = x_[DmpParam::dt];
double t0 = x[DmpParam::t0];
double dt = x[DmpParam::dt];
auto [t_vth, t_vl, ignore1] = gateDelays(ceff_);
double ignore2;
if (dt <= 0.0)
dt = x_[DmpParam::dt] = (t_vl - t_vth) / 100;
dt = x[DmpParam::dt] = (t_vl - t_vth) / 100;
fvec_[DmpFunc::y50] = y(t_vth, t0, dt, ceff_).first - vth_;
fvec_[DmpFunc::y20] = y(t_vl, t0, dt, ceff_).first - vl_;
fvec[DmpFunc::y50] = y(t_vth, t0, dt, ceff_).first - vth_;
fvec[DmpFunc::y20] = y(t_vl, t0, dt, ceff_).first - vl_;
if (debug_->check("dmp_ceff", 4)) {
showX();
showFvec();
showX(x);
showFvec(fvec);
}
std::tie(fjac_[DmpFunc::y20][DmpParam::t0],
fjac_[DmpFunc::y20][DmpParam::dt],
std::tie(fjac(DmpFunc::y20, DmpParam::t0),
fjac(DmpFunc::y20, DmpParam::dt),
ignore2) = dy(t_vl, t0, dt, ceff_);
std::tie(fjac_[DmpFunc::y50][DmpParam::t0],
fjac_[DmpFunc::y50][DmpParam::dt],
std::tie(fjac(DmpFunc::y50, DmpParam::t0),
fjac(DmpFunc::y50, DmpParam::dt),
ignore2) = dy(t_vth, t0, dt, ceff_);
if (debug_->check("dmp_ceff", 4)) {
showJacobian();
showJacobian(fjac);
report_->report(".................");
}
}
@ -870,136 +902,86 @@ DmpZeroC2::voCrossingUpperBound()
// driver_param_tol_ is the scale that all changes in x must be under (1.0 = 100%).
// evalDmpEqns() fills fvec_ and fjac_.
void
DmpAlg::newtonRaphson()
DmpAlg::newtonRaphson(Eigen::Vector3d &x)
{
for (int k = 0; k < newton_raphson_max_iter_; k++) {
evalDmpEqns();
for (int i = 0; i < nr_order_; i++)
// Right-hand side of linear equations.
p_[i] = -fvec_[i];
luDecomp();
luSolve();
Eigen::Vector3d fvec = Eigen::Vector3d::Zero();
Eigen::Matrix3d fjac = Eigen::Matrix3d::Zero();
Eigen::Vector3d p = Eigen::Vector3d::Zero();
for (int k = 0; k < newton_raphson_max_iter_; k++) {
evalDmpEqns(x, fvec, fjac);
p = solveNewtonStep(fjac, fvec);
// Note: 'auto' on Eigen expressions captures the expression template
// and doesn't form a temporary vector/matrix, avoiding extra
// allocations.
auto p_abs = p.head(nr_order_).array().abs();
auto x_tol = x.head(nr_order_).array().abs() * driver_param_tol_;
bool all_under_x_tol = (p_abs <= x_tol).all();
x.head(nr_order_) += p.head(nr_order_);
bool all_under_x_tol = true;
for (int i = 0; i < nr_order_; i++) {
if (std::abs(p_[i]) > std::abs(x_[i]) * driver_param_tol_)
all_under_x_tol = false;
x_[i] += p_[i];
}
if (all_under_x_tol) {
evalDmpEqns();
return;
}
}
throw DmpError("Newton-Raphson max iterations exceeded");
}
// luDecomp, luSolve based on MatClass from C. R. Birchenhall,
// University of Manchester
// ftp://ftp.mcc.ac.uk/pub/matclass/libmat.tar.Z
// Crout's Method of LU decomposition of square matrix, with implicit
// partial pivoting. fjac_ is overwritten. U is explicit in the upper
// triangle and L is in multiplier form in the subdiagionals i.e. subdiag
// a[i,j] is the multiplier used to eliminate the [i,j] term.
// Solves the linear system J * p = -f (Jacobian * step = -residuals) for the Newton step.
//
// Replaces fjac_[0..nr_order_-1][*] by the LU decomposition.
// index_[0..nr_order_-1] is an output vector of the row permutations.
void
DmpAlg::luDecomp()
// This implementation uses a "Determinant Guarded" solver:
// 1. Manually computes/checks the determinant of the Jacobian (safety guard).
// 2. If the determinant is dangerously close to zero (< 1e-12), throws a DmpError.
// 3. Otherwise, uses Eigen's highly optimized analytical inverse (fast path).
//
// Performance Note:
// Analytical solvers are extremely fast for 2x2 and 3x3 matrices because they
// contain no loops or branching, allowing the compiler to unroll them and use
// SIMD instructions. This yields a ~23% speedup over LU decomposition in optimized builds.
//
// Numerical Stability Note:
// If this analytical approach ever causes numerical issues (e.g., in extremely
// ill-conditioned systems where the determinant is > 1e-12 but still causes loss
// of precision), it can be TRIVIALLY swapped back to a robust LU decomposition
// with partial pivoting by replacing the body of this function with:
//
// Eigen::Vector3d p = Eigen::Vector3d::Zero();
// if (nr_order_ == 2) {
// auto lu = fjac.topLeftCorner<2, 2>().partialPivLu();
// if (std::abs(lu.matrixLU().diagonal().prod()) < 1e-12) {
// throw DmpError("Jacobian is singular (order 2)");
// }
// p.head<2>() = lu.solve(-fvec.head<2>());
// return p;
// }
// auto lu = fjac.partialPivLu();
// if (std::abs(lu.matrixLU().diagonal().prod()) < 1e-12) {
// throw DmpError("Jacobian is singular (order 3)");
// }
// p = lu.solve(-fvec);
// return p;
//
Eigen::Vector3d
DmpAlg::solveNewtonStep(const Eigen::Matrix3d &fjac,
const Eigen::Vector3d &fvec)
{
const int size = nr_order_;
Eigen::Vector3d p = Eigen::Vector3d::Zero();
if (nr_order_ == 2) {
double det = fjac.topLeftCorner<2, 2>().determinant();
if (std::abs(det) < 1e-12) {
throw DmpError("Jacobian is singular (order 2)");
}
p.head<2>() = fjac.topLeftCorner<2, 2>().inverse() * -fvec.head<2>();
return p;
}
// Find implicit scaling factors.
for (int i = 0; i < size; i++) {
double big = 0.0;
for (int j = 0; j < size; j++) {
double temp = std::abs(fjac_[i][j]);
big = std::max(temp, big);
}
if (big == 0.0)
throw DmpError("LU decomposition: no non-zero row element");
scale_[i] = 1.0 / big;
}
int size_1 = size - 1;
for (int j = 0; j < size; j++) {
// Run down jth column from top to diag, to form the elements of U.
for (int i = 0; i < j; i++) {
double sum = fjac_[i][j];
for (int k = 0; k < i; k++)
sum -= fjac_[i][k] * fjac_[k][j];
fjac_[i][j] = sum;
}
// Run down jth subdiag to form the residuals after the elimination
// of the first j-1 subdiags. These residuals diviyded by the
// appropriate diagonal term will become the multipliers in the
// elimination of the jth. subdiag. Find index of largest scaled
// term in imax.
double big = 0.0;
int imax = 0;
for (int i = j; i < size; i++) {
double sum = fjac_[i][j];
for (int k = 0; k < j; k++)
sum -= fjac_[i][k] * fjac_[k][j];
fjac_[i][j] = sum;
double dum = scale_[i] * std::abs(sum);
if (dum >= big) {
big = dum;
imax = i;
}
}
// Permute current row with imax.
if (j != imax) {
// Yes, do so...
for (int k = 0; k < size; k++) {
double dum = fjac_[imax][k];
fjac_[imax][k] = fjac_[j][k];
fjac_[j][k] = dum;
}
scale_[imax] = scale_[j];
}
index_[j] = imax;
// If diag term is not zero divide subdiag to form multipliers.
if (fjac_[j][j] == 0.0)
fjac_[j][j] = tiny_double_;
if (j != size_1) {
double pivot = 1.0 / fjac_[j][j];
for (int i = j + 1; i < size; i++)
fjac_[i][j] *= pivot;
}
}
}
// Solves fjac_ * x = p_ for x, assuming fjac_ is LU form from luDecomp.
// Solution overwrites p_.
void
DmpAlg::luSolve()
{
const int size = nr_order_;
// Transform p_ allowing for leading zeros.
int non_zero = -1;
for (int i = 0; i < size; i++) {
int iperm = index_[i];
double sum = p_[iperm];
p_[iperm] = p_[i];
if (non_zero != -1) {
for (int j = non_zero; j <= i - 1; j++)
sum -= fjac_[i][j] * p_[j];
}
else {
if (sum != 0.0)
non_zero = i;
}
p_[i] = sum;
}
// Backsubstitution.
for (int i = size - 1; i >= 0; i--) {
double sum = p_[i];
for (int j = i + 1; j < size; j++)
sum -= fjac_[i][j] * p_[j];
p_[i] = sum / fjac_[i][i];
double det = fjac.determinant();
if (std::abs(det) < 1e-12) {
throw DmpError("Jacobian is singular (order 3)");
}
p = fjac.inverse() * -fvec;
return p;
}
////////////////////////////////////////////////////////////////

View File

@ -27,6 +27,9 @@
#include <optional>
#include <utility>
#include <Eigen/Core>
#include <Eigen/Dense>
#include "LibertyClass.hh"
#include "LumpedCapDelayCalc.hh"
@ -59,18 +62,21 @@ public:
double elmore);
double ceff() { return ceff_; }
// Given x_ as a vector of input parameters, fill fvec_ with the
// equations evaluated at x_ and fjac_ with the jabobian evaluated at x_.
virtual void evalDmpEqns() = 0;
virtual void
evalDmpEqns(Eigen::Vector3d &x,
Eigen::Vector3d &fvec,
Eigen::Matrix3d &fjac) = 0;
// Output response to vs(t) ramp driving pi model load (vo, dvo_dt).
std::pair<double, double> Vo(double t);
// Load response to driver waveform (vl, dvl/dt).
std::pair<double, double> Vl(double t);
protected:
void luDecomp();
void luSolve();
void newtonRaphson();
void newtonRaphson(Eigen::Vector3d &x);
// Solves J * p = -f using a fast analytical solver with singularity checks.
// Can be easily swapped to LU with partial pivoting if needed.
Eigen::Vector3d solveNewtonStep(const Eigen::Matrix3d &fjac,
const Eigen::Vector3d &fvec);
// Find driver parameters t0, delta_t, Ceff.
void findDriverParams(double ceff);
std::pair<double, double> gateCapDelaySlew(double ceff);
@ -84,9 +90,9 @@ protected:
double cl);
double y0dcl(double t,
double cl);
void showX();
void showFvec();
void showJacobian();
void showX(const Eigen::Vector3d &x);
void showFvec(const Eigen::Vector3d &fvec);
void showJacobian(const Eigen::Matrix3d &fjac);
std::pair<double, double> findDriverDelaySlew();
double findVoCrossing(double vth,
double t_lower,
@ -148,12 +154,7 @@ protected:
static constexpr int max_nr_order_ = 3;
std::array<double, max_nr_order_> x_;
std::array<double, max_nr_order_> fvec_;
std::array<std::array<double, max_nr_order_>, max_nr_order_> fjac_;
std::array<double, max_nr_order_> scale_;
std::array<double, max_nr_order_> p_;
std::array<int, max_nr_order_> index_;
// Driver slew used to check load delay.
double drvr_slew_;
@ -194,7 +195,10 @@ public:
std::pair<double, double> gateDelaySlew() override;
std::pair<double, double> loadDelaySlew(const Pin *,
double elmore) override;
void evalDmpEqns() override;
void
evalDmpEqns(Eigen::Vector3d &x,
Eigen::Vector3d &fvec,
Eigen::Matrix3d &fjac) override;
protected:
double voCrossingUpperBound() override;
@ -219,7 +223,10 @@ public:
double rpi,
double c1) override;
std::pair<double, double> gateDelaySlew() override;
void evalDmpEqns() override;
void
evalDmpEqns(Eigen::Vector3d &x,
Eigen::Vector3d &fvec,
Eigen::Matrix3d &fjac) override;
protected:
double voCrossingUpperBound() override;
@ -255,7 +262,10 @@ class DmpOnePole : public DmpAlg
{
public:
DmpOnePole(StaState *sta);
void evalDmpEqns() override;
void
evalDmpEqns(Eigen::Vector3d &x,
Eigen::Vector3d &fvec,
Eigen::Matrix3d &fjac) override;
protected:
double voCrossingUpperBound() override;