diff options
Diffstat (limited to 'src/lib/netlist/solver/nld_ms_gcr.h')
-rw-r--r-- | src/lib/netlist/solver/nld_ms_gcr.h | 562 |
1 files changed, 349 insertions, 213 deletions
diff --git a/src/lib/netlist/solver/nld_ms_gcr.h b/src/lib/netlist/solver/nld_ms_gcr.h index f3e56a342d6..a17ca1ec282 100644 --- a/src/lib/netlist/solver/nld_ms_gcr.h +++ b/src/lib/netlist/solver/nld_ms_gcr.h @@ -10,201 +10,390 @@ #ifndef NLD_MS_GCR_H_ #define NLD_MS_GCR_H_ -#include "plib/mat_cr.h" +#include <algorithm> +#include "../plib/pdynlib.h" +#include "mat_cr.h" #include "nld_ms_direct.h" #include "nld_solver.h" -#include "plib/pdynlib.h" -#include "plib/pstream.h" -#include "plib/vector_ops.h" - -#include <algorithm> +#include "vector_base.h" +#include "../plib/pstream.h" namespace netlist { -namespace devices + namespace devices + { +template <std::size_t m_N, std::size_t storage_N> +class matrix_solver_GCR_t: public matrix_solver_t { +public: + + matrix_solver_GCR_t(netlist_t &anetlist, const pstring &name, + const solver_parameters_t *params, const std::size_t size) + : matrix_solver_t(anetlist, name, matrix_solver_t::ASCENDING, params) + , m_dim(size) + , mat(size) + , m_proc() + { + } - template <typename FT, int SIZE> - class matrix_solver_GCR_t: public matrix_solver_t + virtual ~matrix_solver_GCR_t() override { - public: - - using mat_type = plib::matrix_compressed_rows_t<FT, SIZE>; - // 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 solver_parameters_t *params, const std::size_t size) - : matrix_solver_t(anetlist, name, matrix_solver_t::PREFER_IDENTITY_TOP_LEFT, params) - , m_dim(size) - , RHS(size) - , new_V(size) - , mat(static_cast<typename mat_type::index_type>(size)) - , m_proc() - { - } + } + + constexpr std::size_t N() const { return (m_N == 0) ? m_dim : m_N; } - constexpr std::size_t N() const { return m_dim; } + virtual void vsetup(analog_net_t::list_t &nets) override; + virtual unsigned vsolve_non_dynamic(const bool newton_raphson) override; - void vsetup(analog_net_t::list_t &nets) override; - unsigned vsolve_non_dynamic(const bool newton_raphson) override; + virtual std::pair<pstring, pstring> create_solver_code() override; - std::pair<pstring, pstring> create_solver_code() override; +private: - private: + //typedef typename mat_cr_t<storage_N>::type mattype; + typedef typename mat_cr_t<storage_N>::index_type mattype; - using mat_index_type = typename plib::matrix_compressed_rows_t<FT, SIZE>::index_type; + void csc_private(plib::putf8_fmt_writer &strm); - void csc_private(plib::putf8_fmt_writer &strm); + using extsolver = void (*)(double * RESTRICT m_A, double * RESTRICT RHS, double * RESTRICT V); - using extsolver = void (*)(double * m_A, double * RHS, double * V); + pstring static_compile_name(); - pstring static_compile_name(); + const std::size_t m_dim; + std::vector<unsigned> m_term_cr[storage_N]; + mat_cr_t<storage_N> mat; - const std::size_t m_dim; - plib::parray<FT, SIZE> RHS; - plib::parray<FT, SIZE> new_V; + //extsolver m_proc; + plib::dynproc<void, double * RESTRICT, double * RESTRICT, double * RESTRICT> m_proc; - mat_type mat; +}; - //extsolver m_proc; - plib::dynproc<void, double * , double * , double * > m_proc; +// ---------------------------------------------------------------------------------------- +// matrix_solver - GCR +// ---------------------------------------------------------------------------------------- + +template <std::size_t m_N, std::size_t storage_N> +void matrix_solver_GCR_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets) +{ + setup_base(nets); - }; + mattype nz = 0; + const std::size_t iN = this->N(); - // ---------------------------------------------------------------------------------------- - // matrix_solver - GCR - // ---------------------------------------------------------------------------------------- + /* build the final matrix */ - // FIXME: namespace or static class member - template <typename V> - std::size_t inline get_level(const V &v, std::size_t k) + bool touched[storage_N][storage_N] = { { false } }; + for (std::size_t k = 0; k < iN; k++) { - for (std::size_t i = 0; i < v.size(); i++) - if (plib::container::contains(v[i], k)) - return i; - throw plib::pexception("Error in get_level"); + for (auto &j : this->m_terms[k]->m_nz) + touched[k][j] = true; } - template <typename FT, int SIZE> - void matrix_solver_GCR_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets) - { - setup_base(nets); + unsigned fc = 0; - const std::size_t iN = this->N(); + unsigned ops = 0; - /* build the final matrix */ + for (std::size_t k = 0; k < iN; k++) + { + ops++; // 1/A(k,k) + for (std::size_t row = k + 1; row < iN; row++) + { + if (touched[row][k]) + { + ops++; + fc++; + for (std::size_t col = k + 1; col < iN; col++) + if (touched[k][col]) + { + touched[row][col] = true; + ops += 2; + } + } + } + } - std::vector<std::vector<unsigned>> fill(iN); - std::size_t raw_elements = 0; + for (mattype k=0; k<iN; k++) + { + mat.ia[k] = nz; - for (std::size_t k = 0; k < iN; k++) + for (mattype j=0; j<iN; j++) { - fill[k].resize(iN, decltype(mat)::FILL_INFINITY); - for (auto &j : this->m_terms[k]->m_nz) + if (touched[k][j]) { - fill[k][j] = 0; - raw_elements++; + mat.ja[nz] = j; + if (j == k) + mat.diag[k] = nz; + nz++; } + } + m_term_cr[k].clear(); + /* build pointers into the compressed row format matrix for each terminal */ + for (std::size_t j=0; j< this->m_terms[k]->m_railstart;j++) + { + int other = this->m_terms[k]->connected_net_idx()[j]; + for (auto i = mat.ia[k]; i < nz; i++) + if (other == static_cast<int>(mat.ja[i])) + { + m_term_cr[k].push_back(i); + break; + } } + nl_assert(m_term_cr[k].size() == this->m_terms[k]->m_railstart); + } - auto gr = mat.gaussian_extend_fill_mat(fill); + mat.ia[iN] = nz; + mat.nz_num = nz; - /* FIXME: move this to the cr matrix class and use computed - * parallel ordering once it makes sense. - */ + this->log().verbose("Ops: {1} Occupancy ratio: {2}\n", ops, + static_cast<double>(nz) / static_cast<double>(iN * iN)); + + // FIXME: Move me + + if (netlist().lib().isLoaded()) + { + pstring symname = static_compile_name(); +#if 0 + m_proc = this->netlist().lib().template getsym<extsolver>(symname); + if (m_proc != nullptr) + this->log().verbose("External static solver {1} found ...", symname); + else + this->log().warning("External static solver {1} not found ...", symname); +#else + m_proc.load(this->netlist().lib(), symname); + if (m_proc.resolved()) + this->log().warning("External static solver {1} found ...", symname); + else + this->log().warning("External static solver {1} not found ...", symname); +#endif + } - std::vector<unsigned> levL(iN, 0); - std::vector<unsigned> levU(iN, 0); +} +#if 0 +template <std::size_t m_N, std::size_t storage_N> +void matrix_solver_GCR_t<m_N, storage_N>::csc_private(plib::putf8_fmt_writer &strm) +{ + const std::size_t iN = N(); + for (std::size_t i = 0; i < iN - 1; i++) + { + const auto &nzbd = this->m_terms[i]->m_nzbd; - // parallel scheme for L x = y - for (std::size_t k = 0; k < iN; k++) + if (nzbd.size() > 0) { - unsigned lm=0; - for (std::size_t j = 0; j<k; j++) - if (fill[k][j] < decltype(mat)::FILL_INFINITY) - lm = std::max(lm, levL[j]); - levL[k] = 1+lm; + std::size_t pi = mat.diag[i]; + + //const nl_double f = 1.0 / m_A[pi++]; + strm("const double f{1} = 1.0 / m_A[{2}];\n", i, pi); + pi++; + const std::size_t piie = mat.ia[i+1]; + + //for (auto & j : nzbd) + for (std::size_t j : nzbd) + { + // proceed to column i + std::size_t pj = mat.ia[j]; + + while (mat.ja[pj] < i) + pj++; + + //const nl_double f1 = - m_A[pj++] * f; + strm("\tconst double f{1}_{2} = -f{3} * m_A[{4}];\n", i, j, i, pj); + pj++; + + // subtract row i from j */ + for (std::size_t pii = pi; pii<piie; ) + { + while (mat.ja[pj] < mat.ja[pii]) + pj++; + //m_A[pj++] += m_A[pii++] * f1; + strm("\tm_A[{1}] += m_A[{2}] * f{3}_{4};\n", pj, pii, i, j); + pj++; pii++; + } + //RHS[j] += f1 * RHS[i]; + strm("\tRHS[{1}] += f{2}_{3} * RHS[{4}];\n", j, i, j, i); + } } + } - // parallel scheme for U x = y - for (std::size_t k = iN; k-- > 0; ) + //new_V[iN - 1] = RHS[iN - 1] / mat.A[mat.diag[iN - 1]]; + strm("\tV[{1}] = RHS[{2}] / m_A[{3}];\n", iN - 1, iN - 1, mat.diag[iN - 1]); + for (std::size_t j = iN - 1; j-- > 0;) + { + strm("\tdouble tmp{1} = 0.0;\n", j); + const std::size_t e = mat.ia[j+1]; + for (std::size_t pk = mat.diag[j] + 1; pk < e; pk++) { - unsigned lm=0; - for (std::size_t j = iN; --j > k; ) - if (fill[k][j] < decltype(mat)::FILL_INFINITY) - lm = std::max(lm, levU[j]); - levU[k] = 1+lm; + strm("\ttmp{1} += m_A[{2}] * V[{3}];\n", j, pk, mat.ja[pk]); } + strm("\tV[{1}] = (RHS[{1}] - tmp{1}) / m_A[{4}];\n", j, j, j, mat.diag[j]); + } +} +#else +template <std::size_t m_N, std::size_t storage_N> +void matrix_solver_GCR_t<m_N, storage_N>::csc_private(plib::putf8_fmt_writer &strm) +{ + const std::size_t iN = N(); + for (std::size_t i = 0; i < mat.nz_num; i++) + strm("double m_A{1} = m_A[{2}];\n", i, i); - for (std::size_t k = 0; k < iN; k++) + for (std::size_t i = 0; i < iN - 1; i++) + { + const auto &nzbd = this->m_terms[i]->m_nzbd; + + if (nzbd.size() > 0) { - unsigned fm = 0; - pstring ml = ""; - for (std::size_t j = 0; j < iN; j++) + std::size_t pi = mat.diag[i]; + + //const nl_double f = 1.0 / m_A[pi++]; + strm("const double f{1} = 1.0 / m_A{2};\n", i, pi); + pi++; + const std::size_t piie = mat.ia[i+1]; + + //for (auto & j : nzbd) + for (std::size_t j : nzbd) { - ml += fill[k][j] == 0 ? "X" : fill[k][j] < decltype(mat)::FILL_INFINITY ? "+" : "."; - if (fill[k][j] < decltype(mat)::FILL_INFINITY) - if (fill[k][j] > fm) - fm = fill[k][j]; + // proceed to column i + std::size_t pj = mat.ia[j]; + + while (mat.ja[pj] < i) + pj++; + + //const nl_double f1 = - m_A[pj++] * f; + strm("\tconst double f{1}_{2} = -f{3} * m_A{4};\n", i, j, i, pj); + pj++; + + // subtract row i from j */ + for (std::size_t pii = pi; pii<piie; ) + { + while (mat.ja[pj] < mat.ja[pii]) + pj++; + //m_A[pj++] += m_A[pii++] * f1; + strm("\tm_A{1} += m_A{2} * f{3}_{4};\n", pj, pii, i, j); + pj++; pii++; + } + //RHS[j] += f1 * RHS[i]; + strm("\tRHS[{1}] += f{2}_{3} * RHS[{4}];\n", j, i, j, i); } - this->log().verbose("{1:4} {2} {3:4} {4:4} {5:4} {6:4}", k, ml, levL[k], levU[k], get_level(mat.m_ge_par, k), fm); } + } + + //new_V[iN - 1] = RHS[iN - 1] / mat.A[mat.diag[iN - 1]]; + strm("\tV[{1}] = RHS[{2}] / m_A{3};\n", iN - 1, iN - 1, mat.diag[iN - 1]); + for (std::size_t j = iN - 1; j-- > 0;) + { + strm("\tdouble tmp{1} = 0.0;\n", j); + const std::size_t e = mat.ia[j+1]; + for (std::size_t pk = mat.diag[j] + 1; pk < e; pk++) + { + strm("\ttmp{1} += m_A{2} * V[{3}];\n", j, pk, mat.ja[pk]); + } + strm("\tV[{1}] = (RHS[{1}] - tmp{1}) / m_A{4};\n", j, j, j, mat.diag[j]); + } +} +#endif + +template <std::size_t m_N, std::size_t storage_N> +pstring matrix_solver_GCR_t<m_N, storage_N>::static_compile_name() +{ + plib::postringstream t; + plib::putf8_fmt_writer w(t); + csc_private(w); + std::hash<pstring> h; + + return plib::pfmt("nl_gcr_{1:x}_{2}")(h( t.str() ))(mat.nz_num); +} + +template <std::size_t m_N, std::size_t storage_N> +std::pair<pstring, pstring> matrix_solver_GCR_t<m_N, storage_N>::create_solver_code() +{ + plib::postringstream t; + plib::putf8_fmt_writer strm(t); + pstring name = static_compile_name(); + + strm.writeline(plib::pfmt("extern \"C\" void {1}(double * __restrict m_A, double * __restrict RHS, double * __restrict V)\n")(name)); + strm.writeline("{\n"); + csc_private(strm); + strm.writeline("}\n"); + return std::pair<pstring, pstring>(name, t.str()); +} - mat.build_from_fill_mat(fill); +template <std::size_t m_N, std::size_t storage_N> +unsigned matrix_solver_GCR_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson) +{ + const std::size_t iN = this->N(); + + nl_double RHS[storage_N]; + nl_double new_V[storage_N]; - for (mat_index_type k=0; k<iN; k++) + mat.set_scalar(0.0); + + for (std::size_t k = 0; k < iN; k++) + { + terms_for_net_t *t = this->m_terms[k].get(); + nl_double gtot_t = 0.0; + nl_double RHS_t = 0.0; + + const std::size_t term_count = t->count(); + const std::size_t railstart = t->m_railstart; + const nl_double * const RESTRICT gt = t->gt(); + const nl_double * const RESTRICT go = t->go(); + const nl_double * const RESTRICT Idr = t->Idr(); + const nl_double * const * RESTRICT other_cur_analog = t->connected_net_V(); + const unsigned * const RESTRICT tcr = m_term_cr[k].data(); + +#if 0 + for (std::size_t i = 0; i < term_count; i++) { - 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]->m_railstart;j++) - { - int other = this->m_terms[k]->m_connected_net_idx[j]; - for (auto i = mat.row_idx[k]; i < mat.row_idx[k+1]; i++) - if (other == static_cast<int>(mat.col_idx[i])) - { - m_mat_ptr[k][j] = &mat.A[i]; - cnt++; - break; - } - } - nl_assert(cnt == this->m_terms[k]->m_railstart); - m_mat_ptr[k][this->m_terms[k]->m_railstart] = &mat.A[mat.diag[k]]; + gtot_t += gt[i]; + RHS_t += Idr[i]; } - this->log().verbose("maximum fill: {1}", gr.first); - this->log().verbose("Post elimination occupancy ratio: {2} Ops: {1}", gr.second, - static_cast<double>(mat.nz_num) / static_cast<double>(iN * iN)); - this->log().verbose(" Pre elimination occupancy ratio: {2}", - static_cast<double>(raw_elements) / static_cast<double>(iN * iN)); + for (std::size_t i = railstart; i < term_count; i++) + RHS_t += go[i] * *other_cur_analog[i]; - // FIXME: Move me + RHS[k] = RHS_t; - if (state().lib().isLoaded()) + // add diagonal element + mat.A[mat.diag[k]] = gtot_t; + + for (std::size_t i = 0; i < railstart; i++) + mat.A[tcr[i]] -= go[i]; + } +#else + for (std::size_t i = 0; i < railstart; i++) + mat.A[tcr[i]] -= go[i]; + + for (std::size_t i = 0; i < railstart; i++) { - pstring symname = static_compile_name(); - m_proc.load(this->state().lib(), symname); - if (m_proc.resolved()) - this->log().warning("External static solver {1} found ...", symname); - else - this->log().warning("External static solver {1} not found ...", symname); + gtot_t += gt[i]; + RHS_t += Idr[i]; } - } + for (std::size_t i = railstart; i < term_count; i++) + { + RHS_t += (Idr[i] + go[i] * *other_cur_analog[i]); + gtot_t += gt[i]; + } - template <typename FT, int SIZE> - void matrix_solver_GCR_t<FT, SIZE>::csc_private(plib::putf8_fmt_writer &strm) - { - const std::size_t iN = N(); + RHS[k] = RHS_t; + mat.A[mat.diag[k]] += gtot_t; + } +#endif + mat.ia[iN] = static_cast<mattype>(mat.nz_num); - for (std::size_t i = 0; i < mat.nz_num; i++) - strm("double m_A{1} = m_A[{2}];\n", i, i); + /* now solve it */ + //if (m_proc != nullptr) + if (m_proc.resolved()) + { + //static_solver(m_A, RHS); + m_proc(&mat.A[0], &RHS[0], &new_V[0]); + } + else + { for (std::size_t i = 0; i < iN - 1; i++) { const auto &nzbd = this->m_terms[i]->m_nzbd; @@ -212,114 +401,61 @@ namespace devices if (nzbd.size() > 0) { std::size_t pi = mat.diag[i]; + const nl_double f = 1.0 / mat.A[pi++]; + const std::size_t piie = mat.ia[i+1]; - //const FT f = 1.0 / m_A[pi++]; - strm("const double f{1} = 1.0 / m_A{2};\n", i, pi); - pi++; - const std::size_t piie = mat.row_idx[i+1]; - - //for (auto & j : nzbd) - for (std::size_t j : nzbd) + for (std::size_t j : nzbd) // for (std::size_t j = i + 1; j < iN; j++) { // proceed to column i - std::size_t pj = mat.row_idx[j]; + //__builtin_prefetch(&m_A[mat.diag[j+1]], 1); + std::size_t pj = mat.ia[j]; - while (mat.col_idx[pj] < i) + while (mat.ja[pj] < i) pj++; - //const FT f1 = - m_A[pj++] * f; - strm("\tconst double f{1}_{2} = -f{3} * m_A{4};\n", i, j, i, pj); - pj++; + const nl_double f1 = - mat.A[pj++] * f; // subtract row i from j */ for (std::size_t pii = pi; pii<piie; ) { - while (mat.col_idx[pj] < mat.col_idx[pii]) + while (mat.ja[pj] < mat.ja[pii]) pj++; - //m_A[pj++] += m_A[pii++] * f1; - strm("\tm_A{1} += m_A{2} * f{3}_{4};\n", pj, pii, i, j); - pj++; pii++; + mat.A[pj++] += mat.A[pii++] * f1; } - //RHS[j] += f1 * RHS[i]; - strm("\tRHS[{1}] += f{2}_{3} * RHS[{4}];\n", j, i, j, i); + RHS[j] += f1 * RHS[i]; } } } + /* backward substitution + * + */ + + /* row n-1 */ + new_V[iN - 1] = RHS[iN - 1] / mat.A[mat.diag[iN - 1]]; - //new_V[iN - 1] = RHS[iN - 1] / mat.A[mat.diag[iN - 1]]; - strm("\tV[{1}] = RHS[{2}] / m_A{3};\n", iN - 1, iN - 1, mat.diag[iN - 1]); for (std::size_t j = iN - 1; j-- > 0;) { - strm("\tdouble tmp{1} = 0.0;\n", j); - const std::size_t e = mat.row_idx[j+1]; - for (std::size_t pk = mat.diag[j] + 1; pk < e; pk++) + //__builtin_prefetch(&new_V[j-1], 1); + //if (j>0)__builtin_prefetch(&m_A[mat.diag[j-1]], 0); + double tmp = 0; + auto jdiag = mat.diag[j]; + const std::size_t e = mat.ia[j+1]; + for (std::size_t pk = jdiag + 1; pk < e; pk++) { - strm("\ttmp{1} += m_A{2} * V[{3}];\n", j, pk, mat.col_idx[pk]); + tmp += mat.A[pk] * new_V[mat.ja[pk]]; } - strm("\tV[{1}] = (RHS[{1}] - tmp{1}) / m_A{4};\n", j, j, j, mat.diag[j]); + new_V[j] = (RHS[j] - tmp) / mat.A[jdiag]; } } - template <typename FT, int SIZE> - pstring matrix_solver_GCR_t<FT, SIZE>::static_compile_name() - { - plib::postringstream t; - plib::putf8_fmt_writer w(&t); - csc_private(w); - std::hash<pstring> h; - - return plib::pfmt("nl_gcr_{1:x}_{2}")(h( t.str() ))(mat.nz_num); - } + this->m_stat_calculations++; - template <typename FT, int SIZE> - std::pair<pstring, pstring> matrix_solver_GCR_t<FT, SIZE>::create_solver_code() - { - plib::postringstream t; - plib::putf8_fmt_writer strm(&t); - pstring name = static_compile_name(); - - strm.writeline(plib::pfmt("extern \"C\" void {1}(double * __restrict m_A, double * __restrict RHS, double * __restrict V)\n")(name)); - strm.writeline("{\n"); - csc_private(strm); - strm.writeline("}\n"); - return std::pair<pstring, pstring>(name, t.str()); - } - - template <typename FT, int SIZE> - unsigned matrix_solver_GCR_t<FT, SIZE>::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, m_mat_ptr, RHS); - - /* now solve it */ - - //if (m_proc != nullptr) - if (m_proc.resolved()) - { - //static_solver(m_A, RHS); - m_proc(&mat.A[0], &RHS[0], &new_V[0]); - } - else - { - // mat.gaussian_elimination_parallel(RHS); - mat.gaussian_elimination(RHS); - /* backward substitution */ - mat.gaussian_back_substitution(new_V, RHS); - } - - this->m_stat_calculations++; - - const FT err = (newton_raphson ? delta(new_V) : 0.0); - store(new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; - } + const nl_double err = (newton_raphson ? delta(new_V) : 0.0); + store(new_V); + return (err > this->m_params.m_accuracy) ? 2 : 1; +} -} // namespace devices + } //namespace devices } // namespace netlist #endif /* NLD_MS_GCR_H_ */ |