From 832c62b57952637b5f6959ca9a8c513524481723 Mon Sep 17 00:00:00 2001 From: Drew Lewis Date: Thu, 18 Jun 2026 16:52:14 -0400 Subject: [PATCH] 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% | --- dcalc/DmpCeff.cc | 350 ++++++++++++++++++++++------------------------- dcalc/DmpCeff.hh | 46 ++++--- 2 files changed, 194 insertions(+), 202 deletions(-) diff --git a/dcalc/DmpCeff.cc b/dcalc/DmpCeff.cc index 5d0aef1d..6a37f348 100644 --- a/dcalc/DmpCeff.cc +++ b/dcalc/DmpCeff.cc @@ -31,6 +31,7 @@ // slew voltage is matched instead of y20 in eqn 12. #include "DmpCeff.hh" +#include #include #include @@ -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; } //////////////////////////////////////////////////////////////// diff --git a/dcalc/DmpCeff.hh b/dcalc/DmpCeff.hh index a2e8a144..6db5f338 100644 --- a/dcalc/DmpCeff.hh +++ b/dcalc/DmpCeff.hh @@ -27,6 +27,9 @@ #include #include +#include +#include + #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 Vo(double t); // Load response to driver waveform (vl, dvl/dt). std::pair 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 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 findDriverDelaySlew(); double findVoCrossing(double vth, double t_lower, @@ -148,12 +154,7 @@ protected: static constexpr int max_nr_order_ = 3; - std::array x_; - std::array fvec_; - std::array, max_nr_order_> fjac_; - std::array scale_; - std::array p_; - std::array index_; + // Driver slew used to check load delay. double drvr_slew_; @@ -194,7 +195,10 @@ public: std::pair gateDelaySlew() override; std::pair 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 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;