// license:GPL-2.0+ // copyright-holders:Couriersud #ifndef NLD_MS_GCR_H_ #define NLD_MS_GCR_H_ /// /// \file nld_ms_gcr.h /// /// Gaussian elimination using compressed row format. /// #include "plib/mat_cr.h" //#include "nld_ms_direct.h" #include "nld_matrix_solver_ext.h" #include "nld_solver.h" #include "plib/pdynlib.h" #include "plib/pstream.h" #include "plib/vector_ops.h" #include namespace netlist { namespace solver { template class matrix_solver_GCR_t: public matrix_solver_ext_t { public: using mat_type = plib::pGEmatrix_cr_t>; using base_type = matrix_solver_ext_t; using fptype = typename base_type::fptype; matrix_solver_GCR_t(devices::nld_solver &main_solver, const pstring &name, const matrix_solver_t::net_list_t &nets, const solver::solver_parameters_t *params, const std::size_t size) : matrix_solver_ext_t(main_solver, name, nets, params, size) , mat(static_cast(size)) , m_proc() { const std::size_t iN = this->size(); // build the final matrix std::vector> fill(iN); std::size_t raw_elements = 0; for (std::size_t k = 0; k < iN; k++) { fill[k].resize(iN, decltype(mat)::FILL_INFINITY); for (auto &j : this->m_terms[k].m_nz) { fill[k][j] = 0; raw_elements++; } } auto gr = mat.gaussian_extend_fill_mat(fill); this->log_fill(fill, mat); mat.build_from_fill_mat(fill); for (mat_index_type k=0; km_terms[k].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(mat.col_idx[i])) { this->m_mat_ptr[k][j] = &mat.A[i]; cnt++; break; } } nl_assert(cnt == this->m_terms[k].railstart()); this->m_mat_ptr[k][this->m_terms[k].railstart()] = &mat.A[mat.diag[k]]; } this->state().log().verbose("maximum fill: {1}", gr.first); this->state().log().verbose("Post elimination occupancy ratio: {2} Ops: {1}", gr.second, static_cast(mat.nz_num) / static_cast(iN * iN)); this->state().log().verbose(" Pre elimination occupancy ratio: {1}", static_cast(raw_elements) / static_cast(iN * iN)); // FIXME: Move me // // During extended validation there is no reason to check for // differences in the generated code since during // extended validation this will be different (and non-functional) if (!this->state().is_extended_validation() && this->state().static_solver_lib().isLoaded()) { pstring symname = static_compile_name(); m_proc.load(this->state().static_solver_lib(), symname); if (m_proc.resolved()) { this->state().log().info("External static solver {1} found ...", symname); } else { this->state().log().warning("External static solver {1} not found ...", symname); } } } void vsolve_non_dynamic() override; std::pair create_solver_code(static_compile_target target) override; private: using mat_index_type = typename plib::pmatrix_cr_t::index_type; void generate_code(plib::putf8_fmt_writer &strm); pstring static_compile_name(); mat_type mat; plib::dynproc m_proc; }; // ---------------------------------------------------------------------------------------- // matrix_solver - GCR // ---------------------------------------------------------------------------------------- template void matrix_solver_GCR_t::generate_code(plib::putf8_fmt_writer &strm) { const std::size_t iN = this->size(); pstring fptype(fp_constants::name()); pstring fpsuffix(fp_constants::suffix()); // avoid unused variable warnings strm("\tplib::unused_var({1});\n", "cnV"); for (std::size_t i = 0; i < mat.nz_num; i++) strm("\t{1} m_A{2}(0.0);\n", fptype, i, i); for (std::size_t k = 0; k < iN; k++) { auto &net = this->m_terms[k]; // FIXME: gonn, gtn and Idr - which float types should they have? //auto gtot_t = std::accumulate(gt, gt + term_count, plib::constants::zero()); //*tcr_r[railstart] = static_cast(gtot_t); //mat.A[mat.diag[k]] += gtot_t; auto pd = this->m_mat_ptr[k][net.railstart()] - &this->mat.A[0]; #if 0 pstring terms = plib::pfmt("m_A{1} = gt[{2}]")(pd, this->m_gtn.didx(k,0)); for (std::size_t i=1; i < net.count(); i++) terms += plib::pfmt(" + gt[{1}]")(this->m_gtn.didx(k,i)); strm("\t{1};\n", terms); #else for (std::size_t i=0; i < net.count(); i++) strm("\tm_A{1} += gt[{2}];\n", pd, this->m_gtn.didx(k,i)); #endif //for (std::size_t i = 0; i < railstart; i++) // *tcr_r[i] += static_cast(go[i]); for (std::size_t i = 0; i < net.railstart(); i++) { auto p = this->m_mat_ptr[k][i] - &this->mat.A[0]; strm("\tm_A{1} += go[{2}];\n", p, this->m_gonn.didx(k,i)); } #if 0 //auto RHS_t(std::accumulate(Idr, Idr + term_count, plib::constants::zero())); terms = plib::pfmt("{1} RHS{2} = Idr[{3}]")(fptype, k, this->m_Idrn.didx(k,0)); for (std::size_t i=1; i < net.count(); i++) terms += plib::pfmt(" + Idr[{1}]")(this->m_Idrn.didx(k,i)); //for (std::size_t i = railstart; i < term_count; i++) // RHS_t += (- go[i]) * *cnV[i]; for (std::size_t i = net.railstart(); i < net.count(); i++) terms += plib::pfmt(" - go[{1}] * *cnV[{2}]")(this->m_gonn.didx(k,i), this->m_connected_net_Vn.didx(k,i)); strm("\t{1};\n", terms); #else //auto RHS_t(std::accumulate(Idr, Idr + term_count, plib::constants::zero())); strm("\t{1} RHS{2} = Idr[{3}];\n", fptype, k, this->m_Idrn.didx(k,0)); for (std::size_t i=1; i < net.count(); i++) strm("\tRHS{1} += Idr[{2}];\n", k, this->m_Idrn.didx(k,i)); //for (std::size_t i = railstart; i < term_count; i++) // RHS_t += (- go[i]) * *cnV[i]; for (std::size_t i = net.railstart(); i < net.count(); i++) strm("\tRHS{1} -= go[{2}] * *cnV[{3}];\n", k, this->m_gonn.didx(k,i), this->m_connected_net_Vn.didx(k,i)); #endif } for (std::size_t i = 0; i < iN - 1; i++) { //const auto &nzbd = this->m_terms[i].m_nzbd; const auto *nzbd = mat.nzbd(i); const auto nzbd_count = mat.nzbd_count(i); if (nzbd_count > 0) { std::size_t pi = mat.diag[i]; //const FT f = 1.0 / m_A[pi++]; strm("\tconst {1} f{2} = 1.0{3} / m_A{4};\n", fptype, i, fpsuffix, pi); pi++; const std::size_t piie = mat.row_idx[i+1]; //for (auto & j : nzbd) for (std::size_t jj = 0; jj < nzbd_count; jj++) { std::size_t j = nzbd[jj]; // proceed to column i std::size_t pj = mat.row_idx[j]; while (mat.col_idx[pj] < i) pj++; //const FT f1 = - m_A[pj++] * f; strm("\tconst {1} f{2}_{3} = -f{4} * m_A{5};\n", fptype, i, j, i, pj); pj++; // subtract row i from j for (std::size_t pii = pi; pii 0;) { #if 1 strm("\t{1} tmp{2} = 0.0{3};\n", fptype, j, fpsuffix); const std::size_t e = mat.row_idx[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.col_idx[pk]); } strm("\tV[{1}] = (RHS{1} - tmp{1}) / m_A{4};\n", j, j, j, mat.diag[j]); #else pstring tmp; const std::size_t e = mat.row_idx[j+1]; for (std::size_t pk = mat.diag[j] + 1; pk < e; pk++) { tmp = tmp + plib::pfmt(" + m_A{2} * V[{3}]")(j, pk, mat.col_idx[pk]); } if (tmp.empty()) { strm("\tV[{1}] = RHS{1} / m_A{2};\n", j, mat.diag[j]); } else { //strm("\tconst {1} tmp{2} = {3};\n", fptype, j, tmp.substr(3)); //strm("\tV[{1}] = (RHS{1} - tmp{1}) / m_A{2};\n", j, mat.diag[j]); strm("\tV[{1}] = (RHS{1} - ({2})) / m_A{3};\n", j, tmp.substr(3), mat.diag[j]); } #endif } } template pstring matrix_solver_GCR_t::static_compile_name() { pstring str_floattype(fp_constants::name()); pstring str_fptype(fp_constants::name()); std::stringstream t; t.imbue(std::locale::classic()); plib::putf8_fmt_writer w(&t); generate_code(w); //std::hash::type>::type> h; return plib::pfmt("nl_gcr_{1:x}_{2}_{3}_{4}")(plib::hash( t.str().c_str(), t.str().size() ))(mat.nz_num)(str_fptype)(str_floattype); } template std::pair matrix_solver_GCR_t::create_solver_code(static_compile_target target) { std::stringstream t; t.imbue(std::locale::classic()); plib::putf8_fmt_writer strm(&t); pstring name = static_compile_name(); pstring str_floattype(fp_constants::name()); pstring str_fptype(fp_constants::name()); pstring extqual; if (target == CXX_EXTERNAL_C) extqual = "extern \"C\""; else if (target == CXX_STATIC) extqual = "static"; strm.writeline(plib::pfmt("{1} void {2}({3} * __restrict V, " "const {4} * __restrict go, const {4} * __restrict gt, " "const {4} * __restrict Idr, const {4} * const * __restrict cnV)\n")(extqual, name, str_floattype, str_fptype)); strm.writeline("{\n"); generate_code(strm); strm.writeline("}\n"); // some compilers (_WIN32, _WIN64, mac osx) need an explicit cast return std::pair(name, putf8string(t.str())); } template void matrix_solver_GCR_t::vsolve_non_dynamic() { if (m_proc.resolved()) { m_proc(&this->m_new_V[0], this->m_gonn.data(), this->m_gtn.data(), this->m_Idrn.data(), this->m_connected_net_Vn.data()); } else { // clear matrix mat.set_scalar(plib::constants::zero()); // populate matrix this->fill_matrix_and_rhs(); // now solve it // parallel is slow -- very slow // mat.gaussian_elimination_parallel(RHS); mat.gaussian_elimination(this->m_RHS); // backward substitution mat.gaussian_back_substitution(this->m_new_V, this->m_RHS); } } } // namespace solver } // namespace netlist #endif // NLD_MS_GCR_H_