diff options
Diffstat (limited to 'src/lib/netlist/solver/nld_ms_sm.h')
-rw-r--r-- | src/lib/netlist/solver/nld_ms_sm.h | 416 |
1 files changed, 199 insertions, 217 deletions
diff --git a/src/lib/netlist/solver/nld_ms_sm.h b/src/lib/netlist/solver/nld_ms_sm.h index ed95cd2dc03..d85ab0044f4 100644 --- a/src/lib/netlist/solver/nld_ms_sm.h +++ b/src/lib/netlist/solver/nld_ms_sm.h @@ -33,293 +33,275 @@ #ifndef NLD_MS_SM_H_ #define NLD_MS_SM_H_ -#include <algorithm> - -#include "nld_solver.h" #include "nld_matrix_solver.h" -#include "vector_base.h" +#include "nld_solver.h" +#include "plib/vector_ops.h" + +#include <algorithm> 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 - -template <std::size_t m_N, std::size_t storage_N> -class matrix_solver_sm_t: public matrix_solver_t +namespace devices { - friend class matrix_solver_t; - -public: - matrix_solver_sm_t(netlist_t &anetlist, const pstring &name, - const solver_parameters_t *params, const std::size_t size); - - virtual ~matrix_solver_sm_t() override; - - virtual void vsetup(analog_net_t::list_t &nets) override; - virtual void reset() override { matrix_solver_t::reset(); } - -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; } - - void LE_invert(); - - template <typename T> - void LE_compute_x(T * RESTRICT x); + template <typename FT, int SIZE> + class matrix_solver_sm_t: public matrix_solver_t + { + friend class matrix_solver_t; + public: - template <typename T1, typename T2> - nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r][c]; } - template <typename T1, typename T2> - nl_ext_double &W(const T1 &r, const T2 &c) { return m_W[r][c]; } - template <typename T1, typename T2> - nl_ext_double &Ainv(const T1 &r, const T2 &c) { return m_Ainv[r][c]; } - template <typename T1> - nl_ext_double &RHS(const T1 &r) { return m_RHS[r]; } + using float_ext_type = FT; + using float_type = FT; + // FIXME: dirty hack to make this compile + static constexpr const std::size_t storage_N = 100; + matrix_solver_sm_t(netlist_state_t &anetlist, const pstring &name, + const solver_parameters_t *params, const std::size_t size); - template <typename T1, typename T2> - nl_ext_double &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; } - template <typename T1, typename T2> - nl_ext_double &lAinv(const T1 &r, const T2 &c) { return m_lAinv[r][c]; } + void vsetup(analog_net_t::list_t &nets) override; + void reset() override { matrix_solver_t::reset(); } - nl_double m_last_RHS[storage_N]; // right hand side - contains currents + protected: + unsigned vsolve_non_dynamic(const bool newton_raphson) override; + unsigned solve_non_dynamic(const bool newton_raphson); -private: - static constexpr std::size_t m_pitch = ((( storage_N) + 7) / 8) * 8; - nl_ext_double m_A[storage_N][m_pitch]; - nl_ext_double m_Ainv[storage_N][m_pitch]; - nl_ext_double m_W[storage_N][m_pitch]; - nl_ext_double m_RHS[storage_N]; // right hand side - contains currents + constexpr std::size_t size() const { return m_dim; } - nl_ext_double m_lA[storage_N][m_pitch]; - nl_ext_double m_lAinv[storage_N][m_pitch]; + void LE_invert(); - //nl_ext_double m_RHSx[storage_N]; + template <typename T> + void LE_compute_x(T * x); - const std::size_t m_dim; - std::size_t m_cnt; -}; + template <typename T1, typename T2> + float_ext_type &A(const T1 &r, const T2 &c) { return m_A[r][c]; } + template <typename T1, typename T2> + float_ext_type &W(const T1 &r, const T2 &c) { return m_W[r][c]; } + template <typename T1, typename T2> + float_ext_type &Ainv(const T1 &r, const T2 &c) { return m_Ainv[r][c]; } + template <typename T1> + float_ext_type &RHS(const T1 &r) { return m_RHS[r]; } -// ---------------------------------------------------------------------------------------- -// matrix_solver_direct -// ---------------------------------------------------------------------------------------- -template <std::size_t m_N, std::size_t storage_N> -matrix_solver_sm_t<m_N, storage_N>::~matrix_solver_sm_t() -{ -} + template <typename T1, typename T2> + float_ext_type &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; } + template <typename T1, typename T2> + float_ext_type &lAinv(const T1 &r, const T2 &c) { return m_lAinv[r][c]; } -template <std::size_t m_N, std::size_t storage_N> -void matrix_solver_sm_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets) -{ - matrix_solver_t::setup_base(nets); + private: + static constexpr std::size_t m_pitch = ((( storage_N) + 7) / 8) * 8; + float_ext_type m_A[storage_N][m_pitch]; + float_ext_type m_Ainv[storage_N][m_pitch]; + float_ext_type m_W[storage_N][m_pitch]; + float_ext_type m_RHS[storage_N]; // right hand side - contains currents - netlist().save(*this, m_last_RHS, "m_last_RHS"); + float_ext_type m_lA[storage_N][m_pitch]; + float_ext_type m_lAinv[storage_N][m_pitch]; - for (unsigned k = 0; k < N(); k++) - netlist().save(*this, RHS(k), plib::pfmt("RHS.{1}")(k)); -} + //float_ext_type m_RHSx[storage_N]; + const std::size_t m_dim; + std::size_t m_cnt; + }; -template <std::size_t m_N, std::size_t storage_N> -void matrix_solver_sm_t<m_N, storage_N>::LE_invert() -{ - const std::size_t kN = N(); + // ---------------------------------------------------------------------------------------- + // matrix_solver_direct + // ---------------------------------------------------------------------------------------- - for (std::size_t i = 0; i < kN; i++) + template <typename FT, int SIZE> + void matrix_solver_sm_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets) { - 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,i) = 1.0; + matrix_solver_t::setup_base(nets); + + /* FIXME: 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)); } - /* down */ - for (std::size_t i = 0; i < kN; i++) - { - /* FIXME: Singular matrix? */ - const nl_double f = 1.0 / W(i,i); - const auto * RESTRICT const p = m_terms[i]->m_nzrd.data(); - const std::size_t e = m_terms[i]->m_nzrd.size(); - /* Eliminate column i from row j */ + template <typename FT, int SIZE> + void matrix_solver_sm_t<FT, SIZE>::LE_invert() + { + const std::size_t kN = size(); - const auto * RESTRICT const pb = m_terms[i]->m_nzbd.data(); - const std::size_t eb = m_terms[i]->m_nzbd.size(); - for (std::size_t jb = 0; jb < eb; jb++) + for (std::size_t i = 0; i < kN; i++) { - const unsigned j = pb[jb]; - const nl_double f1 = - W(j,i) * f; - if (f1 != 0.0) + for (std::size_t j = 0; j < kN; j++) { - for (std::size_t k = 0; k < e; k++) - W(j,p[k]) += W(i,p[k]) * f1; - for (std::size_t k = 0; k <= i; k ++) - Ainv(j,k) += Ainv(i,k) * f1; + W(i,j) = lA(i,j) = A(i,j); + Ainv(i,j) = 0.0; } + Ainv(i,i) = 1.0; } - } - /* up */ - for (std::size_t i = kN; i-- > 0; ) - { - /* FIXME: Singular matrix? */ - const nl_double f = 1.0 / W(i,i); - for (std::size_t j = i; j-- > 0; ) + /* down */ + for (std::size_t i = 0; i < kN; i++) { - const nl_double f1 = - W(j,i) * f; - if (f1 != 0.0) + /* FIXME: Singular matrix? */ + const float_type f = 1.0 / W(i,i); + const auto * const p = m_terms[i]->m_nzrd.data(); + const std::size_t e = m_terms[i]->m_nzrd.size(); + + /* Eliminate column i from row j */ + + const auto * const pb = m_terms[i]->m_nzbd.data(); + const std::size_t eb = m_terms[i]->m_nzbd.size(); + for (std::size_t jb = 0; jb < eb; jb++) { - for (std::size_t k = i; k < kN; k++) - W(j,k) += W(i,k) * f1; - for (std::size_t k = 0; k < kN; k++) - Ainv(j,k) += Ainv(i,k) * f1; + const unsigned j = pb[jb]; + const float_type f1 = - W(j,i) * f; + if (f1 != 0.0) + { + for (std::size_t k = 0; k < e; k++) + W(j,p[k]) += W(i,p[k]) * f1; + for (std::size_t k = 0; k <= i; k ++) + Ainv(j,k) += Ainv(i,k) * f1; + } } } - for (std::size_t k = 0; k < kN; k++) + /* up */ + for (std::size_t i = kN; i-- > 0; ) { - Ainv(i,k) *= f; - lAinv(i,k) = Ainv(i,k); + /* FIXME: Singular matrix? */ + const float_type f = 1.0 / W(i,i); + for (std::size_t j = i; j-- > 0; ) + { + const float_type f1 = - W(j,i) * f; + if (f1 != 0.0) + { + for (std::size_t k = i; k < kN; k++) + W(j,k) += W(i,k) * f1; + for (std::size_t k = 0; k < kN; k++) + Ainv(j,k) += Ainv(i,k) * f1; + } + } + for (std::size_t k = 0; k < kN; k++) + { + Ainv(i,k) *= f; + lAinv(i,k) = Ainv(i,k); + } } } -} - -template <std::size_t m_N, std::size_t storage_N> -template <typename T> -void matrix_solver_sm_t<m_N, storage_N>::LE_compute_x( - T * RESTRICT x) -{ - const std::size_t kN = N(); - for (std::size_t i=0; i<kN; i++) - x[i] = 0.0; - - for (std::size_t k=0; k<kN; k++) + template <typename FT, int SIZE> + template <typename T> + void matrix_solver_sm_t<FT, SIZE>::LE_compute_x( + T * x) { - const nl_double f = RHS(k); + const std::size_t kN = size(); for (std::size_t i=0; i<kN; i++) - x[i] += Ainv(i,k) * f; - } -} - - -template <std::size_t m_N, std::size_t storage_N> -unsigned matrix_solver_sm_t<m_N, storage_N>::solve_non_dynamic(const bool newton_raphson) -{ - static constexpr const bool incremental = true; - const std::size_t iN = N(); + x[i] = 0.0; - nl_double new_V[storage_N]; // = { 0.0 }; + for (std::size_t k=0; k<kN; k++) + { + const float_type f = RHS(k); - if ((m_cnt % 50) == 0) - { - /* complete calculation */ - this->LE_invert(); + for (std::size_t i=0; i<kN; i++) + x[i] += Ainv(i,k) * f; + } } - else + + template <typename FT, int SIZE> + unsigned matrix_solver_sm_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphson) { - if (!incremental) + static constexpr const bool incremental = true; + const std::size_t iN = size(); + + float_type new_V[storage_N]; // = { 0.0 }; + + if ((m_cnt % 50) == 0) { - for (std::size_t row = 0; row < iN; row ++) - for (std::size_t k = 0; k < iN; k++) - Ainv(row,k) = lAinv(row, k); + /* complete calculation */ + this->LE_invert(); } - for (std::size_t row = 0; row < iN; row ++) + else { - nl_double v[m_pitch] = {0}; - std::size_t cols[m_pitch]; - std::size_t colcount = 0; - - auto &nz = m_terms[row]->m_nz; - for (unsigned & col : nz) + if (!incremental) { - v[col] = A(row,col) - lA(row,col); - if (incremental) - lA(row,col) = A(row,col); - if (v[col] != 0.0) - cols[colcount++] = col; + for (std::size_t row = 0; row < iN; row ++) + for (std::size_t k = 0; k < iN; k++) + Ainv(row,k) = lAinv(row, k); } - - if (colcount > 0) + for (std::size_t row = 0; row < iN; row ++) { - nl_double lamba = 0.0; - nl_double w[m_pitch] = {0}; - - nl_double z[m_pitch]; - /* compute w and lamba */ - for (std::size_t i = 0; i < iN; i++) - z[i] = Ainv(i, row); /* u is row'th column */ + float_type v[m_pitch] = {0}; + std::size_t cols[m_pitch]; + std::size_t colcount = 0; - for (std::size_t j = 0; j < colcount; j++) - lamba += v[cols[j]] * z[cols[j]]; - - for (std::size_t j=0; j<colcount; j++) + auto &nz = m_terms[row]->m_nz; + for (unsigned & col : nz) { - std::size_t col = cols[j]; - nl_double f = v[col]; - for (std::size_t k = 0; k < iN; k++) - w[k] += Ainv(col,k) * f; /* Transpose(Ainv) * v */ + v[col] = A(row,col) - lA(row,col); + if (incremental) + lA(row,col) = A(row,col); + if (v[col] != 0.0) + cols[colcount++] = col; } - lamba = -1.0 / (1.0 + lamba); - for (std::size_t i=0; i<iN; i++) + if (colcount > 0) { - const nl_double f = lamba * z[i]; - if (f != 0.0) + float_type lamba = 0.0; + float_type w[m_pitch] = {0}; + + float_type z[m_pitch]; + /* compute w and lamba */ + for (std::size_t i = 0; i < iN; i++) + z[i] = Ainv(i, row); /* u is row'th column */ + + for (std::size_t j = 0; j < colcount; j++) + lamba += v[cols[j]] * z[cols[j]]; + + for (std::size_t j=0; j<colcount; j++) + { + std::size_t col = cols[j]; + float_type f = v[col]; for (std::size_t k = 0; k < iN; k++) - Ainv(i,k) += f * w[k]; + w[k] += Ainv(col,k) * f; /* Transpose(Ainv) * v */ + } + + lamba = -1.0 / (1.0 + lamba); + for (std::size_t i=0; i<iN; i++) + { + const float_type f = lamba * z[i]; + if (f != 0.0) + for (std::size_t k = 0; k < iN; k++) + Ainv(i,k) += f * w[k]; + } } - } + } } - } - m_cnt++; + m_cnt++; - this->LE_compute_x(new_V); + this->LE_compute_x(new_V); - const nl_double err = (newton_raphson ? delta(new_V) : 0.0); - store(new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; -} + const float_type err = (newton_raphson ? delta(new_V) : 0.0); + store(new_V); + return (err > this->m_params.m_accuracy) ? 2 : 1; + } -template <std::size_t m_N, std::size_t storage_N> -inline unsigned matrix_solver_sm_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson) -{ - build_LE_A<matrix_solver_sm_t>(); - build_LE_RHS<matrix_solver_sm_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_sm_t<m_N, storage_N>::matrix_solver_sm_t(netlist_t &anetlist, const pstring &name, - const solver_parameters_t *params, const std::size_t size) -: matrix_solver_t(anetlist, name, NOSORT, params) -, m_dim(size) -, m_cnt(0) -{ - for (std::size_t k = 0; k < N(); k++) + template <typename FT, int SIZE> + unsigned matrix_solver_sm_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson) + { + this->build_LE_A(*this); + this->build_LE_RHS(*this); + + this->m_stat_calculations++; + return this->solve_non_dynamic(newton_raphson); + } + + template <typename FT, int SIZE> + matrix_solver_sm_t<FT, SIZE>::matrix_solver_sm_t(netlist_state_t &anetlist, const pstring &name, + const solver_parameters_t *params, const std::size_t size) + : matrix_solver_t(anetlist, name, NOSORT, params) + , m_dim(size) + , m_cnt(0) { - m_last_RHS[k] = 0.0; } -} - } //namespace devices +} // namespace devices } // namespace netlist #endif /* NLD_MS_DIRECT_H_ */ |