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.h386
1 files changed, 222 insertions, 164 deletions
diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h
index 2501742218d..f03cb738028 100644
--- a/src/lib/netlist/solver/nld_ms_direct.h
+++ b/src/lib/netlist/solver/nld_ms_direct.h
@@ -8,228 +8,286 @@
#ifndef NLD_MS_DIRECT_H_
#define NLD_MS_DIRECT_H_
-#include "nld_matrix_solver.h"
+#include <algorithm>
+
#include "nld_solver.h"
-#include "plib/mat_cr.h"
-#include "plib/vector_ops.h"
+#include "nld_matrix_solver.h"
+#include "vector_base.h"
-#include <algorithm>
-#include <cmath>
+/* Disabling dynamic allocation gives a ~10% boost in performance
+ * This flag has been added to support continuous storage for arrays
+ * going forward in case we implement cuda solvers in the future.
+ */
+#define NL_USE_DYNAMIC_ALLOCATION (1)
namespace netlist
{
-namespace devices
-{
-
- template <typename FT, int SIZE>
- class matrix_solver_direct_t: public matrix_solver_t
+ namespace devices
{
- friend class matrix_solver_t;
- public:
-
- using float_type = FT;
+//#define nl_ext_double _float128 // slow, very slow
+//#define nl_ext_double long double // slightly slower
+#define nl_ext_double nl_double
- 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);
- 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);
+template <std::size_t m_N, std::size_t storage_N>
+class matrix_solver_direct_t: public matrix_solver_t
+{
+ friend class matrix_solver_t;
+public:
- constexpr std::size_t size() const { return (SIZE > 0) ? static_cast<std::size_t>(SIZE) : m_dim; }
+ matrix_solver_direct_t(netlist_t &anetlist, const pstring &name, const solver_parameters_t *params, const std::size_t size);
+ matrix_solver_direct_t(netlist_t &anetlist, const pstring &name, const eSortType sort, const solver_parameters_t *params, const std::size_t size);
- void LE_solve();
+ virtual ~matrix_solver_direct_t() override;
- template <typename T>
- void LE_back_subst(T & x);
+ virtual void vsetup(analog_net_t::list_t &nets) override;
+ virtual void reset() override { matrix_solver_t::reset(); }
- 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;
+protected:
+ virtual unsigned vsolve_non_dynamic(const bool newton_raphson) override;
+ unsigned solve_non_dynamic(const bool newton_raphson);
- 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;
+ constexpr std::size_t N() const { return (m_N == 0) ? m_dim : m_N; }
- const std::size_t m_dim;
- const std::size_t m_pitch;
- plib::parray<FT, SIZE * int(m_pitch_ABS)> m_A;
+ void LE_solve();
- };
+ template <typename T>
+ void LE_back_subst(T * RESTRICT x);
+
+#if (NL_USE_DYNAMIC_ALLOCATION)
+ template <typename T1, typename T2>
+ nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r * m_pitch + c]; }
+ template <typename T1>
+ nl_ext_double &RHS(const T1 &r) { return m_A[r * m_pitch + N()]; }
+#else
+ template <typename T1, typename T2>
+ nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r][c]; }
+ template <typename T1>
+ nl_ext_double &RHS(const T1 &r) { return m_A[r][N()]; }
+#endif
+ nl_double m_last_RHS[storage_N]; // right hand side - contains currents
+
+private:
+ //static const std::size_t m_pitch = (((storage_N + 1) + 0) / 1) * 1;
+ static constexpr std::size_t m_pitch = (((storage_N + 1) + 7) / 8) * 8;
+ //static const std::size_t m_pitch = (((storage_N + 1) + 15) / 16) * 16;
+ //static const std::size_t m_pitch = (((storage_N + 1) + 31) / 32) * 32;
+#if (NL_USE_DYNAMIC_ALLOCATION)
+ //nl_ext_double * RESTRICT m_A;
+ std::vector<nl_ext_double> m_A;
+#else
+ nl_ext_double m_A[storage_N][m_pitch];
+#endif
+ //nl_ext_double m_RHSx[storage_N];
+
+ const std::size_t m_dim;
+
+};
+
+// ----------------------------------------------------------------------------------------
+// matrix_solver_direct
+// ----------------------------------------------------------------------------------------
+
+template <std::size_t m_N, std::size_t storage_N>
+matrix_solver_direct_t<m_N, storage_N>::~matrix_solver_direct_t()
+{
+#if (NL_USE_DYNAMIC_ALLOCATION)
+ //plib::pfree_array(m_A);
+#endif
+}
- // ----------------------------------------------------------------------------------------
- // matrix_solver_direct
- // ----------------------------------------------------------------------------------------
+template <std::size_t m_N, std::size_t storage_N>
+void matrix_solver_direct_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets)
+{
+ matrix_solver_t::setup_base(nets);
- template <typename FT, int SIZE>
- void matrix_solver_direct_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets)
+ /* add RHS element */
+ for (std::size_t k = 0; k < N(); k++)
{
- matrix_solver_t::setup_base(nets);
+ terms_for_net_t * t = m_terms[k].get();
- /* 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>(N())))
+ t->m_nzrd.push_back(static_cast<unsigned>(N()));
+ }
- if (!plib::container::contains(t->m_nzrd, static_cast<unsigned>(size())))
- t->m_nzrd.push_back(static_cast<unsigned>(size()));
- }
+ netlist().save(*this, m_last_RHS, "m_last_RHS");
- // 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));
- }
+ for (std::size_t k = 0; k < N(); k++)
+ netlist().save(*this, RHS(k), plib::pfmt("RHS.{1}")(k));
+}
- template <typename FT, int SIZE>
- void matrix_solver_direct_t<FT, SIZE>::LE_solve()
+
+template <std::size_t m_N, std::size_t storage_N>
+void matrix_solver_direct_t<m_N, storage_N>::LE_solve()
+{
+ const std::size_t kN = N();
+ if (!m_params.m_pivot)
{
- const std::size_t kN = size();
- if (!m_params.m_pivot)
+ for (std::size_t i = 0; i < kN; i++)
{
- 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;
- for (std::size_t 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;
- }
+ /* FIXME: Singular matrix? */
+ nl_double *Ai = &A(i, 0);
+ const nl_double f = 1.0 / A(i,i);
+ const auto &nzrd = m_terms[i]->m_nzrd;
+ const auto &nzbd = m_terms[i]->m_nzbd;
+
+ for (std::size_t j : nzbd)
+ {
+ nl_double *Aj = &A(j, 0);
+ const nl_double f1 = -f * Aj[i];
+ for (std::size_t k : nzrd)
+ Aj[k] += Ai[k] * f1;
+ //RHS(j) += RHS(i) * f1;
}
}
- else
+ }
+ else
+ {
+ for (std::size_t i = 0; i < kN; i++)
{
- for (std::size_t i = 0; i < kN; i++)
+ /* Find the row with the largest first value */
+ std::size_t maxrow = i;
+ for (std::size_t j = i + 1; j < kN; j++)
{
- /* Find the row with the largest first value */
- std::size_t maxrow = 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 (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 (maxrow != 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));
- }
- //std::swap(RHS(i), RHS(maxrow));
+ if (maxrow != 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));
}
- /* FIXME: Singular matrix? */
- const FT f = 1.0 / A(i,i);
+ //std::swap(RHS(i), RHS(maxrow));
+ }
+ /* FIXME: Singular matrix? */
+ const nl_double f = 1.0 / A(i,i);
- /* Eliminate column i from row j */
+ /* Eliminate column i from row j */
- for (std::size_t j = i + 1; j < kN; j++)
+ for (std::size_t j = i + 1; j < kN; j++)
+ {
+ const nl_double f1 = - A(j,i) * f;
+ if (f1 != NL_FCONST(0.0))
{
- const FT f1 = - 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);
- //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
- }
+ const nl_double * RESTRICT pi = &A(i,i+1);
+ nl_double * RESTRICT pj = &A(j,i+1);
+#if 1
+ vec_add_mult_scalar_p(kN-i,pi,f1,pj);
+#else
+ vec_add_mult_scalar_p(kN-i-1,pj,f1,pi);
+ //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
}
}
}
}
+}
- template <typename FT, int SIZE>
- template <typename T>
- void matrix_solver_direct_t<FT, SIZE>::LE_back_subst(
- T & x)
- {
- const std::size_t kN = size();
+template <std::size_t m_N, std::size_t storage_N>
+template <typename T>
+void matrix_solver_direct_t<m_N, storage_N>::LE_back_subst(
+ T * RESTRICT x)
+{
+ const std::size_t kN = N();
- /* back substitution */
- if (m_params.m_pivot)
+ /* back substitution */
+ if (m_params.m_pivot)
+ {
+ for (std::size_t j = kN; j-- > 0; )
{
- for (std::size_t j = kN; j-- > 0; )
- {
- 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);
- }
+ T 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);
}
- else
+ }
+ else
+ {
+ for (std::size_t j = kN; j-- > 0; )
{
- for (std::size_t j = kN; j-- > 0; )
+ T tmp = 0;
+
+ const auto *p = m_terms[j]->m_nzrd.data();
+ const auto e = m_terms[j]->m_nzrd.size() - 1; /* exclude RHS element */
+ T * Aj = &A(j,0);
+ for (std::size_t k = 0; k < e; k++)
{
- FT tmp = 0;
- const auto &nzrd = m_terms[j]->m_nzrd;
- 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);
+ const auto pk = p[k];
+ tmp += Aj[pk] * x[pk];
}
+ x[j] = (RHS(j) - tmp) / A(j,j);
}
}
+}
- template <typename FT, int SIZE>
- unsigned matrix_solver_direct_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphson)
- {
- 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;
- }
+template <std::size_t m_N, std::size_t storage_N>
+unsigned matrix_solver_direct_t<m_N, storage_N>::solve_non_dynamic(const bool newton_raphson)
+{
+ nl_double new_V[storage_N]; // = { 0.0 };
- 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);
+ this->LE_solve();
+ this->LE_back_subst(new_V);
- this->m_stat_calculations++;
- return this->solve_non_dynamic(newton_raphson);
- }
+ const nl_double err = (newton_raphson ? delta(new_V) : 0.0);
+ store(new_V);
+ return (err > this->m_params.m_accuracy) ? 2 : 1;
+}
- 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)
+template <std::size_t m_N, std::size_t storage_N>
+inline unsigned matrix_solver_direct_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson)
+{
+ build_LE_A<matrix_solver_direct_t>();
+ build_LE_RHS<matrix_solver_direct_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_direct_t<m_N, storage_N>::matrix_solver_direct_t(netlist_t &anetlist, const pstring &name,
+ const solver_parameters_t *params, const std::size_t size)
+: matrix_solver_t(anetlist, name, ASCENDING, params)
+, m_dim(size)
+{
+#if (NL_USE_DYNAMIC_ALLOCATION)
+ m_A.resize(N() * m_pitch);
+ //m_A = plib::palloc_array<nl_ext_double>(N() * m_pitch);
+#endif
+ for (unsigned k = 0; k < N(); k++)
{
+ m_last_RHS[k] = 0.0;
}
+}
- 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)
+template <std::size_t m_N, std::size_t storage_N>
+matrix_solver_direct_t<m_N, storage_N>::matrix_solver_direct_t(netlist_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_dim(size)
+{
+#if (NL_USE_DYNAMIC_ALLOCATION)
+ m_A.resize(N() * m_pitch);
+ //m_A = plib::palloc_array<nl_ext_double>(N() * m_pitch);
+#endif
+ for (unsigned k = 0; k < N(); k++)
{
+ m_last_RHS[k] = 0.0;
}
+}
-} // namespace devices
+ } //namespace devices
} // namespace netlist
#endif /* NLD_MS_DIRECT_H_ */