diff options
Diffstat (limited to 'src/lib/netlist/solver/nld_ms_direct.h')
-rw-r--r-- | src/lib/netlist/solver/nld_ms_direct.h | 386 |
1 files changed, 222 insertions, 164 deletions
diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h index 2501742218d..f03cb738028 100644 --- a/src/lib/netlist/solver/nld_ms_direct.h +++ b/src/lib/netlist/solver/nld_ms_direct.h @@ -8,228 +8,286 @@ #ifndef NLD_MS_DIRECT_H_ #define NLD_MS_DIRECT_H_ -#include "nld_matrix_solver.h" +#include <algorithm> + #include "nld_solver.h" -#include "plib/mat_cr.h" -#include "plib/vector_ops.h" +#include "nld_matrix_solver.h" +#include "vector_base.h" -#include <algorithm> -#include <cmath> +/* 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) namespace netlist { -namespace devices -{ - - template <typename FT, int SIZE> - class matrix_solver_direct_t: public matrix_solver_t + namespace devices { - friend class matrix_solver_t; - public: - - using float_type = FT; +//#define nl_ext_double _float128 // slow, very slow +//#define nl_ext_double long double // slightly slower +#define nl_ext_double nl_double - matrix_solver_direct_t(netlist_state_t &anetlist, const pstring &name, const solver_parameters_t *params, const std::size_t size); - matrix_solver_direct_t(netlist_state_t &anetlist, const pstring &name, const eSortType sort, const solver_parameters_t *params, const std::size_t size); - void vsetup(analog_net_t::list_t &nets) override; - void reset() override { matrix_solver_t::reset(); } - - protected: - unsigned vsolve_non_dynamic(const bool newton_raphson) override; - unsigned solve_non_dynamic(const bool newton_raphson); +template <std::size_t m_N, std::size_t storage_N> +class matrix_solver_direct_t: public matrix_solver_t +{ + friend class matrix_solver_t; +public: - constexpr std::size_t size() const { return (SIZE > 0) ? static_cast<std::size_t>(SIZE) : m_dim; } + matrix_solver_direct_t(netlist_t &anetlist, const pstring &name, const solver_parameters_t *params, const std::size_t size); + matrix_solver_direct_t(netlist_t &anetlist, const pstring &name, const eSortType sort, const solver_parameters_t *params, const std::size_t size); - void LE_solve(); + virtual ~matrix_solver_direct_t() override; - template <typename T> - void LE_back_subst(T & x); + virtual void vsetup(analog_net_t::list_t &nets) override; + virtual void reset() override { matrix_solver_t::reset(); } - 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 + size()]; } - plib::parray<FT, SIZE> m_new_V; +protected: + virtual unsigned vsolve_non_dynamic(const bool newton_raphson) override; + unsigned solve_non_dynamic(const bool newton_raphson); - private: - 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; + constexpr std::size_t N() const { return (m_N == 0) ? m_dim : m_N; } - const std::size_t m_dim; - const std::size_t m_pitch; - plib::parray<FT, SIZE * int(m_pitch_ABS)> m_A; + void LE_solve(); - }; + template <typename T> + void LE_back_subst(T * RESTRICT x); + +#if (NL_USE_DYNAMIC_ALLOCATION) + template <typename T1, typename T2> + nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r * m_pitch + c]; } + template <typename T1> + nl_ext_double &RHS(const T1 &r) { return m_A[r * m_pitch + N()]; } +#else + template <typename T1, typename T2> + nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r][c]; } + template <typename T1> + nl_ext_double &RHS(const T1 &r) { return m_A[r][N()]; } +#endif + nl_double m_last_RHS[storage_N]; // right hand side - contains currents + +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]; + + const std::size_t m_dim; + +}; + +// ---------------------------------------------------------------------------------------- +// 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() +{ +#if (NL_USE_DYNAMIC_ALLOCATION) + //plib::pfree_array(m_A); +#endif +} - // ---------------------------------------------------------------------------------------- - // matrix_solver_direct - // ---------------------------------------------------------------------------------------- +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) +{ + matrix_solver_t::setup_base(nets); - template <typename FT, int SIZE> - void matrix_solver_direct_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets) + /* add RHS element */ + for (std::size_t k = 0; k < N(); k++) { - matrix_solver_t::setup_base(nets); + terms_for_net_t * t = m_terms[k].get(); - /* add RHS element */ - for (std::size_t k = 0; k < size(); k++) - { - terms_for_net_t * t = m_terms[k].get(); + if (!plib::container::contains(t->m_nzrd, static_cast<unsigned>(N()))) + t->m_nzrd.push_back(static_cast<unsigned>(N())); + } - if (!plib::container::contains(t->m_nzrd, static_cast<unsigned>(size()))) - t->m_nzrd.push_back(static_cast<unsigned>(size())); - } + netlist().save(*this, m_last_RHS, "m_last_RHS"); - // FIXME: This shouldn't be necessary ... - for (std::size_t k = 0; k < size(); k++) - state().save(*this, RHS(k), this->name(), plib::pfmt("RHS.{1}")(k)); - } + for (std::size_t k = 0; k < N(); k++) + netlist().save(*this, RHS(k), plib::pfmt("RHS.{1}")(k)); +} - template <typename FT, int SIZE> - void matrix_solver_direct_t<FT, SIZE>::LE_solve() + +template <std::size_t m_N, std::size_t storage_N> +void matrix_solver_direct_t<m_N, storage_N>::LE_solve() +{ + const std::size_t kN = N(); + if (!m_params.m_pivot) { - const std::size_t kN = size(); - if (!m_params.m_pivot) + for (std::size_t i = 0; i < kN; i++) { - for (std::size_t i = 0; i < kN; i++) - { - /* FIXME: Singular matrix? */ - 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) - { - const FT f1 = -f * A(j, i); - for (std::size_t k : nzrd) - A(j, k) += A(i, k) * f1; - //RHS(j) += RHS(i) * f1; - } + /* FIXME: Singular matrix? */ + nl_double *Ai = &A(i, 0); + const nl_double 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]; + for (std::size_t k : nzrd) + Aj[k] += Ai[k] * f1; + //RHS(j) += RHS(i) * f1; } } - else + } + else + { + for (std::size_t i = 0; i < kN; i++) { - for (std::size_t i = 0; i < kN; i++) + /* Find the row with the largest first value */ + std::size_t maxrow = i; + for (std::size_t j = i + 1; j < kN; j++) { - /* Find the row with the largest first value */ - 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 (A(j,i) * A(j,i) > A(maxrow,i) * A(maxrow,i)) - maxrow = j; - } + //if (std::abs(m_A[j][i]) > std::abs(m_A[maxrow][i])) + if (A(j,i) * A(j,i) > A(maxrow,i) * A(maxrow,i)) + maxrow = j; + } - if (maxrow != i) - { - /* Swap the maxrow and ith row */ - for (std::size_t k = 0; k < kN + 1; k++) { - std::swap(A(i,k), A(maxrow,k)); - } - //std::swap(RHS(i), RHS(maxrow)); + if (maxrow != i) + { + /* Swap the maxrow and ith row */ + for (std::size_t k = 0; k < kN + 1; k++) { + std::swap(A(i,k), A(maxrow,k)); } - /* FIXME: Singular matrix? */ - const FT f = 1.0 / A(i,i); + //std::swap(RHS(i), RHS(maxrow)); + } + /* FIXME: Singular matrix? */ + const nl_double f = 1.0 / A(i,i); - /* Eliminate column i from row j */ + /* Eliminate column i from row j */ - for (std::size_t j = i + 1; j < kN; j++) + for (std::size_t j = i + 1; j < kN; j++) + { + const nl_double f1 = - A(j,i) * f; + if (f1 != NL_FCONST(0.0)) { - const FT f1 = - A(j,i) * f; - if (f1 != plib::constants<FT>::zero()) - { - const FT * pi = &A(i,i+1); - FT * pj = &A(j,i+1); - #if 1 - plib::vec_add_mult_scalar_p(kN-i,pj, pi,f1); - #else - vec_add_mult_scalar_p1(kN-i-1,pj,pi,f1); - //for (unsigned k = i+1; k < kN; k++) - // pj[k] = pj[k] + pi[k] * f1; - //for (unsigned k = i+1; k < kN; k++) - //A(j,k) += A(i,k) * f1; - RHS(j) += RHS(i) * f1; - #endif - } + const nl_double * RESTRICT pi = &A(i,i+1); + nl_double * RESTRICT pj = &A(j,i+1); +#if 1 + vec_add_mult_scalar_p(kN-i,pi,f1,pj); +#else + vec_add_mult_scalar_p(kN-i-1,pj,f1,pi); + //for (unsigned k = i+1; k < kN; k++) + // pj[k] = pj[k] + pi[k] * f1; + //for (unsigned k = i+1; k < kN; k++) + //A(j,k) += A(i,k) * f1; + RHS(j) += RHS(i) * f1; +#endif } } } } +} - template <typename FT, int SIZE> - template <typename T> - void matrix_solver_direct_t<FT, SIZE>::LE_back_subst( - T & x) - { - const std::size_t kN = size(); +template <std::size_t m_N, std::size_t storage_N> +template <typename T> +void matrix_solver_direct_t<m_N, storage_N>::LE_back_subst( + T * RESTRICT x) +{ + const std::size_t kN = N(); - /* back substitution */ - if (m_params.m_pivot) + /* back substitution */ + if (m_params.m_pivot) + { + for (std::size_t j = kN; j-- > 0; ) { - for (std::size_t j = kN; j-- > 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); - } + T 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); } - else + } + else + { + for (std::size_t j = kN; j-- > 0; ) { - 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++) { - 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); + const auto pk = p[k]; + tmp += Aj[pk] * x[pk]; } + x[j] = (RHS(j) - tmp) / A(j,j); } } +} - template <typename FT, int SIZE> - unsigned matrix_solver_direct_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphson) - { - this->LE_solve(); - this->LE_back_subst(m_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>::solve_non_dynamic(const bool newton_raphson) +{ + nl_double new_V[storage_N]; // = { 0.0 }; - template <typename FT, int SIZE> - unsigned matrix_solver_direct_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson) - { - this->build_LE_A(*this); - this->build_LE_RHS(*this); + this->LE_solve(); + this->LE_back_subst(new_V); - this->m_stat_calculations++; - return this->solve_non_dynamic(newton_raphson); - } + const nl_double err = (newton_raphson ? delta(new_V) : 0.0); + store(new_V); + return (err > this->m_params.m_accuracy) ? 2 : 1; +} - template <typename FT, int SIZE> - matrix_solver_direct_t<FT, SIZE>::matrix_solver_direct_t(netlist_state_t &anetlist, const pstring &name, - const solver_parameters_t *params, const std::size_t size) - : matrix_solver_t(anetlist, name, ASCENDING, params) - , 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) +template <std::size_t m_N, std::size_t storage_N> +inline unsigned matrix_solver_direct_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson) +{ + build_LE_A<matrix_solver_direct_t>(); + build_LE_RHS<matrix_solver_direct_t>(); + + for (std::size_t i=0, iN=N(); i < iN; i++) + m_last_RHS[i] = RHS(i); + + this->m_stat_calculations++; + 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_t &anetlist, const pstring &name, + const solver_parameters_t *params, const std::size_t size) +: matrix_solver_t(anetlist, name, ASCENDING, params) +, m_dim(size) +{ +#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 <typename FT, int SIZE> - matrix_solver_direct_t<FT, SIZE>::matrix_solver_direct_t(netlist_state_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_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) +template <std::size_t m_N, std::size_t storage_N> +matrix_solver_direct_t<m_N, storage_N>::matrix_solver_direct_t(netlist_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_dim(size) +{ +#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; } +} -} // namespace devices + } //namespace devices } // namespace netlist #endif /* NLD_MS_DIRECT_H_ */ |