// license:GPL-2.0+ // copyright-holders:Couriersud /* * nld_matrix_solver.h * */ #ifndef NLD_MATRIX_SOLVER_H_ #define NLD_MATRIX_SOLVER_H_ #include //#include "solver/nld_solver.h" #include "nl_base.h" #include "plib/pstream.h" namespace netlist { namespace devices { /* FIXME: these should become proper devices */ struct solver_parameters_t { int m_pivot; nl_double m_accuracy; nl_double m_lte; nl_double m_min_timestep; nl_double m_max_timestep; nl_double m_sor; bool m_dynamic; int m_gs_loops; int m_nr_loops; netlist_time m_nt_sync_delay; bool m_log_stats; }; class terms_t { P_PREVENT_COPYING(terms_t) public: terms_t() : m_railstart(0) , m_last_V(0.0) , m_DD_n_m_1(0.0) , m_h_n_m_1(1e-6) {} void clear() { m_term.clear(); m_net_other.clear(); m_gt.clear(); m_go.clear(); m_Idr.clear(); m_other_curanalog.clear(); } void add(terminal_t *term, int net_other, bool sorted); inline unsigned count() { return m_term.size(); } inline terminal_t **terms() { return m_term.data(); } inline int *net_other() { return m_net_other.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 **other_curanalog() { return m_other_curanalog.data(); } void set_pointers(); unsigned m_railstart; std::vector m_nz; /* all non zero for multiplication */ std::vector m_nzrd; /* non zero right of the diagonal for elimination, may include RHS element */ std::vector 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; private: std::vector m_net_other; std::vector m_go; std::vector m_gt; std::vector m_Idr; std::vector m_other_curanalog; std::vector m_term; }; class proxied_analog_output_t : public analog_output_t { public: proxied_analog_output_t(core_device_t &dev, const pstring &aname) : analog_output_t(dev, aname) , m_proxied_net(nullptr) { } 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; enum eSortType { NOSORT, ASCENDING, DESCENDING }; matrix_solver_t(netlist_t &anetlist, const pstring &name, const eSortType sort, const solver_parameters_t *params) : device_t(anetlist, name) , m_params(*params) , m_stat_calculations(*this, "m_stat_calculations", 0) , m_stat_newton_raphson(*this, "m_stat_newton_raphson", 0) , m_stat_vsolver_calls(*this, "m_stat_vsolver_calls", 0) , m_iterative_fail(*this, "m_iterative_fail", 0) , m_iterative_total(*this, "m_iterative_total", 0) , m_last_step(*this, "m_last_step", netlist_time::quantum()) , m_cur_ts(*this, "m_cur_ts", 0) , m_fb_sync(*this, "FB_sync") , m_Q_sync(*this, "Q_sync") , m_sort(sort) { connect_post_start(m_fb_sync, m_Q_sync); } virtual ~matrix_solver_t(); void setup(analog_net_t::list_t &nets) { vsetup(nets); } const netlist_time solve_base(); const netlist_time solve(); inline bool has_dynamic_devices() const { return m_dynamic_devices.size() > 0; } inline bool has_timestep_devices() const { return m_step_devices.size() > 0; } void update_forced(); void update_after(const netlist_time &after) { m_Q_sync.net().toggle_new_Q(); m_Q_sync.net().reschedule_in_queue(after); } /* netdevice functions */ NETLIB_UPDATEI(); NETLIB_RESETI(); public: int get_net_idx(net_t *net); plib::plog_base &log() { return netlist().log(); } virtual void log_stats(); virtual void create_solver_code(plib::postream &strm) { strm.writeline(plib::pfmt("/* {1} doesn't support static compile */")); } protected: void setup_base(analog_net_t::list_t &nets); void update_dynamic(); 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(); /* virtual */ void add_term(int net_idx, terminal_t *term); template void store(const T * RESTRICT V); template T delta(const T * RESTRICT V); template void build_LE_A(); template void build_LE_RHS(); std::vector m_terms; std::vector m_nets; std::vector> m_inps; std::vector m_rails_temp; const solver_parameters_t &m_params; state_var m_stat_calculations; state_var m_stat_newton_raphson; state_var m_stat_vsolver_calls; state_var m_iterative_fail; state_var m_iterative_total; inline nl_double current_timestep() { return m_cur_ts; } private: state_var m_last_step; state_var m_cur_ts; std::vector m_step_devices; std::vector m_dynamic_devices; logic_input_t m_fb_sync; logic_output_t m_Q_sync; /* calculate matrix */ void setup_matrix(); void step(const netlist_time &delta); void update_inputs(); const eSortType m_sort; }; template T matrix_solver_t::delta(const T * RESTRICT V) { /* FIXME: 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 unsigned iN = this->m_terms.size(); T cerr = 0; for (unsigned i = 0; i < iN; i++) cerr = std::max(cerr, std::abs(V[i] - (T) this->m_nets[i]->m_cur_Analog)); return cerr; } template 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]; } template void matrix_solver_t::build_LE_A() { static_assert(std::is_base_of::value, "T must derive from matrix_solver_t"); T &child = static_cast(*this); const unsigned iN = child.N(); for (unsigned k = 0; k < iN; k++) { for (unsigned i=0; i < iN; i++) child.A(k,i) = 0.0; const unsigned terms_count = m_terms[k]->count(); const unsigned railstart = m_terms[k]->m_railstart; const nl_double * RESTRICT gt = m_terms[k]->gt(); { nl_double akk = 0.0; for (unsigned i = 0; i < terms_count; i++) akk += gt[i]; child.A(k,k) = akk; } const nl_double * RESTRICT go = m_terms[k]->go(); const int * RESTRICT net_other = m_terms[k]->net_other(); for (unsigned i = 0; i < railstart; i++) child.A(k,net_other[i]) -= go[i]; } } template void matrix_solver_t::build_LE_RHS() { static_assert(std::is_base_of::value, "T must derive from matrix_solver_t"); T &child = static_cast(*this); const unsigned iN = child.N(); for (unsigned k = 0; k < iN; k++) { nl_double rhsk_a = 0.0; nl_double rhsk_b = 0.0; 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 (unsigned i = 0; i < terms_count; i++) rhsk_a = rhsk_a + Idr[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]; child.RHS(k) = rhsk_a + rhsk_b; } } } //namespace devices } // namespace netlist #endif /* NLD_MS_DIRECT_H_ */