summaryrefslogtreecommitdiffstatshomepage
path: root/src/lib/netlist/solver/nld_ms_sm.h
diff options
context:
space:
mode:
Diffstat (limited to 'src/lib/netlist/solver/nld_ms_sm.h')
-rw-r--r--src/lib/netlist/solver/nld_ms_sm.h416
1 files changed, 199 insertions, 217 deletions
diff --git a/src/lib/netlist/solver/nld_ms_sm.h b/src/lib/netlist/solver/nld_ms_sm.h
index ed95cd2dc03..d85ab0044f4 100644
--- a/src/lib/netlist/solver/nld_ms_sm.h
+++ b/src/lib/netlist/solver/nld_ms_sm.h
@@ -33,293 +33,275 @@
#ifndef NLD_MS_SM_H_
#define NLD_MS_SM_H_
-#include <algorithm>
-
-#include "nld_solver.h"
#include "nld_matrix_solver.h"
-#include "vector_base.h"
+#include "nld_solver.h"
+#include "plib/vector_ops.h"
+
+#include <algorithm>
namespace netlist
{
- namespace devices
- {
-//#define nl_ext_double _float128 // slow, very slow
-//#define nl_ext_double long double // slightly slower
-#define nl_ext_double nl_double
-
-template <std::size_t m_N, std::size_t storage_N>
-class matrix_solver_sm_t: public matrix_solver_t
+namespace devices
{
- friend class matrix_solver_t;
-
-public:
- matrix_solver_sm_t(netlist_t &anetlist, const pstring &name,
- const solver_parameters_t *params, const std::size_t size);
-
- virtual ~matrix_solver_sm_t() override;
-
- virtual void vsetup(analog_net_t::list_t &nets) override;
- virtual void reset() override { matrix_solver_t::reset(); }
-
-protected:
- virtual unsigned vsolve_non_dynamic(const bool newton_raphson) override;
- unsigned solve_non_dynamic(const bool newton_raphson);
-
- constexpr std::size_t N() const { return (m_N == 0) ? m_dim : m_N; }
-
- void LE_invert();
-
- template <typename T>
- void LE_compute_x(T * RESTRICT x);
+ template <typename FT, int SIZE>
+ class matrix_solver_sm_t: public matrix_solver_t
+ {
+ friend class matrix_solver_t;
+ public:
- template <typename T1, typename T2>
- nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r][c]; }
- template <typename T1, typename T2>
- nl_ext_double &W(const T1 &r, const T2 &c) { return m_W[r][c]; }
- template <typename T1, typename T2>
- nl_ext_double &Ainv(const T1 &r, const T2 &c) { return m_Ainv[r][c]; }
- template <typename T1>
- nl_ext_double &RHS(const T1 &r) { return m_RHS[r]; }
+ using float_ext_type = FT;
+ using float_type = FT;
+ // 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);
- template <typename T1, typename T2>
- nl_ext_double &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; }
- template <typename T1, typename T2>
- nl_ext_double &lAinv(const T1 &r, const T2 &c) { return m_lAinv[r][c]; }
+ void vsetup(analog_net_t::list_t &nets) override;
+ void reset() override { matrix_solver_t::reset(); }
- nl_double m_last_RHS[storage_N]; // right hand side - contains currents
+ protected:
+ unsigned vsolve_non_dynamic(const bool newton_raphson) override;
+ unsigned solve_non_dynamic(const bool newton_raphson);
-private:
- static constexpr std::size_t m_pitch = ((( storage_N) + 7) / 8) * 8;
- nl_ext_double m_A[storage_N][m_pitch];
- nl_ext_double m_Ainv[storage_N][m_pitch];
- nl_ext_double m_W[storage_N][m_pitch];
- nl_ext_double m_RHS[storage_N]; // right hand side - contains currents
+ constexpr std::size_t size() const { return m_dim; }
- nl_ext_double m_lA[storage_N][m_pitch];
- nl_ext_double m_lAinv[storage_N][m_pitch];
+ void LE_invert();
- //nl_ext_double m_RHSx[storage_N];
+ template <typename T>
+ void LE_compute_x(T * x);
- const std::size_t m_dim;
- std::size_t m_cnt;
-};
+ template <typename T1, typename T2>
+ float_ext_type &A(const T1 &r, const T2 &c) { return m_A[r][c]; }
+ template <typename T1, typename T2>
+ float_ext_type &W(const T1 &r, const T2 &c) { return m_W[r][c]; }
+ 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]; }
-// ----------------------------------------------------------------------------------------
-// matrix_solver_direct
-// ----------------------------------------------------------------------------------------
-template <std::size_t m_N, std::size_t storage_N>
-matrix_solver_sm_t<m_N, storage_N>::~matrix_solver_sm_t()
-{
-}
+ template <typename T1, typename T2>
+ float_ext_type &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; }
+ template <typename T1, typename T2>
+ float_ext_type &lAinv(const T1 &r, const T2 &c) { return m_lAinv[r][c]; }
-template <std::size_t m_N, std::size_t storage_N>
-void matrix_solver_sm_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets)
-{
- matrix_solver_t::setup_base(nets);
+ private:
+ 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];
+ float_ext_type m_RHS[storage_N]; // right hand side - contains currents
- netlist().save(*this, m_last_RHS, "m_last_RHS");
+ float_ext_type m_lA[storage_N][m_pitch];
+ float_ext_type m_lAinv[storage_N][m_pitch];
- for (unsigned k = 0; k < N(); k++)
- netlist().save(*this, RHS(k), plib::pfmt("RHS.{1}")(k));
-}
+ //float_ext_type m_RHSx[storage_N];
+ const std::size_t m_dim;
+ std::size_t m_cnt;
+ };
-template <std::size_t m_N, std::size_t storage_N>
-void matrix_solver_sm_t<m_N, storage_N>::LE_invert()
-{
- const std::size_t kN = N();
+ // ----------------------------------------------------------------------------------------
+ // matrix_solver_direct
+ // ----------------------------------------------------------------------------------------
- for (std::size_t i = 0; i < kN; i++)
+ template <typename FT, int SIZE>
+ void matrix_solver_sm_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets)
{
- 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,i) = 1.0;
+ 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));
}
- /* down */
- for (std::size_t i = 0; i < kN; i++)
- {
- /* FIXME: Singular matrix? */
- const nl_double f = 1.0 / W(i,i);
- const auto * RESTRICT const p = m_terms[i]->m_nzrd.data();
- const std::size_t e = m_terms[i]->m_nzrd.size();
- /* Eliminate column i from row j */
+ template <typename FT, int SIZE>
+ void matrix_solver_sm_t<FT, SIZE>::LE_invert()
+ {
+ const std::size_t kN = size();
- const auto * RESTRICT const pb = m_terms[i]->m_nzbd.data();
- const std::size_t eb = m_terms[i]->m_nzbd.size();
- for (std::size_t jb = 0; jb < eb; jb++)
+ for (std::size_t i = 0; i < kN; i++)
{
- const unsigned j = pb[jb];
- const nl_double f1 = - W(j,i) * f;
- if (f1 != 0.0)
+ for (std::size_t j = 0; j < kN; j++)
{
- for (std::size_t k = 0; k < e; k++)
- W(j,p[k]) += W(i,p[k]) * f1;
- for (std::size_t k = 0; k <= i; k ++)
- Ainv(j,k) += Ainv(i,k) * f1;
+ W(i,j) = lA(i,j) = A(i,j);
+ Ainv(i,j) = 0.0;
}
+ Ainv(i,i) = 1.0;
}
- }
- /* up */
- for (std::size_t i = kN; i-- > 0; )
- {
- /* FIXME: Singular matrix? */
- const nl_double f = 1.0 / W(i,i);
- for (std::size_t j = i; j-- > 0; )
+ /* down */
+ for (std::size_t i = 0; i < kN; i++)
{
- const nl_double f1 = - W(j,i) * f;
- if (f1 != 0.0)
+ /* 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();
+
+ /* 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();
+ for (std::size_t jb = 0; jb < eb; jb++)
{
- for (std::size_t k = i; k < kN; k++)
- W(j,k) += W(i,k) * f1;
- for (std::size_t k = 0; k < kN; k++)
- Ainv(j,k) += Ainv(i,k) * f1;
+ const unsigned j = pb[jb];
+ const float_type f1 = - W(j,i) * f;
+ if (f1 != 0.0)
+ {
+ for (std::size_t k = 0; k < e; k++)
+ W(j,p[k]) += W(i,p[k]) * f1;
+ for (std::size_t k = 0; k <= i; k ++)
+ Ainv(j,k) += Ainv(i,k) * f1;
+ }
}
}
- for (std::size_t k = 0; k < kN; k++)
+ /* up */
+ for (std::size_t i = kN; i-- > 0; )
{
- Ainv(i,k) *= f;
- lAinv(i,k) = Ainv(i,k);
+ /* FIXME: Singular matrix? */
+ const float_type f = 1.0 / W(i,i);
+ for (std::size_t j = i; j-- > 0; )
+ {
+ const float_type f1 = - W(j,i) * f;
+ if (f1 != 0.0)
+ {
+ for (std::size_t k = i; k < kN; k++)
+ W(j,k) += W(i,k) * f1;
+ for (std::size_t k = 0; k < kN; k++)
+ Ainv(j,k) += Ainv(i,k) * f1;
+ }
+ }
+ for (std::size_t k = 0; k < kN; k++)
+ {
+ Ainv(i,k) *= f;
+ lAinv(i,k) = Ainv(i,k);
+ }
}
}
-}
-
-template <std::size_t m_N, std::size_t storage_N>
-template <typename T>
-void matrix_solver_sm_t<m_N, storage_N>::LE_compute_x(
- T * RESTRICT x)
-{
- const std::size_t kN = N();
- for (std::size_t i=0; i<kN; i++)
- x[i] = 0.0;
-
- for (std::size_t k=0; k<kN; k++)
+ template <typename FT, int SIZE>
+ template <typename T>
+ void matrix_solver_sm_t<FT, SIZE>::LE_compute_x(
+ T * x)
{
- const nl_double f = RHS(k);
+ const std::size_t kN = size();
for (std::size_t i=0; i<kN; i++)
- x[i] += Ainv(i,k) * f;
- }
-}
-
-
-template <std::size_t m_N, std::size_t storage_N>
-unsigned matrix_solver_sm_t<m_N, storage_N>::solve_non_dynamic(const bool newton_raphson)
-{
- static constexpr const bool incremental = true;
- const std::size_t iN = N();
+ x[i] = 0.0;
- nl_double new_V[storage_N]; // = { 0.0 };
+ for (std::size_t k=0; k<kN; k++)
+ {
+ const float_type f = RHS(k);
- if ((m_cnt % 50) == 0)
- {
- /* complete calculation */
- this->LE_invert();
+ for (std::size_t i=0; i<kN; i++)
+ x[i] += Ainv(i,k) * f;
+ }
}
- else
+
+ template <typename FT, int SIZE>
+ unsigned matrix_solver_sm_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphson)
{
- if (!incremental)
+ static constexpr const bool incremental = true;
+ const std::size_t iN = size();
+
+ float_type new_V[storage_N]; // = { 0.0 };
+
+ if ((m_cnt % 50) == 0)
{
- for (std::size_t row = 0; row < iN; row ++)
- for (std::size_t k = 0; k < iN; k++)
- Ainv(row,k) = lAinv(row, k);
+ /* complete calculation */
+ this->LE_invert();
}
- for (std::size_t row = 0; row < iN; row ++)
+ else
{
- nl_double v[m_pitch] = {0};
- std::size_t cols[m_pitch];
- std::size_t colcount = 0;
-
- auto &nz = m_terms[row]->m_nz;
- for (unsigned & col : nz)
+ if (!incremental)
{
- v[col] = A(row,col) - lA(row,col);
- if (incremental)
- lA(row,col) = A(row,col);
- if (v[col] != 0.0)
- cols[colcount++] = col;
+ for (std::size_t row = 0; row < iN; row ++)
+ for (std::size_t k = 0; k < iN; k++)
+ Ainv(row,k) = lAinv(row, k);
}
-
- if (colcount > 0)
+ for (std::size_t row = 0; row < iN; row ++)
{
- nl_double lamba = 0.0;
- nl_double w[m_pitch] = {0};
-
- nl_double z[m_pitch];
- /* compute w and lamba */
- for (std::size_t i = 0; i < iN; i++)
- z[i] = Ainv(i, row); /* u is row'th column */
+ float_type v[m_pitch] = {0};
+ std::size_t cols[m_pitch];
+ std::size_t colcount = 0;
- for (std::size_t j = 0; j < colcount; j++)
- lamba += v[cols[j]] * z[cols[j]];
-
- for (std::size_t j=0; j<colcount; j++)
+ auto &nz = m_terms[row]->m_nz;
+ for (unsigned & col : nz)
{
- std::size_t col = cols[j];
- nl_double f = v[col];
- for (std::size_t k = 0; k < iN; k++)
- w[k] += Ainv(col,k) * f; /* Transpose(Ainv) * v */
+ v[col] = A(row,col) - lA(row,col);
+ if (incremental)
+ lA(row,col) = A(row,col);
+ if (v[col] != 0.0)
+ cols[colcount++] = col;
}
- lamba = -1.0 / (1.0 + lamba);
- for (std::size_t i=0; i<iN; i++)
+ if (colcount > 0)
{
- const nl_double f = lamba * z[i];
- if (f != 0.0)
+ float_type lamba = 0.0;
+ float_type w[m_pitch] = {0};
+
+ float_type z[m_pitch];
+ /* compute w and lamba */
+ for (std::size_t i = 0; i < iN; i++)
+ 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<colcount; j++)
+ {
+ std::size_t col = cols[j];
+ float_type f = v[col];
for (std::size_t k = 0; k < iN; k++)
- Ainv(i,k) += f * w[k];
+ w[k] += Ainv(col,k) * f; /* Transpose(Ainv) * v */
+ }
+
+ lamba = -1.0 / (1.0 + lamba);
+ for (std::size_t i=0; i<iN; i++)
+ {
+ const float_type f = lamba * z[i];
+ if (f != 0.0)
+ for (std::size_t k = 0; k < iN; k++)
+ Ainv(i,k) += f * w[k];
+ }
}
- }
+ }
}
- }
- m_cnt++;
+ m_cnt++;
- this->LE_compute_x(new_V);
+ this->LE_compute_x(new_V);
- const nl_double err = (newton_raphson ? delta(new_V) : 0.0);
- store(new_V);
- return (err > this->m_params.m_accuracy) ? 2 : 1;
-}
+ const float_type err = (newton_raphson ? delta(new_V) : 0.0);
+ store(new_V);
+ return (err > this->m_params.m_accuracy) ? 2 : 1;
+ }
-template <std::size_t m_N, std::size_t storage_N>
-inline unsigned matrix_solver_sm_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson)
-{
- build_LE_A<matrix_solver_sm_t>();
- build_LE_RHS<matrix_solver_sm_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_sm_t<m_N, storage_N>::matrix_solver_sm_t(netlist_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)
-{
- for (std::size_t k = 0; k < N(); k++)
+ template <typename FT, int SIZE>
+ unsigned matrix_solver_sm_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson)
+ {
+ this->build_LE_A(*this);
+ this->build_LE_RHS(*this);
+
+ this->m_stat_calculations++;
+ return this->solve_non_dynamic(newton_raphson);
+ }
+
+ 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)
{
- m_last_RHS[k] = 0.0;
}
-}
- } //namespace devices
+} // namespace devices
} // namespace netlist
#endif /* NLD_MS_DIRECT_H_ */