diff options
Diffstat (limited to 'src/lib/netlist/solver/nld_matrix_solver.h')
-rw-r--r-- | src/lib/netlist/solver/nld_matrix_solver.h | 322 |
1 files changed, 24 insertions, 298 deletions
diff --git a/src/lib/netlist/solver/nld_matrix_solver.h b/src/lib/netlist/solver/nld_matrix_solver.h index 606ccbb9ee2..6a769775248 100644 --- a/src/lib/netlist/solver/nld_matrix_solver.h +++ b/src/lib/netlist/solver/nld_matrix_solver.h @@ -14,6 +14,7 @@ #include "netlist/plib/palloc.h" #include "netlist/plib/penum.h" #include "netlist/plib/pmatrix2d.h" +#include "netlist/plib/pmempool.h" #include "netlist/plib/putil.h" #include "netlist/plib/vector_ops.h" @@ -183,6 +184,8 @@ namespace solver { public: using list_t = std::vector<matrix_solver_t *>; + using fptype = nl_fptype; + using arena_type = plib::mempool_arena<plib::aligned_arena>; // after every call to solve, update inputs must be called. // this can be done as well as a batch to ease parallel processing. @@ -259,57 +262,31 @@ namespace solver } // return number of floating point operations for solve - std::size_t ops() const { return m_ops; } + constexpr std::size_t ops() const { return m_ops; } protected: - template <typename T> - using aligned_alloc = plib::aligned_allocator<T, PALIGN_VECTOROPT>; - matrix_solver_t(netlist_state_t &anetlist, const pstring &name, const analog_net_t::list_t &nets, const solver_parameters_t *params); virtual void vsolve_non_dynamic() = 0; - virtual netlist_time compute_next_timestep(nl_fptype cur_ts, nl_fptype max_ts) = 0; - virtual bool check_err() = 0; + virtual netlist_time compute_next_timestep(fptype cur_ts, fptype max_ts) = 0; + virtual bool check_err() const = 0; virtual void store() = 0; -#if 0 - plib::pmatrix2d<nl_fptype, aligned_alloc<nl_fptype>> m_gonn; - plib::pmatrix2d<nl_fptype, aligned_alloc<nl_fptype>> m_gtn; - plib::pmatrix2d<nl_fptype, aligned_alloc<nl_fptype>> m_Idrn; - plib::pmatrix2d<nl_fptype *, aligned_alloc<nl_fptype *>> m_connected_net_Vn; -#else - plib::pmatrix2d_vrl<nl_fptype, aligned_alloc<nl_fptype>> m_gonn; - plib::pmatrix2d_vrl<nl_fptype, aligned_alloc<nl_fptype>> m_gtn; - plib::pmatrix2d_vrl<nl_fptype, aligned_alloc<nl_fptype>> m_Idrn; - plib::pmatrix2d_vrl<nl_fptype *, aligned_alloc<nl_fptype *>> m_connected_net_Vn; -#endif - - plib::aligned_vector<terms_for_net_t> m_terms; + plib::pmatrix2d_vrl<fptype, arena_type> m_gonn; + plib::pmatrix2d_vrl<fptype, arena_type> m_gtn; + plib::pmatrix2d_vrl<fptype, arena_type> m_Idrn; + plib::pmatrix2d_vrl<fptype *, arena_type> m_connected_net_Vn; const solver_parameters_t &m_params; state_var<std::size_t> m_iterative_fail; state_var<std::size_t> m_iterative_total; - private: - - plib::aligned_vector<terms_for_net_t> m_rails_temp; - std::vector<device_arena::unique_ptr<proxied_analog_output_t>> m_inps; - - state_var<std::size_t> m_stat_calculations; - state_var<std::size_t> m_stat_newton_raphson; - state_var<std::size_t> m_stat_vsolver_calls; - - state_var<netlist_time_ext> m_last_step; - plib::aligned_vector<nldelegate_ts> m_step_funcs; - plib::aligned_vector<nldelegate_dyn> m_dynamic_funcs; + plib::aligned_vector<terms_for_net_t> m_terms; // setup only - logic_input_t m_fb_sync; - logic_output_t m_Q_sync; - - std::size_t m_ops; + private: // base setup - called from constructor void setup_base(setup_t &setup, const analog_net_t::list_t &nets) noexcept(false); @@ -321,7 +298,7 @@ namespace solver int get_net_idx(const analog_net_t *net) const noexcept; std::pair<int, int> get_left_right_of_diag(std::size_t irow, std::size_t idiag); - nl_fptype get_weight_around_diag(std::size_t row, std::size_t diag); + fptype get_weight_around_diag(std::size_t row, std::size_t diag); void add_term(std::size_t net_idx, terminal_t *term) noexcept(false); @@ -332,272 +309,21 @@ namespace solver analog_net_t *get_connected_net(terminal_t *term); - }; - - template <typename FT, int SIZE> - class matrix_solver_ext_t: public matrix_solver_t - { - public: - - using float_type = FT; - - matrix_solver_ext_t(netlist_state_t &anetlist, const pstring &name, - const analog_net_t::list_t &nets, - const solver_parameters_t *params, const std::size_t size) - : matrix_solver_t(anetlist, name, nets, params) - , m_dim(size) - , m_new_V(size) - , m_RHS(size) - , m_mat_ptr(size, this->max_railstart() + 1) - , m_last_V(size, nlconst::zero()) - , m_DD_n_m_1(size, nlconst::zero()) - , m_h_n_m_1(size, nlconst::magic(1e-6)) // we need a non zero value here - { - // - // save states - // - state().save(*this, m_last_V.as_base(), this->name(), "m_last_V"); - state().save(*this, m_DD_n_m_1.as_base(), this->name(), "m_DD_n_m_1"); - state().save(*this, m_h_n_m_1.as_base(), this->name(), "m_h_n_m_1"); - } - - - private: - const std::size_t m_dim; - - protected: - static constexpr const std::size_t SIZEABS = plib::parray<FT, SIZE>::SIZEABS(); - static constexpr const std::size_t m_pitch_ABS = (((SIZEABS + 0) + 7) / 8) * 8; - - PALIGNAS_VECTOROPT() - plib::parray<float_type, SIZE> m_new_V; - PALIGNAS_VECTOROPT() - plib::parray<float_type, SIZE> m_RHS; - - PALIGNAS_VECTOROPT() - plib::pmatrix2d<float_type *> m_mat_ptr; - - // FIXME: below should be private - // state - variable time_stepping - PALIGNAS_VECTOROPT() - plib::parray<nl_fptype, SIZE> m_last_V; - PALIGNAS_VECTOROPT() - plib::parray<nl_fptype, SIZE> m_DD_n_m_1; - PALIGNAS_VECTOROPT() - plib::parray<nl_fptype, SIZE> m_h_n_m_1; - - std::size_t max_railstart() const noexcept - { - std::size_t max_rail = 0; - for (std::size_t k = 0; k < m_terms.size(); k++) - max_rail = std::max(max_rail, m_terms[k].railstart()); - return max_rail; - } - - - template <typename T, typename M> - void log_fill(const T &fill, M &mat) - { - const std::size_t iN = fill.size(); - - // FIXME: Not yet working, mat_cr.h needs some more work -#if 0 - auto mat_GE = dynamic_cast<plib::pGEmatrix_cr_t<typename M::base> *>(&mat); -#else - plib::unused_var(mat); -#endif - std::vector<unsigned> levL(iN, 0); - std::vector<unsigned> levU(iN, 0); - - // parallel scheme for L x = y - for (std::size_t k = 0; k < iN; k++) - { - unsigned lm=0; - for (std::size_t j = 0; j<k; j++) - if (fill[k][j] < M::FILL_INFINITY) - lm = std::max(lm, levL[j]); - levL[k] = 1+lm; - } - - // parallel scheme for U x = y - for (std::size_t k = iN; k-- > 0; ) - { - unsigned lm=0; - for (std::size_t j = iN; --j > k; ) - if (fill[k][j] < M::FILL_INFINITY) - lm = std::max(lm, levU[j]); - levU[k] = 1+lm; - } - for (std::size_t k = 0; k < iN; k++) - { - unsigned fm = 0; - pstring ml = ""; - for (std::size_t j = 0; j < iN; j++) - { - ml += fill[k][j] == 0 ? 'X' : fill[k][j] < M::FILL_INFINITY ? '+' : '.'; - if (fill[k][j] < M::FILL_INFINITY) - if (fill[k][j] > fm) - fm = fill[k][j]; - } -#if 0 - this->log().verbose("{1:4} {2} {3:4} {4:4} {5:4} {6:4}", k, ml, - levL[k], levU[k], mat_GE ? mat_GE->get_parallel_level(k) : 0, fm); -#else - this->log().verbose("{1:4} {2} {3:4} {4:4} {5:4} {6:4}", k, ml, - levL[k], levU[k], 0, fm); -#endif - } - } - - constexpr std::size_t size() const noexcept - { - return (SIZE > 0) ? static_cast<std::size_t>(SIZE) : m_dim; - } - -#if 1 - void store() override - { - const std::size_t iN = size(); - for (std::size_t i = 0; i < iN; i++) - this->m_terms[i].setV(static_cast<nl_fptype>(m_new_V[i])); - } -#else - // global tanh damping (4.197) - // partially cures the symptoms but not the cause - void store() override - { - const std::size_t iN = size(); - for (std::size_t i = 0; i < iN; i++) - { - auto oldV = this->m_terms[i].template getV<nl_fptype>(); - this->m_terms[i].setV(oldV + 0.02 * plib::tanh((m_new_V[i]-oldV)*50.0)); - } - } -#endif - bool check_err() override - { - // 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 = size(); - const auto reltol(static_cast<float_type>(m_params.m_reltol)); - const auto vntol(static_cast<float_type>(m_params.m_vntol)); - for (std::size_t i = 0; i < iN; i++) - { - const auto vold(static_cast<float_type>(this->m_terms[i].getV())); - const auto vnew(m_new_V[i]); - const auto tol(vntol + reltol * std::max(plib::abs(vnew),plib::abs(vold))); - if (plib::abs(vnew - vold) > tol) - return true; - } - return false; - } - - netlist_time compute_next_timestep(nl_fptype cur_ts, nl_fptype max_ts) override - { - nl_fptype new_solver_timestep(max_ts); - - for (std::size_t k = 0; k < size(); k++) - { - const auto &t = m_terms[k]; - const auto v(static_cast<nl_fptype>(t.getV())); - // avoid floating point exceptions - const nl_fptype DD_n = std::max(-fp_constants<nl_fptype>::TIMESTEP_MAXDIFF(), - std::min(+fp_constants<nl_fptype>::TIMESTEP_MAXDIFF(),(v - m_last_V[k]))); - - m_last_V[k] = v; - const nl_fptype hn = cur_ts; - - nl_fptype DD2 = (DD_n / hn - m_DD_n_m_1[k] / m_h_n_m_1[k]) / (hn + m_h_n_m_1[k]); - nl_fptype new_net_timestep(0); - - m_h_n_m_1[k] = hn; - m_DD_n_m_1[k] = DD_n; - if (plib::abs(DD2) > fp_constants<nl_fptype>::TIMESTEP_MINDIV()) // avoid div-by-zero - new_net_timestep = plib::sqrt(m_params.m_dynamic_lte / plib::abs(nlconst::half()*DD2)); - else - new_net_timestep = m_params.m_max_timestep; - - new_solver_timestep = std::min(new_net_timestep, new_solver_timestep); - } - new_solver_timestep = std::max(new_solver_timestep, m_params.m_min_timestep); - - // FIXME: Factor 2 below is important. Without, we get timing issues. This must be a bug elsewhere. - return std::max(netlist_time::from_fp(new_solver_timestep), netlist_time::quantum() * 2); - } - - template <typename M> - void build_mat_ptr(M &mat) - { - const std::size_t iN = size(); - - for (std::size_t k=0; k<iN; k++) - { - std::size_t cnt(0); - // build pointers into the compressed row format matrix for each terminal - for (std::size_t j=0; j< this->m_terms[k].railstart();j++) - { - int other = this->m_terms[k].m_connected_net_idx[j]; - if (other >= 0) - { - m_mat_ptr[k][j] = &(mat[k][static_cast<std::size_t>(other)]); - cnt++; - } - } - nl_assert_always(cnt == this->m_terms[k].railstart(), "Count and railstart mismatch"); - m_mat_ptr[k][this->m_terms[k].railstart()] = &(mat[k][k]); - } - } - - template <typename M> - void clear_square_mat(M &m) - { - const std::size_t n = size(); - for (std::size_t k=0; k < n; k++) - { - auto *p = &(m[k][0]); - using mat_elem_type = typename std::decay<decltype(*p)>::type; - for (std::size_t i=0; i < n; i++) - p[i] = plib::constants<mat_elem_type>::zero(); - } - } - - void fill_matrix_and_rhs() - { - const std::size_t N = size(); - - for (std::size_t k = 0; k < N; k++) - { - auto &net = m_terms[k]; - auto **tcr_r = &(m_mat_ptr[k][0]); - - using source_type = typename decltype(m_gtn)::value_type; - const std::size_t term_count = net.count(); - const std::size_t railstart = net.railstart(); - const auto &go = m_gonn[k]; - const auto > = m_gtn[k]; - const auto &Idr = m_Idrn[k]; - const auto &cnV = m_connected_net_Vn[k]; - - // FIXME: gonn, gtn and Idr - which float types should they have? - - auto gtot_t = std::accumulate(gt, gt + term_count, plib::constants<source_type>::zero()); - - // update diagonal element ... - *tcr_r[railstart] = static_cast<FT>(gtot_t); //mat.A[mat.diag[k]] += gtot_t; + state_var<std::size_t> m_stat_calculations; + state_var<std::size_t> m_stat_newton_raphson; + state_var<std::size_t> m_stat_vsolver_calls; - for (std::size_t i = 0; i < railstart; i++) - *tcr_r[i] += static_cast<FT>(go[i]); + state_var<netlist_time_ext> m_last_step; + plib::aligned_vector<nldelegate_ts> m_step_funcs; + plib::aligned_vector<nldelegate_dyn> m_dynamic_funcs; + plib::aligned_vector<device_arena::unique_ptr<proxied_analog_output_t>> m_inps; - auto RHS_t(std::accumulate(Idr, Idr + term_count, plib::constants<source_type>::zero())); + logic_input_t m_fb_sync; + logic_output_t m_Q_sync; - for (std::size_t i = railstart; i < term_count; i++) - RHS_t += (- go[i]) * *cnV[i]; + std::size_t m_ops; - m_RHS[k] = static_cast<FT>(RHS_t); - } - } + plib::aligned_vector<terms_for_net_t> m_rails_temp; // setup only }; |