diff options
Diffstat (limited to 'src/lib/netlist/solver/nld_ms_direct.h')
-rw-r--r-- | src/lib/netlist/solver/nld_ms_direct.h | 217 |
1 files changed, 92 insertions, 125 deletions
diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h index 2501742218d..a3e9b6345b5 100644 --- a/src/lib/netlist/solver/nld_ms_direct.h +++ b/src/lib/netlist/solver/nld_ms_direct.h @@ -1,63 +1,56 @@ -// license:GPL-2.0+ +// license:BSD-3-Clause // copyright-holders:Couriersud -/* - * nld_ms_direct.h - * - */ #ifndef NLD_MS_DIRECT_H_ #define NLD_MS_DIRECT_H_ +/// +/// \file nld_ms_direct.h +/// + #include "nld_matrix_solver.h" #include "nld_solver.h" -#include "plib/mat_cr.h" +#include "plib/parray.h" #include "plib/vector_ops.h" +#include "nld_matrix_solver_ext.h" + #include <algorithm> -#include <cmath> -namespace netlist -{ -namespace devices +namespace netlist::solver { template <typename FT, int SIZE> - class matrix_solver_direct_t: public matrix_solver_t + class matrix_solver_direct_t: public matrix_solver_ext_t<FT, SIZE> { - friend class matrix_solver_t; public: using float_type = FT; - 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); + matrix_solver_direct_t(devices::nld_solver &main_solver, const pstring &name, + const matrix_solver_t::net_list_t &nets, + const solver::solver_parameters_t *params, std::size_t size); - void vsetup(analog_net_t::list_t &nets) override; void reset() override { matrix_solver_t::reset(); } + private: + + const std::size_t m_pitch; + protected: - unsigned vsolve_non_dynamic(const bool newton_raphson) override; - unsigned solve_non_dynamic(const bool newton_raphson); + static constexpr const std::size_t SIZEABS = plib::parray<FT, SIZE>::SIZEABS(); + static constexpr const std::size_t m_pitch_ABS = (((SIZEABS + 0) + 7) / 8) * 8; - constexpr std::size_t size() const { return (SIZE > 0) ? static_cast<std::size_t>(SIZE) : m_dim; } + void upstream_solve_non_dynamic() override; + void solve_non_dynamic(); void LE_solve(); template <typename T> void LE_back_subst(T & x); - 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; - - 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; - - const std::size_t m_dim; - const std::size_t m_pitch; - plib::parray<FT, SIZE * int(m_pitch_ABS)> m_A; - + // PALIGNAS_VECTOROPT() `parray` defines alignment already + plib::parray2D<FT, SIZE, m_pitch_ABS> m_A; }; // ---------------------------------------------------------------------------------------- @@ -65,43 +58,26 @@ namespace devices // ---------------------------------------------------------------------------------------- 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); - - /* 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>(size()))) - t->m_nzrd.push_back(static_cast<unsigned>(size())); - } - - // 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)); - } - - template <typename FT, int SIZE> void matrix_solver_direct_t<FT, SIZE>::LE_solve() { - const std::size_t kN = size(); - if (!m_params.m_pivot) + const std::size_t kN = this->size(); + if (!this->m_params.m_pivot) { 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; + // FIXME: Singular matrix? + const auto &Ai = m_A[i]; + const FT f = plib::reciprocal(Ai[i]); + const auto &nzrd = this->m_terms[i].m_nzrd; + const auto &nzbd = this->m_terms[i].m_nzbd; - for (std::size_t j : nzbd) + for (auto &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; + auto &Aj = m_A[j]; + const FT f1 = -f * Aj[i]; + for (auto &k : nzrd) + Aj[k] += Ai[k] * f1; + this->m_RHS[j] += this->m_RHS[i] * f1; } } } @@ -109,45 +85,47 @@ namespace devices { for (std::size_t i = 0; i < kN; i++) { - /* Find the row with the largest first value */ - std::size_t maxrow = i; + // Find the row with the largest first value + std::size_t max_row = 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 (plib::abs(m_A[j][i]) > plib::abs(m_A[max_row][i])) + //#if (m_A[j][i] * m_A[j][i] > m_A[max_row][i] * m_A[max_row][i]) + max_row = j; } - if (maxrow != i) + if (max_row != 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)); +#if 0 + // Swap the max_row and ith row + for (std::size_t k = 0; k < kN; k++) { + std::swap(m_A[i][k], m_A[max_row][k]); } - //std::swap(RHS(i), RHS(maxrow)); +#else + std::swap(m_A[i], m_A[max_row]); +#endif + std::swap(this->m_RHS[i], this->m_RHS[max_row]); } - /* FIXME: Singular matrix? */ - const FT f = 1.0 / A(i,i); + // FIXME: Singular matrix? + const auto &Ai = m_A[i]; + const FT f = plib::reciprocal(Ai[i]); - /* Eliminate column i from row j */ + // Eliminate column i from row j for (std::size_t j = i + 1; j < kN; j++) { - const FT f1 = - A(j,i) * f; + auto &Aj = m_A[j]; + const FT f1 = - m_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); + const FT * pi = &(Ai[i+1]); + FT * pj = &(Aj[i+1]); + plib::vec_add_mult_scalar_p(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 + this->m_RHS[j] += this->m_RHS[i] * f1; } } } @@ -159,77 +137,66 @@ namespace devices void matrix_solver_direct_t<FT, SIZE>::LE_back_subst( T & x) { - const std::size_t kN = size(); + const std::size_t kN = this->size(); - /* back substitution */ - if (m_params.m_pivot) + // back substitution + if (this->m_params.m_pivot) { for (std::size_t j = kN; j-- > 0; ) { - FT tmp = 0; + FT tmp(0); + const auto & Aj(m_A[j]); + for (std::size_t k = j+1; k < kN; k++) - tmp += A(j,k) * x[k]; - x[j] = (RHS(j) - tmp) / A(j,j); + tmp += Aj[k] * x[k]; + x[j] = (this->m_RHS[j] - tmp) / Aj[j]; } } else { for (std::size_t j = kN; j-- > 0; ) { - FT tmp = 0; - const auto &nzrd = m_terms[j]->m_nzrd; - const auto e = nzrd.size() - 1; /* exclude RHS element */ + FT tmp(0); + const auto &nzrd = this->m_terms[j].m_nzrd; + const auto & Aj(m_A[j]); + const auto e = nzrd.size(); + 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); + tmp += Aj[nzrd[k]] * x[nzrd[k]]; + x[j] = (this->m_RHS[j] - tmp) / Aj[j]; } } } template <typename FT, int SIZE> - unsigned matrix_solver_direct_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphson) + void matrix_solver_direct_t<FT, SIZE>::solve_non_dynamic() { 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; + this->LE_back_subst(this->m_new_V); } template <typename FT, int SIZE> - unsigned matrix_solver_direct_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson) + void matrix_solver_direct_t<FT, SIZE>::upstream_solve_non_dynamic() { - this->build_LE_A(*this); - this->build_LE_RHS(*this); + // populate matrix + this->clear_square_mat(m_A); + this->fill_matrix_and_rhs(); - this->m_stat_calculations++; - return this->solve_non_dynamic(newton_raphson); - } - - 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) - { + this->solve_non_dynamic(); } 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) + matrix_solver_direct_t<FT, SIZE>::matrix_solver_direct_t(devices::nld_solver &main_solver, const pstring &name, + const matrix_solver_t::net_list_t &nets, + const solver::solver_parameters_t *params, + std::size_t size) + : matrix_solver_ext_t<FT, SIZE>(main_solver, name, nets, params, size) + , m_pitch(m_pitch_ABS ? m_pitch_ABS : (((size + 0) + 7) / 8) * 8) + , m_A(size, m_pitch) { + this->build_mat_ptr(m_A); } -} // namespace devices -} // namespace netlist +} // namespace netlist::solver -#endif /* NLD_MS_DIRECT_H_ */ +#endif // NLD_MS_DIRECT_H_ |