summaryrefslogtreecommitdiffstatshomepage
path: root/src/lib/netlist/solver/nld_matrix_solver.h
diff options
context:
space:
mode:
Diffstat (limited to 'src/lib/netlist/solver/nld_matrix_solver.h')
-rw-r--r--src/lib/netlist/solver/nld_matrix_solver.h471
1 files changed, 278 insertions, 193 deletions
diff --git a/src/lib/netlist/solver/nld_matrix_solver.h b/src/lib/netlist/solver/nld_matrix_solver.h
index da44370d0d9..f76660e1cb9 100644
--- a/src/lib/netlist/solver/nld_matrix_solver.h
+++ b/src/lib/netlist/solver/nld_matrix_solver.h
@@ -10,282 +10,367 @@
#include "netlist/nl_base.h"
#include "netlist/nl_errstr.h"
-#include "netlist/plib/putil.h"
+#include "plib/palloc.h"
+#include "plib/pmatrix2d.h"
+#include "plib/putil.h"
+#include "plib/vector_ops.h"
+
+#include <cmath>
namespace netlist
{
- namespace devices
- {
+namespace devices
+{
/* FIXME: these should become proper devices */
struct solver_parameters_t
{
- int m_pivot;
+ bool m_pivot;
nl_double m_accuracy;
nl_double m_dynamic_lte;
nl_double m_min_timestep;
nl_double m_max_timestep;
nl_double m_gs_sor;
bool m_dynamic_ts;
- unsigned m_gs_loops;
- unsigned m_nr_loops;
+ std::size_t m_gs_loops;
+ std::size_t m_nr_loops;
netlist_time m_nr_recalc_delay;
+ bool m_use_gabs;
+ bool m_use_linear_prediction;
bool m_log_stats;
};
-class terms_for_net_t : plib::nocopyassignmove
-{
-public:
- terms_for_net_t();
+ class terms_for_net_t : plib::nocopyassignmove
+ {
+ public:
+ terms_for_net_t();
- void clear();
+ void clear();
- void add(terminal_t *term, int net_other, bool sorted);
+ void add(terminal_t *term, int net_other, bool sorted);
- inline std::size_t count() const { return m_terms.size(); }
+ std::size_t count() const { return m_terms.size(); }
- inline terminal_t **terms() { return m_terms.data(); }
- inline int *connected_net_idx() { return m_connected_net_idx.data(); }
- inline nl_double *gt() { return m_gt.data(); }
- inline nl_double *go() { return m_go.data(); }
- inline nl_double *Idr() { return m_Idr.data(); }
- inline nl_double * const *connected_net_V() const { return m_connected_net_V.data(); }
+ terminal_t **terms() { return m_terms.data(); }
- void set_pointers();
+ std::size_t m_railstart;
- std::size_t m_railstart;
+ std::vector<unsigned> m_nz; /* all non zero for multiplication */
+ std::vector<unsigned> m_nzrd; /* non zero right of the diagonal for elimination, may include RHS element */
+ std::vector<unsigned> m_nzbd; /* non zero below of the diagonal for elimination */
- std::vector<unsigned> m_nz; /* all non zero for multiplication */
- std::vector<unsigned> m_nzrd; /* non zero right of the diagonal for elimination, may include RHS element */
- std::vector<unsigned> m_nzbd; /* non zero below of the diagonal for elimination */
+ /* state */
+ nl_double m_last_V;
+ nl_double m_DD_n_m_1;
+ nl_double m_h_n_m_1;
- /* state */
- nl_double m_last_V;
- nl_double m_DD_n_m_1;
- nl_double m_h_n_m_1;
+ std::vector<int> m_connected_net_idx;
+ private:
+ std::vector<terminal_t *> m_terms;
-private:
- std::vector<int> m_connected_net_idx;
- std::vector<nl_double> m_go;
- std::vector<nl_double> m_gt;
- std::vector<nl_double> m_Idr;
- std::vector<nl_double *> m_connected_net_V;
- std::vector<terminal_t *> m_terms;
+ };
-};
+ class proxied_analog_output_t : public analog_output_t
+ {
+ public:
-class proxied_analog_output_t : public analog_output_t
-{
-public:
+ proxied_analog_output_t(core_device_t &dev, const pstring &aname, analog_net_t *pnet)
+ : analog_output_t(dev, aname)
+ , m_proxied_net(pnet)
+ { }
- proxied_analog_output_t(core_device_t &dev, const pstring &aname)
- : analog_output_t(dev, aname)
- , m_proxied_net(nullptr)
- { }
- virtual ~proxied_analog_output_t();
+ analog_net_t *proxied_net() const { return m_proxied_net;}
+ private:
+ analog_net_t *m_proxied_net; // only for proxy nets in analog input logic
+ };
- analog_net_t *m_proxied_net; // only for proxy nets in analog input logic
-};
+ class matrix_solver_t : public device_t
+ {
+ public:
+ using list_t = std::vector<matrix_solver_t *>;
+ enum eSortType
+ {
+ NOSORT,
+ ASCENDING,
+ DESCENDING,
+ PREFER_IDENTITY_TOP_LEFT,
+ PREFER_BAND_MATRIX
+ };
+
+ void setup(analog_net_t::list_t &nets)
+ {
+ vsetup(nets);
+ }
-class matrix_solver_t : public device_t
-{
-public:
- using list_t = std::vector<matrix_solver_t *>;
+ void solve_base();
- enum eSortType
- {
- NOSORT,
- ASCENDING,
- DESCENDING
- };
+ /* after every call to solve, update inputs must be called.
+ * this can be done as well as a batch to ease parallel processing.
+ */
+ const netlist_time solve(netlist_time now);
+ void update_inputs();
- virtual ~matrix_solver_t() override;
+ bool has_dynamic_devices() const { return m_dynamic_devices.size() > 0; }
+ bool has_timestep_devices() const { return m_step_devices.size() > 0; }
- void setup(analog_net_t::list_t &nets)
- {
- vsetup(nets);
- }
+ void update_forced();
+ void update_after(const netlist_time after)
+ {
+ m_Q_sync.net().toggle_and_push_to_queue(after);
+ }
- void solve_base();
+ /* netdevice functions */
+ NETLIB_UPDATEI();
+ NETLIB_RESETI();
- /* after every call to solve, update inputs must be called.
- * this can be done as well as a batch to ease parallel processing.
- */
- const netlist_time solve();
- void update_inputs();
+ public:
+ int get_net_idx(detail::net_t *net);
+ std::pair<int, int> get_left_right_of_diag(std::size_t row, std::size_t diag);
+ double get_weight_around_diag(std::size_t row, std::size_t diag);
- inline bool has_dynamic_devices() const { return m_dynamic_devices.size() > 0; }
- inline bool has_timestep_devices() const { return m_step_devices.size() > 0; }
+ virtual void log_stats();
- void update_forced();
- void update_after(const netlist_time &after)
- {
- m_Q_sync.net().toggle_and_push_to_queue(after);
- }
+ virtual std::pair<pstring, pstring> create_solver_code()
+ {
+ return std::pair<pstring, pstring>("", plib::pfmt("/* solver doesn't support static compile */\n\n"));
+ }
- /* netdevice functions */
- NETLIB_UPDATEI();
- NETLIB_RESETI();
+ /* return number of floating point operations for solve */
+ std::size_t ops() { return m_ops; }
-public:
- int get_net_idx(detail::net_t *net);
+ protected:
- virtual void log_stats();
+ matrix_solver_t(netlist_state_t &anetlist, const pstring &name,
+ eSortType sort, const solver_parameters_t *params);
- virtual std::pair<pstring, pstring> create_solver_code()
- {
- return std::pair<pstring, pstring>("", plib::pfmt("/* solver doesn't support static compile */\n\n"));
- }
+ void sort_terms(eSortType sort);
- /* return number of floating point operations for solve */
- std::size_t ops() { return m_ops; }
+ void setup_base(analog_net_t::list_t &nets);
+ void update_dynamic();
-protected:
+ virtual void vsetup(analog_net_t::list_t &nets) = 0;
+ virtual unsigned vsolve_non_dynamic(const bool newton_raphson) = 0;
- matrix_solver_t(netlist_t &anetlist, const pstring &name,
- const eSortType sort, const solver_parameters_t *params);
+ netlist_time compute_next_timestep(const double cur_ts);
+ /* virtual */ void add_term(std::size_t net_idx, terminal_t *term);
- void setup_base(analog_net_t::list_t &nets);
- void update_dynamic();
+ template <typename T>
+ void store(const T & V);
- virtual void vsetup(analog_net_t::list_t &nets) = 0;
- virtual unsigned vsolve_non_dynamic(const bool newton_raphson) = 0;
+ template <typename T>
+ auto delta(const T & V) -> typename std::decay<decltype(V[0])>::type;
- netlist_time compute_next_timestep(const double cur_ts);
- /* virtual */ void add_term(std::size_t net_idx, terminal_t *term);
+ template <typename T>
+ void build_LE_A(T &child);
+ template <typename T>
+ void build_LE_RHS(T &child);
- template <typename T>
- void store(const T * RESTRICT V);
- template <typename T>
- T delta(const T * RESTRICT V);
+ void set_pointers()
+ {
+ const std::size_t iN = this->m_nets.size();
+
+ std::size_t max_count = 0;
+ std::size_t max_rail = 0;
+ for (std::size_t k = 0; k < iN; k++)
+ {
+ max_count = std::max(max_count, m_terms[k]->count());
+ max_rail = std::max(max_rail, m_terms[k]->m_railstart);
+ }
+
+ m_mat_ptr.resize(iN, max_rail+1);
+ m_gtn.resize(iN, max_count);
+ m_gonn.resize(iN, max_count);
+ m_Idrn.resize(iN, max_count);
+ m_connected_net_Vn.resize(iN, max_count);
+
+ for (std::size_t k = 0; k < iN; k++)
+ {
+ auto count = m_terms[k]->count();
+
+ for (std::size_t i = 0; i < count; i++)
+ {
+ m_terms[k]->terms()[i]->set_ptrs(&m_gtn[k][i], &m_gonn[k][i], &m_Idrn[k][i]);
+ m_connected_net_Vn[k][i] = m_terms[k]->terms()[i]->connected_terminal()->net().Q_Analog_state_ptr();
+ }
+ }
+ }
- template <typename T>
- void build_LE_A();
- template <typename T>
- void build_LE_RHS();
+ template <typename AP, typename FT>
+ void fill_matrix(std::size_t N, AP &tcr, FT &RHS)
+ {
+ for (std::size_t k = 0; k < N; k++)
+ {
+ auto *net = m_terms[k].get();
+ auto **tcr_r = &(tcr[k][0]);
+
+ const std::size_t term_count = net->count();
+ const std::size_t railstart = net->m_railstart;
+ const auto &go = m_gonn[k];
+ const auto &gt = m_gtn[k];
+ const auto &Idr = m_Idrn[k];
+ const auto &cnV = m_connected_net_Vn[k];
+
+ for (std::size_t i = 0; i < railstart; i++)
+ *tcr_r[i] += go[i];
+
+ typename FT::value_type gtot_t = 0.0;
+ typename FT::value_type RHS_t = 0.0;
+
+ for (std::size_t i = 0; i < term_count; i++)
+ {
+ gtot_t += gt[i];
+ RHS_t += Idr[i];
+ }
+ // FIXME: Code above is faster than vec_sum - Check this
+ #if 0
+ auto gtot_t = plib::vec_sum<FT>(term_count, m_gt);
+ auto RHS_t = plib::vec_sum<FT>(term_count, m_Idr);
+ #endif
+
+ for (std::size_t i = railstart; i < term_count; i++)
+ {
+ RHS_t += (/*m_Idr[i]*/ (- go[i]) * *cnV[i]);
+ }
+
+ RHS[k] = RHS_t;
+ // update diagonal element ...
+ *tcr_r[railstart] += gtot_t; //mat.A[mat.diag[k]] += gtot_t;
+ }
- std::vector<std::unique_ptr<terms_for_net_t>> m_terms;
- std::vector<analog_net_t *> m_nets;
- std::vector<std::unique_ptr<proxied_analog_output_t>> m_inps;
+ }
- std::vector<terms_for_net_t *> m_rails_temp;
+ template <typename T>
+ using aligned_alloc = plib::aligned_allocator<T, PALIGN_VECTOROPT>;
- const solver_parameters_t &m_params;
+ plib::pmatrix2d<nl_double, aligned_alloc<nl_double>> m_gonn;
+ plib::pmatrix2d<nl_double, aligned_alloc<nl_double>> m_gtn;
+ plib::pmatrix2d<nl_double, aligned_alloc<nl_double>> m_Idrn;
+ plib::pmatrix2d<nl_double *, aligned_alloc<nl_double *>> m_mat_ptr;
+ plib::pmatrix2d<nl_double *, aligned_alloc<nl_double *>> m_connected_net_Vn;
- state_var<int> m_stat_calculations;
- state_var<int> m_stat_newton_raphson;
- state_var<int> m_stat_vsolver_calls;
- state_var<int> m_iterative_fail;
- state_var<int> m_iterative_total;
+ plib::pmatrix2d<nl_double> m_test;
-private:
+ std::vector<plib::unique_ptr<terms_for_net_t>> m_terms;
+ std::vector<analog_net_t *> m_nets;
+ std::vector<pool_owned_ptr<proxied_analog_output_t>> m_inps;
- state_var<netlist_time> m_last_step;
- std::vector<core_device_t *> m_step_devices;
- std::vector<core_device_t *> m_dynamic_devices;
+ std::vector<plib::unique_ptr<terms_for_net_t>> m_rails_temp;
- logic_input_t m_fb_sync;
- logic_output_t m_Q_sync;
+ const solver_parameters_t &m_params;
- /* calculate matrix */
- void setup_matrix();
+ state_var<int> m_stat_calculations;
+ state_var<int> m_stat_newton_raphson;
+ state_var<int> m_stat_vsolver_calls;
+ state_var<int> m_iterative_fail;
+ state_var<int> m_iterative_total;
- void step(const netlist_time &delta);
+ private:
- std::size_t m_ops;
- const eSortType m_sort;
-};
+ state_var<netlist_time> m_last_step;
+ std::vector<core_device_t *> m_step_devices;
+ std::vector<core_device_t *> m_dynamic_devices;
-template <typename T>
-T matrix_solver_t::delta(const T * RESTRICT V)
-{
- /* NOTE: Ideally we should also include currents (RHS) here. This would
- * need a reevaluation of the right hand side after voltages have been updated
- * and thus belong into a different calculation. This applies to all solvers.
- */
-
- const std::size_t iN = this->m_terms.size();
- T cerr = 0;
- for (std::size_t i = 0; i < iN; i++)
- cerr = std::max(cerr, std::abs(V[i] - static_cast<T>(this->m_nets[i]->Q_Analog())));
- return cerr;
-}
-
-template <typename T>
-void matrix_solver_t::store(const T * RESTRICT V)
-{
- const std::size_t iN = this->m_terms.size();
- for (std::size_t i = 0; i < iN; i++)
- this->m_nets[i]->set_Q_Analog(V[i]);
-}
+ logic_input_t m_fb_sync;
+ logic_output_t m_Q_sync;
-template <typename T>
-void matrix_solver_t::build_LE_A()
-{
- static_assert(std::is_base_of<matrix_solver_t, T>::value, "T must derive from matrix_solver_t");
+ /* calculate matrix */
+ void setup_matrix();
- T &child = static_cast<T &>(*this);
+ void step(const netlist_time &delta);
- const std::size_t iN = child.N();
- for (std::size_t k = 0; k < iN; k++)
+ std::size_t m_ops;
+ const eSortType m_sort;
+ };
+
+ template <typename T>
+ auto matrix_solver_t::delta(const T & V) -> typename std::decay<decltype(V[0])>::type
{
- terms_for_net_t *terms = m_terms[k].get();
- nl_double * Ak = &child.A(k, 0);
+ /* NOTE: Ideally we should also include currents (RHS) here. This would
+ * need a reevaluation of the right hand side after voltages have been updated
+ * and thus belong into a different calculation. This applies to all solvers.
+ */
+
+ const std::size_t iN = this->m_terms.size();
+ typename std::decay<decltype(V[0])>::type cerr = 0;
+ for (std::size_t i = 0; i < iN; i++)
+ cerr = std::max(cerr, std::abs(V[i] - this->m_nets[i]->Q_Analog()));
+ return cerr;
+ }
- for (std::size_t i=0; i < iN; i++)
- Ak[i] = 0.0;
+ template <typename T>
+ void matrix_solver_t::store(const T & V)
+ {
+ const std::size_t iN = this->m_terms.size();
+ for (std::size_t i = 0; i < iN; i++)
+ this->m_nets[i]->set_Q_Analog(V[i]);
+ }
- const std::size_t terms_count = terms->count();
- const std::size_t railstart = terms->m_railstart;
- const nl_double * const RESTRICT gt = terms->gt();
+ template <typename T>
+ void matrix_solver_t::build_LE_A(T &child)
+ {
+ using float_type = typename T::float_type;
+ static_assert(std::is_base_of<matrix_solver_t, T>::value, "T must derive from matrix_solver_t");
+ const std::size_t iN = child.size();
+ for (std::size_t k = 0; k < iN; k++)
{
- nl_double akk = 0.0;
- for (std::size_t i = 0; i < terms_count; i++)
- akk += gt[i];
+ terms_for_net_t *terms = m_terms[k].get();
+ float_type * Ak = &child.A(k, 0ul);
- Ak[k] = akk;
- }
+ for (std::size_t i=0; i < iN; i++)
+ Ak[i] = 0.0;
- const nl_double * const RESTRICT go = terms->go();
- int * RESTRICT net_other = terms->connected_net_idx();
+ const std::size_t terms_count = terms->count();
+ const std::size_t railstart = terms->m_railstart;
+ const float_type * const gt = m_gtn[k];
- for (std::size_t i = 0; i < railstart; i++)
- Ak[net_other[i]] -= go[i];
- }
-}
+ {
+ float_type akk = 0.0;
+ for (std::size_t i = 0; i < terms_count; i++)
+ akk += gt[i];
-template <typename T>
-void matrix_solver_t::build_LE_RHS()
-{
- static_assert(std::is_base_of<matrix_solver_t, T>::value, "T must derive from matrix_solver_t");
- T &child = static_cast<T &>(*this);
+ Ak[k] = akk;
+ }
+
+ const float_type * const go = m_gonn[k];
+ int * net_other = terms->m_connected_net_idx.data();
+
+ for (std::size_t i = 0; i < railstart; i++)
+ Ak[net_other[i]] += go[i];
+ }
+ }
- const std::size_t iN = child.N();
- for (std::size_t k = 0; k < iN; k++)
+ template <typename T>
+ void matrix_solver_t::build_LE_RHS(T &child)
{
- nl_double rhsk_a = 0.0;
- nl_double rhsk_b = 0.0;
+ static_assert(std::is_base_of<matrix_solver_t, T>::value, "T must derive from matrix_solver_t");
+ using float_type = typename T::float_type;
- const std::size_t terms_count = m_terms[k]->count();
- const nl_double * const RESTRICT go = m_terms[k]->go();
- const nl_double * const RESTRICT Idr = m_terms[k]->Idr();
- const nl_double * const * RESTRICT other_cur_analog = m_terms[k]->connected_net_V();
+ const std::size_t iN = child.size();
+ for (std::size_t k = 0; k < iN; k++)
+ {
+ float_type rhsk_a = 0.0;
+ float_type rhsk_b = 0.0;
- for (std::size_t i = 0; i < terms_count; i++)
- rhsk_a = rhsk_a + Idr[i];
+ const std::size_t terms_count = m_terms[k]->count();
+ const float_type * const go = m_gonn[k];
+ const float_type * const Idr = m_Idrn[k];
+ const float_type * const * other_cur_analog = m_connected_net_Vn[k];
- for (std::size_t 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];
+ for (std::size_t i = 0; i < terms_count; i++)
+ rhsk_a = rhsk_a + Idr[i];
+
+ for (std::size_t 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];
- child.RHS(k) = rhsk_a + rhsk_b;
+ child.RHS(k) = rhsk_a + rhsk_b;
+ }
}
-}
- } //namespace devices
+} //namespace devices
} // namespace netlist
#endif /* NLD_MS_DIRECT_H_ */