summaryrefslogtreecommitdiffstatshomepage
path: root/src/lib
diff options
context:
space:
mode:
authorGravatarGravatar couriersud <couriersud@gmx.org>2019-11-02 12:25:50 +0100
committerGravatarGravatar couriersud <couriersud@gmx.org>2019-11-02 12:25:50 +0100
commit6573037934e93af6d103fd532049b28af0d01c8f (patch)
tree8b749a840e9cc34b4fa4e48a7841c6679ba9596c /src/lib
parent248d3629d086f821bd01c6439fb6fb2ba53a68bd (diff)
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.
Diffstat (limited to 'src/lib')
-rw-r--r--src/lib/netlist/analog/nld_bjt.cpp8
-rw-r--r--src/lib/netlist/analog/nld_generic_models.h4
-rw-r--r--src/lib/netlist/analog/nld_mosfet.cpp4
-rw-r--r--src/lib/netlist/analog/nld_switches.cpp2
-rw-r--r--src/lib/netlist/analog/nlid_twoterm.cpp2
-rw-r--r--src/lib/netlist/analog/nlid_twoterm.h4
-rw-r--r--src/lib/netlist/devices/nld_74123.cpp2
-rw-r--r--src/lib/netlist/devices/nld_74ls629.cpp2
-rw-r--r--src/lib/netlist/devices/nld_log.cpp6
-rw-r--r--src/lib/netlist/devices/nld_mm5837.cpp4
-rw-r--r--src/lib/netlist/devices/nld_r2r_dac.cpp2
-rw-r--r--src/lib/netlist/devices/nld_system.cpp2
-rw-r--r--src/lib/netlist/devices/nlid_system.h22
-rw-r--r--src/lib/netlist/nl_config.h26
-rw-r--r--src/lib/netlist/plib/gmres.h25
-rw-r--r--src/lib/netlist/plib/mat_cr.h8
-rw-r--r--src/lib/netlist/plib/pfmtlog.h6
-rw-r--r--src/lib/netlist/plib/ptime.h28
-rw-r--r--src/lib/netlist/plib/putil.h13
-rw-r--r--src/lib/netlist/prg/nltool.cpp15
-rw-r--r--src/lib/netlist/solver/nld_matrix_solver.cpp10
-rw-r--r--src/lib/netlist/solver/nld_matrix_solver.h55
-rw-r--r--src/lib/netlist/solver/nld_ms_direct.h12
-rw-r--r--src/lib/netlist/solver/nld_ms_direct1.h4
-rw-r--r--src/lib/netlist/solver/nld_ms_direct2.h4
-rw-r--r--src/lib/netlist/solver/nld_ms_gcr.h6
-rw-r--r--src/lib/netlist/solver/nld_ms_gmres.h10
-rw-r--r--src/lib/netlist/solver/nld_ms_sm.h30
-rw-r--r--src/lib/netlist/solver/nld_ms_sor.h36
-rw-r--r--src/lib/netlist/solver/nld_ms_sor_mat.h16
-rw-r--r--src/lib/netlist/solver/nld_ms_w.h43
-rw-r--r--src/lib/netlist/solver/nld_solver.cpp212
-rw-r--r--src/lib/netlist/solver/nld_solver.h4
33 files changed, 364 insertions, 263 deletions
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<nl_fptype>(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<nl_fptype>(m_I())));
+ m_writer.writeline(plib::pfmt("{1:.9} {2}").e(exec().time().as_fp<nl_fptype>()).e(static_cast<nl_fptype>(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<nl_fptype>(m_I() - m_I2())));
+ m_writer.writeline(plib::pfmt("{1:.9} {2}").e(exec().time().as_fp<nl_fptype>()).e(static_cast<nl_fptype>(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<nl_fptype>(), 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<nl_fptype>::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<nl_fptype>(1 << m_num())
* static_cast<nl_fptype>(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<decltype(m_offset())>(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<decltype(m_freq())>(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<decltype(m_freq())>(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<decltype(m_freq())>(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<decltype(m_freq())>(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<nl_fptype>();
+ 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<decltype(m_freq())>(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<decltype(m_freq())>(1.0 / (m_freq()*2.0));
std::vector<pstring> pat(plib::psplit(m_pattern(),","));
- m_off = netlist_time::from_double(m_offset());
+ m_off = netlist_time::from_fp(m_offset());
std::array<std::int64_t, 32> 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<long double>
+ {
+ 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<float>
{
- 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<FT>(n, m_v[kp1]));
- if (m_ht[kp1][k] != 0.0)
- vec_scale(n, m_v[kp1], constants<FT>::one() / m_ht[kp1][k]);
+ // FIXME: comparison to zero
+ if (m_ht[kp1][k] != plib::constants<FT>::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<FT>::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<float_type>::zero());
const std::size_t n = size();
@@ -400,7 +401,7 @@ namespace plib
vec_set_scalar(RESTART+1, m_g, +constants<FT>::zero());
m_g[0] = rho;
- vec_mult_scalar(n, m_v[0], residual, constants<FT>::one() / rho);
+ vec_mult_scalar(n, m_v[0], residual, reciprocal(rho));
if (do_k<RESTART-1>(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<T>::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
@@ -128,6 +128,12 @@ struct ptype_traits<double> : ptype_traits_base<double>
};
template<>
+struct ptype_traits<long double> : ptype_traits_base<long double>
+{
+ static char32_t fmt_spec() { return 'f'; }
+};
+
+template<>
struct ptype_traits<char *> : ptype_traits_base<char *>
{
static const char *cast(const char *x) { return x; }
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 <typename FT>
+ constexpr typename std::enable_if<std::is_floating_point<FT>::value, FT>::type
+ as_fp() const noexcept
{
- return static_cast<double>(m_time) * inv_res;
+ return static_cast<FT>(m_time) * inv_res<FT>();
}
+ constexpr double as_double() const noexcept { return as_fp<double>(); }
+ constexpr double as_float() const noexcept { return as_fp<float>(); }
+ constexpr double as_long_double() const noexcept { return as_fp<long double>(); }
+
// 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<internal_type>(std::floor(t * static_cast<double>(RES) + 0.5)), RES); }
+
+ template <typename FT>
+ static constexpr const typename std::enable_if<std::is_floating_point<FT>::value, ptime>::type
+ from_fp(const FT t) noexcept { return ptime(static_cast<internal_type>(std::floor(t * static_cast<FT>(RES) + 0.5)), RES); }
+
+ static constexpr const ptime from_double(const double t) noexcept
+ { return from_fp<double>(t); }
+
+ static constexpr const ptime from_float(const float t) noexcept
+ { return from_fp<float>(t); }
+
+ static constexpr const ptime from_long_double(const long double t) noexcept
+ { return from_fp<long double>(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<double>(RES);
+ template <typename FT>
+ static constexpr FT inv_res() noexcept { return static_cast<FT>(1.0) / static_cast<FT>(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<T>(v); }
};
+ /*! typesafe reciprocal function
+ *
+ * @tparam T type of the argument
+ * @param v argument
+ * @return reciprocal of argument
+ */
+ template <typename T>
+ static inline constexpr typename std::enable_if<std::is_floating_point<T>::value, T>::type
+ reciprocal(T v) noexcept
+ {
+ return constants<T>::one() / v;
+ }
+
static_assert(noexcept(constants<double>::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<pstring>({"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<unsigned> 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<nl_fptype>(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<nl_fptype>();
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<nl_fptype>());
return next_time_step;
}
@@ -558,7 +558,7 @@ namespace solver
static_cast<nl_fptype>(this->m_stat_newton_raphson) / static_cast<nl_fptype>(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<nl_fptype>(this->m_stat_calculations) / this->exec().time().as_double(),
+ static_cast<nl_fptype>(this->m_stat_calculations) / this->exec().time().as_fp<nl_fptype>(),
this->m_iterative_fail,
100.0 * static_cast<nl_fptype>(this->m_iterative_fail)
/ static_cast<nl_fptype>(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<nl_fptype>()) // 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<decltype(m_max_timestep)>();
if (m_dynamic_ts)
{
@@ -87,6 +94,7 @@ namespace solver
param_fp_t m_freq;
param_fp_t m_gs_sor;
param_enum_t<matrix_type_e> m_method;
+ param_enum_t<matrix_fp_type_e> m_fp_type;
param_fp_t m_accuracy;
param_num_t<std::size_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 <typename FT, typename = std::enable_if<std::is_floating_point<FT>::value, void>>
+ FT getV() const noexcept { return static_cast<FT>(m_net->Q_Analog()); }
- void setV(nl_fptype v) noexcept { m_net->set_Q_Analog(v); }
+ template <typename FT, typename = std::enable_if<std::is_floating_point<FT>::value, void>>
+ void setV(FT v) noexcept { m_net->set_Q_Analog(static_cast<nl_fptype>(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<decltype(V[0])>::type;
vtype cerr = 0;
for (std::size_t i = 0; i < iN; i++)
- cerr = std::max(cerr, std::abs(V[i] - static_cast<vtype>(this->m_terms[i].getV())));
+ cerr = std::max(cerr, std::abs(V[i] - this->m_terms[i].getV<vtype>()));
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<float_type>::zero())
- , m_DD_n_m_1(size, plib::constants<float_type>::zero())
- , m_h_n_m_1(size, plib::constants<float_type>::zero())
+ , m_last_V(size, plib::constants<nl_fptype>::zero())
+ , m_DD_n_m_1(size, plib::constants<nl_fptype>::zero())
+ , m_h_n_m_1(size, plib::constants<nl_fptype>::zero())
{
/*
* save states
@@ -407,11 +417,11 @@ namespace solver
plib::parray2D<float_type *, SIZE, 0> m_mat_ptr;
/* state - variable time_stepping */
PALIGNAS_VECTOROPT()
- plib::parray<float_type, SIZE> m_last_V;
+ plib::parray<nl_fptype, SIZE> m_last_V;
PALIGNAS_VECTOROPT()
- plib::parray<float_type, SIZE> m_DD_n_m_1;
+ plib::parray<nl_fptype, SIZE> m_DD_n_m_1;
PALIGNAS_VECTOROPT()
- plib::parray<float_type, SIZE> m_h_n_m_1;
+ plib::parray<nl_fptype, SIZE> m_h_n_m_1;
constexpr std::size_t size() const noexcept { return (SIZE > 0) ? static_cast<std::size_t>(SIZE) : m_dim; }
@@ -428,7 +438,7 @@ namespace solver
// avoid floating point exceptions
const nl_fptype DD_n = std::max(-fp_constants<nl_fptype>::TIMESTEP_MAXDIFF(),
- std::min(+fp_constants<nl_fptype>::TIMESTEP_MAXDIFF(),(t.getV() - m_last_V[k])));
+ std::min(+fp_constants<nl_fptype>::TIMESTEP_MAXDIFF(),(t.template getV<nl_fptype>() - 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<nl_fptype>();
}
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<decltype(*p)>::type;
for (std::size_t i=0; i < n; i++)
- p[i] = 0.0;
+ p[i] = plib::constants<mat_elem_type>::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<FT>(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<nl_fptype>::zero());
+ auto RHS_t (plib::constants<nl_fptype>::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<FT>(RHS_t);
// update diagonal element ...
- *tcr_r[railstart] += gtot_t; //mat.A[mat.diag[k]] += gtot_t;
+ *tcr_r[railstart] += static_cast<FT>(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<FT>::zero());
this->store(m_new_V);
- return (err > this->m_params.m_accuracy) ? 2 : 1;
+ return (err > static_cast<FT>(this->m_params.m_accuracy)) ? 2 : 1;
}
template <typename FT, int SIZE>
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<FT, 1> 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<FT>::zero());
this->store(new_V);
- return (err > this->m_params.m_accuracy) ? 2 : 1;
+ return (err > static_cast<FT>(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<float_type, 2> 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<FT>::zero());
this->store(new_V);
- return (err > this->m_params.m_accuracy) ? 2 : 1;
+ return (err > static_cast<float_type>(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<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson)
{
/* populate matrix */
- mat.set_scalar(0.0);
+ mat.set_scalar(plib::constants<FT>::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<FT>::zero());
this->store(new_V);
- return (err > this->m_params.m_accuracy) ? 2 : 1;
+ return (err > static_cast<FT>(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<FT, SIZE> RHS(iN);
- m_ops.m_mat.set_scalar(0.0);
+ m_ops.m_mat.set_scalar(plib::constants<FT>::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<float_type>();
}
- const float_type accuracy = this->m_params.m_accuracy;
+ const float_type accuracy = static_cast<float_type>(this->m_params.m_accuracy);
auto iter = std::max(plib::constants<std::size_t>::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<FT, SIZE>::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<FT>::zero());
this->store(this->m_new_V);
- return (err > this->m_params.m_accuracy) ? 2 : 1;
+ return (err > static_cast<float_type>(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<FT>::zero();
}
- Ainv(i,i) = 1.0;
+ Ainv(i,i) = plib::constants<FT>::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<float_type>::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<float_type>::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<kN; i++)
- x[i] = 0.0;
+ x[i] = plib::constants<FT>::zero();
for (std::size_t k=0; k<kN; k++)
{
@@ -234,13 +236,14 @@ namespace solver
v[col] = A(row,col) - lA(row,col);
if (incremental)
lA(row,col) = A(row,col);
- if (v[col] != 0.0)
+ // FIXME: comparison to zero
+ if (v[col] != plib::constants<float_type>::zero())
cols[colcount++] = col;
}
if (colcount > 0)
{
- float_type lamba = 0.0;
+ auto lamba(plib::constants<FT>::zero());
std::array<float_type, m_pitch> 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<float_type>::one() + lamba);
for (std::size_t i=0; i<iN; i++)
{
const float_type f = lamba * z[i];
- if (f != 0.0)
+ // FIXME: comparison to zero
+ if (f != plib::constants<float_type>::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<FT>::zero());
this->store(new_V);
- return (err > this->m_params.m_accuracy) ? 2 : 1;
+ return (err > static_cast<FT>(this->m_params.m_accuracy)) ? 2 : 1;
}
template <typename FT, int SIZE>
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<FT, SIZE>(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<FT>::zero())
+ , one_m_w(size, plib::constants<FT>::zero())
+ , RHS(size, plib::constants<FT>::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<float_type>(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<float_type>();
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<float_type>(RHS_t);
if (this->m_params.m_use_gabs)
{
@@ -105,23 +105,23 @@ namespace solver
gabs_t *= plib::constants<nl_fptype>::cast(0.5); // derived by try and error
if (gabs_t <= gtot_t)
{
- w[k] = ws / gtot_t;
+ w[k] = ws / static_cast<float_type>(gtot_t);
one_m_w[k] = plib::constants<FT>::one() - ws;
}
else
{
- w[k] = plib::constants<FT>::one() / (gtot_t + gabs_t);
- one_m_w[k] = plib::constants<FT>::one() - plib::constants<FT>::one() * gtot_t / (gtot_t + gabs_t);
+ w[k] = plib::reciprocal(static_cast<float_type>(gtot_t + gabs_t));
+ one_m_w[k] = plib::constants<FT>::one() - plib::constants<FT>::one() * static_cast<FT>(gtot_t / (gtot_t + gabs_t));
}
}
else
{
- w[k] = ws / gtot_t;
+ w[k] = ws / static_cast<float_type>(gtot_t);
one_m_w[k] = plib::constants<FT>::one() - ws;
}
}
- const float_type accuracy = this->m_params.m_accuracy;
+ const auto accuracy(static_cast<float_type>(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<float_type>::zero();
for (std::size_t i = 0; i < railstart; i++)
- Idrive = Idrive - go[i] * this->m_new_V[static_cast<std::size_t>(net_other[i])];
+ Idrive = Idrive - static_cast<float_type>(go[i]) * this->m_new_V[static_cast<std::size_t>(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<FT, SIZE>::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<FT>::zero());
this->store(this->m_new_V);
- return (err > this->m_params.m_accuracy) ? 2 : 1;
+ return (err > static_cast<float_type>(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<FT, SIZE>(anetlist, name, nets, params, size)
, m_Vdelta(*this, "m_Vdelta", std::vector<float_type>(size))
- , m_omega(*this, "m_omega", params->m_gs_sor)
+ , m_omega(*this, "m_omega", static_cast<float_type>(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<FT>();
do {
resched = false;
- float_type cerr = 0.0;
+ FT cerr = plib::constants<FT>::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<FT>::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<float_type>(this->m_params.m_accuracy))
{
resched = true;
}
@@ -210,9 +210,9 @@ namespace solver
return matrix_solver_direct_t<FT, SIZE>::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<FT>::zero());
this->store(this->m_new_V);
- return (err > this->m_params.m_accuracy) ? 2 : 1;
+ return (err > static_cast<float_type>(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<FT>::zero();
}
- Ainv(i,i) = 1.0;
+ Ainv(i,i) = plib::constants<FT>::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<float_type>::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<float_type>::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<kN; i++)
- x[i] = 0.0;
+ x[i] = plib::constants<FT>::zero();
for (std::size_t k=0; k<kN; k++)
{
@@ -257,7 +259,7 @@ namespace solver
for (unsigned i = 0; i < rowcount; i++)
{
const unsigned r = rows[i];
- nl_fptype tmp = 0.0;
+ FT tmp = plib::constants<FT>::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<FT>::zero();
for (unsigned i = 0; i < rowcount; i++)
- H[i][i] = 1.0;
+ H[i][i] = plib::constants<FT>::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<float_type>::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<float_type>::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<float_type>::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<iN; i++)
{
- float_type tmp = 0.0;
+ float_type tmp = plib::constants<FT>::zero();
for (unsigned j=0; j<rowcount;j++)
{
const unsigned row = rows[j];
@@ -332,18 +337,18 @@ namespace solver
if (false)
for (unsigned i=0; i<iN; i++)
{
- float_type tmp = 0.0;
+ float_type tmp = plib::constants<FT>::zero();
for (unsigned j=0; j<iN; j++)
{
tmp += A(i,j) * new_V[j];
}
- if (std::abs(tmp-RHS(i)) > 1e-6)
+ if (std::abs(tmp-RHS(i)) > static_cast<float_type>(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<FT>::zero());
this->store(new_V);
- return (err > this->m_params.m_accuracy) ? 2 : 1;
+ return (err > static_cast<float_type>(this->m_params.m_accuracy)) ? 2 : 1;
}
template <typename FT, int SIZE>
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<decltype(m_params.m_max_timestep)>(2 * m_params.m_max_timestep));
std::size_t nthreads = std::min(static_cast<std::size_t>(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<solver::matrix_solver_t>();
}
+ template <typename FT>
+ plib::unique_ptr<solver::matrix_solver_t> 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<solver::matrix_solver_direct1_t<FT>>(state(), sname, nets, &m_params);
+ break;
+ case 2:
+ return plib::make_unique<solver::matrix_solver_direct2_t<FT>>(state(), sname, nets, &m_params);
+ break;
+#if 0
+ case 3:
+ return create_solver<FT, 3>(3, sname, nets);
+ break;
+ case 4:
+ return create_solver<FT, 4>(4, sname, nets);
+ break;
+ case 5:
+ return create_solver<FT, 5>(5, sname, nets);
+ break;
+ case 6:
+ return create_solver<FT, 6>(6, sname, nets);
+ break;
+ case 7:
+ return create_solver<FT, 7>(7, sname, nets);
+ break;
+ case 8:
+ return create_solver<FT, 8>(8, sname, nets);
+ break;
+ case 9:
+ return create_solver<FT, 9>(9, sname, nets);
+ break;
+ case 10:
+ return create_solver<FT, 10>(10, sname, nets);
+ break;
+#if 0
+ case 11:
+ return create_solver<FT, 11>(11, sname);
+ break;
+ case 12:
+ return create_solver<FT, 12>(12, sname);
+ break;
+ case 15:
+ return create_solver<FT, 15>(15, sname);
+ break;
+ case 31:
+ return create_solver<FT, 31>(31, sname);
+ break;
+ case 35:
+ return create_solver<FT, 35>(35, sname);
+ break;
+ case 43:
+ return create_solver<FT, 43>(43, sname);
+ break;
+ case 49:
+ return create_solver<FT, 49>(49, sname);
+ break;
+#endif
+#if 1
+ case 87:
+ return create_solver<FT,86>(86, sname, nets);
+ break;
+#endif
+#endif
+ default:
+ log().info(MI_NO_SPECIFIC_SOLVER(net_count));
+ if (net_count <= 8)
+ {
+ return create_solver<FT, -8>(net_count, sname, nets);
+ }
+ else if (net_count <= 16)
+ {
+ return create_solver<FT, -16>(net_count, sname, nets);
+ }
+ else if (net_count <= 32)
+ {
+ return create_solver<FT, -32>(net_count, sname, nets);
+ }
+ else if (net_count <= 64)
+ {
+ return create_solver<FT, -64>(net_count, sname, nets);
+ }
+ else if (net_count <= 128)
+ {
+ return create_solver<FT, -128>(net_count, sname, nets);
+ }
+ else if (net_count <= 256)
+ {
+ return create_solver<FT, -256>(net_count, sname, nets);
+ }
+ else if (net_count <= 512)
+ {
+ return create_solver<FT, -512>(net_count, sname, nets);
+ }
+ else
+ {
+ return create_solver<FT, 0>(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<solver::matrix_solver_t> 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<solver::matrix_solver_direct1_t<nl_mat_fptype>>(state(), sname, grp, &m_params);
- break;
- case 2:
- ms = plib::make_unique<solver::matrix_solver_direct2_t<nl_mat_fptype>>(state(), sname, grp, &m_params);
- break;
-#if 0
- case 3:
- ms = create_solver<nl_mat_fptype, 3>(3, sname, grp);
- break;
- case 4:
- ms = create_solver<nl_mat_fptype, 4>(4, sname, grp);
- break;
- case 5:
- ms = create_solver<nl_mat_fptype, 5>(5, sname, grp);
- break;
- case 6:
- ms = create_solver<nl_mat_fptype, 6>(6, sname, grp);
- break;
- case 7:
- ms = create_solver<nl_mat_fptype, 7>(7, sname, grp);
- break;
- case 8:
- ms = create_solver<nl_mat_fptype, 8>(8, sname, grp);
- break;
- case 9:
- ms = create_solver<nl_mat_fptype, 9>(9, sname, grp);
- break;
- case 10:
- ms = create_solver<nl_mat_fptype, 10>(10, sname, grp);
+ case solver::matrix_fp_type_e::FLOAT:
+ ms = create_solvers<float>(sname, grp);
break;
- #if 0
- case 11:
- ms = create_solver<nl_mat_fptype, 11>(11, sname);
+ case solver::matrix_fp_type_e::DOUBLE:
+ ms = create_solvers<double>(sname, grp);
break;
- case 12:
- ms = create_solver<nl_mat_fptype, 12>(12, sname);
- break;
- case 15:
- ms = create_solver<nl_mat_fptype, 15>(15, sname);
- break;
- case 31:
- ms = create_solver<nl_mat_fptype, 31>(31, sname);
- break;
- case 35:
- ms = create_solver<nl_mat_fptype, 35>(35, sname);
- break;
- case 43:
- ms = create_solver<nl_mat_fptype, 43>(43, sname);
- break;
- case 49:
- ms = create_solver<nl_mat_fptype, 49>(49, sname);
- break;
- #endif
- #if 1
- case 87:
- ms = create_solver<nl_mat_fptype,86>(86, sname, grp);
- break;
- #endif
- #endif
- default:
- log().info(MI_NO_SPECIFIC_SOLVER(net_count));
- if (net_count <= 8)
- {
- ms = create_solver<nl_mat_fptype, -8>(net_count, sname, grp);
- }
- else if (net_count <= 16)
- {
- ms = create_solver<nl_mat_fptype, -16>(net_count, sname, grp);
- }
- else if (net_count <= 32)
- {
- ms = create_solver<nl_mat_fptype, -32>(net_count, sname, grp);
- }
- else if (net_count <= 64)
- {
- ms = create_solver<nl_mat_fptype, -64>(net_count, sname, grp);
- }
- else if (net_count <= 128)
- {
- ms = create_solver<nl_mat_fptype, -128>(net_count, sname, grp);
- }
- else if (net_count <= 256)
- {
- ms = create_solver<nl_mat_fptype, -256>(net_count, sname, grp);
- }
- else if (net_count <= 512)
- {
- ms = create_solver<nl_mat_fptype, -512>(net_count, sname, grp);
- }
- else
- {
- ms = create_solver<nl_mat_fptype, 0>(net_count, sname, grp);
- }
+ case solver::matrix_fp_type_e::LONGDOUBLE:
+ ms = create_solvers<long double>(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 <typename FT, int SIZE>
plib::unique_ptr<solver::matrix_solver_t> create_solver(std::size_t size,
const pstring &solvername, analog_net_t::list_t &nets);
+
+ template <typename FT>
+ plib::unique_ptr<solver::matrix_solver_t> create_solvers(
+ const pstring &sname, analog_net_t::list_t &nets);
};
} // namespace devices