diff options
Diffstat (limited to 'src/lib/netlist/solver/nld_ms_direct.h')
-rw-r--r-- | src/lib/netlist/solver/nld_ms_direct.h | 153 |
1 files changed, 61 insertions, 92 deletions
diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h index 398e20995f8..2c56e9fecb4 100644 --- a/src/lib/netlist/solver/nld_ms_direct.h +++ b/src/lib/netlist/solver/nld_ms_direct.h @@ -14,28 +14,21 @@ #include "nld_solver.h" #include "nld_matrix_solver.h" #include "vector_base.h" - -/* Disabling dynamic allocation gives a ~10% boost in performance - * This flag has been added to support continuous storage for arrays - * going forward in case we implement cuda solvers in the future. - */ -#define NL_USE_DYNAMIC_ALLOCATION (1) +#include "mat_cr.h" namespace netlist { - namespace devices - { -//#define nl_ext_double _float128 // slow, very slow -//#define nl_ext_double long double // slightly slower -#define nl_ext_double nl_double - +namespace devices +{ -template <std::size_t m_N, std::size_t storage_N> +template <typename FT, int SIZE> class matrix_solver_direct_t: public matrix_solver_t { friend class matrix_solver_t; public: + typedef FT float_type; + matrix_solver_direct_t(netlist_base_t &anetlist, const pstring &name, const solver_parameters_t *params, const std::size_t size); matrix_solver_direct_t(netlist_base_t &anetlist, const pstring &name, const eSortType sort, const solver_parameters_t *params, const std::size_t size); @@ -48,36 +41,25 @@ protected: virtual unsigned vsolve_non_dynamic(const bool newton_raphson) override; unsigned solve_non_dynamic(const bool newton_raphson); - constexpr std::size_t N() const { return (m_N == 0) ? m_dim : m_N; } + constexpr std::size_t N() const { return (SIZE > 0) ? static_cast<std::size_t>(SIZE) : m_dim; } void LE_solve(); template <typename T> - void LE_back_subst(T * RESTRICT x); + void LE_back_subst(T & RESTRICT x); -#if (NL_USE_DYNAMIC_ALLOCATION) - nl_ext_double &A(const std::size_t r, const std::size_t c) { return m_A[r * m_pitch + c]; } - nl_ext_double &RHS(const std::size_t r) { return m_A[r * m_pitch + N()]; } -#else - nl_ext_double &A(const std::size_t r, const std::size_t c) { return m_A[r][c]; } - nl_ext_double &RHS(const std::size_t r) { return m_A[r][N()]; } -#endif - nl_double m_last_RHS[storage_N]; // right hand side - contains currents + FT &A(std::size_t r, std::size_t c) { return m_A[r * m_pitch + c]; } + FT &RHS(std::size_t r) { return m_A[r * m_pitch + N()]; } + plib::parray<FT, SIZE> m_last_RHS; + plib::parray<FT, SIZE> m_new_V; private: - //static const std::size_t m_pitch = (((storage_N + 1) + 0) / 1) * 1; - static constexpr std::size_t m_pitch = (((storage_N + 1) + 7) / 8) * 8; - //static const std::size_t m_pitch = (((storage_N + 1) + 15) / 16) * 16; - //static const std::size_t m_pitch = (((storage_N + 1) + 31) / 32) * 32; -#if (NL_USE_DYNAMIC_ALLOCATION) - //nl_ext_double * RESTRICT m_A; - std::vector<nl_ext_double> m_A; -#else - nl_ext_double m_A[storage_N][m_pitch]; -#endif - //nl_ext_double m_RHSx[storage_N]; + static constexpr const std::size_t SIZEABS = plib::parray<FT, SIZE>::SIZEABS(); + static constexpr const std::size_t m_pitch_ABS = (((SIZEABS + 1) + 7) / 8) * 8; const std::size_t m_dim; + const std::size_t m_pitch; + plib::parray<FT, SIZE * m_pitch_ABS> m_A; }; @@ -85,16 +67,13 @@ private: // matrix_solver_direct // ---------------------------------------------------------------------------------------- -template <std::size_t m_N, std::size_t storage_N> -matrix_solver_direct_t<m_N, storage_N>::~matrix_solver_direct_t() +template <typename FT, int SIZE> +matrix_solver_direct_t<FT, SIZE>::~matrix_solver_direct_t() { -#if (NL_USE_DYNAMIC_ALLOCATION) - //plib::pfree_array(m_A); -#endif } -template <std::size_t m_N, std::size_t storage_N> -void matrix_solver_direct_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets) +template <typename FT, int SIZE> +void matrix_solver_direct_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets) { matrix_solver_t::setup_base(nets); @@ -107,34 +86,31 @@ void matrix_solver_direct_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets) t->m_nzrd.push_back(static_cast<unsigned>(N())); } - state().save(*this, m_last_RHS, "m_last_RHS"); + state().save(*this, m_last_RHS.data(), "m_last_RHS", m_last_RHS.size()); for (std::size_t k = 0; k < N(); k++) state().save(*this, RHS(k), plib::pfmt("RHS.{1}")(k)); } -template <std::size_t m_N, std::size_t storage_N> -void matrix_solver_direct_t<m_N, storage_N>::LE_solve() +template <typename FT, int SIZE> +void matrix_solver_direct_t<FT, SIZE>::LE_solve() { const std::size_t kN = N(); if (!m_params.m_pivot) { for (std::size_t i = 0; i < kN; i++) { - /* FIXME: Singular matrix? */ - nl_double *Ai = &A(i, 0); - const nl_double f = 1.0 / A(i,i); + const FT f = 1.0 / A(i,i); const auto &nzrd = m_terms[i]->m_nzrd; const auto &nzbd = m_terms[i]->m_nzbd; for (std::size_t j : nzbd) { - nl_double *Aj = &A(j, 0); - const nl_double f1 = -f * Aj[i]; + const FT f1 = -f * A(j, i); for (std::size_t k : nzrd) - Aj[k] += Ai[k] * f1; + A(j, k) += A(i, k) * f1; //RHS(j) += RHS(i) * f1; } } @@ -161,17 +137,17 @@ void matrix_solver_direct_t<m_N, storage_N>::LE_solve() //std::swap(RHS(i), RHS(maxrow)); } /* FIXME: Singular matrix? */ - const nl_double f = 1.0 / A(i,i); + const FT f = 1.0 / A(i,i); /* Eliminate column i from row j */ for (std::size_t j = i + 1; j < kN; j++) { - const nl_double f1 = - A(j,i) * f; + const FT f1 = - A(j,i) * f; if (f1 != NL_FCONST(0.0)) { - const nl_double * RESTRICT pi = &A(i,i+1); - nl_double * RESTRICT pj = &A(j,i+1); + const FT * RESTRICT pi = &A(i,i+1); + FT * RESTRICT pj = &A(j,i+1); #if 1 vec_add_mult_scalar_p(kN-i,pi,f1,pj); #else @@ -188,10 +164,10 @@ void matrix_solver_direct_t<m_N, storage_N>::LE_solve() } } -template <std::size_t m_N, std::size_t storage_N> +template <typename FT, int SIZE> template <typename T> -void matrix_solver_direct_t<m_N, storage_N>::LE_back_subst( - T * RESTRICT x) +void matrix_solver_direct_t<FT, SIZE>::LE_back_subst( + T & RESTRICT x) { const std::size_t kN = N(); @@ -200,7 +176,7 @@ void matrix_solver_direct_t<m_N, storage_N>::LE_back_subst( { for (std::size_t j = kN; j-- > 0; ) { - T tmp = 0; + FT tmp = 0; for (std::size_t k = j+1; k < kN; k++) tmp += A(j,k) * x[k]; x[j] = (RHS(j) - tmp) / A(j,j); @@ -210,40 +186,33 @@ void matrix_solver_direct_t<m_N, storage_N>::LE_back_subst( { for (std::size_t j = kN; j-- > 0; ) { - T tmp = 0; - - const auto *p = m_terms[j]->m_nzrd.data(); - const auto e = m_terms[j]->m_nzrd.size() - 1; /* exclude RHS element */ - T * Aj = &A(j,0); - for (std::size_t k = 0; k < e; k++) - { - const auto pk = p[k]; - tmp += Aj[pk] * x[pk]; - } + FT tmp = 0; + const auto &nzrd = m_terms[j]->m_nzrd; + const auto e = nzrd.size() - 1; /* exclude RHS element */ + for ( std::size_t k = 0; k < e; k++) + tmp += A(j, nzrd[k]) * x[nzrd[k]]; x[j] = (RHS(j) - tmp) / A(j,j); } } } -template <std::size_t m_N, std::size_t storage_N> -unsigned matrix_solver_direct_t<m_N, storage_N>::solve_non_dynamic(const bool newton_raphson) +template <typename FT, int SIZE> +unsigned matrix_solver_direct_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphson) { - nl_double new_V[storage_N]; // = { 0.0 }; - this->LE_solve(); - this->LE_back_subst(new_V); + this->LE_back_subst(m_new_V); - const nl_double err = (newton_raphson ? delta(new_V) : 0.0); - store(new_V); + const FT err = (newton_raphson ? delta(m_new_V) : 0.0); + store(m_new_V); return (err > this->m_params.m_accuracy) ? 2 : 1; } -template <std::size_t m_N, std::size_t storage_N> -unsigned matrix_solver_direct_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson) +template <typename FT, int SIZE> +unsigned matrix_solver_direct_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson) { - build_LE_A<matrix_solver_direct_t>(); - build_LE_RHS<matrix_solver_direct_t>(); + this->build_LE_A(*this); + this->build_LE_RHS(*this); for (std::size_t i=0, iN=N(); i < iN; i++) m_last_RHS[i] = RHS(i); @@ -252,33 +221,33 @@ unsigned matrix_solver_direct_t<m_N, storage_N>::vsolve_non_dynamic(const bool n return this->solve_non_dynamic(newton_raphson); } -template <std::size_t m_N, std::size_t storage_N> -matrix_solver_direct_t<m_N, storage_N>::matrix_solver_direct_t(netlist_base_t &anetlist, const pstring &name, +template <typename FT, int SIZE> +matrix_solver_direct_t<FT, SIZE>::matrix_solver_direct_t(netlist_base_t &anetlist, const pstring &name, const solver_parameters_t *params, const std::size_t size) : matrix_solver_t(anetlist, name, ASCENDING, params) +, m_last_RHS(size) +, m_new_V(size) , m_dim(size) +, m_pitch(m_pitch_ABS ? m_pitch_ABS : (((m_dim + 1) + 7) / 8) * 8) +, m_A(size * m_pitch) { -#if (NL_USE_DYNAMIC_ALLOCATION) - m_A.resize(N() * m_pitch); - //m_A = plib::palloc_array<nl_ext_double>(N() * m_pitch); -#endif for (unsigned k = 0; k < N(); k++) { m_last_RHS[k] = 0.0; } } -template <std::size_t m_N, std::size_t storage_N> -matrix_solver_direct_t<m_N, storage_N>::matrix_solver_direct_t(netlist_base_t &anetlist, const pstring &name, +template <typename FT, int SIZE> +matrix_solver_direct_t<FT, SIZE>::matrix_solver_direct_t(netlist_base_t &anetlist, const pstring &name, const eSortType sort, const solver_parameters_t *params, const std::size_t size) : matrix_solver_t(anetlist, name, sort, params) +, m_last_RHS(size) +, m_new_V(size) , m_dim(size) +, m_pitch(m_pitch_ABS ? m_pitch_ABS : (((m_dim + 1) + 7) / 8) * 8) +, m_A(size * m_pitch) { -#if (NL_USE_DYNAMIC_ALLOCATION) - m_A.resize(N() * m_pitch); - //m_A = plib::palloc_array<nl_ext_double>(N() * m_pitch); -#endif - for (unsigned k = 0; k < N(); k++) + for (std::size_t k = 0; k < N(); k++) { m_last_RHS[k] = 0.0; } |