From cabfd7188ae3b4442f6d5b6346709ba3dc3be7f9 Mon Sep 17 00:00:00 2001 From: couriersud Date: Sun, 3 Nov 2019 01:33:48 +0100 Subject: netlist: fix regression, preliminary __float128 support. (nw) __float128 is a gnu extension delivering true 128 bit floating point support. Currently not supported by clang. In addition, the quadmath library needs to be linked. For the time being therefore disabled. --- src/lib/netlist/nl_config.h | 47 +++++++++++++++++++-- src/lib/netlist/plib/gmres.h | 12 +++--- src/lib/netlist/plib/pconfig.h | 9 ++++ src/lib/netlist/plib/pfmtlog.h | 23 ++++++++++- src/lib/netlist/plib/putil.h | 66 ++++++++++++++++++++++++++++++ src/lib/netlist/solver/nld_matrix_solver.h | 21 +++++++--- src/lib/netlist/solver/nld_ms_direct.h | 2 +- src/lib/netlist/solver/nld_ms_sor.h | 2 +- src/lib/netlist/solver/nld_ms_sor_mat.h | 4 +- src/lib/netlist/solver/nld_ms_w.h | 4 +- src/lib/netlist/solver/nld_solver.cpp | 5 +++ 11 files changed, 172 insertions(+), 23 deletions(-) diff --git a/src/lib/netlist/nl_config.h b/src/lib/netlist/nl_config.h index 8d1cfa898f5..9baad43ac44 100644 --- a/src/lib/netlist/nl_config.h +++ b/src/lib/netlist/nl_config.h @@ -57,7 +57,8 @@ #define NL_USE_COPY_INSTEAD_OF_REFERENCE (0) #endif -/* +/*! Use the truthtable implementation of 7448 instead of the coded device + * * FIXME: Using truthtable is a lot slower than the explicit device * in breakout. Performance drops by 20%. This can be fixed by * setting param USE_DEACTIVATE for the device. @@ -67,7 +68,8 @@ #define NL_USE_TRUTHTABLE_7448 (0) #endif -/* +/*! Use the truthtable implementation of 74107 instead of the coded device + * * FIXME: The truthtable implementation of 74107 (JK-Flipflop) * is included for educational purposes to demonstrate how * to implement state holding devices as truthtables. @@ -78,6 +80,15 @@ #define NL_USE_TRUTHTABLE_74107 (0) #endif +/*! Use the __float128 type for matrix calculations. + * + * Defaults to PUSE_FLOAT128 + */ + +#ifndef NL_USE_FLOAT128 +#define NL_USE_FLOAT128 PUSE_FLOAT128 +#endif + //============================================================ // DEBUGGING //============================================================ @@ -104,7 +115,7 @@ static constexpr const auto NETLIST_CLOCK = NETLIST_INTERNAL_RES; // due to numeric issues //============================================================ -using nl_fptype = float; +using nl_fptype = double; namespace netlist { @@ -160,6 +171,36 @@ namespace netlist static inline constexpr const char * name() noexcept { return "float"; } static inline constexpr const char * suffix() noexcept { return "f"; } }; + +#if (NL_USE_FLOAT128) + /*! Specific constants for __float128 floating point type + */ + template <> + 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 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 * 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 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"; } + +#endif + }; +#endif } // namespace netlist #endif /* NLCONFIG_H_ */ diff --git a/src/lib/netlist/plib/gmres.h b/src/lib/netlist/plib/gmres.h index 840da0cfb45..686e1444d36 100644 --- a/src/lib/netlist/plib/gmres.h +++ b/src/lib/netlist/plib/gmres.h @@ -259,7 +259,7 @@ namespace plib m_ht[j][k] = vec_mult(n, m_v[kp1], m_v[j]); vec_add_mult_scalar(n, m_v[kp1], m_v[j], -m_ht[j][k]); } - m_ht[kp1][k] = std::sqrt(vec_mult2(n, m_v[kp1])); + m_ht[kp1][k] = plib::sqrt(vec_mult2(n, m_v[kp1])); // FIXME: comparison to zero if (m_ht[kp1][k] != plib::constants::zero()) @@ -268,7 +268,7 @@ namespace plib 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 = reciprocal(std::hypot(m_ht[k][k], m_ht[kp1][k])); + const float_type mu = reciprocal(plib::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; @@ -277,7 +277,7 @@ namespace plib givens_mult(m_c[k], m_s[k], m_g[k], m_g[kp1]); - FT rho = std::abs(m_g[kp1]); + FT rho = plib::abs(m_g[kp1]); // FIXME .. itr_used = itr_used + 1; @@ -359,12 +359,12 @@ namespace plib ops.solve_inplace(Ax); - const float_type rho_to_accuracy = std::sqrt(vec_mult2(n, Ax)) / accuracy; + const float_type rho_to_accuracy = plib::sqrt(vec_mult2(n, Ax)) / accuracy; rho_delta = accuracy * rho_to_accuracy; } else - rho_delta = accuracy * std::sqrt(static_cast(n)); + rho_delta = accuracy * plib::sqrt(static_cast(n)); /* * Using @@ -389,7 +389,7 @@ namespace plib ops.solve_inplace(residual); - rho = std::sqrt(vec_mult2(n, residual)); + rho = plib::sqrt(vec_mult2(n, residual)); if (rho < rho_delta) return itr_used + 1; diff --git a/src/lib/netlist/plib/pconfig.h b/src/lib/netlist/plib/pconfig.h index 102a991e613..58a9387c51c 100644 --- a/src/lib/netlist/plib/pconfig.h +++ b/src/lib/netlist/plib/pconfig.h @@ -33,6 +33,15 @@ #define PHAS_INT128 (0) #endif +/*! Add support for the __float128 floating point type. + * + * + */ + +#ifndef PUSE_FLOAT128 +#define PUSE_FLOAT128 (0) +#endif + /* * OpenMP adds about 10% to 20% performance for analog * netlists like kidniki. diff --git a/src/lib/netlist/plib/pfmtlog.h b/src/lib/netlist/plib/pfmtlog.h index 53ed539078e..fd1f073cb73 100644 --- a/src/lib/netlist/plib/pfmtlog.h +++ b/src/lib/netlist/plib/pfmtlog.h @@ -28,11 +28,22 @@ P_ENUM(plog_level, template struct ptype_traits_base { - static const T cast(const T &x) { return x; } - static const bool is_signed = std::numeric_limits::is_signed; + static constexpr const T cast(const T &x) { return x; } + static constexpr const bool is_signed = std::numeric_limits::is_signed; static char32_t fmt_spec() { return 'u'; } }; +#if (PUSE_FLOAT128) +template <> +struct ptype_traits_base<__float128> +{ + // FIXME: need native support at some time + static constexpr const long double cast(const __float128 &x) { return static_cast(x); } + static constexpr const bool is_signed = true; + static char32_t fmt_spec() { return 'f'; } +}; +#endif + template <> struct ptype_traits_base { @@ -133,6 +144,14 @@ struct ptype_traits : ptype_traits_base static char32_t fmt_spec() { return 'f'; } }; +#if (PUSE_FLOAT128) +template<> +struct ptype_traits<__float128> : ptype_traits_base<__float128> +{ + static char32_t fmt_spec() { return 'f'; } +}; +#endif + template<> struct ptype_traits : ptype_traits_base { diff --git a/src/lib/netlist/plib/putil.h b/src/lib/netlist/plib/putil.h index e72eccfe01f..e575a9adc73 100644 --- a/src/lib/netlist/plib/putil.h +++ b/src/lib/netlist/plib/putil.h @@ -18,6 +18,11 @@ #include #include #include +#include + +#if (PUSE_FLOAT128) +#include +#endif #define PSTRINGIFY_HELP(y) # y #define PSTRINGIFY(x) PSTRINGIFY_HELP(x) @@ -281,6 +286,67 @@ namespace plib return constants::one() / v; } + /*! abs function + * + * @tparam T type of the argument + * @param v argument + * @return absolute value of argument + */ + template + static inline constexpr typename std::enable_if::value, T>::type + abs(T v) noexcept + { + return std::abs(v); + } + + /*! sqrt function + * + * @tparam T type of the argument + * @param v argument + * @return absolute value of argument + */ + template + static inline constexpr typename std::enable_if::value, T>::type + sqrt(T v) noexcept + { + return std::sqrt(v); + } + + /*! hypot function + * + * @tparam T type of the argument + * @param v argument + * @return absolute value of argument + */ + template + static inline constexpr typename std::enable_if::value, T>::type + hypot(T v1, T v2) noexcept + { + return std::hypot(v1, v2); + } + +#if (PUSE_FLOAT128) + static inline constexpr __float128 reciprocal(__float128 v) noexcept + { + return constants<__float128>::one() / v; + } + + static inline __float128 abs(__float128 v) noexcept + { + return fabsq(v); + } + + static inline __float128 sqrt(__float128 v) noexcept + { + return sqrtq(v); + } + + static inline __float128 hypot(__float128 v1, __float128 v2) noexcept + { + return hypotq(v1, v2); + } +#endif + static_assert(noexcept(constants::one()) == true, "Not evaluated as constexpr"); diff --git a/src/lib/netlist/solver/nld_matrix_solver.h b/src/lib/netlist/solver/nld_matrix_solver.h index 09e650be1a3..26f6ea149a7 100644 --- a/src/lib/netlist/solver/nld_matrix_solver.h +++ b/src/lib/netlist/solver/nld_matrix_solver.h @@ -40,11 +40,20 @@ namespace solver GMRES ) +#if (NL_USE_FLOAT128) P_ENUM(matrix_fp_type_e, - FLOAT, - DOUBLE, - LONGDOUBLE + FLOAT + , DOUBLE + , LONGDOUBLE + , FLOAT128 ) +#else + P_ENUM(matrix_fp_type_e, + FLOAT + , DOUBLE + , LONGDOUBLE + ) +#endif struct solver_parameters_t { @@ -56,7 +65,7 @@ namespace solver , m_method(parent, "METHOD", matrix_type_e::MAT_CR) , m_fp_type(parent, "FPTYPE", matrix_fp_type_e::DOUBLE) , m_reltol(parent, "RELTOL", nlconst::magic(1e-3)) ///< SPICE RELTOL parameter - , m_vntol(parent, "VNTOL", nlconst::magic(1e-6)) ///< SPICE VNTOL parameter + , m_vntol(parent, "VNTOL", nlconst::magic(1e-7)) ///< SPICE VNTOL parameter , m_accuracy(parent, "ACCURACY", nlconst::magic(1e-7)) ///< Iterative solver accuracy , m_nr_loops(parent, "NR_LOOPS", 250) ///< Maximum number of Newton-Raphson loops , m_gs_loops(parent, "GS_LOOPS", 9) ///< Maximum number of Gauss-Seidel loops @@ -438,8 +447,8 @@ namespace solver { const auto vold(this->m_terms[i].template getV()); const auto vnew(m_new_V[i]); - const auto tol(vntol + reltol * std::max(std::abs(vnew),std::abs(vold))); - if (std::abs(vnew - vold) > tol) + const auto tol(vntol + reltol * std::max(plib::abs(vnew),plib::abs(vold))); + if (plib::abs(vnew - vold) > tol) return true; } return false; diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h index 20f9fa23964..031cae34a27 100644 --- a/src/lib/netlist/solver/nld_ms_direct.h +++ b/src/lib/netlist/solver/nld_ms_direct.h @@ -89,7 +89,7 @@ namespace solver std::size_t maxrow = i; for (std::size_t j = i + 1; j < kN; j++) { - if (std::abs(m_A[j][i]) > std::abs(m_A[maxrow][i])) + if (plib::abs(m_A[j][i]) > plib::abs(m_A[maxrow][i])) //if (m_A[j][i] * m_A[j][i] > m_A[maxrow][i] * m_A[maxrow][i]) maxrow = j; } diff --git a/src/lib/netlist/solver/nld_ms_sor.h b/src/lib/netlist/solver/nld_ms_sor.h index ee79599a2ff..56075fa090e 100644 --- a/src/lib/netlist/solver/nld_ms_sor.h +++ b/src/lib/netlist/solver/nld_ms_sor.h @@ -134,7 +134,7 @@ namespace solver const float_type new_val = this->m_new_V[k] * one_m_w[k] + (Idrive + this->m_RHS[k]) * w[k]; - err = std::max(std::abs(new_val - this->m_new_V[k]), err); + err = std::max(plib::abs(new_val - this->m_new_V[k]), err); this->m_new_V[k] = new_val; } diff --git a/src/lib/netlist/solver/nld_ms_sor_mat.h b/src/lib/netlist/solver/nld_ms_sor_mat.h index e411b1f1e59..a5ded1a38ba 100644 --- a/src/lib/netlist/solver/nld_ms_sor_mat.h +++ b/src/lib/netlist/solver/nld_ms_sor_mat.h @@ -180,7 +180,7 @@ namespace solver 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]]); + gabs_t = gabs_t + plib::abs(this->m_A[k][p[i]]); gabs_t *= plib::constants::one(); // derived by try and error if (gabs_t > this->m_A[k][k]) @@ -190,7 +190,7 @@ namespace solver } const float_type delta = w * (this->m_RHS[k] - Idrive) ; - cerr = std::max(cerr, std::abs(delta)); + cerr = std::max(cerr, plib::abs(delta)); this->m_new_V[k] += delta; } diff --git a/src/lib/netlist/solver/nld_ms_w.h b/src/lib/netlist/solver/nld_ms_w.h index b0c148bf252..5e47e93d171 100644 --- a/src/lib/netlist/solver/nld_ms_w.h +++ b/src/lib/netlist/solver/nld_ms_w.h @@ -340,8 +340,8 @@ namespace solver { tmp += A(i,j) * this->m_new_V[j]; } - 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)); + if (plib::abs(tmp-RHS(i)) > static_cast(1e-6)) + plib::perrlogger("{} failed on row {}: {} RHS: {}\n", this->name(), i, plib::abs(tmp-RHS(i)), RHS(i)); } bool err(false); diff --git a/src/lib/netlist/solver/nld_solver.cpp b/src/lib/netlist/solver/nld_solver.cpp index 86772fac6a0..ecbf48ab0c2 100644 --- a/src/lib/netlist/solver/nld_solver.cpp +++ b/src/lib/netlist/solver/nld_solver.cpp @@ -343,6 +343,11 @@ namespace devices case solver::matrix_fp_type_e::LONGDOUBLE: ms = create_solvers(sname, grp); break; +#if (NL_USE_FLOAT128) + case solver::matrix_fp_type_e::FLOAT128: + ms = create_solvers<__float128>(sname, grp); + break; +#endif } log().verbose("Solver {1}", ms->name()); -- cgit v1.2.3