From 6573037934e93af6d103fd532049b28af0d01c8f Mon Sep 17 00:00:00 2001 From: couriersud Date: Sat, 2 Nov 2019 12:25:50 +0100 Subject: netlist: Support for float, double and long double solvers. [Couriersud] - Added new solver parameter FPTYPE. This determines in which floating point domain the linear system is solved. May be one of "FLOAT", "DOUBLE" or "LONGDOUBLE" - Added option "--fperr" to nltool. This enables floating point exceptions. This helps debugging the code under gdb. The purpose of this going forward is to have more choice in optimization. Non-dynamic systems should be just fine in the float domain. Dynamic systems (i.e. diodes, bjts, mosfets) should in general work with double. Certain edge cases may require long double resolution. --- src/lib/netlist/analog/nld_bjt.cpp | 8 +- src/lib/netlist/analog/nld_generic_models.h | 4 +- src/lib/netlist/analog/nld_mosfet.cpp | 4 +- src/lib/netlist/analog/nld_switches.cpp | 2 +- src/lib/netlist/analog/nlid_twoterm.cpp | 2 +- src/lib/netlist/analog/nlid_twoterm.h | 4 +- src/lib/netlist/devices/nld_74123.cpp | 2 +- src/lib/netlist/devices/nld_74ls629.cpp | 2 +- src/lib/netlist/devices/nld_log.cpp | 6 +- src/lib/netlist/devices/nld_mm5837.cpp | 4 +- src/lib/netlist/devices/nld_r2r_dac.cpp | 2 +- src/lib/netlist/devices/nld_system.cpp | 2 +- src/lib/netlist/devices/nlid_system.h | 22 +-- src/lib/netlist/nl_config.h | 26 +++- src/lib/netlist/plib/gmres.h | 25 ++-- src/lib/netlist/plib/mat_cr.h | 8 +- src/lib/netlist/plib/pfmtlog.h | 6 + src/lib/netlist/plib/ptime.h | 28 +++- src/lib/netlist/plib/putil.h | 13 ++ src/lib/netlist/prg/nltool.cpp | 15 +- src/lib/netlist/solver/nld_matrix_solver.cpp | 10 +- src/lib/netlist/solver/nld_matrix_solver.h | 55 ++++--- src/lib/netlist/solver/nld_ms_direct.h | 12 +- src/lib/netlist/solver/nld_ms_direct1.h | 4 +- src/lib/netlist/solver/nld_ms_direct2.h | 4 +- src/lib/netlist/solver/nld_ms_gcr.h | 6 +- src/lib/netlist/solver/nld_ms_gmres.h | 10 +- src/lib/netlist/solver/nld_ms_sm.h | 30 ++-- src/lib/netlist/solver/nld_ms_sor.h | 36 ++--- src/lib/netlist/solver/nld_ms_sor_mat.h | 16 +- src/lib/netlist/solver/nld_ms_w.h | 43 +++--- src/lib/netlist/solver/nld_solver.cpp | 212 +++++++++++++++------------ src/lib/netlist/solver/nld_solver.h | 4 + 33 files changed, 364 insertions(+), 263 deletions(-) (limited to 'src/lib') diff --git a/src/lib/netlist/analog/nld_bjt.cpp b/src/lib/netlist/analog/nld_bjt.cpp index 0402f94e8b9..f1d66c2344f 100644 --- a/src/lib/netlist/analog/nld_bjt.cpp +++ b/src/lib/netlist/analog/nld_bjt.cpp @@ -21,18 +21,18 @@ namespace analog class diode { public: - diode() : m_Is(1e-15), m_VT(0.0258), m_VT_inv(1.0 / m_VT) {} + diode() : m_Is(1e-15), m_VT(0.0258), m_VT_inv(plib::reciprocal(m_VT)) {} diode(const nl_fptype Is, const nl_fptype n) { m_Is = Is; m_VT = 0.0258 * n; - m_VT_inv = 1.0 / m_VT; + m_VT_inv = plib::reciprocal(m_VT); } void set(const nl_fptype Is, const nl_fptype n) { m_Is = Is; m_VT = 0.0258 * n; - m_VT_inv = 1.0 / m_VT; + m_VT_inv = plib::reciprocal(m_VT); } nl_fptype I(const nl_fptype V) const { return m_Is * std::exp(V * m_VT_inv) - m_Is; } nl_fptype g(const nl_fptype V) const { return m_Is * m_VT_inv * std::exp(V * m_VT_inv); } @@ -334,7 +334,7 @@ namespace analog /* Base current is 0.005 / beta * as a rough estimate, we just scale the conductance down */ - m_gB = 1.0 / (m_V/(0.005 / BF)); + m_gB = plib::reciprocal((m_V/(0.005 / BF))); //m_gB = d.gI(0.005 / alpha); diff --git a/src/lib/netlist/analog/nld_generic_models.h b/src/lib/netlist/analog/nld_generic_models.h index f33168594fb..fe56e25e674 100644 --- a/src/lib/netlist/analog/nld_generic_models.h +++ b/src/lib/netlist/analog/nld_generic_models.h @@ -71,7 +71,7 @@ namespace analog void timestep(nl_fptype cap, nl_fptype v, nl_fptype step) { - m_h = 1.0 / step; + m_h = plib::reciprocal(step); m_c = cap; m_v = v; } @@ -108,7 +108,7 @@ namespace analog void timestep(nl_fptype cap, nl_fptype v, nl_fptype step) { plib::unused_var(cap); - m_h = 1.0 / step; + m_h = plib::reciprocal(step); m_v = v; } void setparams(nl_fptype gmin) { m_gmin = gmin; } diff --git a/src/lib/netlist/analog/nld_mosfet.cpp b/src/lib/netlist/analog/nld_mosfet.cpp index 3bd4d4fa2c4..e78de46252f 100644 --- a/src/lib/netlist/analog/nld_mosfet.cpp +++ b/src/lib/netlist/analog/nld_mosfet.cpp @@ -449,9 +449,9 @@ namespace analog const nl_fptype k = 3.5; // see "Circuit Simulation", page 185 nl_fptype d = (Vgs - m_Vgs); - Vgs = m_Vgs + 1.0/k * (d < 0 ? -1.0 : 1.0) * std::log1p(k * std::abs(d)); + Vgs = m_Vgs + plib::reciprocal(k) * (d < 0 ? -1.0 : 1.0) * std::log1p(k * std::abs(d)); d = (Vgd - m_Vgd); - Vgd = m_Vgd + 1.0/k * (d < 0 ? -1.0 : 1.0) * std::log1p(k * std::abs(d)); + Vgd = m_Vgd + plib::reciprocal(k) * (d < 0 ? -1.0 : 1.0) * std::log1p(k * std::abs(d)); m_Vgs = Vgs; m_Vgd = Vgd; diff --git a/src/lib/netlist/analog/nld_switches.cpp b/src/lib/netlist/analog/nld_switches.cpp index 0e62fe13952..9b73eb9650f 100644 --- a/src/lib/netlist/analog/nld_switches.cpp +++ b/src/lib/netlist/analog/nld_switches.cpp @@ -12,7 +12,7 @@ /* FIXME : convert to parameters */ -#define R_OFF (1.0 / exec().gmin()) +#define R_OFF (plib::reciprocal(exec().gmin())) #define R_ON 0.01 namespace netlist diff --git a/src/lib/netlist/analog/nlid_twoterm.cpp b/src/lib/netlist/analog/nlid_twoterm.cpp index bc63e625c1a..fe928782810 100644 --- a/src/lib/netlist/analog/nlid_twoterm.cpp +++ b/src/lib/netlist/analog/nlid_twoterm.cpp @@ -53,7 +53,7 @@ NETLIB_UPDATE(twoterm) NETLIB_RESET(R_base) { NETLIB_NAME(twoterm)::reset(); - set_R(1.0 / exec().gmin()); + set_R(plib::reciprocal(exec().gmin())); } // ---------------------------------------------------------------------------------------- diff --git a/src/lib/netlist/analog/nlid_twoterm.h b/src/lib/netlist/analog/nlid_twoterm.h index 2618c3eb20e..cb80e15a34a 100644 --- a/src/lib/netlist/analog/nlid_twoterm.h +++ b/src/lib/netlist/analog/nlid_twoterm.h @@ -411,7 +411,7 @@ namespace analog { m_t += step; m_funcparam[0] = m_t; - this->set_G_V_I(1.0 / m_R(), + this->set_G_V_I(plib::reciprocal(m_R()), m_compiled.evaluate(m_funcparam), 0.0); } @@ -422,7 +422,7 @@ namespace analog NETLIB_RESETI() { NETLIB_NAME(twoterm)::reset(); - this->set_G_V_I(1.0 / m_R(), m_V(), 0.0); + this->set_G_V_I(plib::reciprocal(m_R()), m_V(), 0.0); } private: diff --git a/src/lib/netlist/devices/nld_74123.cpp b/src/lib/netlist/devices/nld_74123.cpp index 213fa9cc396..067e6ca1918 100644 --- a/src/lib/netlist/devices/nld_74123.cpp +++ b/src/lib/netlist/devices/nld_74123.cpp @@ -255,7 +255,7 @@ namespace netlist NETLIB_RESET(74123) { - m_KP = 1.0 / (1.0 + exp(m_K())); + m_KP = plib::reciprocal(1.0 + exp(m_K())); m_RP.reset(); m_RN.reset(); diff --git a/src/lib/netlist/devices/nld_74ls629.cpp b/src/lib/netlist/devices/nld_74ls629.cpp index 0f1c64fc74a..42c9292315c 100644 --- a/src/lib/netlist/devices/nld_74ls629.cpp +++ b/src/lib/netlist/devices/nld_74ls629.cpp @@ -212,7 +212,7 @@ namespace netlist // FIXME: we need a possibility to remove entries from queue ... // or an exact model ... - m_clock.m_inc = netlist_time::from_double(0.5 / freq); + m_clock.m_inc = netlist_time::from_fp(0.5 / freq); //m_clock.update(); //NL_VERBOSE_OUT(("{1} {2} {3} {4}\n", name(), v_freq, v_rng, freq)); diff --git a/src/lib/netlist/devices/nld_log.cpp b/src/lib/netlist/devices/nld_log.cpp index 3b173060986..03165a573ef 100644 --- a/src/lib/netlist/devices/nld_log.cpp +++ b/src/lib/netlist/devices/nld_log.cpp @@ -31,7 +31,7 @@ namespace netlist NETLIB_UPDATEI() { /* use pstring::sprintf, it is a LOT faster */ - m_writer.writeline(plib::pfmt("{1:.9} {2}").e(exec().time().as_double()).e(static_cast(m_I()))); + m_writer.writeline(plib::pfmt("{1:.9} {2}").e(exec().time().as_fp()).e(static_cast(m_I()))); } NETLIB_RESETI() { } @@ -50,7 +50,7 @@ namespace netlist NETLIB_UPDATEI() { - m_writer.writeline(plib::pfmt("{1:.9} {2}").e(exec().time().as_double()).e(static_cast(m_I() - m_I2()))); + m_writer.writeline(plib::pfmt("{1:.9} {2}").e(exec().time().as_fp()).e(static_cast(m_I() - m_I2()))); } NETLIB_RESETI() { } @@ -83,7 +83,7 @@ namespace netlist NETLIB_UPDATE(wav) { - fprintf(m_file, "%e %e\n", netlist().time().as_double(), m_I()); + fprintf(m_file, "%e %e\n", netlist().time().as_fp(), m_I()); } NETLIB_NAME(log)::~NETLIB_NAME(wav)() diff --git a/src/lib/netlist/devices/nld_mm5837.cpp b/src/lib/netlist/devices/nld_mm5837.cpp index 2c58354fd3f..5d32c8ff8cd 100644 --- a/src/lib/netlist/devices/nld_mm5837.cpp +++ b/src/lib/netlist/devices/nld_mm5837.cpp @@ -70,7 +70,7 @@ namespace netlist //m_V0.initial(0.0); //m_RV.do_reset(); m_RV.set_G_V_I(plib::constants::one() / R_LOW, 0.0, 0.0); - m_inc = netlist_time::from_double(1.0 / m_FREQ()); + m_inc = netlist_time::from_fp(plib::reciprocal(m_FREQ())); if (m_FREQ() < 24000 || m_FREQ() > 56000) log().warning(MW_FREQUENCY_OUTSIDE_OF_SPECS_1(m_FREQ())); @@ -80,7 +80,7 @@ namespace netlist NETLIB_UPDATE_PARAM(MM5837_dip) { - m_inc = netlist_time::from_double(1.0 / m_FREQ()); + m_inc = netlist_time::from_fp(plib::reciprocal(m_FREQ())); if (m_FREQ() < 24000 || m_FREQ() > 56000) log().warning(MW_FREQUENCY_OUTSIDE_OF_SPECS_1(m_FREQ())); } diff --git a/src/lib/netlist/devices/nld_r2r_dac.cpp b/src/lib/netlist/devices/nld_r2r_dac.cpp index e2695766332..c1cc4462b7b 100644 --- a/src/lib/netlist/devices/nld_r2r_dac.cpp +++ b/src/lib/netlist/devices/nld_r2r_dac.cpp @@ -44,7 +44,7 @@ namespace netlist nl_fptype V = m_VIN() / static_cast(1 << m_num()) * static_cast(m_val()); - this->set_G_V_I(1.0 / m_R(), V, 0.0); + this->set_G_V_I(plib::reciprocal(m_R()), V, 0.0); } } //namespace analog diff --git a/src/lib/netlist/devices/nld_system.cpp b/src/lib/netlist/devices/nld_system.cpp index 427c2f3dcae..54ef14fa1d0 100644 --- a/src/lib/netlist/devices/nld_system.cpp +++ b/src/lib/netlist/devices/nld_system.cpp @@ -25,7 +25,7 @@ namespace netlist NETLIB_RESET(extclock) { m_cnt = 0; - m_off = netlist_time::from_double(m_offset()); + m_off = netlist_time::from_fp(m_offset()); m_feedback.m_delegate = NETLIB_DELEGATE(extclock, update); //m_feedback.m_delegate .set(&NETLIB_NAME(extclock)::update, this); diff --git a/src/lib/netlist/devices/nlid_system.h b/src/lib/netlist/devices/nlid_system.h index c51f5c36c56..d57de8f4750 100644 --- a/src/lib/netlist/devices/nlid_system.h +++ b/src/lib/netlist/devices/nlid_system.h @@ -54,7 +54,7 @@ namespace devices , m_Q(*this, "Q") , m_freq(*this, "FREQ", 7159000.0 * 5) { - m_inc = netlist_time::from_double(1.0 / (m_freq()*2.0)); + m_inc = netlist_time::from_fp(1.0 / (m_freq()*2.0)); } NETLIB_RESETI() @@ -64,7 +64,7 @@ namespace devices NETLIB_UPDATE_PARAMI() { - m_inc = netlist_time::from_double(1.0 / (m_freq()*2.0)); + m_inc = netlist_time::from_fp(1.0 / (m_freq()*2.0)); } NETLIB_UPDATEI() @@ -93,7 +93,7 @@ namespace devices , m_Q(*this, "Q") , m_freq(*this, "FREQ", 7159000.0 * 5.0) { - m_inc = netlist_time::from_double(1.0 / (m_freq()*2.0)); + m_inc = netlist_time::from_fp(1.0 / (m_freq()*2.0)); connect(m_feedback, m_Q); } @@ -101,7 +101,7 @@ namespace devices NETLIB_UPDATE_PARAMI() { - m_inc = netlist_time::from_double(1.0 / (m_freq() * 2.0)); + m_inc = netlist_time::from_fp(1.0 / (m_freq() * 2.0)); } NETLIB_UPDATEI() @@ -139,8 +139,8 @@ namespace devices NETLIB_UPDATEI() { - m_funcparam[0] = exec().time().as_double(); - const netlist_time m_inc = netlist_time::from_double(m_compiled.evaluate(m_funcparam)); + m_funcparam[0] = exec().time().as_fp(); + const netlist_time m_inc = netlist_time::from_fp(m_compiled.evaluate(m_funcparam)); m_Q.push(!m_feedback(), m_inc); } @@ -168,13 +168,13 @@ namespace devices , m_cnt(*this, "m_cnt", 0) , m_off(*this, "m_off", netlist_time::zero()) { - m_inc[0] = netlist_time::from_double(1.0 / (m_freq() * 2.0)); + m_inc[0] = netlist_time::from_fp(1.0 / (m_freq() * 2.0)); connect(m_feedback, m_Q); - netlist_time base = netlist_time::from_double(1.0 / (m_freq()*2.0)); + netlist_time base = netlist_time::from_fp(1.0 / (m_freq()*2.0)); std::vector pat(plib::psplit(m_pattern(),",")); - m_off = netlist_time::from_double(m_offset()); + m_off = netlist_time::from_fp(m_offset()); std::array pati = { 0 }; @@ -331,8 +331,8 @@ namespace devices NETLIB_RESETI() { - m_RIN.set_G_V_I(1.0 / m_p_RIN(),0,0); - m_ROUT.set_G_V_I(1.0 / m_p_ROUT(),0,0); + m_RIN.set_G_V_I(plib::reciprocal(m_p_RIN()),0,0); + m_ROUT.set_G_V_I(plib::reciprocal(m_p_ROUT()),0,0); } NETLIB_UPDATEI() diff --git a/src/lib/netlist/nl_config.h b/src/lib/netlist/nl_config.h index 80f74cb1582..1633efc523b 100644 --- a/src/lib/netlist/nl_config.h +++ b/src/lib/netlist/nl_config.h @@ -105,10 +105,6 @@ static constexpr const auto NETLIST_CLOCK = NETLIST_INTERNAL_RES; //============================================================ using nl_fptype = double; -//using nl_fptype = float; - -using nl_mat_fptype = nl_fptype; -//using nl_mat_fptype = float; namespace netlist { @@ -120,6 +116,20 @@ namespace netlist struct fp_constants { }; + /*! Specific constants for long double floating point type + */ + template <> + struct fp_constants + { + static inline constexpr long double DIODE_MAXDIFF() noexcept { return 1e100L; } + static inline constexpr long double DIODE_MAXVOLT() noexcept { return 300.0L; } + + static inline constexpr long double TIMESTEP_MAXDIFF() noexcept { return 1e100L; } + static inline constexpr long double TIMESTEP_MINDIV() noexcept { return 1e-60L; } + + static inline constexpr const char * name() noexcept { return "long double"; } + }; + /*! Specific constants for double floating point type */ template <> @@ -130,6 +140,8 @@ namespace netlist static inline constexpr double TIMESTEP_MAXDIFF() noexcept { return 1e100; } static inline constexpr double TIMESTEP_MINDIV() noexcept { return 1e-60; } + + static inline constexpr const char * name() noexcept { return "double"; } }; /*! Specific constants for float floating point type @@ -137,11 +149,13 @@ namespace netlist template <> struct fp_constants { - static inline constexpr float DIODE_MAXDIFF() noexcept { return 1e5; } - static inline constexpr float DIODE_MAXVOLT() noexcept { return 30.0; } + static inline constexpr float DIODE_MAXDIFF() noexcept { return 1e2; } + static inline constexpr float DIODE_MAXVOLT() noexcept { return 20.0; } static inline constexpr float TIMESTEP_MAXDIFF() noexcept { return 1e30f; } static inline constexpr float TIMESTEP_MINDIV() noexcept { return 1e-8f; } + + static inline constexpr const char * name() noexcept { return "float"; } }; } // namespace netlist diff --git a/src/lib/netlist/plib/gmres.h b/src/lib/netlist/plib/gmres.h index 8c5c724c581..840da0cfb45 100644 --- a/src/lib/netlist/plib/gmres.h +++ b/src/lib/netlist/plib/gmres.h @@ -126,7 +126,7 @@ namespace plib // 136% for (std::size_t j = m_mat.row_idx[i]; j< m_mat.row_idx[i+1]; j++) v += m_mat.A[j] * m_mat.A[j]; - m_diag[i] = 1.0 / std::sqrt(v); + m_diag[i] = reciprocal(std::sqrt(v)); #elif 0 // works halfway, i.e. Mame perforamnce 50% // 147% - lowest average solution time with 7.094 @@ -151,14 +151,14 @@ namespace plib // 151% for (std::size_t j = m_mat.row_idx[i]; j< m_mat.row_idx[i+1]; j++) v += std::abs(m_mat.A[j]); - m_diag[i] = 1.0 / v; + m_diag[i] = reciprocal(v); #else // 124% for (std::size_t j = m_mat.row_idx[i]; j< m_mat.row_idx[i+1]; j++) v = std::max(v, std::abs(m_mat.A[j])); - m_diag[i] = 1.0 / v; + m_diag[i] = reciprocal(v); #endif - //m_diag[i] = 1.0 / m_mat.A[m_mat.diag[i]]; + //m_diag[i] = reciprocal(m_mat.A[m_mat.diag[i]]); } } @@ -261,18 +261,19 @@ namespace plib } m_ht[kp1][k] = std::sqrt(vec_mult2(n, m_v[kp1])); - if (m_ht[kp1][k] != 0.0) - vec_scale(n, m_v[kp1], constants::one() / m_ht[kp1][k]); + // FIXME: comparison to zero + if (m_ht[kp1][k] != plib::constants::zero()) + vec_scale(n, m_v[kp1], reciprocal(m_ht[kp1][k])); for (std::size_t j = 0; j < k; j++) givens_mult(m_c[j], m_s[j], m_ht[j][k], m_ht[j+1][k]); - const float_type mu = 1.0 / std::hypot(m_ht[k][k], m_ht[kp1][k]); + const float_type mu = reciprocal(std::hypot(m_ht[k][k], m_ht[kp1][k])); m_c[k] = m_ht[k][k] * mu; m_s[k] = -m_ht[kp1][k] * mu; m_ht[k][k] = m_c[k] * m_ht[k][k] - m_s[k] * m_ht[kp1][k]; - m_ht[kp1][k] = 0.0; + m_ht[kp1][k] = plib::constants::zero(); givens_mult(m_c[k], m_s[k], m_g[k], m_g[kp1]); @@ -287,7 +288,7 @@ namespace plib /* x += m_v[j] * m_y[j] */ for (std::size_t i = k + 1; i-- > 0;) { - double tmp = m_g[i]; + auto tmp(m_g[i]); for (std::size_t j = i + 1; j <= k; j++) tmp -= m_ht[i][j] * m_y[j]; m_y[i] = tmp / m_ht[i][i]; @@ -334,7 +335,7 @@ namespace plib *------------------------------------------------------------------------*/ std::size_t itr_used = 0; - double rho_delta = 0.0; + auto rho_delta(plib::constants::zero()); const std::size_t n = size(); @@ -400,7 +401,7 @@ namespace plib vec_set_scalar(RESTART+1, m_g, +constants::zero()); m_g[0] = rho; - vec_mult_scalar(n, m_v[0], residual, constants::one() / rho); + vec_mult_scalar(n, m_v[0], residual, reciprocal(rho)); if (do_k(ops, x, itr_used, rho_delta, true)) // converged @@ -504,7 +505,7 @@ namespace plib else { beta = alpha * ( c / 2.0)*( c / 2.0); - alpha = 1.0 / (d - beta); + alpha = reciprocal(d - beta); for (std::size_t k = 0; k < size(); k++) p[k] = residual[k] + beta * p[k]; } diff --git a/src/lib/netlist/plib/mat_cr.h b/src/lib/netlist/plib/mat_cr.h index 3a6c59f24ab..132cab5d96f 100644 --- a/src/lib/netlist/plib/mat_cr.h +++ b/src/lib/netlist/plib/mat_cr.h @@ -167,7 +167,7 @@ namespace plib const std::size_t oe = nz_num; while (k < oe) { - T tmp = 0.0; + T tmp = plib::constants::zero(); const std::size_t e = row_idx[row+1]; for (; k < e; k++) tmp += A[k] * x[col_idx[k]]; @@ -294,7 +294,7 @@ namespace plib { std::size_t nzbdp = 0; std::size_t pi = base::diag[i]; - const typename base::value_type f = 1.0 / base::A[pi++]; + auto f = reciprocal(base::A[pi++]); const std::size_t piie = base::row_idx[i+1]; const auto &nz = base::nzbd[i]; @@ -345,7 +345,7 @@ namespace plib { std::size_t nzbdp = 0; std::size_t pi = base::diag[i]; - const typename base::value_type f = 1.0 / base::A[pi++]; + const auto f = reciprocal(base::A[pi++]); const std::size_t piie = base::row_idx[i+1]; const auto &nz = base::nzbd[i]; @@ -358,7 +358,7 @@ namespace plib while (base::col_idx[pj] < i) pj++; - const typename base::value_type f1 = - base::A[pj++] * f; + auto f1 = - base::A[pj++] * f; // subtract row i from j // fill-in available assumed, i.e. matrix was prepared diff --git a/src/lib/netlist/plib/pfmtlog.h b/src/lib/netlist/plib/pfmtlog.h index 13ec2396a8a..17ca3f1d7ad 100644 --- a/src/lib/netlist/plib/pfmtlog.h +++ b/src/lib/netlist/plib/pfmtlog.h @@ -127,6 +127,12 @@ struct ptype_traits : ptype_traits_base static char32_t fmt_spec() { return 'f'; } }; +template<> +struct ptype_traits : ptype_traits_base +{ + static char32_t fmt_spec() { return 'f'; } +}; + template<> struct ptype_traits : ptype_traits_base { diff --git a/src/lib/netlist/plib/ptime.h b/src/lib/netlist/plib/ptime.h index 4044f518e24..c0a6cd5a34b 100644 --- a/src/lib/netlist/plib/ptime.h +++ b/src/lib/netlist/plib/ptime.h @@ -99,11 +99,18 @@ namespace plib } constexpr internal_type as_raw() const noexcept { return m_time; } - constexpr double as_double() const noexcept + + template + constexpr typename std::enable_if::value, FT>::type + as_fp() const noexcept { - return static_cast(m_time) * inv_res; + return static_cast(m_time) * inv_res(); } + constexpr double as_double() const noexcept { return as_fp(); } + constexpr double as_float() const noexcept { return as_fp(); } + constexpr double as_long_double() const noexcept { return as_fp(); } + // for save states .... C14CONSTEXPR internal_type *get_internaltype_ptr() noexcept { return &m_time; } @@ -113,7 +120,19 @@ namespace plib static constexpr const ptime from_sec(const internal_type s) noexcept { return ptime(s, UINT64_C( 1)); } static constexpr const ptime from_hz(const internal_type hz) noexcept { return ptime(1 , hz); } static constexpr const ptime from_raw(const internal_type raw) noexcept { return ptime(raw); } - static constexpr const ptime from_double(const double t) noexcept { return ptime(static_cast(std::floor(t * static_cast(RES) + 0.5)), RES); } + + template + static constexpr const typename std::enable_if::value, ptime>::type + from_fp(const FT t) noexcept { return ptime(static_cast(std::floor(t * static_cast(RES) + 0.5)), RES); } + + static constexpr const ptime from_double(const double t) noexcept + { return from_fp(t); } + + static constexpr const ptime from_float(const float t) noexcept + { return from_fp(t); } + + static constexpr const ptime from_long_double(const long double t) noexcept + { return from_fp(t); } static constexpr const ptime zero() noexcept { return ptime(0, RES); } static constexpr const ptime quantum() noexcept { return ptime(1, RES); } @@ -126,7 +145,8 @@ namespace plib constexpr internal_type in_sec() const noexcept { return m_time / (RES / UINT64_C( 1)); } private: - static constexpr const double inv_res = 1.0 / static_cast(RES); + template + static constexpr FT inv_res() noexcept { return static_cast(1.0) / static_cast(RES); } internal_type m_time; }; diff --git a/src/lib/netlist/plib/putil.h b/src/lib/netlist/plib/putil.h index 9a21a223ba9..ae358b960df 100644 --- a/src/lib/netlist/plib/putil.h +++ b/src/lib/netlist/plib/putil.h @@ -259,6 +259,19 @@ namespace plib static inline constexpr const T cast(V &&v) noexcept { return static_cast(v); } }; + /*! typesafe reciprocal function + * + * @tparam T type of the argument + * @param v argument + * @return reciprocal of argument + */ + template + static inline constexpr typename std::enable_if::value, T>::type + reciprocal(T v) noexcept + { + return constants::one() / v; + } + static_assert(noexcept(constants::one()) == true, "Not evaluated as constexpr"); diff --git a/src/lib/netlist/prg/nltool.cpp b/src/lib/netlist/prg/nltool.cpp index 20a37d24ff1..9ca4579ae4d 100644 --- a/src/lib/netlist/prg/nltool.cpp +++ b/src/lib/netlist/prg/nltool.cpp @@ -53,6 +53,8 @@ public: opt_inp(*this, "i", "input", "", "input file to process (default is none)"), opt_loadstate(*this,"", "loadstate", "", "load state from file and continue from there"), opt_savestate(*this,"", "savestate", "", "save state to file at end of run"), + opt_fperr(*this, "", "fperr", + "raise exception on floating point errors. This is intended to be used during debugging."), opt_grp5(*this, "Options for convert command", "These options are only used by the convert command."), opt_type(*this, "y", "type", 0, std::vector({"spice","eagle","rinf"}), "type of file to be converted: spice,eagle,rinf"), @@ -97,6 +99,7 @@ public: plib::option_str opt_inp; plib::option_str opt_loadstate; plib::option_str opt_savestate; + plib::option_bool opt_fperr; plib::option_group opt_grp5; plib::option_str_limit opt_type; plib::option_group opt_grp6; @@ -319,7 +322,7 @@ struct input_t if (e != 3) throw netlist::nl_exception(plib::pfmt("error {1} scanning line {2}\n")(e)(line)); m_value = static_cast(val); - m_time = netlist::netlist_time::from_double(t); + m_time = netlist::netlist_time::from_fp(t); m_param = setup.find_param(pstring(buf.data()), true); } @@ -395,7 +398,7 @@ void tool_app_t::run() nt.reset(); inps = read_input(nt.setup(), opt_inp()); - ttr = netlist::netlist_time::from_double(opt_ttr()); + ttr = netlist::netlist_time::from_fp(opt_ttr()); } @@ -831,10 +834,6 @@ int tool_app_t::execute() { tool_app_t opts; - /* make SIGFPE actually deliver signals on supoorted platforms */ - plib::fpsignalenabler::global_enable(true); - plib::fpsignalenabler sigen(plib::FP_ALL & ~plib::FP_INEXACT & ~plib::FP_UNDERFLOW); - if (opt_help()) { pout(usage()); @@ -869,6 +868,9 @@ int tool_app_t::execute() try { + plib::fpsignalenabler::global_enable(opt_fperr()); + plib::fpsignalenabler fpprotect(plib::FP_DIVBYZERO | plib::FP_UNDERFLOW | plib::FP_OVERFLOW | plib::FP_INVALID); + pstring cmd = opt_cmd.as_string(); if (cmd == "listdevices") listdevices(); @@ -902,6 +904,7 @@ int tool_app_t::execute() perr("plib exception caught: {}\n", e.text()); return 2; } + #if 0 std::cout.imbue(std::locale("de_DE.utf8")); std::cout.imbue(std::locale("C.UTF-8")); diff --git a/src/lib/netlist/solver/nld_matrix_solver.cpp b/src/lib/netlist/solver/nld_matrix_solver.cpp index afcdec647cd..a1bb4f2f88e 100644 --- a/src/lib/netlist/solver/nld_matrix_solver.cpp +++ b/src/lib/netlist/solver/nld_matrix_solver.cpp @@ -397,13 +397,13 @@ namespace solver if (m_params.m_dynamic_ts && has_timestep_devices()) { - m_Q_sync.net().toggle_and_push_to_queue(netlist_time::from_double(m_params.m_min_timestep)); + m_Q_sync.net().toggle_and_push_to_queue(netlist_time::from_fp(m_params.m_min_timestep)); } } void matrix_solver_t::step(const netlist_time &delta) { - const nl_fptype dd = delta.as_double(); + const nl_fptype dd = delta.as_fp(); for (auto &d : m_step_devices) d->timestep(dd); } @@ -439,7 +439,7 @@ namespace solver if (this_resched > 1 && !m_Q_sync.net().is_queued()) { log().warning(MW_NEWTON_LOOPS_EXCEEDED_ON_NET_1(this->name())); - m_Q_sync.net().toggle_and_push_to_queue(netlist_time::from_double(m_params.m_nr_recalc_delay)); + m_Q_sync.net().toggle_and_push_to_queue(netlist_time::from_fp(m_params.m_nr_recalc_delay())); } } else @@ -447,7 +447,7 @@ namespace solver this->vsolve_non_dynamic(false); } - const netlist_time next_time_step = compute_next_timestep(delta.as_double()); + const netlist_time next_time_step = compute_next_timestep(delta.as_fp()); return next_time_step; } @@ -558,7 +558,7 @@ namespace solver static_cast(this->m_stat_newton_raphson) / static_cast(this->m_stat_vsolver_calls)); log().verbose(" {1:10} invocations ({2:6.0} Hz) {3:10} gs fails ({4:6.2} %) {5:6.3} average", this->m_stat_calculations, - static_cast(this->m_stat_calculations) / this->exec().time().as_double(), + static_cast(this->m_stat_calculations) / this->exec().time().as_fp(), this->m_iterative_fail, 100.0 * static_cast(this->m_iterative_fail) / static_cast(this->m_stat_calculations), diff --git a/src/lib/netlist/solver/nld_matrix_solver.h b/src/lib/netlist/solver/nld_matrix_solver.h index 632d038594f..341374b1df0 100644 --- a/src/lib/netlist/solver/nld_matrix_solver.h +++ b/src/lib/netlist/solver/nld_matrix_solver.h @@ -40,6 +40,12 @@ namespace solver GMRES ) + P_ENUM(matrix_fp_type_e, + FLOAT, + DOUBLE, + LONGDOUBLE + ) + struct solver_parameters_t { solver_parameters_t(device_t &parent) @@ -48,6 +54,7 @@ namespace solver /* iteration parameters */ , m_gs_sor(parent, "SOR_FACTOR", 1.059) , m_method(parent, "METHOD", matrix_type_e::MAT_CR) + , m_fp_type(parent, "FPTYPE", matrix_fp_type_e::DOUBLE) , m_accuracy(parent, "ACCURACY", 1e-7) , m_gs_loops(parent, "GS_LOOPS", 9) // Gauss-Seidel loops @@ -55,7 +62,7 @@ namespace solver , m_gmin(parent, "GMIN", 1e-9) , m_pivot(parent, "PIVOT", false) // use pivoting - on supported solvers , m_nr_loops(parent, "NR_LOOPS", 250) // Newton-Raphson loops - , m_nr_recalc_delay(parent, "NR_RECALC_DELAY", netlist_time::quantum().as_double()) // Delay to next solve attempt if nr loops exceeded + , m_nr_recalc_delay(parent, "NR_RECALC_DELAY", netlist_time::quantum().as_fp()) // Delay to next solve attempt if nr loops exceeded , m_parallel(parent, "PARALLEL", 0) /* automatic time step */ @@ -72,7 +79,7 @@ namespace solver { m_min_timestep = m_dynamic_min_ts(); - m_max_timestep = netlist_time::from_double(1.0 / m_freq()).as_double(); + m_max_timestep = netlist_time::from_fp(plib::reciprocal(m_freq())).as_fp(); if (m_dynamic_ts) { @@ -87,6 +94,7 @@ namespace solver param_fp_t m_freq; param_fp_t m_gs_sor; param_enum_t m_method; + param_enum_t m_fp_type; param_fp_t m_accuracy; param_num_t m_gs_loops; param_fp_t m_gmin; @@ -122,9 +130,11 @@ namespace solver terminal_t **terms() noexcept { return m_terms.data(); } - nl_fptype getV() const noexcept { return m_net->Q_Analog(); } + template ::value, void>> + FT getV() const noexcept { return static_cast(m_net->Q_Analog()); } - void setV(nl_fptype v) noexcept { m_net->set_Q_Analog(v); } + template ::value, void>> + void setV(FT v) noexcept { m_net->set_Q_Analog(static_cast(v)); } bool isNet(const analog_net_t * net) const noexcept { return net == m_net; } @@ -240,7 +250,7 @@ namespace solver using vtype = typename std::decay::type; vtype cerr = 0; for (std::size_t i = 0; i < iN; i++) - cerr = std::max(cerr, std::abs(V[i] - static_cast(this->m_terms[i].getV()))); + cerr = std::max(cerr, std::abs(V[i] - this->m_terms[i].getV())); return cerr; } @@ -383,9 +393,9 @@ namespace solver : matrix_solver_t(anetlist, name, nets, params) , m_dim(size) , m_mat_ptr(size, this->max_railstart() + 1) - , m_last_V(size, plib::constants::zero()) - , m_DD_n_m_1(size, plib::constants::zero()) - , m_h_n_m_1(size, plib::constants::zero()) + , m_last_V(size, plib::constants::zero()) + , m_DD_n_m_1(size, plib::constants::zero()) + , m_h_n_m_1(size, plib::constants::zero()) { /* * save states @@ -407,11 +417,11 @@ namespace solver plib::parray2D m_mat_ptr; /* state - variable time_stepping */ PALIGNAS_VECTOROPT() - plib::parray m_last_V; + plib::parray m_last_V; PALIGNAS_VECTOROPT() - plib::parray m_DD_n_m_1; + plib::parray m_DD_n_m_1; PALIGNAS_VECTOROPT() - plib::parray m_h_n_m_1; + plib::parray m_h_n_m_1; constexpr std::size_t size() const noexcept { return (SIZE > 0) ? static_cast(SIZE) : m_dim; } @@ -428,7 +438,7 @@ namespace solver // avoid floating point exceptions const nl_fptype DD_n = std::max(-fp_constants::TIMESTEP_MAXDIFF(), - std::min(+fp_constants::TIMESTEP_MAXDIFF(),(t.getV() - m_last_V[k]))); + std::min(+fp_constants::TIMESTEP_MAXDIFF(),(t.template getV() - m_last_V[k]))); const nl_fptype hn = cur_ts; //printf("%g %g %g %g\n", DD_n, hn, t.m_DD_n_m_1, t.m_h_n_m_1); @@ -445,7 +455,7 @@ namespace solver if (new_net_timestep < new_solver_timestep) new_solver_timestep = new_net_timestep; - m_last_V[k] = t.getV(); + m_last_V[k] = t.template getV(); } if (new_solver_timestep < m_params.m_min_timestep) { @@ -457,7 +467,7 @@ namespace solver /* * FIXME: Factor 2 below is important. Without, we get timing issues. This must be a bug elsewhere. */ - return std::max(netlist_time::from_double(new_solver_timestep), netlist_time::quantum() * 2); + return std::max(netlist_time::from_fp(new_solver_timestep), netlist_time::quantum() * 2); } @@ -488,12 +498,12 @@ namespace solver void clear_square_mat(M &m) { const std::size_t n = size(); - for (std::size_t k=0; k < n; k++) { auto *p = &(m[k][0]); + using mat_elem_type = typename std::decay::type; for (std::size_t i=0; i < n; i++) - p[i] = 0.0; + p[i] = plib::constants::zero(); } } @@ -514,11 +524,14 @@ namespace solver const auto &Idr = m_Idrn[k]; const auto &cnV = m_connected_net_Vn[k]; + // FIXME: gonn, gtn and Idr - which float types should they have? + for (std::size_t i = 0; i < railstart; i++) - *tcr_r[i] += go[i]; + *tcr_r[i] += static_cast(go[i]); - typename RT::value_type gtot_t = 0.0; - typename RT::value_type RHS_t = 0.0; + /* use native floattype for now */ + auto gtot_t(plib::constants::zero()); + auto RHS_t (plib::constants::zero()); for (std::size_t i = 0; i < term_count; i++) { @@ -536,9 +549,9 @@ namespace solver RHS_t += (/*m_Idr[i]*/ (- go[i]) * *cnV[i]); } - RHS[k] = RHS_t; + RHS[k] = static_cast(RHS_t); // update diagonal element ... - *tcr_r[railstart] += gtot_t; //mat.A[mat.diag[k]] += gtot_t; + *tcr_r[railstart] += static_cast(gtot_t); //mat.A[mat.diag[k]] += gtot_t; } } diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h index 62332dbe6b5..dc7aea3e743 100644 --- a/src/lib/netlist/solver/nld_ms_direct.h +++ b/src/lib/netlist/solver/nld_ms_direct.h @@ -72,7 +72,7 @@ namespace solver for (std::size_t i = 0; i < kN; i++) { /* FIXME: Singular matrix? */ - const FT f = 1.0 / m_A[i][i]; + const FT f = plib::reciprocal(m_A[i][i]); const auto &nzrd = this->m_terms[i].m_nzrd; const auto &nzbd = this->m_terms[i].m_nzbd; @@ -101,13 +101,13 @@ namespace solver if (maxrow != i) { /* Swap the maxrow and ith row */ - for (std::size_t k = 0; k < kN + 1; k++) { + for (std::size_t k = 0; k < kN; k++) { std::swap(m_A[i][k], m_A[maxrow][k]); } - //std::swap(RHS(i), RHS(maxrow)); + std::swap(m_RHS[i], m_RHS[maxrow]); } /* FIXME: Singular matrix? */ - const FT f = 1.0 / m_A[i][i]; + const FT f = plib::reciprocal(m_A[i][i]); /* Eliminate column i from row j */ @@ -168,9 +168,9 @@ namespace solver this->LE_solve(); this->LE_back_subst(m_new_V); - const FT err = (newton_raphson ? this->delta(m_new_V) : 0.0); + const FT err = (newton_raphson ? this->delta(m_new_V) : plib::constants::zero()); this->store(m_new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; + return (err > static_cast(this->m_params.m_accuracy)) ? 2 : 1; } template diff --git a/src/lib/netlist/solver/nld_ms_direct1.h b/src/lib/netlist/solver/nld_ms_direct1.h index 2931aaec29c..fa4d50677fb 100644 --- a/src/lib/netlist/solver/nld_ms_direct1.h +++ b/src/lib/netlist/solver/nld_ms_direct1.h @@ -39,9 +39,9 @@ namespace solver std::array new_V = { this->m_RHS[0] / this->m_A[0][0] }; - const FT err = (newton_raphson ? this->delta(new_V) : 0.0); + const FT err = (newton_raphson ? this->delta(new_V) : plib::constants::zero()); this->store(new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; + return (err > static_cast(this->m_params.m_accuracy)) ? 2 : 1; } }; diff --git a/src/lib/netlist/solver/nld_ms_direct2.h b/src/lib/netlist/solver/nld_ms_direct2.h index 03050083158..a24ceb9ef71 100644 --- a/src/lib/netlist/solver/nld_ms_direct2.h +++ b/src/lib/netlist/solver/nld_ms_direct2.h @@ -47,9 +47,9 @@ namespace solver std::array new_V = {v0, v1}; this->m_stat_calculations++; - const float_type err = (newton_raphson ? this->delta(new_V) : 0.0); + const float_type err = (newton_raphson ? this->delta(new_V) : plib::constants::zero()); this->store(new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; + return (err > static_cast(this->m_params.m_accuracy)) ? 2 : 1; } }; diff --git a/src/lib/netlist/solver/nld_ms_gcr.h b/src/lib/netlist/solver/nld_ms_gcr.h index 66b56397a3d..d7f9716a5d7 100644 --- a/src/lib/netlist/solver/nld_ms_gcr.h +++ b/src/lib/netlist/solver/nld_ms_gcr.h @@ -227,7 +227,7 @@ namespace solver unsigned matrix_solver_GCR_t::vsolve_non_dynamic(const bool newton_raphson) { /* populate matrix */ - mat.set_scalar(0.0); + mat.set_scalar(plib::constants::zero()); this->fill_matrix(RHS); /* now solve it */ @@ -247,9 +247,9 @@ namespace solver this->m_stat_calculations++; - const FT err = (newton_raphson ? this->delta(new_V) : 0.0); + const FT err = (newton_raphson ? this->delta(new_V) : plib::constants::zero()); this->store(new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; + return (err > static_cast(this->m_params.m_accuracy)) ? 2 : 1; } } // namespace solver diff --git a/src/lib/netlist/solver/nld_ms_gmres.h b/src/lib/netlist/solver/nld_ms_gmres.h index 9a3c4e2f1bc..6d9e7929fe3 100644 --- a/src/lib/netlist/solver/nld_ms_gmres.h +++ b/src/lib/netlist/solver/nld_ms_gmres.h @@ -103,17 +103,17 @@ namespace solver plib::parray RHS(iN); - m_ops.m_mat.set_scalar(0.0); + m_ops.m_mat.set_scalar(plib::constants::zero()); /* populate matrix and V for first estimate */ this->fill_matrix(RHS); for (std::size_t k = 0; k < iN; k++) { - this->m_new_V[k] = this->m_terms[k].getV(); + this->m_new_V[k] = this->m_terms[k].template getV(); } - const float_type accuracy = this->m_params.m_accuracy; + const float_type accuracy = static_cast(this->m_params.m_accuracy); auto iter = std::max(plib::constants::one(), this->m_params.m_gs_loops()); auto gsl = m_gmres.solve(m_ops, this->m_new_V, RHS, iter, accuracy); @@ -127,9 +127,9 @@ namespace solver return matrix_solver_direct_t::vsolve_non_dynamic(newton_raphson); } - const float_type err = (newton_raphson ? this->delta(this->m_new_V) : 0.0); + const float_type err = (newton_raphson ? this->delta(this->m_new_V) : plib::constants::zero()); this->store(this->m_new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; + return (err > static_cast(this->m_params.m_accuracy)) ? 2 : 1; } diff --git a/src/lib/netlist/solver/nld_ms_sm.h b/src/lib/netlist/solver/nld_ms_sm.h index 4828dbfe72a..679819a2627 100644 --- a/src/lib/netlist/solver/nld_ms_sm.h +++ b/src/lib/netlist/solver/nld_ms_sm.h @@ -124,15 +124,15 @@ namespace solver for (std::size_t j = 0; j < kN; j++) { W(i,j) = lA(i,j) = A(i,j); - Ainv(i,j) = 0.0; + Ainv(i,j) = plib::constants::zero(); } - Ainv(i,i) = 1.0; + Ainv(i,i) = plib::constants::one(); } /* down */ for (std::size_t i = 0; i < kN; i++) { /* FIXME: Singular matrix? */ - const float_type f = 1.0 / W(i,i); + const float_type f = plib::reciprocal(W(i,i)); const auto * const p = this->m_terms[i].m_nzrd.data(); const std::size_t e = this->m_terms[i].m_nzrd.size(); @@ -144,7 +144,8 @@ namespace solver { const unsigned j = pb[jb]; const float_type f1 = - W(j,i) * f; - if (f1 != 0.0) + // FIXME: comparison to zero + if (f1 != plib::constants::zero()) { for (std::size_t k = 0; k < e; k++) W(j,p[k]) += W(i,p[k]) * f1; @@ -157,11 +158,12 @@ namespace solver for (std::size_t i = kN; i-- > 0; ) { /* FIXME: Singular matrix? */ - const float_type f = 1.0 / W(i,i); + const float_type f = plib::reciprocal(W(i,i)); for (std::size_t j = i; j-- > 0; ) { const float_type f1 = - W(j,i) * f; - if (f1 != 0.0) + // FIXME: comparison to zero + if (f1 != plib::constants::zero()) { for (std::size_t k = i; k < kN; k++) W(j,k) += W(i,k) * f1; @@ -185,7 +187,7 @@ namespace solver const std::size_t kN = this->size(); for (std::size_t i=0; i::zero(); for (std::size_t k=0; k::zero()) cols[colcount++] = col; } if (colcount > 0) { - float_type lamba = 0.0; + auto lamba(plib::constants::zero()); std::array w = {0}; /* compute w and lamba */ @@ -258,11 +261,12 @@ namespace solver w[k] += Ainv(col,k) * f; /* Transpose(Ainv) * v */ } - lamba = -1.0 / (1.0 + lamba); + lamba = -plib::reciprocal(plib::constants::one() + lamba); for (std::size_t i=0; i::zero()) for (std::size_t k = 0; k < iN; k++) Ainv(i,k) += f * w[k]; } @@ -275,9 +279,9 @@ namespace solver this->LE_compute_x(new_V); - const float_type err = (newton_raphson ? this->delta(new_V) : 0.0); + const float_type err = (newton_raphson ? this->delta(new_V) : plib::constants::zero()); this->store(new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; + return (err > static_cast(this->m_params.m_accuracy)) ? 2 : 1; } template diff --git a/src/lib/netlist/solver/nld_ms_sor.h b/src/lib/netlist/solver/nld_ms_sor.h index 8bf4bf48721..08d8a358a45 100644 --- a/src/lib/netlist/solver/nld_ms_sor.h +++ b/src/lib/netlist/solver/nld_ms_sor.h @@ -34,9 +34,9 @@ namespace solver const solver_parameters_t *params, const std::size_t size) : matrix_solver_direct_t(anetlist, name, nets, params, size) , m_lp_fact(*this, "m_lp_fact", 0) - , w(size, 0.0) - , one_m_w(size, 0.0) - , RHS(size, 0.0) + , w(size, plib::constants::zero()) + , one_m_w(size, plib::constants::zero()) + , RHS(size, plib::constants::zero()) //, new_V(size, 0.0) { } @@ -70,13 +70,13 @@ namespace solver * omega = 2.0 / (1.0 + std::sqrt(1-rho)) */ - const float_type ws = this->m_params.m_gs_sor; + const auto ws(static_cast(this->m_params.m_gs_sor)); for (std::size_t k = 0; k < iN; k++) { - float_type gtot_t = 0.0; - float_type gabs_t = 0.0; - float_type RHS_t = 0.0; + nl_fptype gtot_t = 0.0; + nl_fptype gabs_t = 0.0; + nl_fptype RHS_t = 0.0; const std::size_t term_count = this->m_terms[k].count(); const nl_fptype * const gt = this->m_gtn[k]; @@ -84,7 +84,7 @@ namespace solver const nl_fptype * const Idr = this->m_Idrn[k]; auto other_cur_analog = this->m_connected_net_Vn[k]; - this->m_new_V[k] = this->m_terms[k].getV(); + this->m_new_V[k] = this->m_terms[k].template getV(); for (std::size_t i = 0; i < term_count; i++) { @@ -95,7 +95,7 @@ namespace solver for (std::size_t i = this->m_terms[k].railstart(); i < term_count; i++) RHS_t = RHS_t - go[i] * *other_cur_analog[i]; - RHS[k] = RHS_t; + RHS[k] = static_cast(RHS_t); if (this->m_params.m_use_gabs) { @@ -105,23 +105,23 @@ namespace solver gabs_t *= plib::constants::cast(0.5); // derived by try and error if (gabs_t <= gtot_t) { - w[k] = ws / gtot_t; + w[k] = ws / static_cast(gtot_t); one_m_w[k] = plib::constants::one() - ws; } else { - w[k] = plib::constants::one() / (gtot_t + gabs_t); - one_m_w[k] = plib::constants::one() - plib::constants::one() * gtot_t / (gtot_t + gabs_t); + w[k] = plib::reciprocal(static_cast(gtot_t + gabs_t)); + one_m_w[k] = plib::constants::one() - plib::constants::one() * static_cast(gtot_t / (gtot_t + gabs_t)); } } else { - w[k] = ws / gtot_t; + w[k] = ws / static_cast(gtot_t); one_m_w[k] = plib::constants::one() - ws; } } - const float_type accuracy = this->m_params.m_accuracy; + const auto accuracy(static_cast(this->m_params.m_accuracy)); do { resched = false; @@ -132,9 +132,9 @@ namespace solver const std::size_t railstart = this->m_terms[k].railstart(); const nl_fptype * go = this->m_gonn[k]; - float_type Idrive = 0.0; + float_type Idrive = plib::constants::zero(); for (std::size_t i = 0; i < railstart; i++) - Idrive = Idrive - go[i] * this->m_new_V[static_cast(net_other[i])]; + Idrive = Idrive - static_cast(go[i]) * this->m_new_V[static_cast(net_other[i])]; const float_type new_val = this->m_new_V[k] * one_m_w[k] + (Idrive + RHS[k]) * w[k]; @@ -158,9 +158,9 @@ namespace solver return matrix_solver_direct_t::vsolve_non_dynamic(newton_raphson); } - const float_type err = (newton_raphson ? this->delta(this->m_new_V) : 0.0); + const float_type err = (newton_raphson ? this->delta(this->m_new_V) : plib::constants::zero()); this->store(this->m_new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; + return (err > static_cast(this->m_params.m_accuracy)) ? 2 : 1; } } // namespace solver diff --git a/src/lib/netlist/solver/nld_ms_sor_mat.h b/src/lib/netlist/solver/nld_ms_sor_mat.h index 8a1575ebfc7..da862e182ec 100644 --- a/src/lib/netlist/solver/nld_ms_sor_mat.h +++ b/src/lib/netlist/solver/nld_ms_sor_mat.h @@ -37,7 +37,7 @@ namespace solver const solver_parameters_t *params, std::size_t size) : matrix_solver_direct_t(anetlist, name, nets, params, size) , m_Vdelta(*this, "m_Vdelta", std::vector(size)) - , m_omega(*this, "m_omega", params->m_gs_sor) + , m_omega(*this, "m_omega", static_cast(params->m_gs_sor)) , m_lp_fact(*this, "m_lp_fact", 0) { } @@ -84,7 +84,7 @@ namespace solver { float_type sq = 0; float_type sqo = 0; - const float_type rez_cts = 1.0 / this->current_timestep(); + const float_type rez_cts = plib::reciprocal(this->current_timestep()); for (unsigned k = 0; k < this->size(); k++) { const analog_net_t *n = this->m_nets[k]; @@ -158,11 +158,11 @@ namespace solver #endif for (std::size_t k = 0; k < iN; k++) - this->m_new_V[k] = this->m_terms[k].getV(); + this->m_new_V[k] = this->m_terms[k].template getV(); do { resched = false; - float_type cerr = 0.0; + FT cerr = plib::constants::zero(); for (std::size_t k = 0; k < iN; k++) { @@ -177,7 +177,7 @@ namespace solver FT w = m_omega / this->m_A[k][k]; if (this->m_params.m_use_gabs) { - FT gabs_t = 0.0; + FT gabs_t = plib::constants::zero(); for (std::size_t i = 0; i < e; i++) if (p[i] != k) gabs_t = gabs_t + std::abs(this->m_A[k][p[i]]); @@ -194,7 +194,7 @@ namespace solver this->m_new_V[k] += delta; } - if (cerr > this->m_params.m_accuracy) + if (cerr > static_cast(this->m_params.m_accuracy)) { resched = true; } @@ -210,9 +210,9 @@ namespace solver return matrix_solver_direct_t::solve_non_dynamic(newton_raphson); } - const float_type err = (newton_raphson ? this->delta(this->m_new_V) : 0.0); + const float_type err = (newton_raphson ? this->delta(this->m_new_V) : plib::constants::zero()); this->store(this->m_new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; + return (err > static_cast(this->m_params.m_accuracy)) ? 2 : 1; } diff --git a/src/lib/netlist/solver/nld_ms_w.h b/src/lib/netlist/solver/nld_ms_w.h index 3470cd86f22..b89edd70b74 100644 --- a/src/lib/netlist/solver/nld_ms_w.h +++ b/src/lib/netlist/solver/nld_ms_w.h @@ -134,15 +134,15 @@ namespace solver for (std::size_t j = 0; j < kN; j++) { W(i,j) = lA(i,j) = A(i,j); - Ainv(i,j) = 0.0; + Ainv(i,j) = plib::constants::zero(); } - Ainv(i,i) = 1.0; + Ainv(i,i) = plib::constants::one(); } /* down */ for (std::size_t i = 0; i < kN; i++) { /* FIXME: Singular matrix? */ - const float_type f = 1.0 / W(i,i); + const float_type f = plib::reciprocal(W(i,i)); const auto * const p = this->m_terms[i].m_nzrd.data(); const size_t e = this->m_terms[i].m_nzrd.size(); @@ -154,7 +154,8 @@ namespace solver { const auto j = pb[jb]; const float_type f1 = - W(j,i) * f; - if (f1 != 0.0) + // FIXME: comparison to zero + if (f1 != plib::constants::zero()) { for (std::size_t k = 0; k < e; k++) W(j,p[k]) += W(i,p[k]) * f1; @@ -167,11 +168,12 @@ namespace solver for (std::size_t i = kN; i-- > 0; ) { /* FIXME: Singular matrix? */ - const float_type f = 1.0 / W(i,i); + const float_type f = plib::reciprocal(W(i,i)); for (std::size_t j = i; j-- > 0; ) { const float_type f1 = - W(j,i) * f; - if (f1 != 0.0) + // FIXME: comparison to zero + if (f1 != plib::constants::zero()) { for (std::size_t k = i; k < kN; k++) W(j,k) += W(i,k) * f1; @@ -194,7 +196,7 @@ namespace solver const std::size_t kN = this->size(); for (std::size_t i=0; i::zero(); for (std::size_t k=0; k::zero(); for (unsigned k = 0; k < iN; k++) tmp += VT(r,k) * new_V[k]; w[i] = tmp; @@ -265,17 +267,18 @@ namespace solver for (unsigned i = 0; i < rowcount; i++) for (unsigned k=0; k< rowcount; k++) - H[i][k] = 0.0; + H[i][k] = plib::constants::zero(); for (unsigned i = 0; i < rowcount; i++) - H[i][i] = 1.0; + H[i][i] = plib::constants::one(); /* Construct H = (I + VT*Z) */ for (unsigned i = 0; i < rowcount; i++) for (unsigned k=0; k< colcount[i]; k++) { const unsigned col = cols[i][k]; float_type f = VT(rows[i],col); - if (f!=0.0) + // FIXME: comparison to zero + if (f != plib::constants::zero()) for (unsigned j= 0; j < rowcount; j++) H[i][j] += f * Ainv(col,rows[j]); } @@ -283,14 +286,16 @@ namespace solver /* Gaussian elimination of H */ for (unsigned i = 0; i < rowcount; i++) { - if (H[i][i] == 0.0) + // FIXME: comparison to zero + if (H[i][i] == plib::constants::zero()) plib::perrlogger("{} H singular\n", this->name()); - const float_type f = 1.0 / H[i][i]; + const float_type f = plib::reciprocal(H[i][i]); for (unsigned j = i+1; j < rowcount; j++) { const float_type f1 = - f * H[j][i]; - if (f1!=0.0) + // FIXME: comparison to zero + if (f1 != plib::constants::zero()) { float_type *pj = &H[j][i+1]; const float_type *pi = &H[i][i+1]; @@ -317,7 +322,7 @@ namespace solver /* x = y - Zt */ for (unsigned i=0; i::zero(); for (unsigned j=0; j::zero(); for (unsigned j=0; j 1e-6) + if (std::abs(tmp-RHS(i)) > static_cast(1e-6)) plib::perrlogger("{} failed on row {}: {} RHS: {}\n", this->name(), i, std::abs(tmp-RHS(i)), RHS(i)); } - const float_type err = (newton_raphson ? this->delta(new_V) : 0.0); + const float_type err = (newton_raphson ? this->delta(new_V) : plib::constants::zero()); this->store(new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; + return (err > static_cast(this->m_params.m_accuracy)) ? 2 : 1; } template diff --git a/src/lib/netlist/solver/nld_solver.cpp b/src/lib/netlist/solver/nld_solver.cpp index 506bff88dab..86772fac6a0 100644 --- a/src/lib/netlist/solver/nld_solver.cpp +++ b/src/lib/netlist/solver/nld_solver.cpp @@ -77,7 +77,7 @@ namespace devices netlist_time now(exec().time()); /* force solving during start up if there are no time-step devices */ /* FIXME: Needs a more elegant solution */ - bool force_solve = (now < netlist_time::from_double(2 * m_params.m_max_timestep)); + bool force_solve = (now < netlist_time::from_fp(2 * m_params.m_max_timestep)); std::size_t nthreads = std::min(static_cast(m_params.m_parallel()), plib::omp::get_max_threads()); @@ -105,7 +105,7 @@ namespace devices /* step circuit */ if (!m_Q_step.net().is_queued()) { - m_Q_step.net().toggle_and_push_to_queue(netlist_time::from_double(m_params.m_max_timestep)); + m_Q_step.net().toggle_and_push_to_queue(netlist_time::from_fp(m_params.m_max_timestep)); } } @@ -151,6 +151,112 @@ namespace devices return plib::unique_ptr(); } + template + plib::unique_ptr NETLIB_NAME(solver)::create_solvers( + const pstring &sname, + analog_net_t::list_t &nets) + { + std::size_t net_count = nets.size(); + switch (net_count) + { + case 1: + return plib::make_unique>(state(), sname, nets, &m_params); + break; + case 2: + return plib::make_unique>(state(), sname, nets, &m_params); + break; +#if 0 + case 3: + return create_solver(3, sname, nets); + break; + case 4: + return create_solver(4, sname, nets); + break; + case 5: + return create_solver(5, sname, nets); + break; + case 6: + return create_solver(6, sname, nets); + break; + case 7: + return create_solver(7, sname, nets); + break; + case 8: + return create_solver(8, sname, nets); + break; + case 9: + return create_solver(9, sname, nets); + break; + case 10: + return create_solver(10, sname, nets); + break; +#if 0 + case 11: + return create_solver(11, sname); + break; + case 12: + return create_solver(12, sname); + break; + case 15: + return create_solver(15, sname); + break; + case 31: + return create_solver(31, sname); + break; + case 35: + return create_solver(35, sname); + break; + case 43: + return create_solver(43, sname); + break; + case 49: + return create_solver(49, sname); + break; +#endif +#if 1 + case 87: + return create_solver(86, sname, nets); + break; +#endif +#endif + default: + log().info(MI_NO_SPECIFIC_SOLVER(net_count)); + if (net_count <= 8) + { + return create_solver(net_count, sname, nets); + } + else if (net_count <= 16) + { + return create_solver(net_count, sname, nets); + } + else if (net_count <= 32) + { + return create_solver(net_count, sname, nets); + } + else if (net_count <= 64) + { + return create_solver(net_count, sname, nets); + } + else if (net_count <= 128) + { + return create_solver(net_count, sname, nets); + } + else if (net_count <= 256) + { + return create_solver(net_count, sname, nets); + } + else if (net_count <= 512) + { + return create_solver(net_count, sname, nets); + } + else + { + return create_solver(net_count, sname, nets); + } + break; + } + } + struct net_splitter { @@ -212,7 +318,6 @@ namespace devices void NETLIB_NAME(solver)::post_start() { - log().verbose("Scanning net groups ..."); // determine net groups @@ -225,105 +330,18 @@ namespace devices for (auto & grp : splitter.groups) { plib::unique_ptr ms; - std::size_t net_count = grp.size(); pstring sname = plib::pfmt("Solver_{1}")(m_mat_solvers.size()); - switch (net_count) + switch (m_params.m_fp_type()) { - case 1: - ms = plib::make_unique>(state(), sname, grp, &m_params); - break; - case 2: - ms = plib::make_unique>(state(), sname, grp, &m_params); - break; -#if 0 - case 3: - ms = create_solver(3, sname, grp); - break; - case 4: - ms = create_solver(4, sname, grp); - break; - case 5: - ms = create_solver(5, sname, grp); - break; - case 6: - ms = create_solver(6, sname, grp); - break; - case 7: - ms = create_solver(7, sname, grp); - break; - case 8: - ms = create_solver(8, sname, grp); - break; - case 9: - ms = create_solver(9, sname, grp); - break; - case 10: - ms = create_solver(10, sname, grp); + case solver::matrix_fp_type_e::FLOAT: + ms = create_solvers(sname, grp); break; - #if 0 - case 11: - ms = create_solver(11, sname); + case solver::matrix_fp_type_e::DOUBLE: + ms = create_solvers(sname, grp); break; - case 12: - ms = create_solver(12, sname); - break; - case 15: - ms = create_solver(15, sname); - break; - case 31: - ms = create_solver(31, sname); - break; - case 35: - ms = create_solver(35, sname); - break; - case 43: - ms = create_solver(43, sname); - break; - case 49: - ms = create_solver(49, sname); - break; - #endif - #if 1 - case 87: - ms = create_solver(86, sname, grp); - break; - #endif - #endif - default: - log().info(MI_NO_SPECIFIC_SOLVER(net_count)); - if (net_count <= 8) - { - ms = create_solver(net_count, sname, grp); - } - else if (net_count <= 16) - { - ms = create_solver(net_count, sname, grp); - } - else if (net_count <= 32) - { - ms = create_solver(net_count, sname, grp); - } - else if (net_count <= 64) - { - ms = create_solver(net_count, sname, grp); - } - else if (net_count <= 128) - { - ms = create_solver(net_count, sname, grp); - } - else if (net_count <= 256) - { - ms = create_solver(net_count, sname, grp); - } - else if (net_count <= 512) - { - ms = create_solver(net_count, sname, grp); - } - else - { - ms = create_solver(net_count, sname, grp); - } + case solver::matrix_fp_type_e::LONGDOUBLE: + ms = create_solvers(sname, grp); break; } diff --git a/src/lib/netlist/solver/nld_solver.h b/src/lib/netlist/solver/nld_solver.h index dba6b69f358..02b231655ee 100644 --- a/src/lib/netlist/solver/nld_solver.h +++ b/src/lib/netlist/solver/nld_solver.h @@ -60,6 +60,10 @@ namespace devices template plib::unique_ptr create_solver(std::size_t size, const pstring &solvername, analog_net_t::list_t &nets); + + template + plib::unique_ptr create_solvers( + const pstring &sname, analog_net_t::list_t &nets); }; } // namespace devices -- cgit v1.2.3