diff options
Diffstat (limited to 'src/lib/netlist/solver/nld_ms_direct.h')
-rw-r--r-- | src/lib/netlist/solver/nld_ms_direct.h | 88 |
1 files changed, 37 insertions, 51 deletions
diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h index 582a7996947..f3c31dc3371 100644 --- a/src/lib/netlist/solver/nld_ms_direct.h +++ b/src/lib/netlist/solver/nld_ms_direct.h @@ -35,7 +35,15 @@ namespace solver void reset() override { matrix_solver_t::reset(); } + private: + + const std::size_t m_dim; + const std::size_t m_pitch; + protected: + 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; + unsigned vsolve_non_dynamic(const bool newton_raphson) override; unsigned solve_non_dynamic(const bool newton_raphson); @@ -46,24 +54,12 @@ namespace solver template <typename T> void LE_back_subst(T & x); - const FT &A(std::size_t r, std::size_t c) const noexcept { return m_A[r][c]; } - FT &A(std::size_t r, std::size_t c) noexcept { return m_A[r][c]; } - - const FT &RHS(std::size_t r) const noexcept { return m_A[r][size()]; } - FT &RHS(std::size_t r) noexcept { return m_A[r][size()]; } - PALIGNAS_VECTOROPT() 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; - PALIGNAS_VECTOROPT() plib::parray2D<FT, SIZE, m_pitch_ABS> m_A; + PALIGNAS_VECTOROPT() + plib::parray<FT, SIZE> m_RHS; }; // ---------------------------------------------------------------------------------------- @@ -79,16 +75,16 @@ namespace solver for (std::size_t i = 0; i < kN; i++) { /* FIXME: Singular matrix? */ - const FT f = 1.0 / A(i,i); + const FT f = 1.0 / m_A[i][i]; const auto &nzrd = m_terms[i].m_nzrd; const auto &nzbd = m_terms[i].m_nzbd; for (auto &j : nzbd) { - const FT f1 = -f * A(j, i); + const FT f1 = -f * m_A[j][i]; for (auto &k : nzrd) - A(j, k) += A(i, k) * f1; - //RHS(j) += RHS(i) * f1; + m_A[j][k] += m_A[i][k] * f1; + m_RHS[j] += m_RHS[i] * f1; } } } @@ -101,7 +97,7 @@ namespace solver 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)) + if (m_A[j][i] * m_A[j][i] > m_A[maxrow][i] * m_A[maxrow][i]) maxrow = j; } @@ -109,32 +105,28 @@ namespace solver { /* 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(m_A[i][k], m_A[maxrow][k]); } //std::swap(RHS(i), RHS(maxrow)); } /* FIXME: Singular matrix? */ - const FT f = 1.0 / A(i,i); + const FT f = 1.0 / m_A[i][i]; /* Eliminate column i from row j */ for (std::size_t j = i + 1; j < kN; j++) { - const FT f1 = - A(j,i) * f; + 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 = &(m_A[i][i+1]); + FT * pj = &(m_A[j][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 + m_RHS[j] += m_RHS[i] * f1; } } } @@ -155,8 +147,8 @@ namespace solver { 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); + tmp += m_A[j][k] * x[k]; + x[j] = (m_RHS[j] - tmp) / m_A[j][j]; } } else @@ -165,10 +157,10 @@ namespace solver { FT tmp = 0; const auto &nzrd = m_terms[j].m_nzrd; - const auto e = nzrd.size() - 1; /* exclude RHS element */ + 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); + tmp += m_A[j][nzrd[k]] * x[nzrd[k]]; + x[j] = (m_RHS[j] - tmp) / m_A[j][j]; } } } @@ -187,8 +179,12 @@ namespace solver 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); + const std::size_t iN = this->size(); + + /* populate matrix */ + + this->clear_square_mat(iN, m_A); + this->fill_matrix(iN, m_RHS); this->m_stat_calculations++; return this->solve_non_dynamic(newton_raphson); @@ -200,23 +196,13 @@ namespace solver const solver_parameters_t *params, const std::size_t size) : matrix_solver_t(anetlist, name, nets, params) - , m_new_V(size) , m_dim(size) - , m_pitch(m_pitch_ABS ? m_pitch_ABS : (((m_dim + 1) + 7) / 8) * 8) + , m_pitch(m_pitch_ABS ? m_pitch_ABS : (((m_dim + 0) + 7) / 8) * 8) + , m_new_V(size) , m_A(size, m_pitch) + , m_RHS(size) { - /* add RHS element */ - for (std::size_t k = 0; k < this->size(); k++) - { - terms_for_net_t & t = m_terms[k]; - - if (!plib::container::contains(t.m_nzrd, static_cast<unsigned>(this->size()))) - t.m_nzrd.push_back(static_cast<unsigned>(this->size())); - } - - // FIXME: This shouldn't be necessary ... - for (std::size_t k = 0; k < this->size(); k++) - state().save(*this, RHS(k), this->name(), plib::pfmt("RHS.{1}")(k)); + this->build_mat_ptr(size, m_A); } } // namespace solver |