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, 213 insertions, 349 deletions
diff --git a/src/lib/netlist/solver/nld_ms_gcr.h b/src/lib/netlist/solver/nld_ms_gcr.h index a17ca1ec282..f3e56a342d6 100644 --- a/src/lib/netlist/solver/nld_ms_gcr.h +++ b/src/lib/netlist/solver/nld_ms_gcr.h @@ -10,390 +10,201 @@ #ifndef NLD_MS_GCR_H_ #define NLD_MS_GCR_H_ -#include <algorithm> +#include "plib/mat_cr.h" -#include "../plib/pdynlib.h" -#include "mat_cr.h" #include "nld_ms_direct.h" #include "nld_solver.h" -#include "vector_base.h" -#include "../plib/pstream.h" +#include "plib/pdynlib.h" +#include "plib/pstream.h" +#include "plib/vector_ops.h" + +#include <algorithm> namespace netlist { - namespace devices - { -template <std::size_t m_N, std::size_t storage_N> -class matrix_solver_GCR_t: public matrix_solver_t +namespace devices { -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() - { - } - virtual ~matrix_solver_GCR_t() override + template <typename FT, int SIZE> + class matrix_solver_GCR_t: public matrix_solver_t { - } - - constexpr std::size_t N() const { return (m_N == 0) ? m_dim : m_N; } - - virtual void vsetup(analog_net_t::list_t &nets) override; - virtual unsigned vsolve_non_dynamic(const bool newton_raphson) 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() + { + } - virtual std::pair<pstring, pstring> create_solver_code() override; + constexpr std::size_t N() const { return m_dim; } -private: + void vsetup(analog_net_t::list_t &nets) override; + unsigned vsolve_non_dynamic(const bool newton_raphson) override; - //typedef typename mat_cr_t<storage_N>::type mattype; - typedef typename mat_cr_t<storage_N>::index_type mattype; + std::pair<pstring, pstring> create_solver_code() override; - void csc_private(plib::putf8_fmt_writer &strm); + private: - using extsolver = void (*)(double * RESTRICT m_A, double * RESTRICT RHS, double * RESTRICT V); + using mat_index_type = typename plib::matrix_compressed_rows_t<FT, SIZE>::index_type; - pstring static_compile_name(); + void csc_private(plib::putf8_fmt_writer &strm); - const std::size_t m_dim; - std::vector<unsigned> m_term_cr[storage_N]; - mat_cr_t<storage_N> mat; + using extsolver = void (*)(double * m_A, double * RHS, double * V); - //extsolver m_proc; - plib::dynproc<void, double * RESTRICT, double * RESTRICT, double * RESTRICT> m_proc; + pstring static_compile_name(); -}; + const std::size_t m_dim; + plib::parray<FT, SIZE> RHS; + plib::parray<FT, SIZE> new_V; -// ---------------------------------------------------------------------------------------- -// matrix_solver - GCR -// ---------------------------------------------------------------------------------------- + mat_type mat; -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); + //extsolver m_proc; + plib::dynproc<void, double * , double * , double * > m_proc; - mattype nz = 0; - const std::size_t iN = this->N(); + }; - /* build the final matrix */ + // ---------------------------------------------------------------------------------------- + // matrix_solver - GCR + // ---------------------------------------------------------------------------------------- - bool touched[storage_N][storage_N] = { { false } }; - for (std::size_t k = 0; k < iN; k++) + // FIXME: namespace or static class member + template <typename V> + std::size_t inline get_level(const V &v, std::size_t k) { - for (auto &j : this->m_terms[k]->m_nz) - touched[k][j] = true; + 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"); } - unsigned fc = 0; + template <typename FT, int SIZE> + void matrix_solver_GCR_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets) + { + setup_base(nets); - unsigned ops = 0; + const std::size_t iN = this->N(); - 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; - } - } - } - } + /* build the final matrix */ + std::vector<std::vector<unsigned>> fill(iN); - for (mattype k=0; k<iN; k++) - { - mat.ia[k] = nz; + std::size_t raw_elements = 0; - for (mattype j=0; j<iN; j++) + for (std::size_t k = 0; k < iN; k++) { - if (touched[k][j]) + fill[k].resize(iN, decltype(mat)::FILL_INFINITY); + for (auto &j : this->m_terms[k]->m_nz) { - mat.ja[nz] = j; - if (j == k) - mat.diag[k] = nz; - nz++; + fill[k][j] = 0; + raw_elements++; } - } - 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); - } - mat.ia[iN] = nz; - mat.nz_num = nz; + auto gr = mat.gaussian_extend_fill_mat(fill); - 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 - } + /* FIXME: move this to the cr matrix class and use computed + * parallel ordering once it makes sense. + */ -} -#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; + std::vector<unsigned> levL(iN, 0); + std::vector<unsigned> levU(iN, 0); - if (nzbd.size() > 0) + // parallel scheme for L x = y + for (std::size_t k = 0; k < iN; k++) { - 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); - } + 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; } - } - //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++) + // parallel scheme for U x = y + for (std::size_t k = iN; k-- > 0; ) { - strm("\ttmp{1} += m_A[{2}] * V[{3}];\n", j, pk, mat.ja[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("\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 i = 0; i < iN - 1; i++) - { - const auto &nzbd = this->m_terms[i]->m_nzbd; - - if (nzbd.size() > 0) + for (std::size_t k = 0; k < iN; k++) { - 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) + unsigned fm = 0; + pstring ml = ""; + for (std::size_t j = 0; j < iN; 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); + 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]; } + 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()); -} -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]; + mat.build_from_fill_mat(fill); - 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++) + for (mat_index_type k=0; k<iN; k++) { - gtot_t += gt[i]; - RHS_t += Idr[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]]; } - for (std::size_t i = railstart; i < term_count; i++) - RHS_t += go[i] * *other_cur_analog[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)); - RHS[k] = RHS_t; + // FIXME: Move me - // 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++) + if (state().lib().isLoaded()) { - gtot_t += gt[i]; - RHS_t += Idr[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); } - for (std::size_t i = railstart; i < term_count; i++) - { - RHS_t += (Idr[i] + go[i] * *other_cur_analog[i]); - gtot_t += gt[i]; - } - - RHS[k] = RHS_t; - mat.A[mat.diag[k]] += gtot_t; } -#endif - mat.ia[iN] = static_cast<mattype>(mat.nz_num); - - /* 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 + 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(); + + 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 i = 0; i < iN - 1; i++) { const auto &nzbd = this->m_terms[i]->m_nzbd; @@ -401,61 +212,114 @@ unsigned matrix_solver_GCR_t<m_N, storage_N>::vsolve_non_dynamic(const bool newt 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]; - for (std::size_t j : nzbd) // for (std::size_t j = i + 1; j < iN; j++) + //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) { // proceed to column i - //__builtin_prefetch(&m_A[mat.diag[j+1]], 1); - std::size_t pj = mat.ia[j]; + std::size_t pj = mat.row_idx[j]; - while (mat.ja[pj] < i) + while (mat.col_idx[pj] < i) pj++; - const nl_double f1 = - mat.A[pj++] * f; + //const FT 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]) + while (mat.col_idx[pj] < mat.col_idx[pii]) pj++; - mat.A[pj++] += mat.A[pii++] * f1; + //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]; + //RHS[j] += f1 * RHS[i]; + strm("\tRHS[{1}] += f{2}_{3} * RHS[{4}];\n", j, i, j, 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;) { - //__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("\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++) { - tmp += mat.A[pk] * new_V[mat.ja[pk]]; + strm("\ttmp{1} += m_A{2} * V[{3}];\n", j, pk, mat.col_idx[pk]); } - new_V[j] = (RHS[j] - tmp) / mat.A[jdiag]; + strm("\tV[{1}] = (RHS[{1}] - tmp{1}) / m_A{4};\n", j, j, j, mat.diag[j]); } } - this->m_stat_calculations++; + 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); + } - const nl_double err = (newton_raphson ? delta(new_V) : 0.0); - store(new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; -} + 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; + } - } //namespace devices +} // namespace devices } // namespace netlist #endif /* NLD_MS_GCR_H_ */ |