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, 193 insertions, 278 deletions
diff --git a/src/lib/netlist/solver/nld_matrix_solver.h b/src/lib/netlist/solver/nld_matrix_solver.h
index f76660e1cb9..da44370d0d9 100644
--- a/src/lib/netlist/solver/nld_matrix_solver.h
+++ b/src/lib/netlist/solver/nld_matrix_solver.h
@@ -10,367 +10,282 @@
#include "netlist/nl_base.h"
#include "netlist/nl_errstr.h"
-#include "plib/palloc.h"
-#include "plib/pmatrix2d.h"
-#include "plib/putil.h"
-#include "plib/vector_ops.h"
-
-#include <cmath>
+#include "netlist/plib/putil.h"
namespace netlist
{
-namespace devices
-{
+ namespace devices
+ {
/* FIXME: these should become proper devices */
struct solver_parameters_t
{
- bool m_pivot;
+ int 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;
- std::size_t m_gs_loops;
- std::size_t m_nr_loops;
+ unsigned m_gs_loops;
+ unsigned 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();
-
- void clear();
+class terms_for_net_t : plib::nocopyassignmove
+{
+public:
+ terms_for_net_t();
- void add(terminal_t *term, int net_other, bool sorted);
+ void clear();
- std::size_t count() const { return m_terms.size(); }
+ void add(terminal_t *term, int net_other, bool sorted);
- terminal_t **terms() { return m_terms.data(); }
+ inline std::size_t count() const { return m_terms.size(); }
- std::size_t m_railstart;
+ 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(); }
- 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 */
+ void set_pointers();
- /* state */
- nl_double m_last_V;
- nl_double m_DD_n_m_1;
- nl_double m_h_n_m_1;
+ std::size_t m_railstart;
- std::vector<int> m_connected_net_idx;
- private:
- std::vector<terminal_t *> m_terms;
+ 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;
- class proxied_analog_output_t : public analog_output_t
- {
- public:
+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;
- proxied_analog_output_t(core_device_t &dev, const pstring &aname, analog_net_t *pnet)
- : analog_output_t(dev, aname)
- , m_proxied_net(pnet)
- { }
+};
- 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
- };
+class proxied_analog_output_t : public analog_output_t
+{
+public:
- class matrix_solver_t : public device_t
- {
- public:
- using list_t = std::vector<matrix_solver_t *>;
+ 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();
- enum eSortType
- {
- NOSORT,
- ASCENDING,
- DESCENDING,
- PREFER_IDENTITY_TOP_LEFT,
- PREFER_BAND_MATRIX
- };
-
- void setup(analog_net_t::list_t &nets)
- {
- vsetup(nets);
- }
+ analog_net_t *m_proxied_net; // only for proxy nets in analog input logic
+};
- void solve_base();
- /* 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();
+class matrix_solver_t : public device_t
+{
+public:
+ using list_t = std::vector<matrix_solver_t *>;
- bool has_dynamic_devices() const { return m_dynamic_devices.size() > 0; }
- bool has_timestep_devices() const { return m_step_devices.size() > 0; }
+ enum eSortType
+ {
+ NOSORT,
+ ASCENDING,
+ DESCENDING
+ };
- void update_forced();
- void update_after(const netlist_time after)
- {
- m_Q_sync.net().toggle_and_push_to_queue(after);
- }
+ virtual ~matrix_solver_t() override;
- /* netdevice functions */
- NETLIB_UPDATEI();
- NETLIB_RESETI();
+ void setup(analog_net_t::list_t &nets)
+ {
+ vsetup(nets);
+ }
- 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);
+ void solve_base();
- virtual void log_stats();
+ /* 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();
- virtual std::pair<pstring, pstring> create_solver_code()
- {
- return std::pair<pstring, pstring>("", plib::pfmt("/* solver doesn't support static compile */\n\n"));
- }
+ inline bool has_dynamic_devices() const { return m_dynamic_devices.size() > 0; }
+ inline bool has_timestep_devices() const { return m_step_devices.size() > 0; }
- /* return number of floating point operations for solve */
- std::size_t ops() { return m_ops; }
+ void update_forced();
+ void update_after(const netlist_time &after)
+ {
+ m_Q_sync.net().toggle_and_push_to_queue(after);
+ }
- protected:
+ /* netdevice functions */
+ NETLIB_UPDATEI();
+ NETLIB_RESETI();
- matrix_solver_t(netlist_state_t &anetlist, const pstring &name,
- eSortType sort, const solver_parameters_t *params);
+public:
+ int get_net_idx(detail::net_t *net);
- void sort_terms(eSortType sort);
+ virtual void log_stats();
- void setup_base(analog_net_t::list_t &nets);
- void update_dynamic();
+ virtual std::pair<pstring, pstring> create_solver_code()
+ {
+ return std::pair<pstring, pstring>("", plib::pfmt("/* solver doesn't support static compile */\n\n"));
+ }
- virtual void vsetup(analog_net_t::list_t &nets) = 0;
- virtual unsigned vsolve_non_dynamic(const bool newton_raphson) = 0;
+ /* return number of floating point operations for solve */
+ std::size_t ops() { return m_ops; }
- netlist_time compute_next_timestep(const double cur_ts);
- /* virtual */ void add_term(std::size_t net_idx, terminal_t *term);
+protected:
- template <typename T>
- void store(const T & V);
+ matrix_solver_t(netlist_t &anetlist, const pstring &name,
+ const eSortType sort, const solver_parameters_t *params);
- template <typename T>
- auto delta(const T & V) -> typename std::decay<decltype(V[0])>::type;
+ void setup_base(analog_net_t::list_t &nets);
+ void update_dynamic();
- template <typename T>
- void build_LE_A(T &child);
- template <typename T>
- void build_LE_RHS(T &child);
+ virtual void vsetup(analog_net_t::list_t &nets) = 0;
+ virtual unsigned vsolve_non_dynamic(const bool newton_raphson) = 0;
- 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();
- }
- }
- }
+ netlist_time compute_next_timestep(const double cur_ts);
+ /* virtual */ void add_term(std::size_t net_idx, terminal_t *term);
- 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;
- }
+ template <typename T>
+ void store(const T * RESTRICT V);
+ template <typename T>
+ T delta(const T * RESTRICT V);
- }
+ template <typename T>
+ void build_LE_A();
+ template <typename T>
+ void build_LE_RHS();
- template <typename T>
- using aligned_alloc = plib::aligned_allocator<T, PALIGN_VECTOROPT>;
+ 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;
- 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;
+ std::vector<terms_for_net_t *> m_rails_temp;
- plib::pmatrix2d<nl_double> m_test;
+ const solver_parameters_t &m_params;
- 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<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;
- std::vector<plib::unique_ptr<terms_for_net_t>> m_rails_temp;
+private:
- const solver_parameters_t &m_params;
+ state_var<netlist_time> m_last_step;
+ std::vector<core_device_t *> m_step_devices;
+ std::vector<core_device_t *> m_dynamic_devices;
- 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;
+ logic_input_t m_fb_sync;
+ logic_output_t m_Q_sync;
- private:
+ /* calculate matrix */
+ void setup_matrix();
- state_var<netlist_time> m_last_step;
- std::vector<core_device_t *> m_step_devices;
- std::vector<core_device_t *> m_dynamic_devices;
+ void step(const netlist_time &delta);
- logic_input_t m_fb_sync;
- logic_output_t m_Q_sync;
+ std::size_t m_ops;
+ const eSortType m_sort;
+};
- /* calculate matrix */
- void setup_matrix();
+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]);
+}
- void step(const netlist_time &delta);
+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");
- std::size_t m_ops;
- const eSortType m_sort;
- };
+ T &child = static_cast<T &>(*this);
- template <typename T>
- auto matrix_solver_t::delta(const T & V) -> typename std::decay<decltype(V[0])>::type
+ const std::size_t iN = child.N();
+ for (std::size_t k = 0; k < iN; k++)
{
- /* 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;
- }
+ terms_for_net_t *terms = m_terms[k].get();
+ nl_double * Ak = &child.A(k, 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]);
- }
+ for (std::size_t i=0; i < iN; i++)
+ Ak[i] = 0.0;
- 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 terms_count = terms->count();
+ const std::size_t railstart = terms->m_railstart;
+ const nl_double * const RESTRICT gt = terms->gt();
- const std::size_t iN = child.size();
- for (std::size_t k = 0; k < iN; k++)
{
- terms_for_net_t *terms = m_terms[k].get();
- float_type * Ak = &child.A(k, 0ul);
-
- for (std::size_t i=0; i < iN; i++)
- Ak[i] = 0.0;
-
- const std::size_t terms_count = terms->count();
- const std::size_t railstart = terms->m_railstart;
- const float_type * const gt = m_gtn[k];
-
- {
- float_type akk = 0.0;
- for (std::size_t i = 0; i < terms_count; i++)
- akk += gt[i];
+ nl_double akk = 0.0;
+ for (std::size_t i = 0; i < terms_count; i++)
+ akk += gt[i];
- Ak[k] = akk;
- }
+ Ak[k] = akk;
+ }
- const float_type * const go = m_gonn[k];
- int * net_other = terms->m_connected_net_idx.data();
+ const nl_double * const RESTRICT go = terms->go();
+ int * RESTRICT net_other = terms->connected_net_idx();
- for (std::size_t i = 0; i < railstart; i++)
- Ak[net_other[i]] += go[i];
- }
+ for (std::size_t i = 0; i < railstart; i++)
+ Ak[net_other[i]] -= go[i];
}
+}
- template <typename T>
- void matrix_solver_t::build_LE_RHS(T &child)
- {
- 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;
+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);
- 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;
+ const std::size_t iN = child.N();
+ for (std::size_t k = 0; k < iN; k++)
+ {
+ nl_double rhsk_a = 0.0;
+ nl_double rhsk_b = 0.0;
- 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];
+ 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();
- for (std::size_t i = 0; i < terms_count; i++)
- rhsk_a = rhsk_a + Idr[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];
+ 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_ */