diff options
| author | 2016-04-09 20:06:25 +0200 | |
|---|---|---|
| committer | 2016-04-09 23:45:42 +0200 | |
| commit | 0d1e57cc40fb31dcc5f29588b1695d6ce1ab9c8a (patch) | |
| tree | 7b98517a3f134477084b5e063a1c41b8fd20dee1 /src/lib/netlist/solver | |
| parent | db1e7f6b076cd092fccd70f3311e40499725281c (diff) | |
Solver model simplification.
Diffstat (limited to 'src/lib/netlist/solver')
| -rw-r--r-- | src/lib/netlist/solver/nld_matrix_solver.h | 44 | ||||
| -rw-r--r-- | src/lib/netlist/solver/nld_ms_direct.h | 137 | ||||
| -rw-r--r-- | src/lib/netlist/solver/nld_ms_direct_lu.h | 11 | ||||
| -rw-r--r-- | src/lib/netlist/solver/nld_ms_sm.h | 131 | ||||
| -rw-r--r-- | src/lib/netlist/solver/nld_ms_w.h | 136 | ||||
| -rw-r--r-- | src/lib/netlist/solver/nld_solver.cpp | 74 |
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) |
