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.h222
1 files changed, 104 insertions, 118 deletions
diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h
index bdcec31ad44..6911e6aa360 100644
--- a/src/lib/netlist/solver/nld_ms_direct.h
+++ b/src/lib/netlist/solver/nld_ms_direct.h
@@ -13,11 +13,18 @@
#include "solver/nld_solver.h"
#include "solver/vector_base.h"
+/* 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 (0)
+
+
NETLIB_NAMESPACE_DEVICES_START()
//#define nl_ext_double __float128 // slow, very slow
//#define nl_ext_double long double // slightly slower
-#define nl_ext_double double
+#define nl_ext_double nl_double
template <unsigned m_N, unsigned _storage_N>
class matrix_solver_direct_t: public matrix_solver_t
@@ -32,48 +39,57 @@ public:
virtual void vsetup(analog_net_t::list_t &nets) override;
virtual void reset() override { matrix_solver_t::reset(); }
- ATTR_HOT inline unsigned N() const { if (m_N == 0) return m_dim; else return m_N; }
-
- ATTR_HOT inline int vsolve_non_dynamic(const bool newton_raphson);
-
protected:
virtual void add_term(int net_idx, terminal_t *term) override;
+ virtual int vsolve_non_dynamic(const bool newton_raphson) override;
+ int solve_non_dynamic(const bool newton_raphson);
- ATTR_HOT virtual nl_double vsolve() override;
+ inline const unsigned N() const { if (m_N == 0) return m_dim; else return m_N; }
- ATTR_HOT int solve_non_dynamic(const bool newton_raphson);
- ATTR_HOT void build_LE_A();
- ATTR_HOT void build_LE_RHS(nl_double * RESTRICT rhs);
- ATTR_HOT void LE_solve();
- ATTR_HOT void LE_back_subst(nl_double * RESTRICT x);
+ void build_LE_A();
+ void build_LE_RHS();
+ void LE_solve();
- /* Full LU back substitution, not used currently, in for future use */
+ template <typename T>
+ void LE_back_subst(T * RESTRICT x);
- ATTR_HOT void LE_back_subst_full(nl_double * RESTRICT x);
+ template <typename T>
+ T delta(const T * RESTRICT V);
- ATTR_HOT nl_double delta(const nl_double * RESTRICT V);
- ATTR_HOT void store(const nl_double * RESTRICT V);
+ template <typename T>
+ void store(const T * RESTRICT V);
- /* bring the whole system to the current time
- * Don't schedule a new calculation time. The recalculation has to be
- * triggered by the caller after the netlist element was changed.
- */
- ATTR_HOT nl_double compute_next_timestep();
+ virtual netlist_time compute_next_timestep() override;
+#if (NL_USE_DYNAMIC_ALLOCATION)
template <typename T1, typename T2>
- inline nl_ext_double &A(const T1 r, const T2 c) { return m_A[r][c]; }
-
- ATTR_ALIGN nl_double m_RHS[_storage_N];
+ inline nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r * m_pitch + c]; }
+ template <typename T1>
+ inline nl_ext_double &RHS(const T1 &r) { return m_A[r * m_pitch + N()]; }
+#else
+ template <typename T1, typename T2>
+ inline nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r][c]; }
+ template <typename T1>
+ inline nl_ext_double &RHS(const T1 &r) { return m_A[r][N()]; }
+#endif
ATTR_ALIGN nl_double m_last_RHS[_storage_N]; // right hand side - contains currents
ATTR_ALIGN nl_double m_last_V[_storage_N];
- terms_t **m_terms;
+ terms_t * m_terms[_storage_N];
terms_t *m_rails_temp;
private:
- ATTR_ALIGN nl_ext_double m_A[_storage_N][((_storage_N + 7) / 8) * 8];
+ static const std::size_t m_pitch = (((_storage_N + 1) + 7) / 8) * 8;
+#if (NL_USE_DYNAMIC_ALLOCATION)
+ ATTR_ALIGN nl_ext_double * RESTRICT m_A;
+#else
+ ATTR_ALIGN nl_ext_double m_A[_storage_N][m_pitch];
+ ATTR_ALIGN nl_ext_double m_B[_storage_N][m_pitch];
+#endif
+ //ATTR_ALIGN nl_ext_double m_RHSx[_storage_N];
const unsigned m_dim;
+
};
// ----------------------------------------------------------------------------------------
@@ -87,12 +103,14 @@ matrix_solver_direct_t<m_N, _storage_N>::~matrix_solver_direct_t()
{
pfree(m_terms[k]);
}
- pfree_array(m_terms);
pfree_array(m_rails_temp);
+#if (NL_USE_DYNAMIC_ALLOCATION)
+ pfree_array(m_A);
+#endif
}
template <unsigned m_N, unsigned _storage_N>
-ATTR_HOT nl_double matrix_solver_direct_t<m_N, _storage_N>::compute_next_timestep()
+netlist_time matrix_solver_direct_t<m_N, _storage_N>::compute_next_timestep()
{
nl_double new_solver_timestep = m_params.m_max_timestep;
@@ -121,13 +139,15 @@ ATTR_HOT nl_double matrix_solver_direct_t<m_N, _storage_N>::compute_next_timeste
if (new_net_timestep < new_solver_timestep)
new_solver_timestep = new_net_timestep;
+
+ m_last_V[k] = n->Q_Analog();
}
if (new_solver_timestep < m_params.m_min_timestep)
new_solver_timestep = m_params.m_min_timestep;
}
//if (new_solver_timestep > 10.0 * hn)
// new_solver_timestep = 10.0 * hn;
- return new_solver_timestep;
+ return netlist_time::from_double(new_solver_timestep);
}
template <unsigned m_N, unsigned _storage_N>
@@ -166,7 +186,7 @@ ATTR_COLD void matrix_solver_direct_t<m_N, _storage_N>::vsetup(analog_net_t::lis
m_rails_temp[k].clear();
}
- matrix_solver_t::setup(nets);
+ matrix_solver_t::setup_base(nets);
for (unsigned k = 0; k < N(); k++)
{
@@ -249,19 +269,23 @@ ATTR_COLD void matrix_solver_direct_t<m_N, _storage_N>::vsetup(analog_net_t::lis
}
}
- for (unsigned j = 0; j < N(); j++)
+ for (unsigned i = 0; i < t->m_railstart; i++)
{
- for (unsigned i = 0; i < t->m_railstart; i++)
- {
- if (!t->m_nzrd.contains(other[i]) && other[i] >= (int) (k + 1))
- t->m_nzrd.push_back(other[i]);
- if (!t->m_nz.contains(other[i]))
- t->m_nz.push_back(other[i]);
- }
+ if (!t->m_nzrd.contains(other[i]) && other[i] >= (int) (k + 1))
+ t->m_nzrd.push_back(other[i]);
+ if (!t->m_nz.contains(other[i]))
+ t->m_nz.push_back(other[i]);
}
+
+ /* Add RHS element */
+ if (!t->m_nzrd.contains(N()))
+ t->m_nzrd.push_back(N());
+
+ /* and sort */
psort_list(t->m_nzrd);
t->m_nz.push_back(k); // add diagonal
+
psort_list(t->m_nz);
}
@@ -304,7 +328,6 @@ ATTR_COLD void matrix_solver_direct_t<m_N, _storage_N>::vsetup(analog_net_t::lis
/*
* save states
*/
- save(NLNAME(m_RHS));
save(NLNAME(m_last_RHS));
save(NLNAME(m_last_V));
@@ -312,6 +335,8 @@ ATTR_COLD void matrix_solver_direct_t<m_N, _storage_N>::vsetup(analog_net_t::lis
{
pstring num = pfmt("{1}")(k);
+ save(RHS(k), "RHS" + num);
+
save(m_terms[k]->go(),"GO" + num, m_terms[k]->count());
save(m_terms[k]->gt(),"GT" + num, m_terms[k]->count());
save(m_terms[k]->Idr(),"IDR" + num , m_terms[k]->count());
@@ -321,7 +346,7 @@ ATTR_COLD void matrix_solver_direct_t<m_N, _storage_N>::vsetup(analog_net_t::lis
template <unsigned m_N, unsigned _storage_N>
-ATTR_HOT void matrix_solver_direct_t<m_N, _storage_N>::build_LE_A()
+void matrix_solver_direct_t<m_N, _storage_N>::build_LE_A()
{
const unsigned iN = N();
for (unsigned k = 0; k < iN; k++)
@@ -350,7 +375,7 @@ ATTR_HOT void matrix_solver_direct_t<m_N, _storage_N>::build_LE_A()
}
template <unsigned m_N, unsigned _storage_N>
-ATTR_HOT void matrix_solver_direct_t<m_N, _storage_N>::build_LE_RHS(nl_double * RESTRICT rhs)
+void matrix_solver_direct_t<m_N, _storage_N>::build_LE_RHS()
{
const unsigned iN = N();
for (unsigned k = 0; k < iN; k++)
@@ -370,12 +395,12 @@ ATTR_HOT void matrix_solver_direct_t<m_N, _storage_N>::build_LE_RHS(nl_double *
//rhsk = rhsk + go[i] * terms[i]->m_otherterm->net().as_analog().Q_Analog();
rhsk_b = rhsk_b + go[i] * *other_cur_analog[i];
- rhs[k] = rhsk_a + rhsk_b;
+ RHS(k) = rhsk_a + rhsk_b;
}
}
template <unsigned m_N, unsigned _storage_N>
-ATTR_HOT void matrix_solver_direct_t<m_N, _storage_N>::LE_solve()
+void matrix_solver_direct_t<m_N, _storage_N>::LE_solve()
{
const unsigned kN = N();
@@ -395,10 +420,10 @@ ATTR_HOT void matrix_solver_direct_t<m_N, _storage_N>::LE_solve()
if (maxrow != i)
{
/* Swap the maxrow and ith row */
- for (unsigned k = 0; k < kN; k++) {
+ for (unsigned k = 0; k < kN + 1; k++) {
std::swap(A(i,k), A(maxrow,k));
}
- std::swap(m_RHS[i], m_RHS[maxrow]);
+ //std::swap(RHS(i), RHS(maxrow));
}
/* FIXME: Singular matrix? */
const nl_double f = 1.0 / A(i,i);
@@ -410,14 +435,18 @@ ATTR_HOT void matrix_solver_direct_t<m_N, _storage_N>::LE_solve()
const nl_double f1 = - A(j,i) * f;
if (f1 != NL_FCONST(0.0))
{
- nl_double * RESTRICT pi = &m_A[i][i+1];
- nl_double * RESTRICT pj = &m_A[j][i+1];
+ const nl_double * RESTRICT pi = &A(i,i+1);
+ nl_double * RESTRICT pj = &A(j,i+1);
+#if 1
+ vec_add_mult_scalar(kN-i,pi,f1,pj);
+#else
vec_add_mult_scalar(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;
- m_RHS[j] += m_RHS[i] * f1;
+ RHS(j) += RHS(i) * f1;
+#endif
}
}
}
@@ -436,23 +465,18 @@ ATTR_HOT void matrix_solver_direct_t<m_N, _storage_N>::LE_solve()
{
const unsigned j = pb[jb];
const nl_double f1 = - A(j,i) * f;
-#if 0
- nl_double * RESTRICT pi = &m_A[i][i+1];
- nl_double * RESTRICT pj = &m_A[j][i+1];
- vec_add_mult_scalar(kN-i-1,pi,f1,pj);
-#else
for (unsigned k = 0; k < e; k++)
A(j,p[k]) += A(i,p[k]) * f1;
-#endif
- m_RHS[j] += m_RHS[i] * f1;
+ //RHS(j) += RHS(i) * f1;
}
}
}
}
template <unsigned m_N, unsigned _storage_N>
-ATTR_HOT void matrix_solver_direct_t<m_N, _storage_N>::LE_back_subst(
- nl_double * RESTRICT x)
+template <typename T>
+void matrix_solver_direct_t<m_N, _storage_N>::LE_back_subst(
+ T * RESTRICT x)
{
const unsigned kN = N();
@@ -461,68 +485,36 @@ ATTR_HOT void matrix_solver_direct_t<m_N, _storage_N>::LE_back_subst(
{
for (int j = kN - 1; j >= 0; j--)
{
- nl_double tmp = 0;
+ T tmp = 0;
for (unsigned k = j+1; k < kN; k++)
tmp += A(j,k) * x[k];
- x[j] = (m_RHS[j] - tmp) / A(j,j);
+ x[j] = (RHS(j) - tmp) / A(j,j);
}
}
else
{
for (int j = kN - 1; j >= 0; j--)
{
- nl_double tmp = 0;
+ T tmp = 0;
const unsigned *p = m_terms[j]->m_nzrd.data();
- const unsigned e = m_terms[j]->m_nzrd.size();
+ const unsigned e = m_terms[j]->m_nzrd.size() - 1; /* exclude RHS element */
for (unsigned k = 0; k < e; k++)
{
const unsigned pk = p[k];
tmp += A(j,pk) * x[pk];
}
- x[j] = (m_RHS[j] - tmp) / A(j,j);
+ x[j] = (RHS(j) - tmp) / A(j,j);
}
}
}
-template <unsigned m_N, unsigned _storage_N>
-ATTR_HOT void matrix_solver_direct_t<m_N, _storage_N>::LE_back_subst_full(
- nl_double * RESTRICT x)
-{
- const unsigned kN = N();
-
- /* back substitution */
-
- // int ip;
- // ii=-1
-
- //for (int i=0; i < kN; i++)
- // x[i] = m_RHS[i];
-
- for (int i=0; i < kN; i++)
- {
- //ip=indx[i]; USE_PIVOT_SEARCH
- //sum=b[ip];
- //b[ip]=b[i];
- double sum=m_RHS[i];//x[i];
- for (int j=0; j < i; j++)
- sum -= A(i,j) * x[j];
- x[i]=sum;
- }
- for (int i=kN-1; i >= 0; i--)
- {
- double sum=x[i];
- for (int j = i+1; j < kN; j++)
- sum -= A(i,j)*x[j];
- x[i] = sum / A(i,i);
- }
-
-}
template <unsigned m_N, unsigned _storage_N>
-ATTR_HOT nl_double matrix_solver_direct_t<m_N, _storage_N>::delta(
- const nl_double * RESTRICT V)
+template <typename T>
+T matrix_solver_direct_t<m_N, _storage_N>::delta(
+ const T * RESTRICT V)
{
/* FIXME: Ideally we should also include currents (RHS) here. This would
* need a revaluation of the right hand side after voltages have been updated
@@ -530,15 +522,16 @@ ATTR_HOT nl_double matrix_solver_direct_t<m_N, _storage_N>::delta(
*/
const unsigned iN = this->N();
- nl_double cerr = 0;
+ T cerr = 0;
for (unsigned i = 0; i < iN; i++)
- cerr = std::max(cerr, nl_math::abs(V[i] - this->m_nets[i]->m_cur_Analog));
+ cerr = std::fmax(cerr, nl_math::abs(V[i] - (T) this->m_nets[i]->m_cur_Analog));
return cerr;
}
template <unsigned m_N, unsigned _storage_N>
-ATTR_HOT void matrix_solver_direct_t<m_N, _storage_N>::store(
- const nl_double * RESTRICT V)
+template <typename T>
+void matrix_solver_direct_t<m_N, _storage_N>::store(
+ const T * RESTRICT V)
{
for (unsigned i = 0, iN=N(); i < iN; i++)
{
@@ -546,19 +539,13 @@ ATTR_HOT void matrix_solver_direct_t<m_N, _storage_N>::store(
}
}
-template <unsigned m_N, unsigned _storage_N>
-ATTR_HOT nl_double matrix_solver_direct_t<m_N, _storage_N>::vsolve()
-{
- this->solve_base(this);
- return this->compute_next_timestep();
-}
-
template <unsigned m_N, unsigned _storage_N>
-ATTR_HOT int matrix_solver_direct_t<m_N, _storage_N>::solve_non_dynamic(ATTR_UNUSED const bool newton_raphson)
+int matrix_solver_direct_t<m_N, _storage_N>::solve_non_dynamic(ATTR_UNUSED const bool newton_raphson)
{
nl_double new_V[_storage_N]; // = { 0.0 };
+ this->LE_solve();
this->LE_back_subst(new_V);
if (newton_raphson)
@@ -577,18 +564,15 @@ ATTR_HOT int matrix_solver_direct_t<m_N, _storage_N>::solve_non_dynamic(ATTR_UNU
}
template <unsigned m_N, unsigned _storage_N>
-ATTR_HOT inline int matrix_solver_direct_t<m_N, _storage_N>::vsolve_non_dynamic(const bool newton_raphson)
+inline int matrix_solver_direct_t<m_N, _storage_N>::vsolve_non_dynamic(const bool newton_raphson)
{
this->build_LE_A();
- this->build_LE_RHS(m_last_RHS);
+ this->build_LE_RHS();
for (unsigned i=0, iN=N(); i < iN; i++)
- m_RHS[i] = m_last_RHS[i];
-
- this->LE_solve();
+ m_last_RHS[i] = RHS(i);
this->m_stat_calculations++;
-
return this->solve_non_dynamic(newton_raphson);
}
@@ -597,9 +581,10 @@ matrix_solver_direct_t<m_N, _storage_N>::matrix_solver_direct_t(const solver_par
: matrix_solver_t(GAUSSIAN_ELIMINATION, params)
, m_dim(size)
{
- m_terms = palloc_array(terms_t *, N());
m_rails_temp = palloc_array(terms_t, N());
-
+#if (NL_USE_DYNAMIC_ALLOCATION)
+ m_A = palloc_array(nl_ext_double, N() * m_pitch);
+#endif
for (unsigned k = 0; k < N(); k++)
{
m_terms[k] = palloc(terms_t);
@@ -613,9 +598,10 @@ matrix_solver_direct_t<m_N, _storage_N>::matrix_solver_direct_t(const eSolverTyp
: matrix_solver_t(type, params)
, m_dim(size)
{
- m_terms = palloc_array(terms_t *, N());
m_rails_temp = palloc_array(terms_t, N());
-
+#if (NL_USE_DYNAMIC_ALLOCATION)
+ m_A = palloc_array(nl_ext_double, N() * m_pitch);
+#endif
for (unsigned k = 0; k < N(); k++)
{
m_terms[k] = palloc(terms_t);