From 65fb297023d63e41e996452c74db26e27a305b65 Mon Sep 17 00:00:00 2001 From: couriersud Date: Thu, 31 Oct 2019 01:15:43 +0100 Subject: netlist: move solver stuff into separate namespace. (nw) - new namespace "solver" - minor ptime modifications --- src/lib/netlist/nl_base.h | 12 +- src/lib/netlist/plib/ptime.h | 26 +- src/lib/netlist/solver/nld_matrix_solver.cpp | 38 +-- src/lib/netlist/solver/nld_matrix_solver.h | 10 +- src/lib/netlist/solver/nld_ms_direct.h | 12 +- src/lib/netlist/solver/nld_ms_direct1.h | 4 +- src/lib/netlist/solver/nld_ms_direct2.h | 4 +- src/lib/netlist/solver/nld_ms_gcr.h | 4 +- src/lib/netlist/solver/nld_ms_gmres.h | 4 +- src/lib/netlist/solver/nld_ms_sm.h | 4 +- src/lib/netlist/solver/nld_ms_sor.h | 224 ++++++------ src/lib/netlist/solver/nld_ms_sor_mat.h | 4 +- src/lib/netlist/solver/nld_ms_w.h | 488 +++++++++++++-------------- src/lib/netlist/solver/nld_solver.cpp | 46 +-- src/lib/netlist/solver/nld_solver.h | 10 +- 15 files changed, 446 insertions(+), 444 deletions(-) diff --git a/src/lib/netlist/nl_base.h b/src/lib/netlist/nl_base.h index 54e6622799e..1ed7b645c0e 100644 --- a/src/lib/netlist/nl_base.h +++ b/src/lib/netlist/nl_base.h @@ -186,7 +186,6 @@ namespace netlist namespace devices { - class matrix_solver_t; class NETLIB_NAME(solver); class NETLIB_NAME(mainclock); class NETLIB_NAME(base_proxy); @@ -194,6 +193,11 @@ namespace netlist class NETLIB_NAME(base_a_to_d_proxy); } // namespace devices + namespace solver + { + class matrix_solver_t; + } // namespace solver + class logic_output_t; class logic_input_t; class analog_net_t; @@ -902,12 +906,12 @@ namespace netlist nl_double *Q_Analog_state_ptr() NL_NOEXCEPT { return m_cur_Analog.ptr(); } //FIXME: needed by current solver code - devices::matrix_solver_t *solver() const noexcept { return m_solver; } - void set_solver(devices::matrix_solver_t *solver) noexcept { m_solver = solver; } + solver::matrix_solver_t *solver() const noexcept { return m_solver; } + void set_solver(solver::matrix_solver_t *solver) noexcept { m_solver = solver; } private: state_var m_cur_Analog; - devices::matrix_solver_t *m_solver; + solver::matrix_solver_t *m_solver; }; // ----------------------------------------------------------------------------- diff --git a/src/lib/netlist/plib/ptime.h b/src/lib/netlist/plib/ptime.h index 9668e0cf428..4044f518e24 100644 --- a/src/lib/netlist/plib/ptime.h +++ b/src/lib/netlist/plib/ptime.h @@ -58,7 +58,7 @@ namespace plib return ptime(lhs.m_time + rhs.m_time); } - friend constexpr const ptime operator*(ptime lhs, const mult_type &factor) noexcept + friend constexpr const ptime operator*(ptime lhs, const mult_type factor) noexcept { return ptime(lhs.m_time * factor); } @@ -107,18 +107,18 @@ namespace plib // for save states .... C14CONSTEXPR internal_type *get_internaltype_ptr() noexcept { return &m_time; } - static constexpr ptime from_nsec(const internal_type ns) noexcept { return ptime(ns, UINT64_C(1000000000)); } - static constexpr ptime from_usec(const internal_type us) noexcept { return ptime(us, UINT64_C( 1000000)); } - static constexpr ptime from_msec(const internal_type ms) noexcept { return ptime(ms, UINT64_C( 1000)); } - static constexpr ptime from_sec(const internal_type s) noexcept { return ptime(s, UINT64_C( 1)); } - static constexpr ptime from_hz(const internal_type hz) noexcept { return ptime(1 , hz); } - static constexpr ptime from_raw(const internal_type raw) noexcept { return ptime(raw); } - static constexpr ptime from_double(const double t) noexcept { return ptime(static_cast(std::floor(t * static_cast(RES) + 0.5)), RES); } - - static constexpr ptime zero() noexcept { return ptime(0, RES); } - static constexpr ptime quantum() noexcept { return ptime(1, RES); } - static constexpr ptime never() noexcept { return ptime(plib::numeric_limits::max(), RES); } - static constexpr internal_type resolution() noexcept { return RES; } + static constexpr const ptime from_nsec(const internal_type ns) noexcept { return ptime(ns, UINT64_C(1000000000)); } + static constexpr const ptime from_usec(const internal_type us) noexcept { return ptime(us, UINT64_C( 1000000)); } + static constexpr const ptime from_msec(const internal_type ms) noexcept { return ptime(ms, UINT64_C( 1000)); } + static constexpr const ptime from_sec(const internal_type s) noexcept { return ptime(s, UINT64_C( 1)); } + static constexpr const ptime from_hz(const internal_type hz) noexcept { return ptime(1 , hz); } + static constexpr const ptime from_raw(const internal_type raw) noexcept { return ptime(raw); } + static constexpr const ptime from_double(const double t) noexcept { return ptime(static_cast(std::floor(t * static_cast(RES) + 0.5)), RES); } + + static constexpr const ptime zero() noexcept { return ptime(0, RES); } + static constexpr const ptime quantum() noexcept { return ptime(1, RES); } + static constexpr const ptime never() noexcept { return ptime(plib::numeric_limits::max(), RES); } + static constexpr const internal_type resolution() noexcept { return RES; } constexpr internal_type in_nsec() const noexcept { return m_time / (RES / UINT64_C(1000000000)); } constexpr internal_type in_usec() const noexcept { return m_time / (RES / UINT64_C( 1000000)); } diff --git a/src/lib/netlist/solver/nld_matrix_solver.cpp b/src/lib/netlist/solver/nld_matrix_solver.cpp index b762749dc65..c8dfa4ba384 100644 --- a/src/lib/netlist/solver/nld_matrix_solver.cpp +++ b/src/lib/netlist/solver/nld_matrix_solver.cpp @@ -12,7 +12,7 @@ namespace netlist { -namespace devices +namespace solver { terms_for_net_t::terms_for_net_t(analog_net_t * net) @@ -58,11 +58,13 @@ namespace devices { connect_post_start(m_fb_sync, m_Q_sync); setup_base(nets); + + /* now setup the matrix */ + setup_matrix(); } void matrix_solver_t::setup_base(const analog_net_t::list_t &nets) { - log().debug("New solver setup\n"); m_terms.clear(); @@ -130,9 +132,6 @@ namespace devices } log().debug("added net with {1} populated connections\n", net->core_terms().size()); } - - /* now setup the matrix */ - setup_matrix(); } void matrix_solver_t::sort_terms(matrix_sort_type_e sort) @@ -417,8 +416,19 @@ namespace devices d->timestep(dd); } - void matrix_solver_t::solve_base() + const netlist_time matrix_solver_t::solve(netlist_time now) { + const netlist_time delta = now - m_last_step; + + // We are already up to date. Avoid oscillations. + // FIXME: Make this a parameter! + if (delta < netlist_time::quantum()) + return netlist_time::zero(); + + /* update all terminals for new time step */ + m_last_step = now; + step(delta); + ++m_stat_vsolver_calls; if (has_dynamic_devices()) { @@ -444,21 +454,7 @@ namespace devices { this->vsolve_non_dynamic(false); } - } - - const netlist_time matrix_solver_t::solve(netlist_time now) - { - const netlist_time delta = now - m_last_step; - // We are already up to date. Avoid oscillations. - // FIXME: Make this a parameter! - if (delta < netlist_time::quantum()) - return netlist_time::zero(); - - /* update all terminals for new time step */ - m_last_step = now; - step(delta); - solve_base(); const netlist_time next_time_step = compute_next_timestep(delta.as_double()); return next_time_step; @@ -622,6 +618,6 @@ namespace devices } } -} // namespace devices +} // namespace solver } // namespace netlist diff --git a/src/lib/netlist/solver/nld_matrix_solver.h b/src/lib/netlist/solver/nld_matrix_solver.h index cd04d43bdbd..07fafd1ece1 100644 --- a/src/lib/netlist/solver/nld_matrix_solver.h +++ b/src/lib/netlist/solver/nld_matrix_solver.h @@ -20,7 +20,7 @@ namespace netlist { -namespace devices +namespace solver { P_ENUM(matrix_sort_type_e, NOSORT, @@ -162,8 +162,6 @@ namespace devices public: using list_t = std::vector; - 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. */ @@ -481,17 +479,17 @@ namespace devices const float_type * const * other_cur_analog = m_connected_net_Vn[k]; for (std::size_t i = 0; i < terms_count; i++) - rhsk_a = rhsk_a + Idr[i]; + rhsk_a += Idr[i]; for (std::size_t i = m_terms[k].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]; + rhsk_b += - go[i] * *other_cur_analog[i]; child.RHS(k) = rhsk_a + rhsk_b; } } -} //namespace devices +} // namespace solver } // namespace netlist #endif /* NLD_MS_DIRECT_H_ */ diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h index fd745d01d22..582a7996947 100644 --- a/src/lib/netlist/solver/nld_ms_direct.h +++ b/src/lib/netlist/solver/nld_ms_direct.h @@ -18,7 +18,7 @@ namespace netlist { -namespace devices +namespace solver { template @@ -51,6 +51,8 @@ namespace devices const FT &RHS(std::size_t r) const noexcept { return m_A[r][size()]; } FT &RHS(std::size_t r) noexcept { return m_A[r][size()]; } + + PALIGNAS_VECTOROPT() plib::parray m_new_V; private: @@ -59,6 +61,8 @@ namespace devices const std::size_t m_dim; const std::size_t m_pitch; + + PALIGNAS_VECTOROPT() plib::parray2D m_A; }; @@ -79,10 +83,10 @@ namespace devices const auto &nzrd = m_terms[i].m_nzrd; const auto &nzbd = m_terms[i].m_nzbd; - for (auto j : nzbd) + for (auto &j : nzbd) { const FT f1 = -f * A(j, i); - for (auto k : nzrd) + for (auto &k : nzrd) A(j, k) += A(i, k) * f1; //RHS(j) += RHS(i) * f1; } @@ -215,7 +219,7 @@ namespace devices state().save(*this, RHS(k), this->name(), plib::pfmt("RHS.{1}")(k)); } -} // namespace devices +} // namespace solver } // namespace netlist #endif /* NLD_MS_DIRECT_H_ */ diff --git a/src/lib/netlist/solver/nld_ms_direct1.h b/src/lib/netlist/solver/nld_ms_direct1.h index 1fefc9045b7..c6ebd77e198 100644 --- a/src/lib/netlist/solver/nld_ms_direct1.h +++ b/src/lib/netlist/solver/nld_ms_direct1.h @@ -13,7 +13,7 @@ namespace netlist { -namespace devices +namespace solver { template class matrix_solver_direct1_t: public matrix_solver_direct_t @@ -49,7 +49,7 @@ namespace devices -} //namespace devices +} // namespace solver } // namespace netlist diff --git a/src/lib/netlist/solver/nld_ms_direct2.h b/src/lib/netlist/solver/nld_ms_direct2.h index bca20eb583b..20d3bb659e5 100644 --- a/src/lib/netlist/solver/nld_ms_direct2.h +++ b/src/lib/netlist/solver/nld_ms_direct2.h @@ -13,7 +13,7 @@ namespace netlist { -namespace devices +namespace solver { // ---------------------------------------------------------------------------------------- @@ -54,7 +54,7 @@ namespace devices }; -} //namespace devices +} // namespace solver } // namespace netlist #endif /* NLD_MS_DIRECT2_H_ */ diff --git a/src/lib/netlist/solver/nld_ms_gcr.h b/src/lib/netlist/solver/nld_ms_gcr.h index edb17ad9e73..3aa39b7bb54 100644 --- a/src/lib/netlist/solver/nld_ms_gcr.h +++ b/src/lib/netlist/solver/nld_ms_gcr.h @@ -22,7 +22,7 @@ namespace netlist { -namespace devices +namespace solver { template @@ -262,7 +262,7 @@ namespace devices return (err > this->m_params.m_accuracy) ? 2 : 1; } -} // namespace devices +} // namespace solver } // namespace netlist #endif /* NLD_MS_GCR_H_ */ diff --git a/src/lib/netlist/solver/nld_ms_gmres.h b/src/lib/netlist/solver/nld_ms_gmres.h index 7a12737297c..b5e3f605cc0 100644 --- a/src/lib/netlist/solver/nld_ms_gmres.h +++ b/src/lib/netlist/solver/nld_ms_gmres.h @@ -21,7 +21,7 @@ namespace netlist { -namespace devices +namespace solver { template @@ -135,7 +135,7 @@ namespace devices -} // namespace devices +} // namespace solver } // namespace netlist #endif /* NLD_MS_GMRES_H_ */ diff --git a/src/lib/netlist/solver/nld_ms_sm.h b/src/lib/netlist/solver/nld_ms_sm.h index a6bc83e27a6..30ad2031628 100644 --- a/src/lib/netlist/solver/nld_ms_sm.h +++ b/src/lib/netlist/solver/nld_ms_sm.h @@ -41,7 +41,7 @@ namespace netlist { -namespace devices +namespace solver { template @@ -297,7 +297,7 @@ namespace devices } -} // namespace devices +} // namespace solver } // namespace netlist #endif /* NLD_MS_DIRECT_H_ */ diff --git a/src/lib/netlist/solver/nld_ms_sor.h b/src/lib/netlist/solver/nld_ms_sor.h index 41ab4460ced..98e2c8ddc6d 100644 --- a/src/lib/netlist/solver/nld_ms_sor.h +++ b/src/lib/netlist/solver/nld_ms_sor.h @@ -19,151 +19,151 @@ namespace netlist { - namespace devices +namespace solver { -template -class matrix_solver_SOR_t: public matrix_solver_direct_t -{ -public: - - using float_type = FT; - - matrix_solver_SOR_t(netlist_state_t &anetlist, const pstring &name, - analog_net_t::list_t &nets, - const solver_parameters_t *params, const std::size_t size) - : matrix_solver_direct_t(anetlist, name, nets, params, size) - , m_lp_fact(*this, "m_lp_fact", 0) - , w(size, 0.0) - , one_m_w(size, 0.0) - , RHS(size, 0.0) - //, new_V(size, 0.0) - { - } - - unsigned vsolve_non_dynamic(const bool newton_raphson) override; - -private: - state_var m_lp_fact; - std::vector w; - std::vector one_m_w; - std::vector RHS; - //std::vector new_V; -}; + template + class matrix_solver_SOR_t: public matrix_solver_direct_t + { + public: + + using float_type = FT; + + matrix_solver_SOR_t(netlist_state_t &anetlist, const pstring &name, + analog_net_t::list_t &nets, + const solver_parameters_t *params, const std::size_t size) + : matrix_solver_direct_t(anetlist, name, nets, params, size) + , m_lp_fact(*this, "m_lp_fact", 0) + , w(size, 0.0) + , one_m_w(size, 0.0) + , RHS(size, 0.0) + //, new_V(size, 0.0) + { + } -// ---------------------------------------------------------------------------------------- -// matrix_solver - Gauss - Seidel -// ---------------------------------------------------------------------------------------- + unsigned vsolve_non_dynamic(const bool newton_raphson) override; -template -unsigned matrix_solver_SOR_t::vsolve_non_dynamic(const bool newton_raphson) -{ - const std::size_t iN = this->size(); - bool resched = false; - unsigned resched_cnt = 0; + private: + state_var m_lp_fact; + std::vector w; + std::vector one_m_w; + std::vector RHS; + //std::vector new_V; + }; - /* ideally, we could get an estimate for the spectral radius of - * Inv(D - L) * U - * - * and estimate using - * - * omega = 2.0 / (1.0 + std::sqrt(1-rho)) - */ + // ---------------------------------------------------------------------------------------- + // matrix_solver - Gauss - Seidel + // ---------------------------------------------------------------------------------------- - const float_type ws = this->m_params.m_gs_sor; - - for (std::size_t k = 0; k < iN; k++) + template + unsigned matrix_solver_SOR_t::vsolve_non_dynamic(const bool newton_raphson) { - float_type gtot_t = 0.0; - float_type gabs_t = 0.0; - float_type RHS_t = 0.0; + const std::size_t iN = this->size(); + bool resched = false; + unsigned resched_cnt = 0; - const std::size_t term_count = this->m_terms[k].count(); - const float_type * const gt = this->m_gtn[k]; - const float_type * const go = this->m_gonn[k]; - const float_type * const Idr = this->m_Idrn[k]; - auto other_cur_analog = this->m_connected_net_Vn[k]; + /* ideally, we could get an estimate for the spectral radius of + * Inv(D - L) * U + * + * and estimate using + * + * omega = 2.0 / (1.0 + std::sqrt(1-rho)) + */ - this->m_new_V[k] = this->m_terms[k].getV(); + const float_type ws = this->m_params.m_gs_sor; - for (std::size_t i = 0; i < term_count; i++) + for (std::size_t k = 0; k < iN; k++) { - gtot_t = gtot_t + gt[i]; - RHS_t = RHS_t + Idr[i]; - } + float_type gtot_t = 0.0; + float_type gabs_t = 0.0; + float_type RHS_t = 0.0; - for (std::size_t i = this->m_terms[k].railstart(); i < term_count; i++) - RHS_t = RHS_t - go[i] * *other_cur_analog[i]; + const std::size_t term_count = this->m_terms[k].count(); + const float_type * const gt = this->m_gtn[k]; + const float_type * const go = this->m_gonn[k]; + const float_type * const Idr = this->m_Idrn[k]; + auto other_cur_analog = this->m_connected_net_Vn[k]; - RHS[k] = RHS_t; + this->m_new_V[k] = this->m_terms[k].getV(); - if (this->m_params.m_use_gabs) - { for (std::size_t i = 0; i < term_count; i++) - gabs_t = gabs_t + std::abs(go[i]); + { + gtot_t = gtot_t + gt[i]; + RHS_t = RHS_t + Idr[i]; + } + + for (std::size_t i = this->m_terms[k].railstart(); i < term_count; i++) + RHS_t = RHS_t - go[i] * *other_cur_analog[i]; + + RHS[k] = RHS_t; - gabs_t *= plib::constants::cast(0.5); // derived by try and error - if (gabs_t <= gtot_t) + if (this->m_params.m_use_gabs) { - w[k] = ws / gtot_t; - one_m_w[k] = plib::constants::one() - ws; + for (std::size_t i = 0; i < term_count; i++) + gabs_t = gabs_t + std::abs(go[i]); + + gabs_t *= plib::constants::cast(0.5); // derived by try and error + if (gabs_t <= gtot_t) + { + w[k] = ws / gtot_t; + one_m_w[k] = plib::constants::one() - ws; + } + else + { + w[k] = plib::constants::one() / (gtot_t + gabs_t); + one_m_w[k] = plib::constants::one() - plib::constants::one() * gtot_t / (gtot_t + gabs_t); + } } else { - w[k] = plib::constants::one() / (gtot_t + gabs_t); - one_m_w[k] = plib::constants::one() - plib::constants::one() * gtot_t / (gtot_t + gabs_t); + w[k] = ws / gtot_t; + one_m_w[k] = plib::constants::one() - ws; } } - else - { - w[k] = ws / gtot_t; - one_m_w[k] = plib::constants::one() - ws; - } - } - const float_type accuracy = this->m_params.m_accuracy; + const float_type accuracy = this->m_params.m_accuracy; - do { - resched = false; - float_type err = 0; - for (std::size_t k = 0; k < iN; k++) - { - const int * net_other = this->m_terms[k].m_connected_net_idx.data(); - const std::size_t railstart = this->m_terms[k].railstart(); - const float_type * go = this->m_gonn[k]; + do { + resched = false; + float_type err = 0; + for (std::size_t k = 0; k < iN; k++) + { + const int * net_other = this->m_terms[k].m_connected_net_idx.data(); + const std::size_t railstart = this->m_terms[k].railstart(); + const float_type * go = this->m_gonn[k]; - float_type Idrive = 0.0; - for (std::size_t i = 0; i < railstart; i++) - Idrive = Idrive - go[i] * this->m_new_V[static_cast(net_other[i])]; + float_type Idrive = 0.0; + for (std::size_t i = 0; i < railstart; i++) + Idrive = Idrive - go[i] * this->m_new_V[static_cast(net_other[i])]; - const float_type new_val = this->m_new_V[k] * one_m_w[k] + (Idrive + RHS[k]) * w[k]; + const float_type new_val = this->m_new_V[k] * one_m_w[k] + (Idrive + RHS[k]) * w[k]; - err = std::max(std::abs(new_val - this->m_new_V[k]), err); - this->m_new_V[k] = new_val; - } + err = std::max(std::abs(new_val - this->m_new_V[k]), err); + this->m_new_V[k] = new_val; + } - if (err > accuracy) - resched = true; + if (err > accuracy) + resched = true; - resched_cnt++; - } while (resched && (resched_cnt < this->m_params.m_gs_loops)); + resched_cnt++; + } while (resched && (resched_cnt < this->m_params.m_gs_loops)); - this->m_iterative_total += resched_cnt; - this->m_stat_calculations++; + this->m_iterative_total += resched_cnt; + this->m_stat_calculations++; - if (resched) - { - // Fallback to direct solver ... - this->m_iterative_fail++; - return matrix_solver_direct_t::vsolve_non_dynamic(newton_raphson); - } + if (resched) + { + // Fallback to direct solver ... + this->m_iterative_fail++; + return matrix_solver_direct_t::vsolve_non_dynamic(newton_raphson); + } - const float_type err = (newton_raphson ? this->delta(this->m_new_V) : 0.0); - this->store(this->m_new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; -} + const float_type err = (newton_raphson ? this->delta(this->m_new_V) : 0.0); + this->store(this->m_new_V); + return (err > this->m_params.m_accuracy) ? 2 : 1; + } - } //namespace devices +} // namespace solver } // namespace netlist #endif /* NLD_MS_SOR_H_ */ diff --git a/src/lib/netlist/solver/nld_ms_sor_mat.h b/src/lib/netlist/solver/nld_ms_sor_mat.h index 943ae1ad989..2b705858854 100644 --- a/src/lib/netlist/solver/nld_ms_sor_mat.h +++ b/src/lib/netlist/solver/nld_ms_sor_mat.h @@ -20,7 +20,7 @@ namespace netlist { -namespace devices +namespace solver { template @@ -218,7 +218,7 @@ namespace devices } -} // namespace devices +} // namespace solver } // namespace netlist #endif /* NLD_MS_GAUSS_SEIDEL_H_ */ diff --git a/src/lib/netlist/solver/nld_ms_w.h b/src/lib/netlist/solver/nld_ms_w.h index 33e15db1a49..0e627c86a72 100644 --- a/src/lib/netlist/solver/nld_ms_w.h +++ b/src/lib/netlist/solver/nld_ms_w.h @@ -48,325 +48,325 @@ namespace netlist { - namespace devices - { - -template -class matrix_solver_w_t: public matrix_solver_t +namespace solver { - friend class matrix_solver_t; -public: - using float_ext_type = FT; - using float_type = FT; + template + class matrix_solver_w_t: public matrix_solver_t + { + friend class matrix_solver_t; - // FIXME: dirty hack to make this compile - static constexpr const std::size_t storage_N = 100; + public: + using float_ext_type = FT; + using float_type = FT; - matrix_solver_w_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_cnt(0) - , m_dim(size) - { - // FIXME: This shouldn't be necessary, recalculate on each entry ... - for (std::size_t k = 0; k < this->size(); k++) - state().save(*this, RHS(k), this->name(), plib::pfmt("RHS.{1}")(k)); - } + // FIXME: dirty hack to make this compile + static constexpr const std::size_t storage_N = 100; - void reset() override { matrix_solver_t::reset(); } + matrix_solver_w_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_cnt(0) + , m_dim(size) + { + // FIXME: This shouldn't be necessary, recalculate on each entry ... + for (std::size_t k = 0; k < this->size(); k++) + state().save(*this, RHS(k), this->name(), plib::pfmt("RHS.{1}")(k)); + } -protected: - unsigned vsolve_non_dynamic(const bool newton_raphson) override; - unsigned solve_non_dynamic(const bool newton_raphson); + void reset() override { matrix_solver_t::reset(); } - constexpr std::size_t size() const { return m_dim; } + protected: + unsigned vsolve_non_dynamic(const bool newton_raphson) override; + unsigned solve_non_dynamic(const bool newton_raphson); - void LE_invert(); + constexpr std::size_t size() const { return m_dim; } - template - void LE_compute_x(T & x); + void LE_invert(); + template + void LE_compute_x(T & x); - template - float_ext_type &A(const T1 &r, const T2 &c) { return m_A[r][c]; } - template - float_ext_type &W(const T1 &r, const T2 &c) { return m_W[r][c]; } - /* access to Ainv for fixed columns over row, there store transposed */ - template - float_ext_type &Ainv(const T1 &r, const T2 &c) { return m_Ainv[c][r]; } - template - float_ext_type &RHS(const T1 &r) { return m_RHS[r]; } + template + float_ext_type &A(const T1 &r, const T2 &c) { return m_A[r][c]; } + template + float_ext_type &W(const T1 &r, const T2 &c) { return m_W[r][c]; } + /* access to Ainv for fixed columns over row, there store transposed */ + template + float_ext_type &Ainv(const T1 &r, const T2 &c) { return m_Ainv[c][r]; } + template + float_ext_type &RHS(const T1 &r) { return m_RHS[r]; } - template - float_ext_type &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; } + template + float_ext_type &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; } -private: - template - using array2D = std::array, N>; - static constexpr std::size_t m_pitch = ((( storage_N) + 7) / 8) * 8; - array2D m_A; - array2D m_Ainv; - array2D m_W; - std::array m_RHS; // right hand side - contains currents - array2D m_lA; + private: + template + using array2D = std::array, N>; + static constexpr std::size_t m_pitch = ((( storage_N) + 7) / 8) * 8; + array2D m_A; + array2D m_Ainv; + array2D m_W; + std::array m_RHS; // right hand side - contains currents - /* temporary */ - array2D H; - std::array rows; - array2D cols; - std::array colcount; + array2D m_lA; - unsigned m_cnt; + /* temporary */ + array2D H; + std::array rows; + array2D cols; + std::array colcount; - //float_ext_type m_RHSx[storage_N]; + unsigned m_cnt; - const std::size_t m_dim; + //float_ext_type m_RHSx[storage_N]; -}; + const std::size_t m_dim; -// ---------------------------------------------------------------------------------------- -// matrix_solver_direct -// ---------------------------------------------------------------------------------------- + }; -template -void matrix_solver_w_t::LE_invert() -{ - const std::size_t kN = size(); + // ---------------------------------------------------------------------------------------- + // matrix_solver_direct + // ---------------------------------------------------------------------------------------- - for (std::size_t i = 0; i < kN; i++) + template + void matrix_solver_w_t::LE_invert() { - for (std::size_t j = 0; j < kN; j++) - { - W(i,j) = lA(i,j) = A(i,j); - Ainv(i,j) = 0.0; - } - Ainv(i,i) = 1.0; - } - /* down */ - for (std::size_t i = 0; i < kN; i++) - { - /* FIXME: Singular matrix? */ - const float_type f = 1.0 / W(i,i); - const auto * const p = m_terms[i].m_nzrd.data(); - const size_t e = m_terms[i].m_nzrd.size(); + const std::size_t kN = size(); - /* Eliminate column i from row j */ - - const auto * const pb = m_terms[i].m_nzbd.data(); - const size_t eb = m_terms[i].m_nzbd.size(); - for (std::size_t jb = 0; jb < eb; jb++) + for (std::size_t i = 0; i < kN; i++) { - const auto j = pb[jb]; - const float_type f1 = - W(j,i) * f; - if (f1 != 0.0) + for (std::size_t j = 0; j < kN; j++) { - for (std::size_t k = 0; k < e; k++) - W(j,p[k]) += W(i,p[k]) * f1; - for (std::size_t k = 0; k <= i; k ++) - Ainv(j,k) += Ainv(i,k) * f1; + W(i,j) = lA(i,j) = A(i,j); + Ainv(i,j) = 0.0; } + Ainv(i,i) = 1.0; } - } - /* up */ - for (std::size_t i = kN; i-- > 0; ) - { - /* FIXME: Singular matrix? */ - const float_type f = 1.0 / W(i,i); - for (std::size_t j = i; j-- > 0; ) + /* down */ + for (std::size_t i = 0; i < kN; i++) { - const float_type f1 = - W(j,i) * f; - if (f1 != 0.0) + /* FIXME: Singular matrix? */ + const float_type f = 1.0 / W(i,i); + const auto * const p = m_terms[i].m_nzrd.data(); + const size_t e = m_terms[i].m_nzrd.size(); + + /* Eliminate column i from row j */ + + const auto * const pb = m_terms[i].m_nzbd.data(); + const size_t eb = m_terms[i].m_nzbd.size(); + for (std::size_t jb = 0; jb < eb; jb++) { - for (std::size_t k = i; k < kN; k++) - W(j,k) += W(i,k) * f1; - for (std::size_t k = 0; k < kN; k++) - Ainv(j,k) += Ainv(i,k) * f1; + const auto j = pb[jb]; + const float_type f1 = - W(j,i) * f; + if (f1 != 0.0) + { + for (std::size_t k = 0; k < e; k++) + W(j,p[k]) += W(i,p[k]) * f1; + for (std::size_t k = 0; k <= i; k ++) + Ainv(j,k) += Ainv(i,k) * f1; + } } } - for (std::size_t k = 0; k < kN; k++) + /* up */ + for (std::size_t i = kN; i-- > 0; ) { - Ainv(i,k) *= f; + /* FIXME: Singular matrix? */ + const float_type f = 1.0 / W(i,i); + for (std::size_t j = i; j-- > 0; ) + { + const float_type f1 = - W(j,i) * f; + if (f1 != 0.0) + { + for (std::size_t k = i; k < kN; k++) + W(j,k) += W(i,k) * f1; + for (std::size_t k = 0; k < kN; k++) + Ainv(j,k) += Ainv(i,k) * f1; + } + } + for (std::size_t k = 0; k < kN; k++) + { + Ainv(i,k) *= f; + } } } -} - -template -template -void matrix_solver_w_t::LE_compute_x( - T & x) -{ - const std::size_t kN = size(); - - for (std::size_t i=0; i + template + void matrix_solver_w_t::LE_compute_x( + T & x) { - const float_type f = RHS(k); + const std::size_t kN = size(); for (std::size_t i=0; i -unsigned matrix_solver_w_t::solve_non_dynamic(const bool newton_raphson) -{ - const auto iN = size(); + for (std::size_t i=0; i new_V; - // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) - std::array t; // FIXME: convert to member - // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) - std::array w; - if ((m_cnt % 50) == 0) - { - /* complete calculation */ - this->LE_invert(); - this->LE_compute_x(new_V); - } - else + template + unsigned matrix_solver_w_t::solve_non_dynamic(const bool newton_raphson) { - /* Solve Ay = b for y */ - this->LE_compute_x(new_V); - - /* determine changed rows */ + const auto iN = size(); - unsigned rowcount=0; - #define VT(r,c) (A(r,c) - lA(r,c)) + // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) + std::array new_V; + // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) + std::array t; // FIXME: convert to member + // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) + std::array w; - for (unsigned row = 0; row < iN; row ++) + if ((m_cnt % 50) == 0) { - unsigned cc=0; - auto &nz = m_terms[row].m_nz; - for (auto & col : nz) - { - if (A(row,col) != lA(row,col)) - cols[rowcount][cc++] = col; - } - if (cc > 0) - { - colcount[rowcount] = cc; - rows[rowcount++] = row; - } + /* complete calculation */ + this->LE_invert(); + this->LE_compute_x(new_V); } - if (rowcount > 0) + else { - /* construct w = transform(V) * y - * dim: rowcount x iN - * */ - for (unsigned i = 0; i < rowcount; i++) + /* Solve Ay = b for y */ + this->LE_compute_x(new_V); + + /* determine changed rows */ + + unsigned rowcount=0; + #define VT(r,c) (A(r,c) - lA(r,c)) + + for (unsigned row = 0; row < iN; row ++) { - const unsigned r = rows[i]; - double tmp = 0.0; - for (unsigned k = 0; k < iN; k++) - tmp += VT(r,k) * new_V[k]; - w[i] = tmp; + unsigned cc=0; + auto &nz = m_terms[row].m_nz; + for (auto & col : nz) + { + if (A(row,col) != lA(row,col)) + cols[rowcount][cc++] = col; + } + if (cc > 0) + { + colcount[rowcount] = cc; + rows[rowcount++] = row; + } } + if (rowcount > 0) + { + /* construct w = transform(V) * y + * dim: rowcount x iN + * */ + for (unsigned i = 0; i < rowcount; i++) + { + const unsigned r = rows[i]; + double tmp = 0.0; + for (unsigned k = 0; k < iN; k++) + tmp += VT(r,k) * new_V[k]; + w[i] = tmp; + } - for (unsigned i = 0; i < rowcount; i++) - for (unsigned k=0; k< rowcount; k++) - H[i][k] = 0.0; + for (unsigned i = 0; i < rowcount; i++) + for (unsigned k=0; k< rowcount; k++) + H[i][k] = 0.0; - for (unsigned i = 0; i < rowcount; i++) - H[i][i] = 1.0; - /* Construct H = (I + VT*Z) */ - for (unsigned i = 0; i < rowcount; i++) - for (unsigned k=0; k< colcount[i]; k++) + for (unsigned i = 0; i < rowcount; i++) + H[i][i] = 1.0; + /* Construct H = (I + VT*Z) */ + for (unsigned i = 0; i < rowcount; i++) + for (unsigned k=0; k< colcount[i]; k++) + { + const unsigned col = cols[i][k]; + float_type f = VT(rows[i],col); + if (f!=0.0) + for (unsigned j= 0; j < rowcount; j++) + H[i][j] += f * Ainv(col,rows[j]); + } + + /* Gaussian elimination of H */ + for (unsigned i = 0; i < rowcount; i++) { - const unsigned col = cols[i][k]; - float_type f = VT(rows[i],col); - if (f!=0.0) - for (unsigned j= 0; j < rowcount; j++) - H[i][j] += f * Ainv(col,rows[j]); + if (H[i][i] == 0.0) + plib::perrlogger("{} H singular\n", this->name()); + const float_type f = 1.0 / H[i][i]; + for (unsigned j = i+1; j < rowcount; j++) + { + const float_type f1 = - f * H[j][i]; + + if (f1!=0.0) + { + float_type *pj = &H[j][i+1]; + const float_type *pi = &H[i][i+1]; + for (unsigned k = 0; k < rowcount-i-1; k++) + pj[k] += f1 * pi[k]; + //H[j][k] += f1 * H[i][k]; + w[j] += f1 * w[i]; + } + } } - - /* Gaussian elimination of H */ - for (unsigned i = 0; i < rowcount; i++) - { - if (H[i][i] == 0.0) - plib::perrlogger("{} H singular\n", this->name()); - const float_type f = 1.0 / H[i][i]; - for (unsigned j = i+1; j < rowcount; j++) + /* Back substitution */ + //inv(H) w = t w = H t + for (unsigned j = rowcount; j-- > 0; ) { - const float_type f1 = - f * H[j][i]; + float_type tmp = 0; + const float_type *pj = &H[j][j+1]; + const float_type *tj = &t[j+1]; + for (unsigned k = 0; k < rowcount-j-1; k++) + tmp += pj[k] * tj[k]; + //tmp += H[j][k] * t[k]; + t[j] = (w[j] - tmp) / H[j][j]; + } - if (f1!=0.0) + /* x = y - Zt */ + for (unsigned i=0; i 0; ) - { - float_type tmp = 0; - const float_type *pj = &H[j][j+1]; - const float_type *tj = &t[j+1]; - for (unsigned k = 0; k < rowcount-j-1; k++) - tmp += pj[k] * tj[k]; - //tmp += H[j][k] * t[k]; - t[j] = (w[j] - tmp) / H[j][j]; - } + } + m_cnt++; - /* x = y - Zt */ + if (false) for (unsigned i=0; i 1e-6) + plib::perrlogger("{} failed on row {}: {} RHS: {}\n", this->name(), i, std::abs(tmp-RHS(i)), RHS(i)); } - if (std::abs(tmp-RHS(i)) > 1e-6) - plib::perrlogger("{} failed on row {}: {} RHS: {}\n", this->name(), i, std::abs(tmp-RHS(i)), RHS(i)); - } - const float_type err = (newton_raphson ? delta(new_V) : 0.0); - store(new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; -} + const float_type err = (newton_raphson ? delta(new_V) : 0.0); + store(new_V); + return (err > this->m_params.m_accuracy) ? 2 : 1; + } -template -unsigned matrix_solver_w_t::vsolve_non_dynamic(const bool newton_raphson) -{ - this->build_LE_A(*this); - this->build_LE_RHS(*this); + template + unsigned matrix_solver_w_t::vsolve_non_dynamic(const bool newton_raphson) + { + this->build_LE_A(*this); + this->build_LE_RHS(*this); - this->m_stat_calculations++; - return this->solve_non_dynamic(newton_raphson); -} + this->m_stat_calculations++; + return this->solve_non_dynamic(newton_raphson); + } - } //namespace devices +} // namespace solver } // namespace netlist #endif /* NLD_MS_DIRECT_H_ */ diff --git a/src/lib/netlist/solver/nld_solver.cpp b/src/lib/netlist/solver/nld_solver.cpp index b0dcbec2523..e19e216a09f 100644 --- a/src/lib/netlist/solver/nld_solver.cpp +++ b/src/lib/netlist/solver/nld_solver.cpp @@ -81,7 +81,7 @@ namespace devices std::size_t nthreads = std::min(static_cast(m_params.m_parallel()), plib::omp::get_max_threads()); - std::vector &solvers = (force_solve ? m_mat_solvers_all : m_mat_solvers_timestepping); + std::vector &solvers = (force_solve ? m_mat_solvers_all : m_mat_solvers_timestepping); if (nthreads > 1 && solvers.size() > 1) { @@ -110,45 +110,45 @@ namespace devices } template - plib::unique_ptr create_it(netlist_state_t &nl, pstring name, + plib::unique_ptr create_it(netlist_state_t &nl, pstring name, analog_net_t::list_t &nets, - solver_parameters_t ¶ms, std::size_t size) + solver::solver_parameters_t ¶ms, std::size_t size) { return plib::make_unique(nl, name, nets, ¶ms, size); } template - plib::unique_ptr NETLIB_NAME(solver)::create_solver(std::size_t size, + plib::unique_ptr NETLIB_NAME(solver)::create_solver(std::size_t size, const pstring &solvername, analog_net_t::list_t &nets) { switch (m_params.m_method()) { - case matrix_type_e::MAT_CR: + case solver::matrix_type_e::MAT_CR: if (size > 0) // GCR always outperforms MAT solver { - return create_it>(state(), solvername, nets, m_params, size); + return create_it>(state(), solvername, nets, m_params, size); } else { - return create_it>(state(), solvername, nets, m_params, size); + return create_it>(state(), solvername, nets, m_params, size); } - case matrix_type_e::SOR_MAT: - return create_it>(state(), solvername, nets, m_params, size); - case matrix_type_e::MAT: - return create_it>(state(), solvername, nets, m_params, size); - case matrix_type_e::SM: + case solver::matrix_type_e::SOR_MAT: + return create_it>(state(), solvername, nets, m_params, size); + case solver::matrix_type_e::MAT: + return create_it>(state(), solvername, nets, m_params, size); + case solver::matrix_type_e::SM: /* Sherman-Morrison Formula */ - return create_it>(state(), solvername, nets, m_params, size); - case matrix_type_e::W: + return create_it>(state(), solvername, nets, m_params, size); + case solver::matrix_type_e::W: /* Woodbury Formula */ - return create_it>(state(), solvername, nets, m_params, size); - case matrix_type_e::SOR: - return create_it>(state(), solvername, nets, m_params, size); - case matrix_type_e::GMRES: - return create_it>(state(), solvername, nets, m_params, size); + return create_it>(state(), solvername, nets, m_params, size); + case solver::matrix_type_e::SOR: + return create_it>(state(), solvername, nets, m_params, size); + case solver::matrix_type_e::GMRES: + return create_it>(state(), solvername, nets, m_params, size); } - return plib::unique_ptr(); + return plib::unique_ptr(); } struct net_splitter @@ -224,17 +224,17 @@ namespace devices log().verbose("Found {1} net groups in {2} nets\n", splitter.groups.size(), state().nets().size()); for (auto & grp : splitter.groups) { - plib::unique_ptr ms; + plib::unique_ptr ms; std::size_t net_count = grp.size(); pstring sname = plib::pfmt("Solver_{1}")(m_mat_solvers.size()); switch (net_count) { case 1: - ms = plib::make_unique>(state(), sname, grp, &m_params); + ms = plib::make_unique>(state(), sname, grp, &m_params); break; case 2: - ms = plib::make_unique>(state(), sname, grp, &m_params); + ms = plib::make_unique>(state(), sname, grp, &m_params); break; #if 0 case 3: diff --git a/src/lib/netlist/solver/nld_solver.h b/src/lib/netlist/solver/nld_solver.h index 1153b1a6011..d4d506d3683 100644 --- a/src/lib/netlist/solver/nld_solver.h +++ b/src/lib/netlist/solver/nld_solver.h @@ -51,14 +51,14 @@ namespace devices logic_input_t m_fb_step; logic_output_t m_Q_step; - std::vector> m_mat_solvers; - std::vector m_mat_solvers_all; - std::vector m_mat_solvers_timestepping; + std::vector> m_mat_solvers; + std::vector m_mat_solvers_all; + std::vector m_mat_solvers_timestepping; - solver_parameters_t m_params; + solver::solver_parameters_t m_params; template - plib::unique_ptr create_solver(std::size_t size, + plib::unique_ptr create_solver(std::size_t size, const pstring &solvername, analog_net_t::list_t &nets); }; -- cgit v1.2.3