From a7a81862837cf8df3bc19c601b9d06ffb16262c9 Mon Sep 17 00:00:00 2001 From: couriersud Date: Sat, 23 Feb 2019 20:34:03 +0100 Subject: Move fillmatrix to to matrix solver base class. (nw) --- src/lib/netlist/plib/mat_cr.h | 13 +++-- src/lib/netlist/plib/parray.h | 1 + src/lib/netlist/solver/nld_matrix_solver.h | 76 ++++++++++++++++-------------- src/lib/netlist/solver/nld_ms_gcr.h | 6 +-- src/lib/netlist/solver/nld_ms_gmres.h | 3 +- 5 files changed, 56 insertions(+), 43 deletions(-) diff --git a/src/lib/netlist/plib/mat_cr.h b/src/lib/netlist/plib/mat_cr.h index 3cf2570b725..7fc5d46b85d 100644 --- a/src/lib/netlist/plib/mat_cr.h +++ b/src/lib/netlist/plib/mat_cr.h @@ -67,13 +67,20 @@ namespace plib , m_size(n) { for (index_type i=0; i((N>0) ? N : m_size); } + void clear() + { + nz_num = 0; + for (index_type i=0; i < size() + 1; i++) + row_idx[i] = 0; + } + void set_scalar(const T scalar) { for (index_type i=0, e=nz_num; i 0 && col_idx[ri] == c) + if (ri < row_idx[r+1] && col_idx[ri] == c) A[ri] = val; else { @@ -97,7 +104,7 @@ namespace plib } A[ri] = val; col_idx[ri] = c; - for (C i = row_idx[r]; i < size()+1;i++) + for (C i = r + 1; i < size() + 1; i++) row_idx[i]++; nz_num++; if (c==r) diff --git a/src/lib/netlist/plib/parray.h b/src/lib/netlist/plib/parray.h index 5ab52524a07..9b64f5343d2 100644 --- a/src/lib/netlist/plib/parray.h +++ b/src/lib/netlist/plib/parray.h @@ -60,6 +60,7 @@ namespace plib { using size_type = typename base_type::size_type; using reference = typename base_type::reference; using const_reference = typename base_type::const_reference; + using value_type = typename base_type::value_type; template parray(size_type size, typename std::enable_if::type = 0) diff --git a/src/lib/netlist/solver/nld_matrix_solver.h b/src/lib/netlist/solver/nld_matrix_solver.h index cd5a4f3bd0c..903c39257d0 100644 --- a/src/lib/netlist/solver/nld_matrix_solver.h +++ b/src/lib/netlist/solver/nld_matrix_solver.h @@ -53,42 +53,6 @@ namespace devices void set_pointers(); - /* FIXME: this works a bit better for larger matrices */ - template - void fill_matrix/*_larger*/(AP &tcr, FT &RHS) - { - - const std::size_t term_count = this->count(); - const std::size_t railstart = this->m_railstart; - - for (std::size_t i = 0; i < railstart; i++) - *tcr[i] -= m_go[i]; - - #if 1 - FT gtot_t = 0.0; - FT RHS_t = 0.0; - - for (std::size_t i = 0; i < term_count; i++) - { - gtot_t += m_gt[i]; - RHS_t += m_Idr[i]; - } - // FIXME: Code above is faster than vec_sum - Check this - #else - 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]*/ + m_go[i] * *m_connected_net_V[i]); - } - - RHS = RHS_t; - // update diagonal element ... - *tcr[railstart] += gtot_t; //mat.A[mat.diag[k]] += gtot_t; - } - std::size_t m_railstart; std::vector m_nz; /* all non zero for multiplication */ @@ -207,6 +171,46 @@ namespace devices template void build_LE_RHS(T &child); + template + 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].data(); + + const std::size_t term_count = net->count(); + const std::size_t railstart = net->m_railstart; + + for (std::size_t i = 0; i < railstart; i++) + *tcr_r[i] -= net->m_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 += net->m_gt[i]; + RHS_t += net->m_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]*/ + net->m_go[i] * *net->m_connected_net_V[i]); + } + + RHS[k] = RHS_t; + // update diagonal element ... + *tcr_r[railstart] += gtot_t; //mat.A[mat.diag[k]] += gtot_t; + } + + } + std::vector> m_terms; std::vector m_nets; std::vector> m_inps; diff --git a/src/lib/netlist/solver/nld_ms_gcr.h b/src/lib/netlist/solver/nld_ms_gcr.h index 481d82a0e6b..fce7b5bae2c 100644 --- a/src/lib/netlist/solver/nld_ms_gcr.h +++ b/src/lib/netlist/solver/nld_ms_gcr.h @@ -67,7 +67,7 @@ private: plib::parray new_V; std::array, storage_N> m_term_cr; -// std::array, storage_N> m_term_cr; + // std::array, storage_N> m_term_cr; mat_type mat; @@ -296,8 +296,8 @@ unsigned matrix_solver_GCR_t::vsolve_non_dynamic(const bool newton_rap mat.set_scalar(0.0); /* populate matrix */ - for (std::size_t k = 0; k < iN; k++) - this->m_terms[k]->fill_matrix(m_term_cr[k], RHS[k]); + + this->fill_matrix(iN, m_term_cr, RHS); /* now solve it */ diff --git a/src/lib/netlist/solver/nld_ms_gmres.h b/src/lib/netlist/solver/nld_ms_gmres.h index 59c6259aafe..95de1bbae31 100644 --- a/src/lib/netlist/solver/nld_ms_gmres.h +++ b/src/lib/netlist/solver/nld_ms_gmres.h @@ -111,9 +111,10 @@ namespace devices m_ops.m_mat.set_scalar(0.0); /* populate matrix and V for first estimate */ + this->fill_matrix(iN, m_term_cr, RHS); + for (std::size_t k = 0; k < iN; k++) { - this->m_terms[k]->fill_matrix(m_term_cr[k], RHS[k]); this->m_new_V[k] = this->m_nets[k]->Q_Analog(); } -- cgit v1.2.3