diff options
Diffstat (limited to 'src/lib/netlist/solver/nld_matrix_solver.cpp')
-rw-r--r-- | src/lib/netlist/solver/nld_matrix_solver.cpp | 839 |
1 files changed, 486 insertions, 353 deletions
diff --git a/src/lib/netlist/solver/nld_matrix_solver.cpp b/src/lib/netlist/solver/nld_matrix_solver.cpp index 24acc8db576..d382da0f7d5 100644 --- a/src/lib/netlist/solver/nld_matrix_solver.cpp +++ b/src/lib/netlist/solver/nld_matrix_solver.cpp @@ -1,37 +1,46 @@ -// license:GPL-2.0+ +// license:BSD-3-Clause // copyright-holders:Couriersud -/* - * nld_matrix_solver.cpp - * - */ + +// Names +// spell-checker: words Raphson, Seidel +// +// Specific technical terms +// spell-checker: words vsolver #include "nld_matrix_solver.h" -#include "plib/putil.h" -#include <cmath> // <<= needed by windows build +#include "nl_setup.h" +#include "nld_solver.h" -namespace netlist -{ -namespace devices +#include "core/setup.h" + +#include "plib/putil.h" + +namespace netlist::solver { - terms_for_net_t::terms_for_net_t() - : m_railstart(0) - , m_last_V(0.0) - , m_DD_n_m_1(0.0) - , m_h_n_m_1(1e-12) + terms_for_net_t::terms_for_net_t(arena_type &arena, analog_net_t *net) + : m_nz(arena) + , m_nzrd(arena) + , m_nzbd(arena) + , m_connected_net_idx(arena) + , m_terms(arena) + , m_net(net) + , m_rail_start(0) { } - void terms_for_net_t::add(terminal_t *term, int net_other, bool sorted) + 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++) + 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); + plib::container::insert_at(m_connected_net_idx, i, + net_other); return; } } @@ -43,265 +52,326 @@ namespace devices // matrix_solver // ---------------------------------------------------------------------------------------- - matrix_solver_t::matrix_solver_t(netlist_state_t &anetlist, const pstring &name, - const eSortType sort, const solver_parameters_t *params) - : device_t(anetlist, name) - , m_params(*params) - , 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_iterative_fail(*this, "m_iterative_fail", 0) - , m_iterative_total(*this, "m_iterative_total", 0) - , m_last_step(*this, "m_last_step", netlist_time::zero()) - , m_fb_sync(*this, "FB_sync") - , m_Q_sync(*this, "Q_sync") - , m_ops(0) - , m_sort(sort) + matrix_solver_t::matrix_solver_t(devices::nld_solver &main_solver, + const pstring &name, + const net_list_t &nets, + const solver::solver_parameters_t *params) + //: device_t(static_cast<device_t &>(main_solver), name) + : device_t( + device_data_t{main_solver.state(), main_solver.name() + "." + name}) + , m_params(*params) + , m_gonn(m_arena) + , m_gtn(m_arena) + , m_Idrn(m_arena) + , m_connected_net_Vn(m_arena) + , m_iterative_fail(*this, "m_iterative_fail", 0) + , m_iterative_total(*this, "m_iterative_total", 0) + , m_main_solver(main_solver) + , m_stat_calculations(*this, "m_stat_calculations", 0) + , m_stat_newton_raphson(*this, "m_stat_newton_raphson", 0) + , m_stat_newton_raphson_fail(*this, "m_stat_newton_raphson_fail", 0) + , m_stat_vsolver_calls(*this, "m_stat_vsolver_calls", 0) + , m_last_step(*this, "m_last_step", netlist_time_ext::zero()) + , m_step_funcs(m_arena) + , m_dynamic_funcs(m_arena) + , m_inputs(m_arena) + , m_ops(0) + { + setup_base(this->state().setup(), nets); + + // now setup the matrix + setup_matrix(); + // printf("Freq: %f\n", m_params.m_freq()); + } + + analog_net_t *matrix_solver_t::get_connected_net(terminal_t *term) { - connect_post_start(m_fb_sync, m_Q_sync); + return &state().setup().get_connected_terminal(*term)->net(); } - void matrix_solver_t::setup_base(analog_net_t::list_t &nets) + void matrix_solver_t::reschedule(netlist_time ts) { + m_main_solver.reschedule(this, ts); + } + void matrix_solver_t::setup_base([[maybe_unused]] setup_t &setup, + const net_list_t &nets) + { log().debug("New solver setup\n"); + std::vector<core_device_t *> step_devices; + std::vector<core_device_t *> dynamic_devices; - m_nets.clear(); m_terms.clear(); - for (auto & net : nets) + for (const auto &net : nets) { - m_nets.push_back(net); - m_terms.push_back(plib::make_unique<terms_for_net_t>()); - m_rails_temp.push_back(plib::make_unique<terms_for_net_t>()); + m_terms.emplace_back(m_arena, net); + m_rails_temp.emplace_back(m_arena); } for (std::size_t k = 0; k < nets.size(); k++) { - analog_net_t *net = nets[k]; + std::vector<detail::core_terminal_t *> temp; + + analog_net_t &net = *nets[k]; - log().debug("setting up net\n"); + // FIXME: add size() to list + // log().debug("adding net with {1} populated connections\n", + // net.core_terms().size()); - net->set_solver(this); + net.set_solver(this); - for (auto &p : net->core_terms()) + for (detail::core_terminal_t *p : net.core_terms_copy()) { - log().debug("{1} {2} {3}\n", p->name(), net->name(), net->isRailNet()); + nl_assert_always(&p->net() == &net, "Net integrity violated"); + + 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(m_step_devices, &p->device())) - m_step_devices.push_back(&p->device()); + if (p->device().is_time_step()) + if (!plib::container::contains(step_devices, + &p->device())) + step_devices.push_back(&p->device()); if (p->device().is_dynamic()) - if (!plib::container::contains(m_dynamic_devices, &p->device())) - m_dynamic_devices.push_back(&p->device()); + if (!plib::container::contains(dynamic_devices, + &p->device())) + dynamic_devices.push_back(&p->device()); { - auto *pterm = dynamic_cast<terminal_t *>(p); - add_term(k, pterm); + auto pterm = plib::dynamic_downcast<terminal_t *>( + p); + nl_assert_always(bool(pterm), + "cast to terminal_t * failed"); + 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) + { + proxied_analog_output_t *net_proxy_output = nullptr; + for (auto &input : m_inputs) + if (input->proxied_net() == &p->net()) { - pstring nname = this->name() + "." + pstring(plib::pfmt("m{1}")(m_inps.size())); - nl_assert(p->net().is_analog()); - auto net_proxy_output_u = pool().make_poolptr<proxied_analog_output_t>(*this, nname, static_cast<analog_net_t *>(&p->net())); - net_proxy_output = net_proxy_output_u.get(); - m_inps.push_back(std::move(net_proxy_output_u)); + net_proxy_output = input.get(); + break; } - net_proxy_output->net().add_terminal(*p); - // FIXME: repeated calling - kind of brute force - net_proxy_output->net().rebuild_list(); - log().debug("Added input\n"); + + if (net_proxy_output == nullptr) + { + pstring new_name( + this->name() + "." + + pstring(plib::pfmt("m{1}")(m_inputs.size()))); + auto proxied_net = plib::dynamic_downcast< + analog_net_t *>(p->net()); + nl_assert_always(proxied_net, + "Net is not an analog net"); + auto net_proxy_output_u + = state() + .make_pool_object< + proxied_analog_output_t>( + *this, new_name, *proxied_net); + net_proxy_output = net_proxy_output_u.get(); + m_inputs.emplace_back( + std::move(net_proxy_output_u)); } - break; + net.remove_terminal(*p); + 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())); - break; + throw nl_exception( + MF_UNHANDLED_ELEMENT_1_FOUND(p->name())); } } - log().debug("added net with {1} populated connections\n", net->core_terms().size()); + net.rebuild_list(); } - - /* now setup the matrix */ - setup_matrix(); + for (auto &d : step_devices) + m_step_funcs.emplace_back( + nl_delegate_ts(&core_device_t::time_step, d)); + for (auto &d : dynamic_devices) + m_dynamic_funcs.emplace_back( + nl_delegate_dyn(&core_device_t::update_terminals, d)); } - void matrix_solver_t::sort_terms(eSortType sort) + /// \brief Sort terminals + /// + /// @param sort Sort algorithm to use. + /// + /// Sort in descending order by number of connected matrix voltages. + /// The idea is, that for Gauss-Seidel algorithm the first voltage computed + /// depends on the greatest number of previous voltages thus taking into + /// account the maximum amount 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. + /// FIXME: This entry needs an update. + /// + 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_nets.size(); + const std::size_t iN = m_terms.size(); switch (sort) { - case PREFER_BAND_MATRIX: + case matrix_sort_type_e::PREFER_BAND_MATRIX: + { + for (std::size_t k = 0; k < iN - 1; k++) { - for (std::size_t k = 0; k < iN - 1; k++) + auto pk = get_weight_around_diagonal(k, k); + for (std::size_t i = k + 1; i < iN; i++) { - auto pk = get_weight_around_diag(k,k); - for (std::size_t i = k+1; i < iN; i++) + auto pi = get_weight_around_diagonal(i, k); + if (pi < pk) { - auto pi = get_weight_around_diag(i,k); - if (pi < pk) - { - std::swap(m_terms[i], m_terms[k]); - std::swap(m_nets[i], m_nets[k]); - pk = get_weight_around_diag(k,k); - } + std::swap(m_terms[i], m_terms[k]); + pk = get_weight_around_diagonal(k, k); } } } - break; - case PREFER_IDENTITY_TOP_LEFT: + } + break; + case matrix_sort_type_e::PREFER_IDENTITY_TOP_LEFT: + { + for (std::size_t k = 0; k < iN - 1; k++) { - for (std::size_t k = 0; k < iN - 1; k++) + auto pk = get_left_right_of_diagonal(k, k); + for (std::size_t i = k + 1; i < iN; i++) { - 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_diagonal(i, k); + if (pi.first <= pk.first && pi.second >= pk.second) { - 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]); - std::swap(m_nets[i], m_nets[k]); - pk = get_left_right_of_diag(k,k); - } + std::swap(m_terms[i], m_terms[k]); + pk = get_left_right_of_diagonal(k, k); } } } - break; - case ASCENDING: - case DESCENDING: - { - int sort_order = (m_sort == DESCENDING ? 1 : -1); + } + 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++) + for (std::size_t k = 0; k < iN - 1; k++) + for (std::size_t i = k + 1; i < iN; i++) + { + if ((static_cast<int>(m_terms[k].rail_start()) + - static_cast<int>(m_terms[i].rail_start())) + * sort_order + < 0) { - if ((static_cast<int>(m_terms[k]->m_railstart) - static_cast<int>(m_terms[i]->m_railstart)) * sort_order < 0) - { - std::swap(m_terms[i], m_terms[k]); - std::swap(m_nets[i], m_nets[k]); - } + std::swap(m_terms[i], m_terms[k]); } - } - break; - case NOSORT: - break; + } + } + break; + case matrix_sort_type_e::NOSORT: break; } - /* rebuild */ + // 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 (other[i] != -1) - other[i] = get_net_idx(&term->terms()[i]->connected_terminal()->net()); + // 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_nets.size(); + const std::size_t iN = m_terms.size(); for (std::size_t k = 0; k < iN; k++) { - m_terms[k]->m_railstart = m_terms[k]->count(); - for (std::size_t i = 0; i < m_rails_temp[k]->count(); i++) - this->m_terms[k]->add(m_rails_temp[k]->terms()[i], m_rails_temp[k]->m_connected_net_idx.data()[i], false); + m_terms[k].set_rail_start(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_sort); + sort_terms(m_params.m_sort_type); this->set_pointers(); - /* create a list of non zero elements. */ + // create a list of non zero elements. for (unsigned k = 0; k < iN; k++) { - terms_for_net_t * t = m_terms[k].get(); - /* pretty brutal */ - int *other = t->m_connected_net_idx.data(); + terms_for_net_t &t = m_terms[k]; + // pretty brutal + int *other = t.m_connected_net_idx.data(); - t->m_nz.clear(); + t.m_nz.clear(); - for (std::size_t i = 0; i < t->m_railstart; i++) - if (!plib::container::contains(t->m_nz, static_cast<unsigned>(other[i]))) - t->m_nz.push_back(static_cast<unsigned>(other[i])); + for (std::size_t i = 0; i < t.rail_start(); i++) + if (!plib::container::contains(t.m_nz, + static_cast<unsigned>(other[i]))) + t.m_nz.push_back(static_cast<unsigned>(other[i])); - t->m_nz.push_back(k); // add diagonal + t.m_nz.push_back(k); // add diagonal - /* and sort */ - std::sort(t->m_nz.begin(), t->m_nz.end()); + // 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. - */ + // 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].get(); - /* pretty brutal */ - int *other = t->m_connected_net_idx.data(); + 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(); + 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(); ) + 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); + j = t.m_nzrd.erase(j); else ++j; } } - for (std::size_t i = 0; i < t->m_railstart; i++) - if (!plib::container::contains(t->m_nzrd, static_cast<unsigned>(other[i])) && other[i] >= static_cast<int>(k + 1)) - t->m_nzrd.push_back(static_cast<unsigned>(other[i])); + for (std::size_t i = 0; i < t.rail_start(); i++) + if (!plib::container::contains(t.m_nzrd, + static_cast<unsigned>(other[i])) + && other[i] >= static_cast<int>(k + 1)) + t.m_nzrd.push_back(static_cast<unsigned>(other[i])); - /* and sort */ - std::sort(t->m_nzrd.begin(), t->m_nzrd.end()); + // 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 ... - */ + // create a list of non zero elements below diagonal k + // This should reduce cache misses ... std::vector<std::vector<bool>> touched(iN, std::vector<bool>(iN)); @@ -309,8 +379,8 @@ namespace devices { 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; + 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; @@ -322,8 +392,8 @@ namespace devices if (touched[row][k]) { m_ops++; - if (!plib::container::contains(m_terms[k]->m_nzbd, row)) - m_terms[k]->m_nzbd.push_back(row); + 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]) { @@ -333,152 +403,246 @@ namespace devices } } } - log().verbose("Number of mults/adds for {1}: {2}", name(), m_ops); + log().verbose("Number of multiplications/additions for {1}: {2}", + name(), m_ops); + // Dumps non zero elements right of diagonal -> to much output, disabled + // NOLINTNEXTLINE(readability-simplify-boolean-expr) 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) + for (const auto &nzrd : m_terms[k].m_nzrd) line += plib::pfmt(" {1:3}")(nzrd); log().verbose("{1}", line); } - /* - * save states - */ + // + // save states + // + for (std::size_t k = 0; k < iN; k++) { pstring num = plib::pfmt("{1}")(k); - state().save(*this, m_terms[k]->m_last_V, this->name(), "lastV." + num); - state().save(*this, m_terms[k]->m_DD_n_m_1, this->name(), "m_DD_n_m_1." + num); - state().save(*this, m_terms[k]->m_h_n_m_1, this->name(), "m_h_n_m_1." + num); + 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].rail_start()); + } + + 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); - // FIXME: This shouldn't be necessary, recalculate on each entry ... - 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()); + // Initialize arrays to 0 (in case the vrl one is used + for (std::size_t k = 0; k < iN; k++) + for (std::size_t j = 0; j < m_terms[k].count(); j++) + { + m_gtn.set(k, j, nlconst::zero()); + m_gonn.set(k, j, nlconst::zero()); + m_Idrn.set(k, j, nlconst::zero()); + m_connected_net_Vn.set(k, j, nullptr); + } + + 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) + for (auto &inp : m_inputs) inp->push(inp->proxied_net()->Q_Analog()); } - void matrix_solver_t::update_dynamic() + void matrix_solver_t::update_dynamic() noexcept { - /* update all non-linear devices */ - for (auto &dyn : m_dynamic_devices) - dyn->update_terminals(); + // update all non-linear devices + for (auto &dyn : m_dynamic_funcs) + dyn(); } void matrix_solver_t::reset() { - m_last_step = netlist_time::zero(); + // m_last_step = netlist_time_ext::zero(); } - void matrix_solver_t::update() NL_NOEXCEPT + void matrix_solver_t::step(detail::time_step_type ts_type, + netlist_time delta) noexcept { - const netlist_time new_timestep = solve(exec().time()); - update_inputs(); + const auto dd(delta.as_fp<fptype>()); + for (auto &d : m_step_funcs) + d(ts_type, dd); + } - if (m_params.m_dynamic_ts && has_timestep_devices() && new_timestep > netlist_time::zero()) + bool matrix_solver_t::solve_nr_base() + { + bool this_resched(false); + std::size_t newton_loops = 0; + do { - m_Q_sync.net().toggle_and_push_to_queue(new_timestep); - } + update_dynamic(); + // Gauss-Seidel will revert to Gaussian elimination if steps + // exceeded. + this->m_stat_calculations++; + this->upstream_solve_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; + if (this_resched) + m_stat_newton_raphson_fail++; + return this_resched; } - /* update_forced is called from within param_update - * - * this should only occur outside of execution and thus - * using time should be safe. - * - */ - void matrix_solver_t::update_forced() + netlist_time matrix_solver_t::newton_loops_exceeded(netlist_time delta) { - const netlist_time new_timestep = solve(exec().time()); - plib::unused_var(new_timestep); + netlist_time next_time_step; + bool resched(false); - update_inputs(); + restore(); + step(detail::time_step_type::RESTORE, delta); - if (m_params.m_dynamic_ts && has_timestep_devices()) + for (std::size_t i = 0; i < 10; i++) { - m_Q_sync.net().toggle_and_push_to_queue(netlist_time::from_double(m_params.m_min_timestep)); + backup(); + step(detail::time_step_type::FORWARD, + netlist_time::from_fp(m_params.m_min_ts_ts())); + resched = solve_nr_base(); + // update time step calculation + next_time_step = compute_next_time_step(m_params.m_min_ts_ts(), + m_params.m_min_ts_ts(), + m_params.m_max_time_step); + delta -= netlist_time::from_fp(m_params.m_min_ts_ts()); + } + // try remaining time using compute_next_time step + while (delta > netlist_time::zero()) + { + if (next_time_step > delta) + next_time_step = delta; + backup(); + step(detail::time_step_type::FORWARD, next_time_step); + delta -= next_time_step; + resched = solve_nr_base(); + next_time_step = compute_next_time_step( + next_time_step.as_fp<nl_fptype>(), m_params.m_min_ts_ts(), + m_params.m_max_time_step); } - } - void matrix_solver_t::step(const netlist_time &delta) - { - const nl_double dd = delta.as_double(); - for (auto &d : m_step_devices) - d->timestep(dd); - } + if (m_stat_newton_raphson % 100 == 0) + log().warning(MW_NEWTON_LOOPS_EXCEEDED_INVOCATION_3( + 100, this->name(), exec().time().as_double() * 1e6)); - void matrix_solver_t::solve_base() - { - ++m_stat_vsolver_calls; - if (has_dynamic_devices()) + if (resched) { - std::size_t this_resched; - std::size_t newton_loops = 0; - do - { - update_dynamic(); - // Gauss-Seidel will revert to Gaussian elemination if steps exceeded. - this_resched = this->vsolve_non_dynamic(true); - newton_loops++; - } while (this_resched > 1 && newton_loops < m_params.m_nr_loops); - - m_stat_newton_raphson += newton_loops; // reschedule .... - if (this_resched > 1 && !m_Q_sync.net().is_queued()) - { - log().warning(MW_NEWTON_LOOPS_EXCEEDED_ON_NET_1(this->name())); - m_Q_sync.net().toggle_and_push_to_queue(m_params.m_nr_recalc_delay); - } - } - else - { - this->vsolve_non_dynamic(false); + log().warning(MW_NEWTON_LOOPS_EXCEEDED_ON_NET_2( + this->name(), exec().time().as_double() * 1e6)); + return netlist_time::from_fp(m_params.m_nr_recalc_delay()); } + if (m_params.m_dynamic_ts) + return next_time_step; + + return netlist_time::from_fp(m_params.m_max_time_step); } - const netlist_time matrix_solver_t::solve(netlist_time now) + netlist_time matrix_solver_t::solve(netlist_time_ext now, + [[maybe_unused]] const char *source) { - const netlist_time delta = now - m_last_step; + auto delta = static_cast<netlist_time>(now - m_last_step()); + PFDEBUG(printf("solve %.10f\n", delta.as_double());) // We are already up to date. Avoid oscillations. // FIXME: Make this a parameter! if (delta < netlist_time::quantum()) - return netlist_time::zero(); + { + // printf("solve return %s at %f\n", source, now.as_double()); + return time_step_device_count() > 0 + ? netlist_time::from_fp(m_params.m_min_time_step) + : netlist_time::zero(); + } - /* update all terminals for new time step */ + backup(); // save voltages for backup and time step calculation + // update all terminals for new time step m_last_step = now; - step(delta); - solve_base(); - const netlist_time next_time_step = compute_next_timestep(delta.as_double()); - return next_time_step; + ++m_stat_vsolver_calls; + if (dynamic_device_count() != 0) + { + step(detail::time_step_type::FORWARD, delta); + const auto resched = solve_nr_base(); + + if (resched) + return newton_loops_exceeded(delta); + } + else + { + step(detail::time_step_type::FORWARD, delta); + this->m_stat_calculations++; + this->upstream_solve_non_dynamic(); + this->store(); + } + + if (m_params.m_dynamic_ts) + { + if (time_step_device_count() > 0) + return compute_next_time_step(delta.as_fp<nl_fptype>(), + m_params.m_min_time_step, + m_params.m_max_time_step); + } + + if (time_step_device_count() > 0) + return netlist_time::from_fp(m_params.m_max_time_step); + + return netlist_time::zero(); } - int matrix_solver_t::get_net_idx(detail::net_t *net) + int matrix_solver_t::get_net_idx(const analog_net_t *net) const noexcept { - for (std::size_t k = 0; k < m_nets.size(); k++) - if (m_nets[k] == net) + for (std::size_t k = 0; k < m_terms.size(); k++) + if (m_terms[k].is_net(net)) return static_cast<int>(k); return -1; } - std::pair<int, int> matrix_solver_t::get_left_right_of_diag(std::size_t irow, std::size_t idiag) + std::pair<int, int> + matrix_solver_t::get_left_right_of_diagonal(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) - */ + // + // 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<int>(irow); const auto diag = static_cast<int>(idiag); @@ -488,13 +652,15 @@ namespace devices auto &term = m_terms[irow]; - for (std::size_t i = 0; i < term->count(); i++) + for (std::size_t i = 0; i < term.count(); i++) { - auto col = get_net_idx(&term->terms()[i]->connected_terminal()->net()); + 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 == row) + col = diag; + else if (col == diag) + col = row; if (col > diag && col < colmin) colmin = col; @@ -505,125 +671,92 @@ namespace devices return {colmax, colmin}; } - double matrix_solver_t::get_weight_around_diag(std::size_t row, std::size_t diag) + matrix_solver_t::fptype + matrix_solver_t::get_weight_around_diagonal(std::size_t row, + std::size_t diag) { { - /* - * return average absolute distance - */ + // + // return average absolute distance + // std::vector<bool> touched(1024, false); // FIXME! - double weight = 0.0; - auto &term = m_terms[row]; - for (std::size_t i = 0; i < term->count(); i++) + fptype weight = nlconst::zero(); + auto &term = m_terms[row]; + for (std::size_t i = 0; i < term.count(); i++) { - auto col = get_net_idx(&term->terms()[i]->connected_terminal()->net()); + auto col = get_net_idx(get_connected_net(term.terms()[i])); if (col >= 0) { auto colu = static_cast<std::size_t>(col); if (!touched[colu]) { - if (colu==row) colu = static_cast<unsigned>(diag); - else if (colu==diag) colu = static_cast<unsigned>(row); - - weight = weight + std::abs(static_cast<double>(colu) - static_cast<double>(diag)); + if (colu == row) + colu = static_cast<unsigned>(diag); + else if (colu == diag) + colu = static_cast<unsigned>(row); + + weight = weight + + plib::abs(static_cast<fptype>(colu) + - static_cast<fptype>(diag)); touched[colu] = true; } } } - return weight; // / static_cast<double>(term->m_railstart); + return weight; } } - void matrix_solver_t::add_term(std::size_t k, terminal_t *term) + void matrix_solver_t::add_term(std::size_t net_idx, terminal_t *term) { - if (term->connected_terminal()->net().isRailNet()) + if (get_connected_net(term)->is_rail_net()) { - m_rails_temp[k]->add(term, -1, false); + m_rails_temp[net_idx].add_terminal(term, -1, false); } else { - int ot = get_net_idx(&term->connected_terminal()->net()); - if (ot>=0) + int ot = get_net_idx(get_connected_net(term)); + if (ot >= 0) { - m_terms[k]->add(term, ot, true); + m_terms[net_idx].add_terminal(term, ot, true); } - /* Should this be allowed ? */ - else // if (ot<0) + else { - m_rails_temp[k]->add(term, ot, true); log().fatal(MF_FOUND_TERM_WITH_MISSING_OTHERNET(term->name())); + throw nl_exception( + MF_FOUND_TERM_WITH_MISSING_OTHERNET(term->name())); } } } - netlist_time matrix_solver_t::compute_next_timestep(const double cur_ts) - { - nl_double new_solver_timestep = m_params.m_max_timestep; - - if (m_params.m_dynamic_ts) - { - for (std::size_t k = 0, iN=m_terms.size(); k < iN; k++) - { - analog_net_t *n = m_nets[k]; - terms_for_net_t *t = m_terms[k].get(); - - //const nl_double DD_n = (n->Q_Analog() - t->m_last_V); - // avoid floating point exceptions - const nl_double DD_n = std::max(-1e100, std::min(1e100,(n->Q_Analog() - t->m_last_V))); - const nl_double hn = cur_ts; - - //printf("%g %g %g %g\n", DD_n, hn, t->m_DD_n_m_1, t->m_h_n_m_1); - nl_double DD2 = (DD_n / hn - t->m_DD_n_m_1 / t->m_h_n_m_1) / (hn + t->m_h_n_m_1); - nl_double new_net_timestep; - - t->m_h_n_m_1 = hn; - t->m_DD_n_m_1 = DD_n; - if (std::fabs(DD2) > plib::constants<nl_double>::cast(1e-60)) // avoid div-by-zero - new_net_timestep = std::sqrt(m_params.m_dynamic_lte / std::fabs(plib::constants<nl_double>::cast(0.5)*DD2)); - else - new_net_timestep = m_params.m_max_timestep; - - if (new_net_timestep < new_solver_timestep) - new_solver_timestep = new_net_timestep; - - t->m_last_V = n->Q_Analog(); - } - if (new_solver_timestep < m_params.m_min_timestep) - { - new_solver_timestep = m_params.m_min_timestep; - } - } - //if (new_solver_timestep > 10.0 * hn) - // new_solver_timestep = 10.0 * hn; - /* - * FIXME: Factor 2 below is important. Without, we get timing issues. This must be a bug elsewhere. - */ - return std::max(netlist_time::from_double(new_solver_timestep), netlist_time::quantum() * 2); - } - void matrix_solver_t::log_stats() { - if (this->m_stat_calculations != 0 && this->m_stat_vsolver_calls && log().verbose.is_enabled()) + 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_nets.size()); //, (*(*groups[i].first())->m_core_terms.first())->name()); - log().verbose(" has {1} elements", this->has_dynamic_devices() ? "dynamic" : "no dynamic"); - log().verbose(" has {1} elements", this->has_timestep_devices() ? "timestep" : "no timestep"); - log().verbose(" {1:6.3} average newton raphson loops", - static_cast<double>(this->m_stat_newton_raphson) / static_cast<double>(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<double>(this->m_stat_calculations) / this->exec().time().as_double(), - this->m_iterative_fail, - 100.0 * static_cast<double>(this->m_iterative_fail) - / static_cast<double>(this->m_stat_calculations), - static_cast<double>(this->m_iterative_total) / static_cast<double>(this->m_stat_calculations)); + log().verbose(" ==> {1} nets", this->m_terms.size()); + log().verbose(" has {1} dynamic elements", + this->dynamic_device_count()); + log().verbose(" has {1} time step elements", + this->time_step_device_count()); + log().verbose( + " {1:6.3} average newton raphson loops", + static_cast<fptype>(this->m_stat_newton_raphson) + / static_cast<fptype>(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<fptype>(this->m_stat_calculations) + / this->exec().time().as_fp<fptype>(), + this->m_iterative_fail, + nlconst::hundred() * static_cast<fptype>(this->m_iterative_fail) + / static_cast<fptype>(this->m_stat_calculations), + static_cast<fptype>(this->m_iterative_total) + / static_cast<fptype>(this->m_stat_calculations)); } } -} // namespace devices -} // namespace netlist - +} // namespace netlist::solver |