// license:BSD-3-Clause // copyright-holders:Couriersud // Names // spell-checker: words Woodbury, #include "nld_solver.h" #include "nl_errstr.h" #include "nl_factory.h" #include "nl_setup.h" // FIXME: only needed for splitter code #include "nld_matrix_solver.h" #include "nld_ms_direct.h" #include "nld_ms_direct1.h" #include "nld_ms_direct2.h" #include "nld_ms_gcr.h" #include "nld_ms_gmres.h" #include "nld_ms_sm.h" #include "nld_ms_sor.h" #include "nld_ms_sor_mat.h" #include "nld_ms_w.h" #include "core/setup.h" #include "plib/pomp.h" #include "plib/ptimed_queue.h" #include #include namespace netlist::devices { // ------------------------------------------------------------------------- // solver // ------------------------------------------------------------------------- nld_solver::nld_solver(constructor_param_t data) : device_t(data) , m_fb_step(*this, "FB_step", NETLIB_DELEGATE(fb_step)) , m_Q_step(*this, "Q_step") , m_params(*this, "", solver::solver_parameter_defaults::get_instance()) , m_queue( this->state().pool(), config::max_solver_queue_size(), queue_type::id_delegate(&NETLIB_NAME(solver)::get_solver_id, this), queue_type::obj_delegate(&NETLIB_NAME(solver)::solver_by_id, this)) { // internal stuff state().save(*this, static_cast(m_queue), this->name(), "m_queue"); connect("FB_step", "Q_step"); } NETLIB_RESET(solver) { if (exec().stats_enabled()) m_fb_step.set_delegate(NETLIB_DELEGATE(fb_step)); for (auto &s : m_mat_solvers) s->reset(); for (auto &s : m_mat_solvers) m_queue.push({netlist_time_ext::zero(), s.get()}); } void NETLIB_NAME(solver)::stop() { for (auto &s : m_mat_solvers) s->log_stats(); } #if 1 template NETLIB_HANDLER(solver, fb_step) { const netlist_time_ext now(exec().time()); const std::size_t nthreads = m_params.m_parallel() < 2 ? 1 : std::min( static_cast( m_params.m_parallel()), plib::omp::get_max_threads()); const netlist_time_ext sched( now + (nthreads <= 1 ? netlist_time_ext::zero() : netlist_time_ext::from_nsec(100))); plib::uninitialised_array tmp; // NOLINT plib::uninitialised_array nt; // NOLINT std::size_t p = 0; while (!m_queue.empty()) { const auto t = m_queue.top().exec_time(); auto * o = m_queue.top().object(); if (t != now) if (t > sched) break; tmp[p++] = o; m_queue.pop(); } // FIXME: Disabled for now since parallel processing will decrease // performance // for tested applications. More testing required here if (true || nthreads < 2) { if (!KEEP_STATS) { for (std::size_t i = 0; i < p; i++) nt[i] = tmp[i]->solve(now, "no-parallel"); } else { stats()->m_stat_total_time.stop(); for (std::size_t i = 0; i < p; i++) { tmp[i]->stats()->m_stat_call_count.inc(); auto g(tmp[i]->stats()->m_stat_total_time.guard()); nt[i] = tmp[i]->solve(now, "no-parallel"); } stats()->m_stat_total_time.start(); } for (std::size_t i = 0; i < p; i++) { if (nt[i] != netlist_time::zero()) m_queue.push({now + nt[i], tmp[i]}); tmp[i]->update_inputs(); } } else { plib::omp::set_num_threads(nthreads); plib::omp::for_static(static_cast(0), p, [&tmp, &nt, now](std::size_t i) { nt[i] = tmp[i]->solve(now, "parallel"); }); for (std::size_t i = 0; i < p; i++) { if (nt[i] != netlist_time::zero()) m_queue.push({now + nt[i], tmp[i]}); tmp[i]->update_inputs(); } } if (!m_queue.empty()) m_Q_step.net().toggle_and_push_to_queue( static_cast(m_queue.top().exec_time() - now)); } void NETLIB_NAME(solver)::reschedule(solver::matrix_solver_t *solv, netlist_time ts) { const netlist_time_ext now(exec().time()); const netlist_time_ext sched(now + ts); m_queue.remove(solv); m_queue.push({sched, solv}); if (m_Q_step.net().is_queued()) { if (m_Q_step.net().next_scheduled_time() > sched) m_Q_step.net().toggle_and_push_to_queue(ts); } else m_Q_step.net().toggle_and_push_to_queue(ts); } #else NETLIB_HANDLER(solver, fb_step) { if (m_params.m_dynamic_ts) return; netlist_time_ext now(exec().time()); // force solving during start up if there are no time-step devices // FIXME: Needs a more elegant solution bool force_solve = (now < netlist_time_ext::from_fp< decltype(m_params.m_max_time_step)>( 2 * m_params.m_max_time_step)); std::size_t nthreads = std::min( static_cast(m_params.m_parallel()), plib::omp::get_max_threads()); std::vector &solvers = (force_solve ? m_mat_solvers_all : m_mat_solvers_time_stepping); if (nthreads > 1 && solvers.size() > 1) { plib::omp::set_num_threads(nthreads); plib::omp::for_static( static_cast(0), solvers.size(), [&solvers, now](std::size_t i) { [[maybe_unused]] const netlist_time ts = solvers[i] ->ptr->solve( now); }); } else for (auto &solver : solvers) { [[maybe_unused]] const netlist_time ts = solver->ptr->solve( now); } for (auto &solver : solvers) solver->ptr->update_inputs(); // step circuit if (!m_Q_step.net().is_queued()) { m_Q_step.net().toggle_and_push_to_queue( netlist_time::from_fp(m_params.m_max_time_step)); } } #endif // FIXME: should be created in device space template NETLIB_NAME(solver)::solver_ptr create_it(A &arena, NETLIB_NAME(solver) &main_solver, pstring name, NETLIB_NAME(solver)::net_list_t & nets, const solver::solver_parameters_t *params, std::size_t size) { return plib::make_unique(arena, main_solver, name, nets, params, size); } template NETLIB_NAME(solver)::solver_ptr NETLIB_NAME(solver)::create_solver( std::size_t size, const pstring &solver_name, const solver::solver_parameters_t *params, NETLIB_NAME(solver)::net_list_t & nets) { switch (params->m_method()) { case solver::matrix_type_e::MAT_CR: return create_it>( state().pool(), *this, solver_name, nets, params, size); case solver::matrix_type_e::MAT: return create_it>( state().pool(), *this, solver_name, nets, params, size); case solver::matrix_type_e::GMRES: return create_it>( state().pool(), *this, solver_name, nets, params, size); #if (NL_USE_ACADEMIC_SOLVERS) case solver::matrix_type_e::SOR: return create_it>( state().pool(), *this, solver_name, nets, params, size); case solver::matrix_type_e::SOR_MAT: return create_it>( state().pool(), *this, solver_name, nets, params, size); case solver::matrix_type_e::SM: // Sherman-Morrison Formula return create_it>( state().pool(), *this, solver_name, nets, params, size); case solver::matrix_type_e::W: // Woodbury Formula return create_it>( state().pool(), *this, solver_name, nets, params, size); #else // case solver::matrix_type_e::GMRES: case solver::matrix_type_e::SOR: case solver::matrix_type_e::SOR_MAT: case solver::matrix_type_e::SM: case solver::matrix_type_e::W: state().log().warning(MW_SOLVER_METHOD_NOT_SUPPORTED( params->m_method().name(), "MAT_CR")); return create_it>( state().pool(), *this, solver_name, nets, params, size); #endif } return solver_ptr(); } template NETLIB_NAME(solver)::solver_ptr NETLIB_NAME(solver)::create_solvers( const pstring &sname, const solver::solver_parameters_t *params, net_list_t &nets) { std::size_t net_count = nets.size(); switch (net_count) { #if !defined(__EMSCRIPTEN__) case 1: return plib::make_unique>( state().pool(), *this, sname, nets, params); case 2: return plib::make_unique>( state().pool(), *this, sname, nets, params); case 3: return create_solver(3, sname, params, nets); case 4: return create_solver(4, sname, params, nets); case 5: return create_solver(5, sname, params, nets); case 6: return create_solver(6, sname, params, nets); case 7: return create_solver(7, sname, params, nets); case 8: return create_solver(8, sname, params, nets); #endif default: log().info(MI_NO_SPECIFIC_SOLVER(net_count)); if (net_count <= 16) { return create_solver(net_count, sname, params, nets); } if (net_count <= 32) { return create_solver(net_count, sname, params, nets); } if (net_count <= 64) { return create_solver(net_count, sname, params, nets); } if (net_count <= 128) { return create_solver(net_count, sname, params, nets); } if (net_count <= 256) { return create_solver(net_count, sname, params, nets); } if (net_count <= 512) { return create_solver(net_count, sname, params, nets); } return create_solver(net_count, sname, params, nets); } } struct net_splitter { void run(netlist_state_t &nl_state) { for (auto &net : nl_state.nets()) { nl_state.log().verbose("processing {1}", net->name()); if (!net->is_rail_net() && !net->core_terms_empty()) { nl_state.log().verbose(" ==> not a rail net"); // Must be an analog net auto n = plib::dynamic_downcast(net.get()); nl_assert_always(bool(n), "Unable to cast to analog_net_t &"); if (!already_processed(*(*n))) { groupspre.emplace_back( NETLIB_NAME(solver)::net_list_t()); process_net(nl_state, *(*n)); } } } for (auto &g : groupspre) if (!g.empty()) groups.push_back(g); } std::vector groups; private: bool already_processed(const analog_net_t &n) const { // no need to process rail nets - these are known variables if (n.is_rail_net()) return true; // if it's already processed - no need to continue for (const auto &grp : groups) if (plib::container::contains(grp, &n)) return true; return false; } bool check_if_processed_and_join(const analog_net_t &n) { // no need to process rail nets - these are known variables if (n.is_rail_net()) return true; // First check if it is in a previous group. // In this case we need to merge this group into the current group if (groupspre.size() > 1) { for (std::size_t i = 0; i < groupspre.size() - 1; i++) if (plib::container::contains(groupspre[i], &n)) { // copy all nets for (auto &cn : groupspre[i]) if (!plib::container::contains(groupspre.back(), cn)) groupspre.back().push_back(cn); // clear groupspre[i].clear(); return true; } } // if it's already processed - no need to continue if (!groupspre.empty() && plib::container::contains(groupspre.back(), &n)) return true; return false; } // NOLINTNEXTLINE(misc-no-recursion) void process_net(netlist_state_t &nl_state, analog_net_t &n) { // ignore empty nets. FIXME: print a warning message nl_state.log().verbose("Net {}", n.name()); auto terminals(n.core_terms_copy()); if (!terminals.empty()) { // add the net groupspre.back().push_back(&n); // process all terminals connected to this net for (detail::core_terminal_t *term : terminals) { nl_state.log().verbose("Term {} {}", term->name(), static_cast(term->type())); // only process analog terminals if (term->is_type(detail::terminal_type::TERMINAL)) { auto pt = plib::dynamic_downcast(term); nl_assert_always(bool(pt), "Error casting *term to terminal_t &"); // check the connected terminal const auto *const connected_terminals = nl_state.setup().get_connected_terminals(*(*pt)); // NOLINTNEXTLINE proposal does not work for VS for (auto ct = connected_terminals->begin(); *ct != nullptr; ct++) { analog_net_t &connected_net = (*ct)->net(); nl_state.log().verbose(" Connected net {}", connected_net.name()); if (!check_if_processed_and_join(connected_net)) process_net(nl_state, connected_net); } } } } } std::vector groupspre; }; void NETLIB_NAME(solver)::post_start() { log().verbose("Scanning net groups ..."); // determine net groups net_splitter splitter; splitter.run(state()); log().verbose("Found {1} net groups in {2} nets\n", splitter.groups.size(), state().nets().size()); int num_errors = 0; log().verbose("checking net consistency ..."); for (const auto &grp : splitter.groups) { int rail_terminals = 0; pstring nets_in_grp; for (const auto &n : grp) { nets_in_grp += (n->name() + " "); if (!n->is_analog()) { state().log().error( ME_SOLVER_CONSISTENCY_NOT_ANALOG_NET(n->name())); num_errors++; } if (n->is_rail_net()) { state().log().error( ME_SOLVER_CONSISTENCY_RAIL_NET(n->name())); num_errors++; } for (detail::core_terminal_t *t : n->core_terms_copy()) { if (!t->has_net()) { state().log().error( ME_SOLVER_TERMINAL_NO_NET(t->name())); num_errors++; } else { if (auto other_terminal = plib::dynamic_downcast< terminal_t *>(t)) if (state() .setup() .get_connected_terminal(*(*other_terminal)) ->net() .is_rail_net()) rail_terminals++; } } } if (rail_terminals == 0) { state().log().error(ME_SOLVER_NO_RAIL_TERMINAL(nets_in_grp)); num_errors++; } } if (num_errors > 0) throw nl_exception(MF_SOLVER_CONSISTENCY_ERRORS(num_errors)); // setup the solvers for (auto &grp : splitter.groups) { solver_ptr ms; pstring sname = plib::pfmt("Solver_{1}")(m_mat_solvers.size()); params_uptr params = plib::make_unique( state().pool(), *this, sname + ".", m_params); switch (params->m_fp_type()) { case solver::matrix_fp_type_e::FLOAT: if (!config::use_float_matrix::value) log().info("FPTYPE {1} not supported. Using DOUBLE", params->m_fp_type().name()); ms = create_solvers>( sname, params.get(), grp); break; case solver::matrix_fp_type_e::DOUBLE: ms = create_solvers(sname, params.get(), grp); break; case solver::matrix_fp_type_e::LONGDOUBLE: if (!config::use_long_double_matrix::value) log().info("FPTYPE {1} not supported. Using DOUBLE", params->m_fp_type().name()); ms = create_solvers>(sname, params.get(), grp); break; case solver::matrix_fp_type_e::FLOATQ128: #if (NL_USE_FLOAT128) ms = create_solvers(sname, params.get(), grp); #else log().info("FPTYPE {1} not supported. Using DOUBLE", params->m_fp_type().name()); ms = create_solvers(sname, params.get(), grp); #endif break; } state().register_device( ms->name(), device_arena::owned_ptr(ms.get(), false)); log().verbose("Solver {1}", ms->name()); log().verbose(" ==> {1} nets", grp.size()); log().verbose(" has {1} dynamic elements", ms->dynamic_device_count()); log().verbose(" has {1} time step elements", ms->time_step_device_count()); for (auto &n : grp) { log().verbose("Net {1}", n->name()); for (const detail::core_terminal_t *t : n->core_terms_copy()) { log().verbose(" {1}", t->name()); } } m_mat_params.push_back(std::move(params)); m_mat_solvers.push_back(std::move(ms)); } } solver::static_compile_container NETLIB_NAME(solver)::create_solver_code( solver::static_compile_target target) { solver::static_compile_container mp; for (auto &s : m_mat_solvers) { auto r = s->create_solver_code(target); if (!r.first.empty()) // ignore solvers not supporting static // compile mp.push_back(r); } return mp; } std::size_t NETLIB_NAME(solver)::get_solver_id( const solver::matrix_solver_t *net) const { for (std::size_t i = 0; i < m_mat_solvers.size(); i++) if (m_mat_solvers[i].get() == net) return i; return std::numeric_limits::max(); } solver::matrix_solver_t *NETLIB_NAME(solver)::solver_by_id( std::size_t id) const { return m_mat_solvers[id].get(); } NETLIB_DEVICE_IMPL(solver, "SOLVER", "FREQ") } // namespace netlist::devices