// license:GPL-2.0+ // copyright-holders:Couriersud #include "nld_matrix_solver.h" #include "nl_setup.h" #include "plib/putil.h" namespace netlist { namespace solver { terms_for_net_t::terms_for_net_t(analog_net_t * net) : m_net(net) , m_railstart(0) { } void terms_for_net_t::add_terminal(terminal_t *term, int net_other, bool sorted) { if (sorted) for (std::size_t i=0; i < m_connected_net_idx.size(); i++) { if (m_connected_net_idx[i] > net_other) { plib::container::insert_at(m_terms, i, term); plib::container::insert_at(m_connected_net_idx, i, net_other); return; } } m_terms.push_back(term); m_connected_net_idx.push_back(net_other); } // ---------------------------------------------------------------------------------------- // matrix_solver // ---------------------------------------------------------------------------------------- matrix_solver_t::matrix_solver_t(netlist_state_t &anetlist, const pstring &name, const analog_net_t::list_t &nets, const solver_parameters_t *params) : device_t(anetlist, name) , m_params(*params) , m_iterative_fail(*this, "m_iterative_fail", 0) , m_iterative_total(*this, "m_iterative_total", 0) , m_stat_calculations(*this, "m_stat_calculations", 0) , m_stat_newton_raphson(*this, "m_stat_newton_raphson", 0) , m_stat_vsolver_calls(*this, "m_stat_vsolver_calls", 0) , m_last_step(*this, "m_last_step", netlist_time_ext::zero()) , m_fb_sync(*this, "FB_sync") , m_Q_sync(*this, "Q_sync") , m_ops(0) { connect_post_start(m_fb_sync, m_Q_sync); setup_base(nets); // now setup the matrix setup_matrix(); } analog_net_t *matrix_solver_t::get_connected_net(terminal_t *term) { return &state().setup().get_connected_terminal(*term)->net(); } void matrix_solver_t::setup_base(const analog_net_t::list_t &nets) { log().debug("New solver setup\n"); std::vector step_devices; std::vector dynamic_devices; m_terms.clear(); for (const auto & net : nets) { m_terms.emplace_back(net); m_rails_temp.emplace_back(); } for (std::size_t k = 0; k < nets.size(); k++) { analog_net_t *net = nets[k]; log().debug("adding net with {1} populated connections\n", net->core_terms().size()); net->set_solver(this); for (auto &p : net->core_terms()) { log().debug("{1} {2} {3}\n", p->name(), net->name(), net->is_rail_net()); switch (p->type()) { case detail::terminal_type::TERMINAL: if (p->device().is_timestep()) if (!plib::container::contains(step_devices, &p->device())) step_devices.push_back(&p->device()); if (p->device().is_dynamic()) if (!plib::container::contains(dynamic_devices, &p->device())) dynamic_devices.push_back(&p->device()); { auto *pterm = dynamic_cast(p); add_term(k, pterm); } log().debug("Added terminal {1}\n", p->name()); break; case detail::terminal_type::INPUT: { proxied_analog_output_t *net_proxy_output = nullptr; for (auto & input : m_inps) if (input->proxied_net() == &p->net()) { net_proxy_output = input.get(); break; } if (net_proxy_output == nullptr) { pstring nname(this->name() + "." + pstring(plib::pfmt("m{1}")(m_inps.size()))); nl_assert(p->net().is_analog()); auto net_proxy_output_u = state().make_object(*this, nname, static_cast(&p->net())); net_proxy_output = net_proxy_output_u.get(); m_inps.emplace_back(std::move(net_proxy_output_u)); } net_proxy_output->net().add_terminal(*p); // FIXME: repeated calling - kind of brute force net_proxy_output->net().rebuild_list(); log().debug("Added input {1}", net_proxy_output->name()); } break; case detail::terminal_type::OUTPUT: log().fatal(MF_UNHANDLED_ELEMENT_1_FOUND(p->name())); throw nl_exception(MF_UNHANDLED_ELEMENT_1_FOUND(p->name())); } } } for (auto &d : step_devices) m_step_funcs.emplace_back(nldelegate_ts(&core_device_t::timestep, d)); for (auto &d : dynamic_devices) m_dynamic_funcs.emplace_back(nldelegate_dyn(&core_device_t::update_terminals, d)); } void matrix_solver_t::sort_terms(matrix_sort_type_e sort) { // Sort in descending order by number of connected matrix voltages. // The idea is, that for Gauss-Seidel algo the first voltage computed // depends on the greatest number of previous voltages thus taking into // account the maximum amout of information. // // This actually improves performance on popeye slightly. Average // GS computations reduce from 2.509 to 2.370 // // Smallest to largest : 2.613 // Unsorted : 2.509 // Largest to smallest : 2.370 // // Sorting as a general matrix pre-conditioning is mentioned in // literature but I have found no articles about Gauss Seidel. // // For Gaussian Elimination however increasing order is better suited. // NOTE: Even better would be to sort on elements right of the matrix diagonal. // const std::size_t iN = m_terms.size(); switch (sort) { case matrix_sort_type_e::PREFER_BAND_MATRIX: { for (std::size_t k = 0; k < iN - 1; k++) { auto pk = get_weight_around_diag(k,k); for (std::size_t i = k+1; i < iN; i++) { auto pi = get_weight_around_diag(i,k); if (pi < pk) { std::swap(m_terms[i], m_terms[k]); pk = get_weight_around_diag(k,k); } } } } break; case matrix_sort_type_e::PREFER_IDENTITY_TOP_LEFT: { for (std::size_t k = 0; k < iN - 1; k++) { auto pk = get_left_right_of_diag(k,k); for (std::size_t i = k+1; i < iN; i++) { auto pi = get_left_right_of_diag(i,k); if (pi.first <= pk.first && pi.second >= pk.second) { std::swap(m_terms[i], m_terms[k]); pk = get_left_right_of_diag(k,k); } } } } break; case matrix_sort_type_e::ASCENDING: case matrix_sort_type_e::DESCENDING: { int sort_order = (sort == matrix_sort_type_e::DESCENDING ? 1 : -1); for (std::size_t k = 0; k < iN - 1; k++) for (std::size_t i = k+1; i < iN; i++) { if ((static_cast(m_terms[k].railstart()) - static_cast(m_terms[i].railstart())) * sort_order < 0) { std::swap(m_terms[i], m_terms[k]); } } } break; case matrix_sort_type_e::NOSORT: break; } // rebuild for (auto &term : m_terms) { //int *other = term.m_connected_net_idx.data(); for (std::size_t i = 0; i < term.count(); i++) //FIXME: this is weird if (term.m_connected_net_idx[i] != -1) term.m_connected_net_idx[i] = get_net_idx(get_connected_net(term.terms()[i])); } } void matrix_solver_t::setup_matrix() { const std::size_t iN = m_terms.size(); for (std::size_t k = 0; k < iN; k++) { m_terms[k].set_railstart(m_terms[k].count()); for (std::size_t i = 0; i < m_rails_temp[k].count(); i++) this->m_terms[k].add_terminal(m_rails_temp[k].terms()[i], m_rails_temp[k].m_connected_net_idx.data()[i], false); } // free all - no longer needed m_rails_temp.clear(); sort_terms(m_params.m_sort_type); this->set_pointers(); // create a list of non zero elements. for (unsigned k = 0; k < iN; k++) { terms_for_net_t & t = m_terms[k]; // pretty brutal int *other = t.m_connected_net_idx.data(); t.m_nz.clear(); for (std::size_t i = 0; i < t.railstart(); i++) if (!plib::container::contains(t.m_nz, static_cast(other[i]))) t.m_nz.push_back(static_cast(other[i])); t.m_nz.push_back(k); // add diagonal // and sort std::sort(t.m_nz.begin(), t.m_nz.end()); } // create a list of non zero elements right of the diagonal // These list anticipate the population of array elements by // Gaussian elimination. for (std::size_t k = 0; k < iN; k++) { terms_for_net_t & t = m_terms[k]; // pretty brutal int *other = t.m_connected_net_idx.data(); if (k==0) t.m_nzrd.clear(); else { t.m_nzrd = m_terms[k-1].m_nzrd; for (auto j = t.m_nzrd.begin(); j != t.m_nzrd.end(); ) { if (*j < k + 1) j = t.m_nzrd.erase(j); else ++j; } } for (std::size_t i = 0; i < t.railstart(); i++) if (!plib::container::contains(t.m_nzrd, static_cast(other[i])) && other[i] >= static_cast(k + 1)) t.m_nzrd.push_back(static_cast(other[i])); // and sort std::sort(t.m_nzrd.begin(), t.m_nzrd.end()); } // create a list of non zero elements below diagonal k // This should reduce cache misses ... std::vector> touched(iN, std::vector(iN)); for (std::size_t k = 0; k < iN; k++) { for (std::size_t j = 0; j < iN; j++) touched[k][j] = false; for (std::size_t j = 0; j < m_terms[k].m_nz.size(); j++) touched[k][m_terms[k].m_nz[j]] = true; } m_ops = 0; for (unsigned k = 0; k < iN; k++) { m_ops++; // 1/A(k,k) for (unsigned row = k + 1; row < iN; row++) { if (touched[row][k]) { m_ops++; if (!plib::container::contains(m_terms[k].m_nzbd, row)) m_terms[k].m_nzbd.push_back(row); for (std::size_t col = k + 1; col < iN; col++) if (touched[k][col]) { touched[row][col] = true; m_ops += 2; } } } } log().verbose("Number of mults/adds for {1}: {2}", name(), m_ops); if ((false)) for (std::size_t k = 0; k < iN; k++) { pstring line = plib::pfmt("{1:3}")(k); for (const auto & nzrd : m_terms[k].m_nzrd) line += plib::pfmt(" {1:3}")(nzrd); log().verbose("{1}", line); } // // save states // for (std::size_t k = 0; k < iN; k++) { pstring num = plib::pfmt("{1}")(k); state().save(*this, m_gonn[k],"GO" + num, this->name(), m_terms[k].count()); state().save(*this, m_gtn[k],"GT" + num, this->name(), m_terms[k].count()); state().save(*this, m_Idrn[k],"IDR" + num, this->name(), m_terms[k].count()); } } void matrix_solver_t::set_pointers() { const std::size_t iN = this->m_terms.size(); std::size_t max_count = 0; std::size_t max_rail = 0; for (std::size_t k = 0; k < iN; k++) { max_count = std::max(max_count, m_terms[k].count()); max_rail = std::max(max_rail, m_terms[k].railstart()); } m_gtn.resize(iN, max_count); m_gonn.resize(iN, max_count); m_Idrn.resize(iN, max_count); m_connected_net_Vn.resize(iN, max_count); for (std::size_t k = 0; k < iN; k++) { auto count = m_terms[k].count(); for (std::size_t i = 0; i < count; i++) { m_terms[k].terms()[i]->set_ptrs(&m_gtn[k][i], &m_gonn[k][i], &m_Idrn[k][i]); m_connected_net_Vn[k][i] = get_connected_net(m_terms[k].terms()[i])->Q_Analog_state_ptr(); } } } void matrix_solver_t::update_inputs() { // avoid recursive calls. Inputs are updated outside this call for (auto &inp : m_inps) inp->push(inp->proxied_net()->Q_Analog()); } bool matrix_solver_t::updates_net(const analog_net_t *net) const noexcept { if (net != nullptr) { for (const auto &t : m_terms ) if (t.is_net(net)) return true; for (const auto &inp : m_inps) if (&inp->net() == net) return true; } return false; } void matrix_solver_t::update_dynamic() noexcept { // update all non-linear devices for (auto &dyn : m_dynamic_funcs) dyn(); } void matrix_solver_t::reset() { m_last_step = netlist_time_ext::zero(); } void matrix_solver_t::update() noexcept { const netlist_time new_timestep = solve(exec().time()); update_inputs(); if (m_params.m_dynamic_ts && (timestep_device_count() != 0) && new_timestep > netlist_time::zero()) { m_Q_sync.net().toggle_and_push_to_queue(new_timestep); } } void matrix_solver_t::step(netlist_time delta) noexcept { const auto dd(delta.as_fp()); for (auto &d : m_step_funcs) d(dd); } netlist_time matrix_solver_t::solve(netlist_time_ext now) { const netlist_time_ext delta = now - m_last_step(); // We are already up to date. Avoid oscillations. // FIXME: Make this a parameter! if (delta < netlist_time_ext::quantum()) return netlist_time::zero(); // update all terminals for new time step m_last_step = now; step(static_cast(delta)); ++m_stat_vsolver_calls; if (dynamic_device_count() != 0) { bool this_resched(false); std::size_t newton_loops = 0; do { update_dynamic(); // Gauss-Seidel will revert to Gaussian elemination if steps exceeded. this->m_stat_calculations++; this->vsolve_non_dynamic(); this_resched = this->check_err(); this->store(); newton_loops++; } while (this_resched && newton_loops < m_params.m_nr_loops); m_stat_newton_raphson += newton_loops; // reschedule .... if (this_resched && !m_Q_sync.net().is_queued()) { log().warning(MW_NEWTON_LOOPS_EXCEEDED_ON_NET_1(this->name())); // FIXME: test and enable - this is working better, though not optimal yet #if 0 // Don't store, the result can not be used return netlist_time::from_fp(m_params.m_nr_recalc_delay()); #else m_Q_sync.net().toggle_and_push_to_queue(netlist_time::from_fp(m_params.m_nr_recalc_delay())); #endif } } else { this->m_stat_calculations++; this->vsolve_non_dynamic(); this->store(); } if (m_params.m_dynamic_ts) return compute_next_timestep(delta.as_fp(), m_params.m_max_timestep); return netlist_time::from_fp(m_params.m_max_timestep); } int matrix_solver_t::get_net_idx(const analog_net_t *net) const noexcept { for (std::size_t k = 0; k < m_terms.size(); k++) if (m_terms[k].is_net(net)) return static_cast(k); return -1; } std::pair matrix_solver_t::get_left_right_of_diag(std::size_t irow, std::size_t idiag) { // // return the maximum column left of the diagonal (-1 if no cols found) // return the minimum column right of the diagonal (999999 if no cols found) // const auto row = static_cast(irow); const auto diag = static_cast(idiag); int colmax = -1; int colmin = 999999; auto &term = m_terms[irow]; for (std::size_t i = 0; i < term.count(); i++) { auto col = get_net_idx(get_connected_net(term.terms()[i])); if (col != -1) { if (col==row) col = diag; else if (col==diag) col = row; if (col > diag && col < colmin) colmin = col; else if (col < diag && col > colmax) colmax = col; } } return {colmax, colmin}; } nl_fptype matrix_solver_t::get_weight_around_diag(std::size_t row, std::size_t diag) { { // // return average absolute distance // std::vector touched(1024, false); // FIXME! nl_fptype weight = nlconst::zero(); auto &term = m_terms[row]; for (std::size_t i = 0; i < term.count(); i++) { auto col = get_net_idx(get_connected_net(term.terms()[i])); if (col >= 0) { auto colu = static_cast(col); if (!touched[colu]) { if (colu==row) colu = static_cast(diag); else if (colu==diag) colu = static_cast(row); weight = weight + plib::abs(static_cast(colu) - static_cast(diag)); touched[colu] = true; } } } return weight; } } void matrix_solver_t::add_term(std::size_t net_idx, terminal_t *term) { if (get_connected_net(term)->is_rail_net()) { m_rails_temp[net_idx].add_terminal(term, -1, false); } else { int ot = get_net_idx(get_connected_net(term)); if (ot>=0) { m_terms[net_idx].add_terminal(term, ot, true); } else { log().fatal(MF_FOUND_TERM_WITH_MISSING_OTHERNET(term->name())); throw nl_exception(MF_FOUND_TERM_WITH_MISSING_OTHERNET(term->name())); } } } void matrix_solver_t::log_stats() { if (this->m_stat_calculations != 0 && this->m_stat_vsolver_calls && log().verbose.is_enabled()) { log().verbose("=============================================="); log().verbose("Solver {1}", this->name()); log().verbose(" ==> {1} nets", this->m_terms.size()); //, (*(*groups[i].first())->m_core_terms.first())->name()); log().verbose(" has {1} dynamic elements", this->dynamic_device_count()); log().verbose(" has {1} timestep elements", this->timestep_device_count()); log().verbose(" {1:6.3} average newton raphson loops", static_cast(this->m_stat_newton_raphson) / static_cast(this->m_stat_vsolver_calls)); log().verbose(" {1:10} invocations ({2:6.0} Hz) {3:10} gs fails ({4:6.2} %) {5:6.3} average", this->m_stat_calculations, static_cast(this->m_stat_calculations) / this->exec().time().as_fp(), this->m_iterative_fail, nlconst::magic(100.0) * static_cast(this->m_iterative_fail) / static_cast(this->m_stat_calculations), static_cast(this->m_iterative_total) / static_cast(this->m_stat_calculations)); } } } // namespace solver } // namespace netlist