summaryrefslogtreecommitdiffstatshomepage
path: root/src/lib/netlist/solver
diff options
context:
space:
mode:
author couriersud <couriersud@arcor.de>2016-04-09 20:06:25 +0200
committer couriersud <couriersud@arcor.de>2016-04-09 23:45:42 +0200
commit0d1e57cc40fb31dcc5f29588b1695d6ce1ab9c8a (patch)
tree7b98517a3f134477084b5e063a1c41b8fd20dee1 /src/lib/netlist/solver
parentdb1e7f6b076cd092fccd70f3311e40499725281c (diff)
Solver model simplification.
Diffstat (limited to 'src/lib/netlist/solver')
-rw-r--r--src/lib/netlist/solver/nld_matrix_solver.h44
-rw-r--r--src/lib/netlist/solver/nld_ms_direct.h137
-rw-r--r--src/lib/netlist/solver/nld_ms_direct_lu.h11
-rw-r--r--src/lib/netlist/solver/nld_ms_sm.h131
-rw-r--r--src/lib/netlist/solver/nld_ms_w.h136
-rw-r--r--src/lib/netlist/solver/nld_solver.cpp74
6 files changed, 133 insertions, 400 deletions
diff --git a/src/lib/netlist/solver/nld_matrix_solver.h b/src/lib/netlist/solver/nld_matrix_solver.h
index 715eab6756a..264f5826b87 100644
--- a/src/lib/netlist/solver/nld_matrix_solver.h
+++ b/src/lib/netlist/solver/nld_matrix_solver.h
@@ -16,8 +16,8 @@ class terms_t
{
P_PREVENT_COPYING(terms_t)
- public:
- ATTR_COLD terms_t() : m_railstart(0)
+public:
+ ATTR_COLD terms_t() : m_railstart(0), m_last_V(0)
{}
ATTR_COLD void clear()
@@ -48,6 +48,10 @@ class terms_t
pvector_t<unsigned> m_nz; /* all non zero for multiplication */
pvector_t<unsigned> m_nzrd; /* non zero right of the diagonal for elimination, may include RHS element */
pvector_t<unsigned> m_nzbd; /* non zero below of the diagonal for elimination */
+
+ /* state */
+ nl_double m_last_V;
+
private:
pvector_t<terminal_t *> m_term;
pvector_t<int> m_net_other;
@@ -55,6 +59,7 @@ private:
pvector_t<nl_double> m_gt;
pvector_t<nl_double> m_Idr;
pvector_t<nl_double *> m_other_curanalog;
+
};
class matrix_solver_t : public device_t
@@ -104,14 +109,23 @@ protected:
ATTR_COLD void setup_base(analog_net_t::list_t &nets);
void update_dynamic();
- virtual void add_term(int net_idx, terminal_t *term) = 0;
virtual void vsetup(analog_net_t::list_t &nets) = 0;
virtual int vsolve_non_dynamic(const bool newton_raphson) = 0;
- virtual netlist_time compute_next_timestep() = 0;
+ /* virtual */ netlist_time compute_next_timestep();
+ /* virtual */ void add_term(int net_idx, terminal_t *term);
+
+ template <typename T>
+ void store(const T * RESTRICT V);
+ template <typename T>
+ T delta(const T * RESTRICT V);
+
+ pvector_t<terms_t *> m_terms;
pvector_t<analog_net_t *> m_nets;
pvector_t<analog_output_t *> m_inps;
+ pvector_t<terms_t *> m_rails_temp;
+
int m_stat_calculations;
int m_stat_newton_raphson;
int m_stat_vsolver_calls;
@@ -138,6 +152,28 @@ private:
const eSolverType m_type;
};
+template <typename T>
+T matrix_solver_t::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
+ * and thus belong into a different calculation. This applies to all solvers.
+ */
+
+ const unsigned iN = this->m_terms.size();
+ T cerr = 0;
+ for (unsigned i = 0; i < iN; i++)
+ cerr = nl_math::max(cerr, nl_math::abs(V[i] - (T) this->m_nets[i]->m_cur_Analog));
+ return cerr;
+}
+
+template <typename T>
+void matrix_solver_t::store(const T * RESTRICT V)
+{
+ for (unsigned i = 0, iN=m_terms.size(); i < iN; i++)
+ this->m_nets[i]->m_cur_Analog = V[i];
+}
+
NETLIB_NAMESPACE_DEVICES_END()
diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h
index 77deb330544..b360c748019 100644
--- a/src/lib/netlist/solver/nld_ms_direct.h
+++ b/src/lib/netlist/solver/nld_ms_direct.h
@@ -129,7 +129,6 @@ public:
virtual void reset() override { matrix_solver_t::reset(); }
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);
@@ -142,14 +141,6 @@ protected:
template <typename T>
void LE_back_subst(T * RESTRICT x);
- template <typename T>
- T delta(const T * RESTRICT V);
-
- template <typename T>
- void store(const T * RESTRICT V);
-
- virtual netlist_time compute_next_timestep() override;
-
#if TEST_PARALLEL
int x_i[10];
int x_start[10];
@@ -169,10 +160,6 @@ protected:
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[_storage_N];
- terms_t *m_rails_temp;
private:
static const std::size_t m_pitch = (((_storage_N + 1) + 7) / 8) * 8;
@@ -196,109 +183,36 @@ private:
template <unsigned m_N, unsigned _storage_N>
matrix_solver_direct_t<m_N, _storage_N>::~matrix_solver_direct_t()
{
- for (unsigned k = 0; k < N(); k++)
- {
- pfree(m_terms[k]);
- }
- pfree_array(m_rails_temp);
#if (NL_USE_DYNAMIC_ALLOCATION)
pfree_array(m_A);
#endif
#if TEST_PARALLEL
thr_dispose();
#endif
-
-}
-
-template <unsigned m_N, unsigned _storage_N>
-netlist_time matrix_solver_direct_t<m_N, _storage_N>::compute_next_timestep()
-{
- nl_double new_solver_timestep = m_params.m_max_timestep;
-
- if (m_params.m_dynamic)
- {
- /*
- * FIXME: We should extend the logic to use either all nets or
- * only output nets.
- */
- for (unsigned k = 0, iN=N(); k < iN; k++)
- {
- analog_net_t *n = m_nets[k];
-
- const nl_double DD_n = (n->Q_Analog() - m_last_V[k]);
- const nl_double hn = current_timestep();
-
- nl_double DD2 = (DD_n / hn - n->m_DD_n_m_1 / n->m_h_n_m_1) / (hn + n->m_h_n_m_1);
- nl_double new_net_timestep;
-
- n->m_h_n_m_1 = hn;
- n->m_DD_n_m_1 = DD_n;
- if (nl_math::abs(DD2) > NL_FCONST(1e-30)) // avoid div-by-zero
- new_net_timestep = nl_math::sqrt(m_params.m_lte / nl_math::abs(NL_FCONST(0.5)*DD2));
- else
- new_net_timestep = m_params.m_max_timestep;
-
- 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 netlist_time::from_double(new_solver_timestep);
}
template <unsigned m_N, unsigned _storage_N>
-ATTR_COLD void matrix_solver_direct_t<m_N, _storage_N>::add_term(int k, terminal_t *term)
-{
- if (term->m_otherterm->net().isRailNet())
- {
- m_rails_temp[k].add(term, -1, false);
- }
- else
- {
- int ot = get_net_idx(&term->m_otherterm->net());
- if (ot>=0)
- {
- m_terms[k]->add(term, ot, true);
- }
- /* Should this be allowed ? */
- else // if (ot<0)
- {
- m_rails_temp[k].add(term, ot, true);
- log().fatal("found term with missing othernet {1}\n", term->name());
- }
- }
-}
-
-
-template <unsigned m_N, unsigned _storage_N>
ATTR_COLD void matrix_solver_direct_t<m_N, _storage_N>::vsetup(analog_net_t::list_t &nets)
{
if (m_dim < nets.size())
log().fatal("Dimension {1} less than {2}", m_dim, nets.size());
- for (unsigned k = 0; k < N(); k++)
- {
- m_terms[k]->clear();
- m_rails_temp[k].clear();
- }
-
matrix_solver_t::setup_base(nets);
for (unsigned k = 0; k < N(); k++)
{
m_terms[k]->m_railstart = m_terms[k]->count();
- for (unsigned i = 0; i < m_rails_temp[k].count(); i++)
- this->m_terms[k]->add(m_rails_temp[k].terms()[i], m_rails_temp[k].net_other()[i], false);
+ for (unsigned i = 0; i < m_rails_temp[k]->count(); i++)
+ this->m_terms[k]->add(m_rails_temp[k]->terms()[i], m_rails_temp[k]->net_other()[i], false);
- m_rails_temp[k].clear(); // no longer needed
+ m_rails_temp[k]->clear(); // no longer needed
m_terms[k]->set_pointers();
}
+ for (unsigned k = 0; k < N(); k++)
+ pfree(m_rails_temp[k]); // no longer needed
+
+ m_rails_temp.clear();
#if 1
/* Sort in descending order by number of connected matrix voltages.
@@ -430,13 +344,13 @@ ATTR_COLD void matrix_solver_direct_t<m_N, _storage_N>::vsetup(analog_net_t::lis
* save states
*/
save(NLNAME(m_last_RHS));
- save(NLNAME(m_last_V));
for (unsigned k = 0; k < N(); k++)
{
pstring num = pfmt("{1}")(k);
save(RHS(k), "RHS" + num);
+ save(m_terms[k]->m_last_V, "lastV" + num);
save(m_terms[k]->go(),"GO" + num, m_terms[k]->count());
save(m_terms[k]->gt(),"GT" + num, m_terms[k]->count());
@@ -683,35 +597,6 @@ void matrix_solver_direct_t<m_N, _storage_N>::LE_back_subst(
template <unsigned m_N, unsigned _storage_N>
-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
- * and thus belong into a different calculation. This applies to all solvers.
- */
-
- const unsigned iN = this->N();
- T cerr = 0;
- for (unsigned i = 0; i < iN; i++)
- 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>
-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++)
- {
- this->m_nets[i]->m_cur_Analog = V[i];
- }
-}
-
-
-template <unsigned m_N, unsigned _storage_N>
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 };
@@ -752,15 +637,12 @@ 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_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);
m_last_RHS[k] = 0.0;
- m_last_V[k] = 0.0;
}
#if TEST_PARALLEL
thr_initialize();
@@ -772,15 +654,12 @@ matrix_solver_direct_t<m_N, _storage_N>::matrix_solver_direct_t(const eSolverTyp
: matrix_solver_t(type, params)
, m_dim(size)
{
- 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);
m_last_RHS[k] = 0.0;
- m_last_V[k] = 0.0;
}
#if TEST_PARALLEL
thr_initialize();
diff --git a/src/lib/netlist/solver/nld_ms_direct_lu.h b/src/lib/netlist/solver/nld_ms_direct_lu.h
index 757970e5c07..e6ccf819653 100644
--- a/src/lib/netlist/solver/nld_ms_direct_lu.h
+++ b/src/lib/netlist/solver/nld_ms_direct_lu.h
@@ -144,7 +144,6 @@ protected:
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_rails_temp;
private:
@@ -161,11 +160,6 @@ private:
template <unsigned m_N, unsigned _storage_N>
matrix_solver_direct_t<m_N, _storage_N>::~matrix_solver_direct_t()
{
- for (unsigned k = 0; k < N(); k++)
- {
- pfree(m_terms[k]);
- }
- pfree_array(m_terms);
pfree_array(m_rails_temp);
}
@@ -608,14 +602,11 @@ matrix_solver_direct_t<m_N, _storage_N>::matrix_solver_direct_t(const solver_par
, m_dim(size)
, m_lp_fact(0)
{
- m_terms = palloc_array(terms_t *, N());
m_rails_temp = palloc_array(terms_t, N());
for (unsigned k = 0; k < N(); k++)
{
- m_terms[k] = palloc(terms_t);
m_last_RHS[k] = 0.0;
- m_last_V[k] = 0.0;
}
}
@@ -625,14 +616,12 @@ matrix_solver_direct_t<m_N, _storage_N>::matrix_solver_direct_t(const eSolverTyp
, m_dim(size)
, m_lp_fact(0)
{
- m_terms = palloc_array(terms_t *, N());
m_rails_temp = palloc_array(terms_t, N());
for (unsigned k = 0; k < N(); k++)
{
m_terms[k] = palloc(terms_t);
m_last_RHS[k] = 0.0;
- m_last_V[k] = 0.0;
}
}
diff --git a/src/lib/netlist/solver/nld_ms_sm.h b/src/lib/netlist/solver/nld_ms_sm.h
index c515fd1e80f..f9585b441f2 100644
--- a/src/lib/netlist/solver/nld_ms_sm.h
+++ b/src/lib/netlist/solver/nld_ms_sm.h
@@ -59,7 +59,6 @@ public:
virtual void reset() override { matrix_solver_t::reset(); }
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);
@@ -72,13 +71,6 @@ protected:
template <typename T>
void LE_compute_x(T * RESTRICT x);
- template <typename T>
- T delta(const T * RESTRICT V);
-
- template <typename T>
- void store(const T * RESTRICT V);
-
- virtual netlist_time compute_next_timestep() override;
template <typename T1, typename T2>
inline nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r][c]; }
@@ -96,10 +88,6 @@ protected:
inline nl_ext_double &lAinv(const T1 &r, const T2 &c) { return m_lAinv[r][c]; }
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[_storage_N];
- terms_t *m_rails_temp;
private:
static const std::size_t m_pitch = ((( _storage_N) + 7) / 8) * 8;
@@ -124,102 +112,26 @@ private:
template <unsigned m_N, unsigned _storage_N>
matrix_solver_sm_t<m_N, _storage_N>::~matrix_solver_sm_t()
{
- for (unsigned k = 0; k < N(); k++)
- {
- pfree(m_terms[k]);
- }
- pfree_array(m_rails_temp);
#if (NL_USE_DYNAMIC_ALLOCATION)
pfree_array(m_A);
#endif
}
template <unsigned m_N, unsigned _storage_N>
-netlist_time matrix_solver_sm_t<m_N, _storage_N>::compute_next_timestep()
-{
- nl_double new_solver_timestep = m_params.m_max_timestep;
-
- if (m_params.m_dynamic)
- {
- /*
- * FIXME: We should extend the logic to use either all nets or
- * only output nets.
- */
- for (unsigned k = 0, iN=N(); k < iN; k++)
- {
- analog_net_t *n = m_nets[k];
-
- const nl_double DD_n = (n->Q_Analog() - m_last_V[k]);
- const nl_double hn = current_timestep();
-
- nl_double DD2 = (DD_n / hn - n->m_DD_n_m_1 / n->m_h_n_m_1) / (hn + n->m_h_n_m_1);
- nl_double new_net_timestep;
-
- n->m_h_n_m_1 = hn;
- n->m_DD_n_m_1 = DD_n;
- if (nl_math::abs(DD2) > NL_FCONST(1e-30)) // avoid div-by-zero
- new_net_timestep = nl_math::sqrt(m_params.m_lte / nl_math::abs(NL_FCONST(0.5)*DD2));
- else
- new_net_timestep = m_params.m_max_timestep;
-
- 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 netlist_time::from_double(new_solver_timestep);
-}
-
-template <unsigned m_N, unsigned _storage_N>
-ATTR_COLD void matrix_solver_sm_t<m_N, _storage_N>::add_term(int k, terminal_t *term)
-{
- if (term->m_otherterm->net().isRailNet())
- {
- m_rails_temp[k].add(term, -1, false);
- }
- else
- {
- int ot = get_net_idx(&term->m_otherterm->net());
- if (ot>=0)
- {
- m_terms[k]->add(term, ot, true);
- }
- /* Should this be allowed ? */
- else // if (ot<0)
- {
- m_rails_temp[k].add(term, ot, true);
- log().fatal("found term with missing othernet {1}\n", term->name());
- }
- }
-}
-
-
-template <unsigned m_N, unsigned _storage_N>
ATTR_COLD void matrix_solver_sm_t<m_N, _storage_N>::vsetup(analog_net_t::list_t &nets)
{
if (m_dim < nets.size())
log().fatal("Dimension {1} less than {2}", m_dim, nets.size());
- for (unsigned k = 0; k < N(); k++)
- {
- m_terms[k]->clear();
- m_rails_temp[k].clear();
- }
-
matrix_solver_t::setup_base(nets);
for (unsigned k = 0; k < N(); k++)
{
m_terms[k]->m_railstart = m_terms[k]->count();
- for (unsigned i = 0; i < m_rails_temp[k].count(); i++)
- this->m_terms[k]->add(m_rails_temp[k].terms()[i], m_rails_temp[k].net_other()[i], false);
+ for (unsigned i = 0; i < m_rails_temp[k]->count(); i++)
+ this->m_terms[k]->add(m_rails_temp[k]->terms()[i], m_rails_temp[k]->net_other()[i], false);
- m_rails_temp[k].clear(); // no longer needed
+ m_rails_temp[k]->clear(); // no longer needed
m_terms[k]->set_pointers();
}
@@ -315,13 +227,13 @@ ATTR_COLD void matrix_solver_sm_t<m_N, _storage_N>::vsetup(analog_net_t::list_t
* save states
*/
save(NLNAME(m_last_RHS));
- save(NLNAME(m_last_V));
for (unsigned k = 0; k < N(); k++)
{
pstring num = pfmt("{1}")(k);
save(RHS(k), "RHS" + num);
+ save(m_terms[k]->m_last_V, "lastV" + num);
save(m_terms[k]->go(),"GO" + num, m_terms[k]->count());
save(m_terms[k]->gt(),"GT" + num, m_terms[k]->count());
@@ -469,35 +381,6 @@ void matrix_solver_sm_t<m_N, _storage_N>::LE_compute_x(
template <unsigned m_N, unsigned _storage_N>
-template <typename T>
-T matrix_solver_sm_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
- * and thus belong into a different calculation. This applies to all solvers.
- */
-
- const unsigned iN = this->N();
- T cerr = 0;
- for (unsigned i = 0; i < iN; i++)
- 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>
-template <typename T>
-void matrix_solver_sm_t<m_N, _storage_N>::store(
- const T * RESTRICT V)
-{
- for (unsigned i = 0, iN=N(); i < iN; i++)
- {
- this->m_nets[i]->m_cur_Analog = V[i];
- }
-}
-
-
-template <unsigned m_N, unsigned _storage_N>
int matrix_solver_sm_t<m_N, _storage_N>::solve_non_dynamic(ATTR_UNUSED const bool newton_raphson)
{
static const bool incremental = true;
@@ -605,15 +488,12 @@ matrix_solver_sm_t<m_N, _storage_N>::matrix_solver_sm_t(const solver_parameters_
: matrix_solver_t(GAUSSIAN_ELIMINATION, params)
, m_dim(size)
{
- 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);
m_last_RHS[k] = 0.0;
- m_last_V[k] = 0.0;
}
}
@@ -622,15 +502,12 @@ matrix_solver_sm_t<m_N, _storage_N>::matrix_solver_sm_t(const eSolverType type,
: matrix_solver_t(type, params)
, m_dim(size)
{
- 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);
m_last_RHS[k] = 0.0;
- m_last_V[k] = 0.0;
}
}
diff --git a/src/lib/netlist/solver/nld_ms_w.h b/src/lib/netlist/solver/nld_ms_w.h
index 154f618a697..7a4ee61dd8c 100644
--- a/src/lib/netlist/solver/nld_ms_w.h
+++ b/src/lib/netlist/solver/nld_ms_w.h
@@ -66,7 +66,6 @@ public:
virtual void reset() override { matrix_solver_t::reset(); }
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);
@@ -79,13 +78,6 @@ protected:
template <typename T>
void LE_compute_x(T * RESTRICT x);
- template <typename T>
- T delta(const T * RESTRICT V);
-
- template <typename T>
- void store(const T * RESTRICT V);
-
- virtual netlist_time compute_next_timestep() override;
template <typename T1, typename T2>
inline nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r][c]; }
@@ -103,10 +95,6 @@ protected:
inline nl_ext_double &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; }
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[_storage_N];
- terms_t *m_rails_temp;
private:
static const std::size_t m_pitch = ((( _storage_N) + 7) / 8) * 8;
@@ -138,102 +126,26 @@ private:
template <unsigned m_N, unsigned _storage_N>
matrix_solver_w_t<m_N, _storage_N>::~matrix_solver_w_t()
{
- for (unsigned k = 0; k < N(); k++)
- {
- pfree(m_terms[k]);
- }
- pfree_array(m_rails_temp);
#if (NL_USE_DYNAMIC_ALLOCATION)
pfree_array(m_A);
#endif
}
template <unsigned m_N, unsigned _storage_N>
-netlist_time matrix_solver_w_t<m_N, _storage_N>::compute_next_timestep()
-{
- nl_double new_solver_timestep = m_params.m_max_timestep;
-
- if (m_params.m_dynamic)
- {
- /*
- * FIXME: We should extend the logic to use either all nets or
- * only output nets.
- */
- for (unsigned k = 0, iN=N(); k < iN; k++)
- {
- analog_net_t *n = m_nets[k];
-
- const nl_double DD_n = (n->Q_Analog() - m_last_V[k]);
- const nl_double hn = current_timestep();
-
- nl_double DD2 = (DD_n / hn - n->m_DD_n_m_1 / n->m_h_n_m_1) / (hn + n->m_h_n_m_1);
- nl_double new_net_timestep;
-
- n->m_h_n_m_1 = hn;
- n->m_DD_n_m_1 = DD_n;
- if (nl_math::abs(DD2) > NL_FCONST(1e-30)) // avoid div-by-zero
- new_net_timestep = nl_math::sqrt(m_params.m_lte / nl_math::abs(NL_FCONST(0.5)*DD2));
- else
- new_net_timestep = m_params.m_max_timestep;
-
- 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 netlist_time::from_double(new_solver_timestep);
-}
-
-template <unsigned m_N, unsigned _storage_N>
-ATTR_COLD void matrix_solver_w_t<m_N, _storage_N>::add_term(int k, terminal_t *term)
-{
- if (term->m_otherterm->net().isRailNet())
- {
- m_rails_temp[k].add(term, -1, false);
- }
- else
- {
- int ot = get_net_idx(&term->m_otherterm->net());
- if (ot>=0)
- {
- m_terms[k]->add(term, ot, true);
- }
- /* Should this be allowed ? */
- else // if (ot<0)
- {
- m_rails_temp[k].add(term, ot, true);
- log().fatal("found term with missing othernet {1}\n", term->name());
- }
- }
-}
-
-
-template <unsigned m_N, unsigned _storage_N>
ATTR_COLD void matrix_solver_w_t<m_N, _storage_N>::vsetup(analog_net_t::list_t &nets)
{
if (m_dim < nets.size())
log().fatal("Dimension {1} less than {2}", m_dim, nets.size());
- for (unsigned k = 0; k < N(); k++)
- {
- m_terms[k]->clear();
- m_rails_temp[k].clear();
- }
-
matrix_solver_t::setup_base(nets);
for (unsigned k = 0; k < N(); k++)
{
m_terms[k]->m_railstart = m_terms[k]->count();
- for (unsigned i = 0; i < m_rails_temp[k].count(); i++)
- this->m_terms[k]->add(m_rails_temp[k].terms()[i], m_rails_temp[k].net_other()[i], false);
+ for (unsigned i = 0; i < m_rails_temp[k]->count(); i++)
+ this->m_terms[k]->add(m_rails_temp[k]->terms()[i], m_rails_temp[k]->net_other()[i], false);
- m_rails_temp[k].clear(); // no longer needed
+ m_rails_temp[k]->clear(); // no longer needed
m_terms[k]->set_pointers();
}
@@ -329,13 +241,13 @@ ATTR_COLD void matrix_solver_w_t<m_N, _storage_N>::vsetup(analog_net_t::list_t &
* save states
*/
save(NLNAME(m_last_RHS));
- save(NLNAME(m_last_V));
for (unsigned k = 0; k < N(); k++)
{
pstring num = pfmt("{1}")(k);
save(RHS(k), "RHS" + num);
+ save(m_terms[k]->m_last_V, "lastV" + num);
save(m_terms[k]->go(),"GO" + num, m_terms[k]->count());
save(m_terms[k]->gt(),"GT" + num, m_terms[k]->count());
@@ -383,15 +295,15 @@ void matrix_solver_w_t<m_N, _storage_N>::build_LE_RHS()
nl_double rhsk_a = 0.0;
nl_double rhsk_b = 0.0;
- const int terms_count = m_terms[k]->count();
+ const unsigned terms_count = m_terms[k]->count();
const nl_double * RESTRICT go = m_terms[k]->go();
const nl_double * RESTRICT Idr = m_terms[k]->Idr();
const nl_double * const * RESTRICT other_cur_analog = m_terms[k]->other_curanalog();
- for (int i = 0; i < terms_count; i++)
+ for (unsigned i = 0; i < terms_count; i++)
rhsk_a = rhsk_a + Idr[i];
- for (int i = m_terms[k]->m_railstart; i < terms_count; i++)
+ for (unsigned i = m_terms[k]->m_railstart; i < terms_count; i++)
//rhsk = rhsk + go[i] * terms[i]->m_otherterm->net().as_analog().Q_Analog();
rhsk_b = rhsk_b + go[i] * *other_cur_analog[i];
@@ -482,34 +394,6 @@ void matrix_solver_w_t<m_N, _storage_N>::LE_compute_x(
template <unsigned m_N, unsigned _storage_N>
-template <typename T>
-T matrix_solver_w_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
- * and thus belong into a different calculation. This applies to all solvers.
- */
-
- const unsigned iN = this->N();
- T cerr = 0;
- for (unsigned i = 0; i < iN; i++)
- 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>
-template <typename T>
-void matrix_solver_w_t<m_N, _storage_N>::store(
- const T * RESTRICT V)
-{
- for (unsigned i = 0, iN=N(); i < iN; i++)
- {
- this->m_nets[i]->m_cur_Analog = V[i];
- }
-}
-
-template <unsigned m_N, unsigned _storage_N>
int matrix_solver_w_t<m_N, _storage_N>::solve_non_dynamic(ATTR_UNUSED const bool newton_raphson)
{
const auto iN = N();
@@ -668,15 +552,12 @@ matrix_solver_w_t<m_N, _storage_N>::matrix_solver_w_t(const solver_parameters_t
,m_cnt(0)
, m_dim(size)
{
- 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);
m_last_RHS[k] = 0.0;
- m_last_V[k] = 0.0;
}
}
@@ -686,12 +567,9 @@ matrix_solver_w_t<m_N, _storage_N>::matrix_solver_w_t(const eSolverType type, co
,m_cnt(0)
, m_dim(size)
{
- m_rails_temp = palloc_array(terms_t, N());
for (unsigned k = 0; k < N(); k++)
{
- m_terms[k] = palloc(terms_t);
m_last_RHS[k] = 0.0;
- m_last_V[k] = 0.0;
}
}
diff --git a/src/lib/netlist/solver/nld_solver.cpp b/src/lib/netlist/solver/nld_solver.cpp
index ed011ff254d..bcd7feb7c69 100644
--- a/src/lib/netlist/solver/nld_solver.cpp
+++ b/src/lib/netlist/solver/nld_solver.cpp
@@ -108,6 +108,11 @@ ATTR_COLD matrix_solver_t::matrix_solver_t(const eSolverType type, const solver_
ATTR_COLD matrix_solver_t::~matrix_solver_t()
{
m_inps.clear_and_free();
+ for (unsigned k = 0; k < m_terms.size(); k++)
+ {
+ pfree(m_terms[k]);
+ }
+
}
ATTR_COLD void matrix_solver_t::setup_base(analog_net_t::list_t &nets)
@@ -115,9 +120,14 @@ ATTR_COLD void matrix_solver_t::setup_base(analog_net_t::list_t &nets)
log().debug("New solver setup\n");
m_nets.clear();
+ m_terms.clear();
for (auto & net : nets)
+ {
m_nets.push_back(net);
+ m_terms.push_back(palloc(terms_t));
+ m_rails_temp.push_back(palloc(terms_t));
+ }
for (std::size_t k = 0; k < nets.size(); k++)
{
@@ -308,6 +318,70 @@ ATTR_COLD int matrix_solver_t::get_net_idx(net_t *net)
return -1;
}
+ATTR_COLD void matrix_solver_t::add_term(int k, terminal_t *term)
+{
+ if (term->m_otherterm->net().isRailNet())
+ {
+ m_rails_temp[k]->add(term, -1, false);
+ }
+ else
+ {
+ int ot = get_net_idx(&term->m_otherterm->net());
+ if (ot>=0)
+ {
+ m_terms[k]->add(term, ot, true);
+ }
+ /* Should this be allowed ? */
+ else // if (ot<0)
+ {
+ m_rails_temp[k]->add(term, ot, true);
+ log().fatal("found term with missing othernet {1}\n", term->name());
+ }
+ }
+}
+
+netlist_time matrix_solver_t::compute_next_timestep()
+{
+ nl_double new_solver_timestep = m_params.m_max_timestep;
+
+ if (m_params.m_dynamic)
+ {
+ /*
+ * FIXME: We should extend the logic to use either all nets or
+ * only output nets.
+ */
+ for (unsigned k = 0, iN=m_terms.size(); k < iN; k++)
+ {
+ analog_net_t *n = m_nets[k];
+
+ const nl_double DD_n = (n->Q_Analog() - m_terms[k]->m_last_V);
+ const nl_double hn = current_timestep();
+
+ nl_double DD2 = (DD_n / hn - n->m_DD_n_m_1 / n->m_h_n_m_1) / (hn + n->m_h_n_m_1);
+ nl_double new_net_timestep;
+
+ n->m_h_n_m_1 = hn;
+ n->m_DD_n_m_1 = DD_n;
+ if (nl_math::abs(DD2) > NL_FCONST(1e-30)) // avoid div-by-zero
+ new_net_timestep = nl_math::sqrt(m_params.m_lte / nl_math::abs(NL_FCONST(0.5)*DD2));
+ else
+ new_net_timestep = m_params.m_max_timestep;
+
+ if (new_net_timestep < new_solver_timestep)
+ new_solver_timestep = new_net_timestep;
+
+ m_terms[k]->m_last_V = 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 netlist_time::from_double(new_solver_timestep);
+}
+
+
+
void matrix_solver_t::log_stats()
{
if (this->m_stat_calculations != 0 && this->m_stat_vsolver_calls && this->m_params.m_log_stats)