From 47d938b149bdbe1f7ea0a875f619a88928e4186e Mon Sep 17 00:00:00 2001 From: couriersud Date: Fri, 1 Nov 2019 01:21:43 +0100 Subject: netlist: further solver refactoring. (nw) --- src/lib/netlist/plib/parray.h | 2 + src/lib/netlist/solver/nld_matrix_solver.cpp | 52 ----- src/lib/netlist/solver/nld_matrix_solver.h | 316 +++++++++++++++++---------- src/lib/netlist/solver/nld_ms_direct.h | 37 ++-- 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 | 34 +-- src/lib/netlist/solver/nld_ms_gmres.h | 2 +- src/lib/netlist/solver/nld_ms_sm.h | 35 ++- src/lib/netlist/solver/nld_ms_sor_mat.h | 5 +- src/lib/netlist/solver/nld_ms_w.h | 39 ++-- 11 files changed, 270 insertions(+), 260 deletions(-) diff --git a/src/lib/netlist/plib/parray.h b/src/lib/netlist/plib/parray.h index c1196e1e7d4..a1e3c2c504a 100644 --- a/src/lib/netlist/plib/parray.h +++ b/src/lib/netlist/plib/parray.h @@ -105,6 +105,8 @@ namespace plib { throw plib::pexception("parray: size error " + plib::to_string(size) + ">" + plib::to_string(SIZE)); } + base_type &as_base() noexcept { return m_a; } + inline size_type size() const noexcept { return SIZE <= 0 ? m_size : SIZEABS(); } constexpr size_type max_size() const noexcept { return base_type::max_size(); } diff --git a/src/lib/netlist/solver/nld_matrix_solver.cpp b/src/lib/netlist/solver/nld_matrix_solver.cpp index 18f0b842f01..afcdec647cd 100644 --- a/src/lib/netlist/solver/nld_matrix_solver.cpp +++ b/src/lib/netlist/solver/nld_matrix_solver.cpp @@ -341,14 +341,6 @@ namespace solver * save states */ - m_last_V.resize(iN, plib::constants::zero()); - m_DD_n_m_1.resize(iN, plib::constants::zero()); - m_h_n_m_1.resize(iN, plib::constants::zero()); - - 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"); - for (std::size_t k = 0; k < iN; k++) { pstring num = plib::pfmt("{1}")(k); @@ -553,50 +545,6 @@ namespace solver } } - netlist_time matrix_solver_t::compute_next_timestep(const nl_fptype cur_ts) - { - nl_fptype new_solver_timestep = m_params.m_max_timestep; - - if (m_params.m_dynamic_ts) - { - for (std::size_t k = 0; k < m_terms.size(); k++) - { - auto &t = m_terms[k]; - //const nl_fptype DD_n = (n->Q_Analog() - t->m_last_V); - // avoid floating point exceptions - - const nl_fptype DD_n = std::max(-fp_constants::TIMESTEP_MAXDIFF, std::min(+fp_constants::TIMESTEP_MAXDIFF,(t.getV() - m_last_V[k]))); - const nl_fptype hn = cur_ts; - - //printf("%g %g %g %g\n", DD_n, hn, t.m_DD_n_m_1, t.m_h_n_m_1); - 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 (std::fabs(DD2) > fp_constants::TIMESTEP_MINDIV) // avoid div-by-zero - new_net_timestep = std::sqrt(m_params.m_dynamic_lte / std::fabs(plib::constants::cast(0.5)*DD2)); - else - new_net_timestep = m_params.m_max_timestep; - - if (new_net_timestep < new_solver_timestep) - new_solver_timestep = new_net_timestep; - - m_last_V[k] = t.getV(); - } - if (new_solver_timestep < m_params.m_min_timestep) - { - new_solver_timestep = m_params.m_min_timestep; - } - } - //if (new_solver_timestep > 10.0 * hn) - // new_solver_timestep = 10.0 * hn; - /* - * FIXME: Factor 2 below is important. Without, we get timing issues. This must be a bug elsewhere. - */ - return std::max(netlist_time::from_double(new_solver_timestep), netlist_time::quantum() * 2); - } - void matrix_solver_t::log_stats() { if (this->m_stat_calculations != 0 && this->m_stat_vsolver_calls && log().verbose.is_enabled()) diff --git a/src/lib/netlist/solver/nld_matrix_solver.h b/src/lib/netlist/solver/nld_matrix_solver.h index fdfaff7a399..d639218681c 100644 --- a/src/lib/netlist/solver/nld_matrix_solver.h +++ b/src/lib/netlist/solver/nld_matrix_solver.h @@ -207,15 +207,42 @@ namespace solver void update_dynamic(); virtual unsigned vsolve_non_dynamic(const bool newton_raphson) = 0; + virtual netlist_time compute_next_timestep(const nl_fptype cur_ts) = 0; - netlist_time compute_next_timestep(const nl_fptype cur_ts); - /* virtual */ void add_term(std::size_t net_idx, terminal_t *term); + void add_term(std::size_t net_idx, terminal_t *term); + + 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 - void store(const T & V); + void store(const T & V) + { + const std::size_t iN = this->m_terms.size(); + for (std::size_t i = 0; i < iN; i++) + this->m_terms[i].setV(V[i]); + } + template - auto delta(const T & V) -> typename std::decay::type; + auto delta(const T & V) -> typename std::decay::type + { + /* 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(); + using vtype = typename std::decay::type; + vtype cerr = 0; + for (std::size_t i = 0; i < iN; i++) + cerr = std::max(cerr, std::abs(V[i] - static_cast(this->m_terms[i].getV()))); + return cerr; + } void set_pointers() { @@ -229,7 +256,6 @@ namespace solver max_rail = std::max(max_rail, m_terms[k].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); @@ -247,49 +273,6 @@ namespace solver } } - template - void fill_matrix(std::size_t N, FT &RHS) - { - for (std::size_t k = 0; k < N; k++) - { - auto &net = m_terms[k]; - auto **tcr_r = &(m_mat_ptr[k][0]); - - 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]; - - 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(term_count, m_gt); - auto RHS_t = plib::vec_sum(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 void log_fill(const T &fill, M &mat) @@ -345,61 +328,17 @@ namespace solver } } - template - void clear_square_mat(std::size_t n, M &m) - { - for (std::size_t k=0; k < n; k++) - { - auto *p = &(m[k][0]); - for (std::size_t i=0; i < n; i++) - p[i] = 0.0; - } - } - - template - void build_mat_ptr(std::size_t iN, M &mat) - { - for (std::size_t k=0; km_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(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 using aligned_alloc = plib::aligned_allocator; plib::pmatrix2d> m_gonn; plib::pmatrix2d> m_gtn; plib::pmatrix2d> m_Idrn; - plib::pmatrix2d> m_mat_ptr; plib::pmatrix2d> m_connected_net_Vn; plib::aligned_vector m_terms; plib::aligned_vector m_rails_temp; - /* state - variable time_stepping */ - plib::aligned_vector m_last_V; - plib::aligned_vector m_DD_n_m_1; - plib::aligned_vector m_h_n_m_1; - - // FIXME: it should be like this, however dimensions are determined - // in vsetup. - //state_container> m_last_V; - //state_container> m_DD_n_m_1; - //state_container> m_h_n_m_1; - std::vector> m_inps; const solver_parameters_t &m_params; @@ -430,29 +369,182 @@ namespace solver std::size_t m_ops; }; - template - auto matrix_solver_t::delta(const T & V) -> typename std::decay::type + template + class matrix_solver_ext_t: public matrix_solver_t { - /* 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. - */ + friend class matrix_solver_t; + public: - const std::size_t iN = this->m_terms.size(); - using vtype = typename std::decay::type; - vtype cerr = 0; - for (std::size_t i = 0; i < iN; i++) - cerr = std::max(cerr, std::abs(V[i] - static_cast(this->m_terms[i].getV()))); - return cerr; - } + using float_type = FT; - template - 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_terms[i].setV(V[i]); - } + 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_mat_ptr(size, this->max_railstart() + 1) + , m_last_V(size, plib::constants::zero()) + , m_DD_n_m_1(size, plib::constants::zero()) + , m_h_n_m_1(size, plib::constants::zero()) + { + /* + * 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::SIZEABS(); + static constexpr const std::size_t m_pitch_ABS = (((SIZEABS + 0) + 7) / 8) * 8; + + plib::parray2D m_mat_ptr; + /* state - variable time_stepping */ + plib::parray m_last_V; + plib::parray m_DD_n_m_1; + plib::parray m_h_n_m_1; + + // FIXME: it should be like this, however dimensions are determined + // in vsetup. + //state_container> m_last_V; + //state_container> m_DD_n_m_1; + //state_container> m_h_n_m_1; + + constexpr std::size_t size() const noexcept { return (SIZE > 0) ? static_cast(SIZE) : m_dim; } + + netlist_time compute_next_timestep(const nl_fptype cur_ts) override + { + nl_fptype new_solver_timestep = m_params.m_max_timestep; + + if (m_params.m_dynamic_ts) + { + for (std::size_t k = 0; k < m_terms.size(); k++) + { + auto &t = m_terms[k]; + //const nl_fptype DD_n = (n->Q_Analog() - t->m_last_V); + // avoid floating point exceptions + + const nl_fptype DD_n = std::max(-fp_constants::TIMESTEP_MAXDIFF, std::min(+fp_constants::TIMESTEP_MAXDIFF,(t.getV() - m_last_V[k]))); + const nl_fptype hn = cur_ts; + + //printf("%g %g %g %g\n", DD_n, hn, t.m_DD_n_m_1, t.m_h_n_m_1); + 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 (std::fabs(DD2) > fp_constants::TIMESTEP_MINDIV) // avoid div-by-zero + new_net_timestep = std::sqrt(m_params.m_dynamic_lte / std::fabs(plib::constants::cast(0.5)*DD2)); + else + new_net_timestep = m_params.m_max_timestep; + + if (new_net_timestep < new_solver_timestep) + new_solver_timestep = new_net_timestep; + + m_last_V[k] = t.getV(); + } + if (new_solver_timestep < m_params.m_min_timestep) + { + new_solver_timestep = m_params.m_min_timestep; + } + } + //if (new_solver_timestep > 10.0 * hn) + // new_solver_timestep = 10.0 * hn; + /* + * FIXME: Factor 2 below is important. Without, we get timing issues. This must be a bug elsewhere. + */ + return std::max(netlist_time::from_double(new_solver_timestep), netlist_time::quantum() * 2); + } + + + template + void build_mat_ptr(M &mat) + { + const std::size_t iN = size(); + + for (std::size_t k=0; km_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(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 + 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]); + for (std::size_t i=0; i < n; i++) + p[i] = 0.0; + } + } + + template + void fill_matrix(RT &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]); + + 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]; + + for (std::size_t i = 0; i < railstart; i++) + *tcr_r[i] += go[i]; + + typename RT::value_type gtot_t = 0.0; + typename RT::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(term_count, m_gt); + auto RHS_t = plib::vec_sum(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; + } + + } + + }; } // namespace solver } // namespace netlist diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h index f3c31dc3371..62332dbe6b5 100644 --- a/src/lib/netlist/solver/nld_ms_direct.h +++ b/src/lib/netlist/solver/nld_ms_direct.h @@ -22,7 +22,7 @@ namespace solver { template - class matrix_solver_direct_t: public matrix_solver_t + class matrix_solver_direct_t: public matrix_solver_ext_t { friend class matrix_solver_t; public: @@ -37,7 +37,6 @@ namespace solver private: - const std::size_t m_dim; const std::size_t m_pitch; protected: @@ -47,8 +46,6 @@ namespace solver unsigned vsolve_non_dynamic(const bool newton_raphson) override; unsigned solve_non_dynamic(const bool newton_raphson); - constexpr std::size_t size() const noexcept { return (SIZE > 0) ? static_cast(SIZE) : m_dim; } - void LE_solve(); template @@ -69,15 +66,15 @@ namespace solver template void matrix_solver_direct_t::LE_solve() { - const std::size_t kN = size(); - if (!m_params.m_pivot) + const std::size_t kN = this->size(); + if (!this->m_params.m_pivot) { for (std::size_t i = 0; i < kN; i++) { /* FIXME: Singular matrix? */ const FT f = 1.0 / m_A[i][i]; - const auto &nzrd = m_terms[i].m_nzrd; - const auto &nzbd = m_terms[i].m_nzbd; + const auto &nzrd = this->m_terms[i].m_nzrd; + const auto &nzbd = this->m_terms[i].m_nzbd; for (auto &j : nzbd) { @@ -138,10 +135,10 @@ namespace solver void matrix_solver_direct_t::LE_back_subst( T & x) { - const std::size_t kN = size(); + const std::size_t kN = this->size(); /* back substitution */ - if (m_params.m_pivot) + if (this->m_params.m_pivot) { for (std::size_t j = kN; j-- > 0; ) { @@ -156,7 +153,7 @@ namespace solver for (std::size_t j = kN; j-- > 0; ) { FT tmp = 0; - const auto &nzrd = m_terms[j].m_nzrd; + const auto &nzrd = this->m_terms[j].m_nzrd; const auto e = nzrd.size(); // - 1; /* exclude RHS element */ for ( std::size_t k = 0; k < e; k++) tmp += m_A[j][nzrd[k]] * x[nzrd[k]]; @@ -171,20 +168,17 @@ namespace solver this->LE_solve(); this->LE_back_subst(m_new_V); - const FT err = (newton_raphson ? delta(m_new_V) : 0.0); - store(m_new_V); + const FT err = (newton_raphson ? this->delta(m_new_V) : 0.0); + this->store(m_new_V); return (err > this->m_params.m_accuracy) ? 2 : 1; } template unsigned matrix_solver_direct_t::vsolve_non_dynamic(const bool newton_raphson) { - const std::size_t iN = this->size(); - /* populate matrix */ - - this->clear_square_mat(iN, m_A); - this->fill_matrix(iN, m_RHS); + this->clear_square_mat(m_A); + this->fill_matrix(m_RHS); this->m_stat_calculations++; return this->solve_non_dynamic(newton_raphson); @@ -195,14 +189,13 @@ namespace solver 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_pitch(m_pitch_ABS ? m_pitch_ABS : (((m_dim + 0) + 7) / 8) * 8) + : matrix_solver_ext_t(anetlist, name, nets, params, size) + , m_pitch(m_pitch_ABS ? m_pitch_ABS : (((size + 0) + 7) / 8) * 8) , m_new_V(size) , m_A(size, m_pitch) , m_RHS(size) { - this->build_mat_ptr(size, m_A); + this->build_mat_ptr(m_A); } } // namespace solver diff --git a/src/lib/netlist/solver/nld_ms_direct1.h b/src/lib/netlist/solver/nld_ms_direct1.h index 1eb992104a8..2931aaec29c 100644 --- a/src/lib/netlist/solver/nld_ms_direct1.h +++ b/src/lib/netlist/solver/nld_ms_direct1.h @@ -34,8 +34,8 @@ namespace solver // ---------------------------------------------------------------------------------------- unsigned vsolve_non_dynamic(const bool newton_raphson) override { - this->clear_square_mat(1, this->m_A); - this->fill_matrix(1, this->m_RHS); + this->clear_square_mat(this->m_A); + this->fill_matrix(this->m_RHS); std::array new_V = { this->m_RHS[0] / this->m_A[0][0] }; diff --git a/src/lib/netlist/solver/nld_ms_direct2.h b/src/lib/netlist/solver/nld_ms_direct2.h index 4479ce8003e..03050083158 100644 --- a/src/lib/netlist/solver/nld_ms_direct2.h +++ b/src/lib/netlist/solver/nld_ms_direct2.h @@ -34,8 +34,8 @@ namespace solver {} unsigned vsolve_non_dynamic(const bool newton_raphson) override { - this->clear_square_mat(2, this->m_A); - this->fill_matrix(2, this->m_RHS); + this->clear_square_mat(this->m_A); + this->fill_matrix(this->m_RHS); const float_type a = this->m_A[0][0]; const float_type b = this->m_A[0][1]; diff --git a/src/lib/netlist/solver/nld_ms_gcr.h b/src/lib/netlist/solver/nld_ms_gcr.h index 6c942f3b1d6..66b56397a3d 100644 --- a/src/lib/netlist/solver/nld_ms_gcr.h +++ b/src/lib/netlist/solver/nld_ms_gcr.h @@ -26,25 +26,22 @@ namespace solver { template - class matrix_solver_GCR_t: public matrix_solver_t + class matrix_solver_GCR_t: public matrix_solver_ext_t { public: using mat_type = plib::pGEmatrix_cr_t>; - // FIXME: dirty hack to make this compile - static constexpr const std::size_t storage_N = 100; matrix_solver_GCR_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) + : matrix_solver_ext_t(anetlist, name, nets, params, size) , RHS(size) , new_V(size) , mat(static_cast(size)) , m_proc() { - const std::size_t iN = this->N(); + const std::size_t iN = this->size(); /* build the final matrix */ @@ -65,7 +62,7 @@ namespace solver auto gr = mat.gaussian_extend_fill_mat(fill); - log_fill(fill, mat); + this->log_fill(fill, mat); mat.build_from_fill_mat(fill); @@ -79,13 +76,13 @@ namespace solver for (auto i = mat.row_idx[k]; i < mat.row_idx[k+1]; i++) if (other == static_cast(mat.col_idx[i])) { - m_mat_ptr[k][j] = &mat.A[i]; + this->m_mat_ptr[k][j] = &mat.A[i]; cnt++; break; } } nl_assert(cnt == this->m_terms[k].railstart()); - m_mat_ptr[k][this->m_terms[k].railstart()] = &mat.A[mat.diag[k]]; + this->m_mat_ptr[k][this->m_terms[k].railstart()] = &mat.A[mat.diag[k]]; } this->log().verbose("maximum fill: {1}", gr.first); @@ -96,7 +93,7 @@ namespace solver // FIXME: Move me - if (state().lib().isLoaded()) + if (this->state().lib().isLoaded()) { pstring symname = static_compile_name(); m_proc.load(this->state().lib(), symname); @@ -111,8 +108,6 @@ namespace solver } } - constexpr std::size_t N() const { return m_dim; } - unsigned vsolve_non_dynamic(const bool newton_raphson) override; std::pair create_solver_code() override; @@ -125,7 +120,6 @@ namespace solver pstring static_compile_name(); - const std::size_t m_dim; plib::parray RHS; plib::parray new_V; @@ -142,7 +136,7 @@ namespace solver template void matrix_solver_GCR_t::generate_code(plib::putf8_fmt_writer &strm) { - const std::size_t iN = N(); + const std::size_t iN = this->size(); for (std::size_t i = 0; i < mat.nz_num; i++) strm("double m_A{1} = m_A[{2}];\n", i, i); @@ -232,13 +226,9 @@ namespace solver template unsigned matrix_solver_GCR_t::vsolve_non_dynamic(const bool newton_raphson) { - const std::size_t iN = this->N(); - - mat.set_scalar(0.0); - /* populate matrix */ - - this->fill_matrix(iN, RHS); + mat.set_scalar(0.0); + this->fill_matrix(RHS); /* now solve it */ @@ -257,8 +247,8 @@ namespace solver this->m_stat_calculations++; - const FT err = (newton_raphson ? delta(new_V) : 0.0); - store(new_V); + const FT err = (newton_raphson ? this->delta(new_V) : 0.0); + this->store(new_V); return (err > this->m_params.m_accuracy) ? 2 : 1; } diff --git a/src/lib/netlist/solver/nld_ms_gmres.h b/src/lib/netlist/solver/nld_ms_gmres.h index a4926e2fc1d..9a3c4e2f1bc 100644 --- a/src/lib/netlist/solver/nld_ms_gmres.h +++ b/src/lib/netlist/solver/nld_ms_gmres.h @@ -106,7 +106,7 @@ namespace solver m_ops.m_mat.set_scalar(0.0); /* populate matrix and V for first estimate */ - this->fill_matrix(iN, RHS); + this->fill_matrix(RHS); for (std::size_t k = 0; k < iN; k++) { diff --git a/src/lib/netlist/solver/nld_ms_sm.h b/src/lib/netlist/solver/nld_ms_sm.h index 83e309407f7..4828dbfe72a 100644 --- a/src/lib/netlist/solver/nld_ms_sm.h +++ b/src/lib/netlist/solver/nld_ms_sm.h @@ -45,7 +45,7 @@ namespace solver { template - class matrix_solver_sm_t: public matrix_solver_t + class matrix_solver_sm_t: public matrix_solver_ext_t { friend class matrix_solver_t; @@ -59,11 +59,10 @@ namespace solver matrix_solver_sm_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) + : matrix_solver_ext_t(anetlist, name, nets, params, size) , m_cnt(0) { - this->build_mat_ptr(this->size(), m_A); + this->build_mat_ptr(m_A); } void reset() override { matrix_solver_t::reset(); } @@ -72,8 +71,6 @@ namespace solver unsigned vsolve_non_dynamic(const bool newton_raphson) override; unsigned solve_non_dynamic(const bool newton_raphson); - constexpr std::size_t size() const { return m_dim; } - void LE_invert(); template @@ -109,7 +106,6 @@ namespace solver //float_ext_type m_RHSx[storage_N]; - const std::size_t m_dim; std::size_t m_cnt; }; @@ -121,7 +117,7 @@ namespace solver template void matrix_solver_sm_t::LE_invert() { - const std::size_t kN = size(); + const std::size_t kN = this->size(); for (std::size_t i = 0; i < kN; i++) { @@ -137,13 +133,13 @@ namespace solver { /* FIXME: Singular matrix? */ const float_type f = 1.0 / W(i,i); - const auto * const p = m_terms[i].m_nzrd.data(); - const std::size_t e = m_terms[i].m_nzrd.size(); + const auto * const p = this->m_terms[i].m_nzrd.data(); + const std::size_t e = this->m_terms[i].m_nzrd.size(); /* Eliminate column i from row j */ - const auto * const pb = m_terms[i].m_nzbd.data(); - const std::size_t eb = m_terms[i].m_nzbd.size(); + const auto * const pb = this->m_terms[i].m_nzbd.data(); + const std::size_t eb = this->m_terms[i].m_nzbd.size(); for (std::size_t jb = 0; jb < eb; jb++) { const unsigned j = pb[jb]; @@ -186,7 +182,7 @@ namespace solver void matrix_solver_sm_t::LE_compute_x( T & x) { - const std::size_t kN = size(); + const std::size_t kN = this->size(); for (std::size_t i=0; i::solve_non_dynamic(const bool newton_raphson) { static constexpr const bool incremental = true; - const std::size_t iN = size(); + const std::size_t iN = this->size(); // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) std::array new_V; @@ -232,7 +228,7 @@ namespace solver { std::size_t colcount = 0; - auto &nz = m_terms[row].m_nz; + auto &nz = this->m_terms[row].m_nz; for (unsigned & col : nz) { v[col] = A(row,col) - lA(row,col); @@ -279,8 +275,8 @@ namespace solver this->LE_compute_x(new_V); - const float_type err = (newton_raphson ? delta(new_V) : 0.0); - store(new_V); + const float_type err = (newton_raphson ? this->delta(new_V) : 0.0); + this->store(new_V); return (err > this->m_params.m_accuracy) ? 2 : 1; } @@ -288,9 +284,8 @@ namespace solver unsigned matrix_solver_sm_t::vsolve_non_dynamic(const bool newton_raphson) { - const std::size_t iN = this->size(); - this->clear_square_mat(iN, this->m_A); - this->fill_matrix(iN, this->m_RHS); + this->clear_square_mat(this->m_A); + this->fill_matrix(this->m_RHS); this->m_stat_calculations++; return this->solve_non_dynamic(newton_raphson); diff --git a/src/lib/netlist/solver/nld_ms_sor_mat.h b/src/lib/netlist/solver/nld_ms_sor_mat.h index f73e88a4d02..8a1575ebfc7 100644 --- a/src/lib/netlist/solver/nld_ms_sor_mat.h +++ b/src/lib/netlist/solver/nld_ms_sor_mat.h @@ -114,11 +114,10 @@ namespace solver * Need something like that for gaussian elimination as well. */ - const std::size_t iN = this->size(); - this->clear_square_mat(iN, this->m_A); - this->fill_matrix(iN, this->m_RHS); + this->clear_square_mat(this->m_A); + this->fill_matrix(this->m_RHS); bool resched = false; diff --git a/src/lib/netlist/solver/nld_ms_w.h b/src/lib/netlist/solver/nld_ms_w.h index 2fd62250b3c..3470cd86f22 100644 --- a/src/lib/netlist/solver/nld_ms_w.h +++ b/src/lib/netlist/solver/nld_ms_w.h @@ -52,7 +52,7 @@ namespace solver { template - class matrix_solver_w_t: public matrix_solver_t + class matrix_solver_w_t: public matrix_solver_ext_t { friend class matrix_solver_t; @@ -66,11 +66,10 @@ namespace solver 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) + : matrix_solver_ext_t(anetlist, name, nets, params, size) , m_cnt(0) - , m_dim(size) { - this->build_mat_ptr(this->size(), m_A); + this->build_mat_ptr(m_A); } void reset() override { matrix_solver_t::reset(); } @@ -79,8 +78,6 @@ namespace solver unsigned vsolve_non_dynamic(const bool newton_raphson) override; unsigned solve_non_dynamic(const bool newton_raphson); - constexpr std::size_t size() const { return m_dim; } - void LE_invert(); template @@ -121,11 +118,6 @@ namespace solver std::array colcount; unsigned m_cnt; - - //float_ext_type m_RHSx[storage_N]; - - const std::size_t m_dim; - }; // ---------------------------------------------------------------------------------------- @@ -135,7 +127,7 @@ namespace solver template void matrix_solver_w_t::LE_invert() { - const std::size_t kN = size(); + const std::size_t kN = this->size(); for (std::size_t i = 0; i < kN; i++) { @@ -151,13 +143,13 @@ namespace solver { /* 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 auto * const p = this->m_terms[i].m_nzrd.data(); + const size_t e = this->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(); + const auto * const pb = this->m_terms[i].m_nzbd.data(); + const size_t eb = this->m_terms[i].m_nzbd.size(); for (std::size_t jb = 0; jb < eb; jb++) { const auto j = pb[jb]; @@ -199,7 +191,7 @@ namespace solver void matrix_solver_w_t::LE_compute_x( T & x) { - const std::size_t kN = size(); + const std::size_t kN = this->size(); for (std::size_t i=0; i unsigned matrix_solver_w_t::solve_non_dynamic(const bool newton_raphson) { - const auto iN = size(); + const auto iN = this->size(); // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) std::array new_V; @@ -245,7 +237,7 @@ namespace solver for (unsigned row = 0; row < iN; row ++) { unsigned cc=0; - auto &nz = m_terms[row].m_nz; + auto &nz = this->m_terms[row].m_nz; for (auto & col : nz) { if (A(row,col) != lA(row,col)) @@ -349,17 +341,16 @@ namespace solver 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); + const float_type err = (newton_raphson ? this->delta(new_V) : 0.0); + this->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) { - const std::size_t iN = this->size(); - this->clear_square_mat(iN, this->m_A); - this->fill_matrix(iN, this->m_RHS); + this->clear_square_mat(this->m_A); + this->fill_matrix(this->m_RHS); this->m_stat_calculations++; return this->solve_non_dynamic(newton_raphson); -- cgit v1.2.3