diff options
Diffstat (limited to 'src/lib/netlist/solver/nld_ms_sm.h')
-rw-r--r-- | src/lib/netlist/solver/nld_ms_sm.h | 226 |
1 files changed, 104 insertions, 122 deletions
diff --git a/src/lib/netlist/solver/nld_ms_sm.h b/src/lib/netlist/solver/nld_ms_sm.h index 5b59effecf8..8a0b67f18ce 100644 --- a/src/lib/netlist/solver/nld_ms_sm.h +++ b/src/lib/netlist/solver/nld_ms_sm.h @@ -1,54 +1,50 @@ -// license:GPL-2.0+ +// license:BSD-3-Clause // copyright-holders:Couriersud -/* - * nld_ms_direct.h - * - * - * Sherman-Morrison Solver - * - * Computes the updated inverse of A given that the change in A is - * - * A <- A + (u x v) u,v vectors - * - * In this specific implementation, u is a unit vector specifying the row which - * changed. Thus v contains the changed column. - * - * Than z = A^-1 u , w = transpose(A^-1) v , lambda = v z - * - * A^-1 <- 1.0 / (1.0 + lambda) * (z x w) - * - * The approach is iterative and applied for each row changed. - * - * The performance for a typical circuit like kidniki compared to Gaussian - * elimination is poor: - * - * a) The code needs to be run for each row change. - * b) The inverse of A typically is fully occupied. - * - * It may have advantages for circuits with a high number of elements and only - * few dynamic/active components. - * - */ #ifndef NLD_MS_SM_H_ #define NLD_MS_SM_H_ +/// +/// \file nld_ms_sm.h +/// +/// Sherman-Morrison Solver +/// +/// Computes the updated inverse of A given that the change in A is +/// +/// A <- A + (u x v) u,v vectors +/// +/// In this specific implementation, u is a unit vector specifying the row which +/// changed. Thus v contains the changed column. +/// +/// Than z = A^-1 u , w = transpose(A^-1) v , lambda = v z +/// +/// A^-1 <- 1.0 / (1.0 + lambda) * (z x w) +/// +/// The approach is iterative and applied for each row changed. +/// +/// The performance for a typical circuit like `kidniki` compared to Gaussian +/// elimination is poor: +/// +/// a) The code needs to be run for each row change. +/// b) The inverse of A typically is fully occupied. +/// +/// It may have advantages for circuits with a high number of elements and only +/// few dynamic/active components. +/// + #include "nld_matrix_solver.h" +#include "nld_matrix_solver_ext.h" #include "nld_solver.h" #include "plib/vector_ops.h" #include <algorithm> -namespace netlist -{ -namespace devices +namespace netlist::solver { template <typename FT, int SIZE> - class matrix_solver_sm_t: public matrix_solver_t + class matrix_solver_sm_t: public matrix_solver_ext_t<FT, SIZE> { - friend class matrix_solver_t; - public: using float_ext_type = FT; @@ -56,17 +52,20 @@ namespace devices // 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); + matrix_solver_sm_t(devices::nld_solver &main_solver, const pstring &name, + const matrix_solver_t::net_list_t &nets, + const solver_parameters_t *params, const std::size_t size) + : matrix_solver_ext_t<FT, SIZE>(main_solver, name, nets, params, size) + , m_cnt(0) + { + this->build_mat_ptr(m_A); + } - 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); - - constexpr std::size_t size() const { return m_dim; } + void upstream_solve_non_dynamic() override; + void solve_non_dynamic(); void LE_invert(); @@ -81,7 +80,7 @@ namespace devices 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]; } + float_ext_type &RHS(const T1 &r) { return this->m_RHS[r]; } template <typename T1, typename T2> @@ -90,18 +89,18 @@ namespace devices float_ext_type &lAinv(const T1 &r, const T2 &c) { return m_lAinv[r][c]; } private: + template <typename T, std::size_t N, std::size_t M> + using array2D = std::array<std::array<T, M>, N>; 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]; - std::array<float_ext_type, storage_N> m_RHS; // right hand side - contains currents + array2D<float_ext_type, storage_N, m_pitch> m_A; + array2D<float_ext_type, storage_N, m_pitch> m_Ainv; + array2D<float_ext_type, storage_N, m_pitch> m_W; - float_ext_type m_lA[storage_N][m_pitch]; - float_ext_type m_lAinv[storage_N][m_pitch]; + array2D<float_ext_type, storage_N, m_pitch> m_lA; + array2D<float_ext_type, storage_N, m_pitch> m_lAinv; //float_ext_type m_RHSx[storage_N]; - const std::size_t m_dim; std::size_t m_cnt; }; @@ -111,46 +110,37 @@ namespace devices // ---------------------------------------------------------------------------------------- template <typename FT, int SIZE> - void matrix_solver_sm_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets) - { - 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)); - } - - template <typename FT, int SIZE> void matrix_solver_sm_t<FT, SIZE>::LE_invert() { - const std::size_t kN = size(); + const std::size_t kN = this->size(); for (std::size_t i = 0; i < kN; i++) { 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,j) = plib::constants<FT>::zero(); } - Ainv(i,i) = 1.0; + Ainv(i,i) = plib::constants<FT>::one(); } - /* down */ + // down for (std::size_t i = 0; i < kN; i++) { - /* 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(); + // FIXME: Singular matrix? + const float_type f = plib::reciprocal(W(i,i)); + const auto * const p = this->m_terms[i].m_nzrd.data(); + const std::size_t e = this->m_terms[i].m_nzrd.size(); - /* Eliminate column i from row j */ + // 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(); + const auto * const pb = this->m_terms[i].m_nzbd.data(); + const std::size_t eb = this->m_terms[i].m_nzbd.size(); for (std::size_t jb = 0; jb < eb; jb++) { const unsigned j = pb[jb]; const float_type f1 = - W(j,i) * f; - if (f1 != 0.0) + // FIXME: comparison to zero + if (f1 != plib::constants<float_type>::zero()) { for (std::size_t k = 0; k < e; k++) W(j,p[k]) += W(i,p[k]) * f1; @@ -159,15 +149,16 @@ namespace devices } } } - /* up */ + // up for (std::size_t i = kN; i-- > 0; ) { - /* FIXME: Singular matrix? */ - const float_type f = 1.0 / W(i,i); + // FIXME: Singular matrix? + const float_type f = plib::reciprocal(W(i,i)); for (std::size_t j = i; j-- > 0; ) { const float_type f1 = - W(j,i) * f; - if (f1 != 0.0) + // FIXME: comparison to zero + if (f1 != plib::constants<float_type>::zero()) { for (std::size_t k = i; k < kN; k++) W(j,k) += W(i,k) * f1; @@ -188,10 +179,10 @@ namespace devices void matrix_solver_sm_t<FT, SIZE>::LE_compute_x( T & x) { - const std::size_t kN = size(); + const std::size_t kN = this->size(); for (std::size_t i=0; i<kN; i++) - x[i] = 0.0; + x[i] = plib::constants<FT>::zero(); for (std::size_t k=0; k<kN; k++) { @@ -203,16 +194,21 @@ namespace devices } template <typename FT, int SIZE> - unsigned matrix_solver_sm_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphson) + void matrix_solver_sm_t<FT, SIZE>::solve_non_dynamic() { static constexpr const bool incremental = true; - const std::size_t iN = size(); + const std::size_t iN = this->size(); - std::array<float_type, storage_N> new_V; // = { 0.0 }; + // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) + std::array<float_type, m_pitch> v; + // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) + std::array<std::size_t, m_pitch> cols; + // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) + std::array<float_type, m_pitch> z; if ((m_cnt % 50) == 0) { - /* complete calculation */ + // complete calculation this->LE_invert(); } else @@ -225,46 +221,45 @@ namespace devices } for (std::size_t row = 0; row < iN; row ++) { - std::array<float_type, m_pitch> v = {0}; - std::array<std::size_t, m_pitch> cols; - std::size_t colcount = 0; + std::size_t col_count = 0; - auto &nz = m_terms[row]->m_nz; + auto &nz = this->m_terms[row].m_nz; for (unsigned & col : nz) { v[col] = A(row,col) - lA(row,col); if (incremental) lA(row,col) = A(row,col); - if (v[col] != 0.0) - cols[colcount++] = col; + // FIXME: comparison to zero + if (v[col] != plib::constants<float_type>::zero()) + cols[col_count++] = col; } - if (colcount > 0) + if (col_count > 0) { - float_type lamba = 0.0; + auto lambda(plib::constants<FT>::zero()); std::array<float_type, m_pitch> w = {0}; - std::array<float_type, m_pitch> z; - /* compute w and lamba */ + // compute w and lambda for (std::size_t i = 0; i < iN; i++) - z[i] = Ainv(i, row); /* u is row'th column */ + 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 < col_count; j++) + lambda += v[cols[j]] * z[cols[j]]; - for (std::size_t j=0; j<colcount; j++) + for (std::size_t j=0; j<col_count; j++) { std::size_t col = cols[j]; float_type f = v[col]; for (std::size_t k = 0; k < iN; k++) - w[k] += Ainv(col,k) * f; /* Transpose(Ainv) * v */ + w[k] += Ainv(col,k) * f; //# Transpose(Ainv) * v } - lamba = -1.0 / (1.0 + lamba); + lambda = -plib::reciprocal(plib::constants<float_type>::one() + lambda); for (std::size_t i=0; i<iN; i++) { - const float_type f = lamba * z[i]; - if (f != 0.0) + const float_type f = lambda * z[i]; + // FIXME: comparison to zero + if (f != plib::constants<float_type>::zero()) for (std::size_t k = 0; k < iN; k++) Ainv(i,k) += f * w[k]; } @@ -275,33 +270,20 @@ namespace devices m_cnt++; - this->LE_compute_x(new_V); - - const float_type err = (newton_raphson ? delta(new_V) : 0.0); - store(new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; + this->LE_compute_x(this->m_new_V); } template <typename FT, int SIZE> - unsigned matrix_solver_sm_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson) + void matrix_solver_sm_t<FT, SIZE>::upstream_solve_non_dynamic() { - this->build_LE_A(*this); - this->build_LE_RHS(*this); - this->m_stat_calculations++; - return this->solve_non_dynamic(newton_raphson); - } + this->clear_square_mat(this->m_A); + this->fill_matrix_and_rhs(); - 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) - { + this->solve_non_dynamic(); } -} // namespace devices -} // namespace netlist -#endif /* NLD_MS_DIRECT_H_ */ +} // namespace netlist::solver + +#endif // NLD_MS_SM_H_ |