summaryrefslogtreecommitdiffstatshomepage
path: root/src/lib/netlist/solver/nld_ms_direct.h
diff options
context:
space:
mode:
Diffstat (limited to 'src/lib/netlist/solver/nld_ms_direct.h')
-rw-r--r--src/lib/netlist/solver/nld_ms_direct.h217
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_