netlist: code maintenance and edge case fixes. (nw)

"Edge cases" are compiles were the floting point type used for the core
may be different from the one used by solvers, e.g. double/float.
This does not happen in MAME, here currently this is always
double/double. Since floats are better suited for GPUs whose double
performance is limited it was time to review the code.
This commit fixes warnings about type conversions and precision loss.
In addition there is better support for INT128 and FLOAT128.

The MAT solver has seen some minor performance increases.
This commit is contained in:
couriersud 2020-04-20 20:52:21 +02:00
parent 4c9e3dc171
commit 3243efee70
19 changed files with 175 additions and 132 deletions

View File

@ -239,7 +239,7 @@ native:
$(MAKE) CEXTRAFLAGS="-march=native -msse4.2 -Wall -Wpedantic -Wsign-compare -Wextra "
gcc9:
$(MAKE) CC=g++-9 LD=g++-9 CEXTRAFLAGS="-march=native -msse4.2 -Wall -pedantic -Wpedantic -Wsign-compare -Wextra" EXTRALIBS="-lquadmath"
$(MAKE) CC=g++-9 LD=g++-9 CEXTRAFLAGS="-march=native -fext-numeric-literals -msse4.2 -Wall -pedantic -Wpedantic -Wsign-compare -Wextra" EXTRALIBS="-lquadmath"
clang:
$(MAKE) CC=clang++-11 LD=clang++-11 OBJ=obj/clang CEXTRAFLAGS="-march=native -msse4.2 -Weverything -Wall -pedantic -Wpedantic -Werror -Wno-padded -Wno-weak-vtables -Wno-unused-template -Wno-missing-variable-declarations -Wno-float-equal -Wconversion -Wno-c++98-compat -Wno-c++98-compat-pedantic -Wno-format-nonliteral -Wno-exit-time-destructors"

View File

@ -117,18 +117,19 @@ namespace netlist
nl_fptype sup = (m_supply.VCC().Q_Analog() - m_supply.GND().Q_Analog());
nl_fptype in = m_DUM1.m_P.net().Q_Analog() - m_supply.GND().Q_Analog();
nl_fptype rON = m_base_r() * nlconst::magic(5.0) / sup;
nl_fptype R = std::exp(-(in / sup - 0.55) * 25.0) + rON;
nl_fptype R = std::exp(-(in / sup - nlconst::magic(0.55)) * nlconst::magic(25.0)) + rON;
nl_fptype G = plib::reciprocal(R);
// dI/dVin = (VR1-VR2)*(1.0/sup*b) * exp((Vin/sup-a) * b)
const auto dfdz = 25.0/(R*sup) * m_R.deltaV();
const auto dfdz = nlconst::magic(25.0)/(R*sup) * m_R.deltaV();
const auto Ieq = dfdz * in;
m_R.set_mat( G, -G, 0.0,
-G, G, 0.0);
const auto zero(nlconst::zero());
m_R.set_mat( G, -G, zero,
-G, G, zero);
//VIN VR1
m_DUM1.set_mat( 0.0, 0.0, 0.0, // IIN
dfdz, 0.0, Ieq); // IR1
m_DUM2.set_mat( 0.0, 0.0, 0.0, // IIN
-dfdz, 0.0, -Ieq); // IR2
m_DUM1.set_mat( zero, zero, zero, // IIN
dfdz, zero, Ieq); // IR1
m_DUM2.set_mat( zero, zero, zero, // IIN
-dfdz, zero, -Ieq); // IR2
}
NETLIB_IS_DYNAMIC(true)

View File

@ -1789,7 +1789,7 @@ namespace netlist
plib::pfunction<nl_fptype> func;
func.compile_infix(p, {});
auto valx = func.evaluate();
if (std::is_integral<T>::value)
if (plib::is_integral<T>::value)
if (plib::abs(valx - plib::trunc(valx)) > nlconst::magic(1e-6))
throw nl_exception(MF_INVALID_NUMBER_CONVERSION_1_2(device.name() + "." + name, p));
m_param = static_cast<T>(valx);

View File

@ -55,10 +55,12 @@
/// \brief Compile in academic solvers
///
/// Set to 0 to disable compiling in the following solvers:
/// Set to 0 to disable compiling the following solvers:
///
/// Sherman-Morrison, Woodbury, SOR and GMRES
///
/// In addition, compilation of FLOAT, LONGDOUBLE and FLOATQ128 matrix
/// solvers will be disabled.
/// GMRES may be added to productive solvers in the future
/// again. Compiling in all solvers may increase compile
/// time significantly.
@ -113,18 +115,20 @@
/// \brief Support float type for matrix calculations.
///
/// Defaults to off to provide faster build times
/// Defaults to NL_USE_ACADEMIC_SOLVERS to provide faster build times
#ifndef NL_USE_FLOAT_MATRIX
#define NL_USE_FLOAT_MATRIX (0)
#define NL_USE_FLOAT_MATRIX (NL_USE_ACADEMIC_SOLVERS)
//#define NL_USE_FLOAT_MATRIX 1
#endif
/// \brief Support long double type for matrix calculations.
///
/// Defaults to off to provide faster build times
/// Defaults to NL_USE_ACADEMIC_SOLVERS to provide faster build times
#ifndef NL_USE_LONG_DOUBLE_MATRIX
#define NL_USE_LONG_DOUBLE_MATRIX (0)
#define NL_USE_LONG_DOUBLE_MATRIX (NL_USE_ACADEMIC_SOLVERS)
//#define NL_USE_LONG_DOUBLE_MATRIX 1
#endif
//============================================================
@ -187,6 +191,8 @@ static constexpr const int NETLIST_CLOCK = 1'000'000'000;
///
using nl_fptype = double;
//using nl_fptype = long double;
//using nl_fptype = float;
namespace netlist
{
@ -244,27 +250,27 @@ namespace netlist
};
#if (NL_USE_FLOAT128)
/// \brief Specific constants for __float128 floating point type
/// \brief Specific constants for FLOAT128 floating point type
///
template <>
struct fp_constants<__float128>
struct fp_constants<FLOAT128>
{
#if 0
// MAME compile doesn't support Q
static inline constexpr __float128 DIODE_MAXDIFF() noexcept { return 1e100Q; }
static inline constexpr __float128 DIODE_MAXVOLT() noexcept { return 300.0Q; }
static inline constexpr FLOAT128 DIODE_MAXDIFF() noexcept { return 1e100Q; }
static inline constexpr FLOAT128 DIODE_MAXVOLT() noexcept { return 300.0Q; }
static inline constexpr __float128 TIMESTEP_MAXDIFF() noexcept { return 1e100Q; }
static inline constexpr __float128 TIMESTEP_MINDIV() noexcept { return 1e-60Q; }
static inline constexpr FLOAT128 TIMESTEP_MAXDIFF() noexcept { return 1e100Q; }
static inline constexpr FLOAT128 TIMESTEP_MINDIV() noexcept { return 1e-60Q; }
static inline constexpr const char * name() noexcept { return "__float128"; }
static inline constexpr const char * name() noexcept { return "FLOAT128"; }
static inline constexpr const char * suffix() noexcept { return "Q"; }
#else
static inline constexpr __float128 DIODE_MAXDIFF() noexcept { return static_cast<__float128>(1e100L); }
static inline constexpr __float128 DIODE_MAXVOLT() noexcept { return static_cast<__float128>(300.0L); }
static inline constexpr FLOAT128 DIODE_MAXDIFF() noexcept { return static_cast<FLOAT128>(1e100L); }
static inline constexpr FLOAT128 DIODE_MAXVOLT() noexcept { return static_cast<FLOAT128>(300.0L); }
static inline constexpr __float128 TIMESTEP_MAXDIFF() noexcept { return static_cast<__float128>(1e100L); }
static inline constexpr __float128 TIMESTEP_MINDIV() noexcept { return static_cast<__float128>(1e-60L); }
static inline constexpr FLOAT128 TIMESTEP_MAXDIFF() noexcept { return static_cast<FLOAT128>(1e100L); }
static inline constexpr FLOAT128 TIMESTEP_MINDIV() noexcept { return static_cast<FLOAT128>(1e-60L); }
static inline constexpr const char * name() noexcept { return "__float128"; }
static inline constexpr const char * suffix() noexcept { return "Q"; }

View File

@ -254,25 +254,19 @@ namespace netlist
void register_link(const pstring &sin, const pstring &sout);
void register_link_arr(const pstring &terms);
void register_param(const pstring &param, const pstring &value);
// FIXME: quick hack
void register_param_x(const pstring &param, nl_fptype value);
template <typename T>
typename std::enable_if<std::is_floating_point<T>::value || std::is_integral<T>::value>::type
typename std::enable_if<plib::is_floating_point<T>::value || plib::is_integral<T>::value>::type
register_param_val(const pstring &param, T value)
{
register_param_x(param, static_cast<nl_fptype>(value));
}
#if PUSE_FLOAT128
void register_param(const pstring &param, __float128 value)
{
register_param_x(param, static_cast<nl_fptype>(value));
}
#endif
void register_lib_entry(const pstring &name, const pstring &sourcefile);
void register_frontier(const pstring &attach, const pstring &r_IN, const pstring &r_OUT);

View File

@ -134,6 +134,10 @@ typedef __uint128_t UINT128;
typedef __int128_t INT128;
#endif
#if (PUSE_FLOAT128)
typedef __float128 FLOAT128;
#endif
//============================================================
// Check for OpenMP
//============================================================

View File

@ -47,12 +47,12 @@ namespace plib {
#if (PUSE_FLOAT128)
template <>
struct ptype_traits_base<__float128>
struct ptype_traits_base<FLOAT128>
{
// FIXME: need native support at some time
static constexpr const bool is_signed = true;
static char32_t fmt_spec() { return 'f'; }
static inline void streamify(std::ostream &s, const __float128 &v)
static inline void streamify(std::ostream &s, const FLOAT128 &v)
{
s << static_cast<long double>(v);
}
@ -153,7 +153,7 @@ namespace plib {
#if (PUSE_FLOAT128)
template<>
struct ptype_traits<__float128> : ptype_traits_base<__float128>
struct ptype_traits<FLOAT128> : ptype_traits_base<FLOAT128>
{
static char32_t fmt_spec() { return 'f'; }
};
@ -202,20 +202,20 @@ namespace plib {
operator pstring() const { return m_str; }
template <typename T>
typename std::enable_if<std::is_floating_point<T>::value, pfmt &>::type
typename std::enable_if<plib::is_floating_point<T>::value, pfmt &>::type
f(const T &x) {return format_element('f', x); }
template <typename T>
typename std::enable_if<std::is_floating_point<T>::value, pfmt &>::type
typename std::enable_if<plib::is_floating_point<T>::value, pfmt &>::type
e(const T &x) {return format_element('e', x); }
#if 0
#if PUSE_FLOAT128
// FIXME: should use quadmath_snprintf
pfmt & e(const __float128 &x) {return format_element('e', static_cast<long double>(x)); }
pfmt & e(const FLOAT128 &x) {return format_element('e', static_cast<long double>(x)); }
#endif
#endif
template <typename T>
typename std::enable_if<std::is_floating_point<T>::value, pfmt &>::type
typename std::enable_if<plib::is_floating_point<T>::value, pfmt &>::type
g(const T &x) {return format_element('g', x); }
pfmt &operator ()(const void *x) {return format_element('p', x); }
@ -245,14 +245,14 @@ namespace plib {
}
template<typename T>
typename std::enable_if<std::is_integral<T>::value, pfmt &>::type
typename std::enable_if<plib::is_integral<T>::value, pfmt &>::type
x(const T &x)
{
return format_element('x', x);
}
template<typename T>
typename std::enable_if<std::is_integral<T>::value, pfmt &>::type
typename std::enable_if<plib::is_integral<T>::value, pfmt &>::type
o(const T &x)
{
return format_element('o', x);

View File

@ -236,7 +236,7 @@ namespace plib {
}
template <typename NT>
static inline typename std::enable_if<std::is_floating_point<NT>::value, NT>::type
static inline typename std::enable_if<plib::is_floating_point<NT>::value, NT>::type
lfsr_random(std::uint16_t &lfsr) noexcept
{
std::uint16_t lsb = lfsr & 1;
@ -247,7 +247,7 @@ namespace plib {
}
template <typename NT>
static inline typename std::enable_if<std::is_integral<NT>::value, NT>::type
static inline typename std::enable_if<plib::is_integral<NT>::value, NT>::type
lfsr_random(std::uint16_t &lfsr) noexcept
{
std::uint16_t lsb = lfsr & 1;
@ -304,7 +304,7 @@ namespace plib {
template class pfunction<double>;
template class pfunction<long double>;
#if (PUSE_FLOAT128)
template class pfunction<__float128>;
template class pfunction<FLOAT128>;
#endif
} // namespace plib

View File

@ -120,8 +120,9 @@ namespace plib {
extern template class pfunction<double>;
extern template class pfunction<long double>;
#if (PUSE_FLOAT128)
extern template class pfunction<__float128>;
extern template class pfunction<FLOAT128>;
#endif
} // namespace plib
#endif // PEXCEPTION_H_

View File

@ -247,85 +247,85 @@ namespace plib
/// FIXME: limited implementation
///
template <typename T1, typename T2>
static inline T1
pow(T1 v, T2 p) noexcept
static inline
auto pow(T1 v, T2 p) noexcept -> decltype(std::pow(v, p))
{
return std::pow(v, p);
}
#if (PUSE_FLOAT128)
static inline constexpr __float128 reciprocal(__float128 v) noexcept
static inline constexpr FLOAT128 reciprocal(FLOAT128 v) noexcept
{
return constants<__float128>::one() / v;
return constants<FLOAT128>::one() / v;
}
static inline __float128 abs(__float128 v) noexcept
static inline FLOAT128 abs(FLOAT128 v) noexcept
{
return fabsq(v);
}
static inline __float128 sqrt(__float128 v) noexcept
static inline FLOAT128 sqrt(FLOAT128 v) noexcept
{
return sqrtq(v);
}
static inline __float128 hypot(__float128 v1, __float128 v2) noexcept
static inline FLOAT128 hypot(FLOAT128 v1, FLOAT128 v2) noexcept
{
return hypotq(v1, v2);
}
static inline __float128 exp(__float128 v) noexcept
static inline FLOAT128 exp(FLOAT128 v) noexcept
{
return expq(v);
}
static inline __float128 log(__float128 v) noexcept
static inline FLOAT128 log(FLOAT128 v) noexcept
{
return logq(v);
}
static inline __float128 tanh(__float128 v) noexcept
static inline FLOAT128 tanh(FLOAT128 v) noexcept
{
return tanhq(v);
}
static inline __float128 floor(__float128 v) noexcept
static inline FLOAT128 floor(FLOAT128 v) noexcept
{
return floorq(v);
}
static inline __float128 log1p(__float128 v) noexcept
static inline FLOAT128 log1p(FLOAT128 v) noexcept
{
return log1pq(v);
}
static inline __float128 sin(__float128 v) noexcept
static inline FLOAT128 sin(FLOAT128 v) noexcept
{
return sinq(v);
}
static inline __float128 cos(__float128 v) noexcept
static inline FLOAT128 cos(FLOAT128 v) noexcept
{
return cosq(v);
}
static inline __float128 trunc(__float128 v) noexcept
static inline FLOAT128 trunc(FLOAT128 v) noexcept
{
return truncq(v);
}
template <typename T>
static inline __float128 pow(__float128 v, T p) noexcept
static inline FLOAT128 pow(FLOAT128 v, T p) noexcept
{
return powq(v, static_cast<__float128>(p));
return powq(v, static_cast<FLOAT128>(p));
}
static inline __float128 pow(__float128 v, int p) noexcept
static inline FLOAT128 pow(FLOAT128 v, int p) noexcept
{
if (p==2)
return v*v;
else
return powq(v, static_cast<__float128>(p));
return powq(v, static_cast<FLOAT128>(p));
}
#endif
@ -351,7 +351,7 @@ namespace plib
///
template<typename T>
constexpr
typename std::enable_if<std::is_integral<T>::value && std::is_signed<T>::value, T>::type
typename std::enable_if<plib::is_integral<T>::value && plib::is_signed<T>::value, T>::type
abs(T v) noexcept
{
return v < 0 ? -v : v;
@ -365,7 +365,7 @@ namespace plib
///
template<typename T>
constexpr
typename std::enable_if<std::is_integral<T>::value && std::is_unsigned<T>::value, T>::type
typename std::enable_if<plib::is_integral<T>::value && plib::is_unsigned<T>::value, T>::type
abs(T v) noexcept
{
return v;
@ -386,8 +386,8 @@ namespace plib
constexpr typename std::common_type<M, N>::type
gcd(M m, N n) noexcept
{
static_assert(std::is_integral<M>::value, "gcd: M must be an integer");
static_assert(std::is_integral<N>::value, "gcd: N must be an integer");
static_assert(plib::is_integral<M>::value, "gcd: M must be an integer");
static_assert(plib::is_integral<N>::value, "gcd: N must be an integer");
return m == 0 ? plib::abs(n)
: n == 0 ? plib::abs(m)
@ -409,8 +409,8 @@ namespace plib
constexpr typename std::common_type<M, N>::type
lcm(M m, N n) noexcept
{
static_assert(std::is_integral<M>::value, "lcm: M must be an integer");
static_assert(std::is_integral<N>::value, "lcm: N must be an integer");
static_assert(plib::is_integral<M>::value, "lcm: M must be an integer");
static_assert(plib::is_integral<N>::value, "lcm: N must be an integer");
return (m != 0 && n != 0) ? (plib::abs(m) / gcd(m, n)) * plib::abs(n) : 0;
}

View File

@ -51,7 +51,7 @@ public:
{
return datatype_t(sizeof(T),
plib::is_integral<T>::value || std::is_enum<T>::value,
std::is_floating_point<T>::value);
plib::is_floating_point<T>::value);
}
class callback_t
@ -97,7 +97,7 @@ public:
state_manager_t() = default;
template<typename C> //, typename std::enable_if<std::is_integral<C>::value || std::is_floating_point<C>::value>::type>
template<typename C>
void save_item(const void *owner, C &state, const pstring &stname)
{
save_state_ptr( owner, stname, dtype<C>(), 1, &state);

View File

@ -50,7 +50,7 @@ namespace plib
struct pstonum_helper;
template<typename T>
struct pstonum_helper<T, typename std::enable_if<std::is_integral<T>::value && std::is_signed<T>::value>::type>
struct pstonum_helper<T, typename std::enable_if<plib::is_integral<T>::value && plib::is_signed<T>::value>::type>
{
template <typename S>
long long operator()(std::locale loc, const S &arg, std::size_t *idx)
@ -61,7 +61,7 @@ namespace plib
};
template<typename T>
struct pstonum_helper<T, typename std::enable_if<std::is_integral<T>::value && !std::is_signed<T>::value>::type>
struct pstonum_helper<T, typename std::enable_if<plib::is_integral<T>::value && !plib::is_signed<T>::value>::type>
{
template <typename S>
unsigned long long operator()(std::locale loc, const S &arg, std::size_t *idx)
@ -83,13 +83,13 @@ namespace plib
#if PUSE_FLOAT128
template<>
struct pstonum_helper<__float128>
struct pstonum_helper<FLOAT128>
{
// FIXME: use strtoflt128 from quadmath.h
template <typename S>
__float128 operator()(std::locale loc, const S &arg, std::size_t *idx)
FLOAT128 operator()(std::locale loc, const S &arg, std::size_t *idx)
{
return static_cast<__float128>(pstonum_locale<long double>(loc, arg, idx));
return static_cast<FLOAT128>(pstonum_locale<long double>(loc, arg, idx));
}
};
#endif
@ -101,9 +101,8 @@ namespace plib
std::size_t idx(0);
auto ret = pstonum_helper<T>()(loc, cstr, &idx);
using ret_type = decltype(ret);
if (ret >= static_cast<ret_type>(std::numeric_limits<T>::lowest())
&& ret <= static_cast<ret_type>(std::numeric_limits<T>::max()))
//&& (ret == T(0) || plib::abs(ret) >= std::numeric_limits<T>::min() ))
if (ret >= static_cast<ret_type>(plib::numeric_limits<T>::lowest())
&& ret <= static_cast<ret_type>(plib::numeric_limits<T>::max()))
{
if (cstr[idx] != 0)
throw pexception(pstring("Continuation after numeric value ends: ") + pstring(cstr));

View File

@ -57,14 +57,14 @@ namespace plib
C14CONSTEXPR ptime &operator=(ptime &&rhs) noexcept = default;
constexpr explicit ptime(const double t) = delete;
//: m_time((internal_type) ( t * (double) resolution)) { }
constexpr explicit ptime(const internal_type nom, const internal_type den) noexcept
: m_time(nom * (RES / den)) { }
// FIXME: check for overflow
template <typename O>
constexpr explicit ptime(const ptime<O, RES> &rhs) noexcept
: m_time(rhs.m_time) { }
: m_time(static_cast<TYPE>(rhs.m_time))
{
}
template <typename O>
C14CONSTEXPR ptime &operator+=(const ptime<O, RES> &rhs) noexcept
@ -149,21 +149,13 @@ namespace plib
constexpr internal_type as_raw() const noexcept { return m_time; }
template <typename FT, typename = std::enable_if<std::is_floating_point<FT>::value, FT>>
template <typename FT, typename = std::enable_if<plib::is_floating_point<FT>::value, FT>>
constexpr FT
as_fp() const noexcept
{
return static_cast<FT>(m_time) * inv_res<FT>();
}
#if PUSE_FLOAT128
constexpr __float128
as_fp() const noexcept
{
return static_cast<__float128>(m_time) * inv_res<__float128>();
}
#endif
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>(); }
@ -183,11 +175,7 @@ namespace plib
static constexpr ptime from_raw(internal_type raw) noexcept { return ptime(raw); }
template <typename FT>
static constexpr typename std::enable_if<std::is_floating_point<FT>::value
#if PUSE_FLOAT128
|| std::is_same<FT, __float128>::value
#endif
, ptime>::type
static constexpr typename std::enable_if<plib::is_floating_point<FT>::value, ptime>::type
from_fp(FT t) noexcept { return ptime(static_cast<internal_type>(plib::floor(t * static_cast<FT>(RES) + static_cast<FT>(0.5))), RES); }
static constexpr ptime from_double(double t) noexcept

View File

@ -14,6 +14,10 @@
#include <string>
#include <type_traits>
#if (PUSE_FLOAT128)
#include <quadmath.h>
#endif
// noexcept on move operator -> issue with macosx clang
#define COPYASSIGNMOVE(name, def) \
name(const name &) = def; \
@ -28,24 +32,50 @@
namespace plib
{
template<typename T> struct is_integral : public std::is_integral<T> { };
template<typename T> struct is_signed : public std::is_signed<T> { };
template<typename T> struct is_unsigned : public std::is_unsigned<T> { };
template<typename T> struct numeric_limits : public std::numeric_limits<T> { };
// 128 bit support at least on GCC is not fully supported
#if PHAS_INT128
template<> struct is_integral<UINT128> { static constexpr bool value = true; };
template<> struct is_integral<INT128> { static constexpr bool value = true; };
template<> struct is_signed<UINT128> { static constexpr bool value = false; };
template<> struct is_signed<INT128> { static constexpr bool value = true; };
template<> struct is_unsigned<UINT128> { static constexpr bool value = true; };
template<> struct is_unsigned<INT128> { static constexpr bool value = false; };
template<> struct numeric_limits<UINT128>
{
static constexpr UINT128 max() noexcept
{
return ~((UINT128)0);
return ~(static_cast<UINT128>(0));
}
};
template<> struct numeric_limits<INT128>
{
static constexpr INT128 max() noexcept
{
return (~((UINT128)0)) >> 1;
return (~static_cast<UINT128>(0)) >> 1;
}
};
#endif
template<typename T> struct is_floating_point : public std::is_floating_point<T> { };
#if PUSE_FLOAT128
template<> struct is_floating_point<FLOAT128> { static constexpr bool value = true; };
template<> struct numeric_limits<FLOAT128>
{
static constexpr FLOAT128 max() noexcept
{
return FLT128_MAX;
}
static constexpr FLOAT128 lowest() noexcept
{
return -FLT128_MAX;
}
};
#endif

View File

@ -471,8 +471,13 @@ namespace solver
if (this_resched && !m_Q_sync.net().is_queued())
{
log().warning(MW_NEWTON_LOOPS_EXCEEDED_ON_NET_1(this->name()));
// FIXME: may collide with compute_next_timestep
// FIXME: test and enable - this is working better, though not optimal yet
#if 0
// Don't store, the result can not be used
return netlist_time::from_fp(m_params.m_nr_recalc_delay());
#else
m_Q_sync.net().toggle_and_push_to_queue(netlist_time::from_fp(m_params.m_nr_recalc_delay()));
#endif
}
}
else

View File

@ -51,7 +51,7 @@ namespace solver
FLOAT
, DOUBLE
, LONGDOUBLE
, FLOAT128
, FLOATQ128
)
using static_compile_container = std::vector<std::pair<pstring, pstring>>;
@ -143,10 +143,10 @@ namespace solver
terminal_t **terms() noexcept { return m_terms.data(); }
template <typename FT, typename = std::enable_if<std::is_floating_point<FT>::value, void>>
template <typename FT, typename = std::enable_if<plib::is_floating_point<FT>::value, void>>
FT getV() const noexcept { return static_cast<FT>(m_net->Q_Analog()); }
template <typename FT, typename = std::enable_if<std::is_floating_point<FT>::value, void>>
template <typename FT, typename = std::enable_if<plib::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; }
@ -299,7 +299,7 @@ namespace solver
, m_mat_ptr(size, this->max_railstart() + 1)
, m_last_V(size, nlconst::zero())
, m_DD_n_m_1(size, nlconst::zero())
, m_h_n_m_1(size, nlconst::zero())
, m_h_n_m_1(size, nlconst::magic(1e-6)) // we need a non zero value here
{
//
// save states
@ -318,9 +318,9 @@ namespace solver
static constexpr const std::size_t m_pitch_ABS = (((SIZEABS + 0) + 7) / 8) * 8;
PALIGNAS_VECTOROPT()
plib::parray<FT, SIZE> m_new_V;
plib::parray<float_type, SIZE> m_new_V;
PALIGNAS_VECTOROPT()
plib::parray<FT, SIZE> m_RHS;
plib::parray<float_type, SIZE> m_RHS;
PALIGNAS_VECTOROPT()
plib::pmatrix2d<float_type *> m_mat_ptr;
@ -416,11 +416,11 @@ namespace solver
// and thus belong into a different calculation. This applies to all solvers.
const std::size_t iN = size();
const auto reltol(static_cast<FT>(m_params.m_reltol));
const auto vntol(static_cast<FT>(m_params.m_vntol));
const auto reltol(static_cast<float_type>(m_params.m_reltol));
const auto vntol(static_cast<float_type>(m_params.m_vntol));
for (std::size_t i = 0; i < iN; i++)
{
const auto vold(this->m_terms[i].template getV<FT>());
const auto vold(this->m_terms[i].template getV<float_type>());
const auto vnew(m_new_V[i]);
const auto tol(vntol + reltol * std::max(plib::abs(vnew),plib::abs(vold)));
if (plib::abs(vnew - vold) > tol)
@ -511,6 +511,7 @@ namespace solver
auto &net = m_terms[k];
auto **tcr_r = &(m_mat_ptr[k][0]);
using source_type = typename decltype(m_gtn)::value_type;
const std::size_t term_count = net.count();
const std::size_t railstart = net.railstart();
const auto &go = m_gonn[k];
@ -520,7 +521,7 @@ namespace solver
// FIXME: gonn, gtn and Idr - which float types should they have?
auto gtot_t = std::accumulate(gt, gt + term_count, plib::constants<FT>::zero());
auto gtot_t = std::accumulate(gt, gt + term_count, plib::constants<source_type>::zero());
// update diagonal element ...
*tcr_r[railstart] = static_cast<FT>(gtot_t); //mat.A[mat.diag[k]] += gtot_t;
@ -528,7 +529,7 @@ namespace solver
for (std::size_t i = 0; i < railstart; i++)
*tcr_r[i] += static_cast<FT>(go[i]);
auto RHS_t(std::accumulate(Idr, Idr + term_count, plib::constants<FT>::zero()));
auto RHS_t(std::accumulate(Idr, Idr + term_count, plib::constants<source_type>::zero()));
for (std::size_t i = railstart; i < term_count; i++)
RHS_t += (- go[i]) * *cnV[i];

View File

@ -67,15 +67,17 @@ namespace solver
for (std::size_t i = 0; i < kN; i++)
{
// FIXME: Singular matrix?
const FT f = plib::reciprocal(m_A[i][i]);
const auto &Ai = m_A[i];
const FT f = plib::reciprocal(Ai[i]);
const auto &nzrd = this->m_terms[i].m_nzrd;
const auto &nzbd = this->m_terms[i].m_nzbd;
for (auto &j : nzbd)
{
const FT f1 = -f * m_A[j][i];
auto &Aj = m_A[j];
const FT f1 = -f * Aj[i];
for (auto &k : nzrd)
m_A[j][k] += m_A[i][k] * f1;
Aj[k] += Ai[k] * f1;
this->m_RHS[j] += this->m_RHS[i] * f1;
}
}
@ -95,24 +97,30 @@ namespace solver
if (maxrow != i)
{
#if 0
// Swap the maxrow and ith row
for (std::size_t k = 0; k < kN; k++) {
std::swap(m_A[i][k], m_A[maxrow][k]);
}
#else
std::swap(m_A[i], m_A[maxrow]);
#endif
std::swap(this->m_RHS[i], this->m_RHS[maxrow]);
}
// FIXME: Singular matrix?
const FT f = plib::reciprocal(m_A[i][i]);
const auto &Ai = m_A[i];
const FT f = plib::reciprocal(Ai[i]);
// Eliminate column i from row j
for (std::size_t j = i + 1; j < kN; j++)
{
auto &Aj = m_A[j];
const FT f1 = - m_A[j][i] * f;
if (f1 != plib::constants<FT>::zero())
{
const FT * pi = &(m_A[i][i+1]);
FT * pj = &(m_A[j][i+1]);
const FT * pi = &(Ai[i+1]);
FT * pj = &(Aj[i+1]);
plib::vec_add_mult_scalar_p(kN-i-1,pj,pi,f1);
//for (unsigned k = i+1; k < kN; k++)
// pj[k] = pj[k] + pi[k] * f1;
@ -137,22 +145,26 @@ namespace solver
{
for (std::size_t j = kN; j-- > 0; )
{
FT tmp = 0;
FT tmp(0);
const auto & Aj(m_A[j]);
for (std::size_t k = j+1; k < kN; k++)
tmp += m_A[j][k] * x[k];
x[j] = (this->m_RHS[j] - tmp) / m_A[j][j];
tmp += Aj[k] * x[k];
x[j] = (this->m_RHS[j] - tmp) / Aj[j];
}
}
else
{
for (std::size_t j = kN; j-- > 0; )
{
FT tmp = 0;
FT tmp(0);
const auto &nzrd = this->m_terms[j].m_nzrd;
const auto e = nzrd.size(); // - 1; // exclude RHS element
const auto & Aj(m_A[j]);
const auto e = nzrd.size();
for ( std::size_t k = 0; k < e; k++)
tmp += m_A[j][nzrd[k]] * x[nzrd[k]];
x[j] = (this->m_RHS[j] - tmp) / m_A[j][j];
tmp += Aj[nzrd[k]] * x[nzrd[k]];
x[j] = (this->m_RHS[j] - tmp) / Aj[j];
}
}
}

View File

@ -123,8 +123,7 @@ namespace solver
pstring static_compile_name();
mat_type mat;
plib::dynproc<void, FT *, FT *, FT *, FT *, FT ** > m_proc;
plib::dynproc<void, FT *, nl_fptype *, nl_fptype *, nl_fptype *, nl_fptype ** > m_proc;
};

View File

@ -329,6 +329,7 @@ namespace devices
#if (NL_USE_FLOAT_MATRIX)
ms = create_solvers<float>(sname, grp);
#else
log().info("FPTYPE {1} not supported. Using DOUBLE", m_params.m_fp_type().name());
ms = create_solvers<double>(sname, grp);
#endif
break;
@ -339,13 +340,15 @@ namespace devices
#if (NL_USE_LONG_DOUBLE_MATRIX)
ms = create_solvers<long double>(sname, grp);
#else
log().info("FPTYPE {1} not supported. Using DOUBLE", m_params.m_fp_type().name());
ms = create_solvers<double>(sname, grp);
#endif
break;
case solver::matrix_fp_type_e::FLOAT128:
case solver::matrix_fp_type_e::FLOATQ128:
#if (NL_USE_FLOAT128)
ms = create_solvers<__float128>(sname, grp);
ms = create_solvers<FLOAT128>(sname, grp);
#else
log().info("FPTYPE {1} not supported. Using DOUBLE", m_params.m_fp_type().name());
ms = create_solvers<double>(sname, grp);
#endif
break;