diff options
Diffstat (limited to 'src/lib/netlist/solver')
-rw-r--r-- | src/lib/netlist/solver/nld_matrix_solver.cpp | 839 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_matrix_solver.h | 659 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_matrix_solver_ext.h | 288 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_ms_direct.h | 217 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_ms_direct1.h | 41 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_ms_direct2.h | 58 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_ms_direct_lu.h | 170 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_ms_gcr.h | 453 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_ms_gmres.h | 152 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_ms_sm.h | 226 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_ms_sor.h | 249 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_ms_sor_mat.h | 165 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_ms_w.h | 585 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_solver.cpp | 762 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_solver.h | 143 |
15 files changed, 2781 insertions, 2226 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 diff --git a/src/lib/netlist/solver/nld_matrix_solver.h b/src/lib/netlist/solver/nld_matrix_solver.h index 29017add8a9..c0b792fd597 100644 --- a/src/lib/netlist/solver/nld_matrix_solver.h +++ b/src/lib/netlist/solver/nld_matrix_solver.h @@ -1,375 +1,464 @@ -// license:GPL-2.0+ +// license:BSD-3-Clause // copyright-holders:Couriersud -/* - * nld_matrix_solver.h - * - */ #ifndef NLD_MATRIX_SOLVER_H_ #define NLD_MATRIX_SOLVER_H_ -#include "netlist/nl_base.h" -#include "netlist/nl_errstr.h" +// Names +// spell-checker: words Raphson, Seidel + +/// +/// \file nld_matrix_solver.h +/// + +#include "nl_errstr.h" +#include "nltypes.h" + +#include "../core/analog.h" +#include "../core/device.h" +#include "../core/device_macros.h" +#include "../core/param.h" + #include "plib/palloc.h" +#include "plib/penum.h" #include "plib/pmatrix2d.h" +#include "plib/pmatrix_cr.h" +#include "plib/pmempool.h" #include "plib/putil.h" #include "plib/vector_ops.h" -#include <cmath> +#include <numeric> -namespace netlist -{ -namespace devices +// FIXME: remove again + +#define PFDEBUG(x) + +namespace netlist::solver { - /* FIXME: these should become proper devices */ - struct solver_parameters_t + enum static_compile_target { - bool m_pivot; - nl_double m_accuracy; - nl_double m_dynamic_lte; - nl_double m_min_timestep; - nl_double m_max_timestep; - nl_double m_gs_sor; - bool m_dynamic_ts; - std::size_t m_gs_loops; - std::size_t m_nr_loops; - netlist_time m_nr_recalc_delay; - bool m_use_gabs; - bool m_use_linear_prediction; + CXX_EXTERNAL_C, + CXX_STATIC }; - - class terms_for_net_t : plib::nocopyassignmove + // clang-format off + + PENUM(matrix_sort_type_e, + NOSORT, + ASCENDING, + DESCENDING, + PREFER_IDENTITY_TOP_LEFT, + PREFER_BAND_MATRIX + ) + + PENUM(matrix_type_e, + SOR_MAT, + MAT_CR, + MAT, + SM, + W, + SOR, + GMRES + ) + + PENUM(matrix_fp_type_e, + FLOAT + , DOUBLE + , LONGDOUBLE + , FLOATQ128 + ) + + // clang-format on + + using arena_type = plib::mempool_arena<plib::aligned_arena<>, 1024>; + using static_compile_container = std::vector<std::pair<pstring, pstring>>; + + struct solver_parameter_defaults { - public: - terms_for_net_t(); - - void clear(); - - void add(terminal_t *term, int net_other, bool sorted); - - std::size_t count() const { return m_terms.size(); } - - terminal_t **terms() { return m_terms.data(); } + static constexpr nl_fptype m_freq() { return nlconst::magic(48000.0); } - std::size_t m_railstart; + // iteration parameters + static constexpr nl_fptype m_gs_sor() { return nlconst::magic(1.059); } + static constexpr matrix_type_e m_method() + { + return matrix_type_e::MAT_CR; + } + static constexpr matrix_fp_type_e m_fp_type() + { + return matrix_fp_type_e::DOUBLE; + } + static constexpr nl_fptype m_reltol() { return nlconst::magic(1e-3); } + static constexpr nl_fptype m_vntol() { return nlconst::magic(1e-7); } + static constexpr nl_fptype m_accuracy() { return nlconst::magic(1e-7); } + static constexpr std::size_t m_nr_loops() { return 250; } + static constexpr std::size_t m_gs_loops() { return 50; } + + // general parameters + static constexpr nl_fptype m_gmin() { return nlconst::magic(1e-9); } + static constexpr bool m_pivot() { return false; } + static constexpr nl_fptype m_nr_recalc_delay() + { + return netlist_time::quantum().as_fp<nl_fptype>(); + } + static constexpr int m_parallel() { return 0; } - std::vector<unsigned> m_nz; /* all non zero for multiplication */ - std::vector<unsigned> m_nzrd; /* non zero right of the diagonal for elimination, may include RHS element */ - std::vector<unsigned> m_nzbd; /* non zero below of the diagonal for elimination */ + static constexpr nl_fptype m_min_ts_ts() + { + return nlconst::magic(1e-9); + } + // automatic time step + static constexpr bool m_dynamic_ts() { return false; } + static constexpr nl_fptype m_dynamic_lte() + { + return nlconst::magic(1e-5); + } + static constexpr nl_fptype m_dynamic_min_ts() + { + return nlconst::magic(1e-6); + } - /* state */ - nl_double m_last_V; - nl_double m_DD_n_m_1; - nl_double m_h_n_m_1; + // matrix sorting + static constexpr matrix_sort_type_e m_sort_type() + { + return matrix_sort_type_e::PREFER_IDENTITY_TOP_LEFT; + } - std::vector<int> m_connected_net_idx; - private: - std::vector<terminal_t *> m_terms; + // special + static constexpr bool m_use_gabs() { return true; } + static solver_parameter_defaults &get_instance() + { + static solver_parameter_defaults s; + return s; + } }; - class proxied_analog_output_t : public analog_output_t + struct solver_parameters_t { - public: + template <typename D> + solver_parameters_t(device_t &parent, const pstring &prefix, + D &defaults) + : m_freq(parent, prefix + "FREQ", defaults.m_freq()) + + // iteration parameters + , m_gs_sor(parent, prefix + "SOR_FACTOR", defaults.m_gs_sor()) + , m_method(parent, prefix + "METHOD", defaults.m_method()) + , m_fp_type(parent, prefix + "FPTYPE", defaults.m_fp_type()) + , m_reltol(parent, prefix + "RELTOL", + defaults.m_reltol()) //!< SPICE RELTOL parameter + , m_vntol(parent, prefix + "VNTOL", defaults.m_vntol()) //!< SPICE VNTOL + //!< parameter + , m_accuracy(parent, prefix + "ACCURACY", + defaults.m_accuracy()) //!< Iterative solver accuracy + , m_nr_loops(parent, prefix + "NR_LOOPS", + defaults.m_nr_loops()) //!< Maximum number of + //!< Newton-Raphson loops + , m_gs_loops(parent, prefix + "GS_LOOPS", + defaults.m_gs_loops()) //!< Maximum number of Gauss-Seidel + //!< loops + + // general parameters + , m_gmin(parent, prefix + "GMIN", defaults.m_gmin()) + , m_pivot(parent, prefix + "PIVOT", defaults.m_pivot()) //!< use + //!< pivoting on + //!< supported + //!< solvers + , m_nr_recalc_delay(parent, prefix + "NR_RECALC_DELAY", + defaults.m_nr_recalc_delay()) //!< Delay to next + //!< solve attempt if + //!< nr loops exceeded + , m_parallel(parent, prefix + "PARALLEL", defaults.m_parallel()) + , m_min_ts_ts(parent, prefix + "MIN_TS_TS", + defaults.m_min_ts_ts()) //!< The minimum time step for + //!< solvers with time stepping + //!< devices. + + // automatic time step + , m_dynamic_ts(parent, prefix + "DYNAMIC_TS", + defaults.m_dynamic_ts()) //!< Use dynamic time stepping + , m_dynamic_lte(parent, prefix + "DYNAMIC_LTE", + defaults.m_dynamic_lte()) //!< dynamic time stepping + //!< slope + , m_dynamic_min_ts(parent, prefix + "DYNAMIC_MIN_TIMESTEP", + defaults.m_dynamic_min_ts()) //!< smallest time step + //!< allowed + + // matrix sorting + , m_sort_type(parent, prefix + "SORT_TYPE", defaults.m_sort_type()) + + // special + , m_use_gabs(parent, prefix + "USE_GABS", defaults.m_use_gabs()) + , m_min_time_step(m_dynamic_min_ts()) + { + m_max_time_step = netlist_time::from_fp(plib::reciprocal(m_freq())) + .as_fp<decltype(m_max_time_step)>(); - proxied_analog_output_t(core_device_t &dev, const pstring &aname, analog_net_t *pnet) - : analog_output_t(dev, aname) - , m_proxied_net(pnet) - { } + if (m_dynamic_ts) + { + m_max_time_step *= 1; // NL_FCONST(1000.0); + } + else + { + m_min_time_step = m_max_time_step; + } + } - analog_net_t *proxied_net() const { return m_proxied_net;} - private: - analog_net_t *m_proxied_net; // only for proxy nets in analog input logic + param_fp_t m_freq; + param_fp_t m_gs_sor; + param_enum_t<matrix_type_e> m_method; + param_enum_t<matrix_fp_type_e> m_fp_type; + param_fp_t m_reltol; + param_fp_t m_vntol; + param_fp_t m_accuracy; + param_num_t<std::size_t> m_nr_loops; + param_num_t<std::size_t> m_gs_loops; + param_fp_t m_gmin; + param_logic_t m_pivot; + param_fp_t m_nr_recalc_delay; + param_int_t m_parallel; + param_fp_t m_min_ts_ts; + param_logic_t m_dynamic_ts; + param_fp_t m_dynamic_lte; + param_fp_t m_dynamic_min_ts; + param_enum_t<matrix_sort_type_e> m_sort_type; + + param_logic_t m_use_gabs; + + nl_fptype m_min_time_step; + nl_fptype m_max_time_step; }; - class matrix_solver_t : public device_t + class terms_for_net_t { public: - using list_t = std::vector<matrix_solver_t *>; + terms_for_net_t(arena_type &arena, analog_net_t *net = nullptr); - enum eSortType - { - NOSORT, - ASCENDING, - DESCENDING, - PREFER_IDENTITY_TOP_LEFT, - PREFER_BAND_MATRIX - }; - - void setup(analog_net_t::list_t &nets) - { - vsetup(nets); - } - - void solve_base(); + void clear(); - /* after every call to solve, update inputs must be called. - * this can be done as well as a batch to ease parallel processing. - */ - const netlist_time solve(netlist_time now); - void update_inputs(); + void add_terminal(terminal_t *term, int net_other, bool sorted); - bool has_dynamic_devices() const { return m_dynamic_devices.size() > 0; } - bool has_timestep_devices() const { return m_step_devices.size() > 0; } + std::size_t count() const noexcept { return m_terms.size(); } - void update_forced(); - void update_after(const netlist_time after) - { - m_Q_sync.net().toggle_and_push_to_queue(after); - } + std::size_t rail_start() const noexcept { return m_rail_start; } - /* netdevice functions */ - NETLIB_UPDATEI(); - NETLIB_RESETI(); + terminal_t **terms() noexcept { return m_terms.data(); } - public: - int get_net_idx(detail::net_t *net); - std::pair<int, int> get_left_right_of_diag(std::size_t row, std::size_t diag); - double get_weight_around_diag(std::size_t row, std::size_t diag); + nl_fptype getV() const noexcept { return m_net->Q_Analog(); } - virtual void log_stats(); + void setV(nl_fptype v) noexcept { m_net->set_Q_Analog(v); } - virtual std::pair<pstring, pstring> create_solver_code() + bool is_net(const analog_net_t *net) const noexcept { - return std::pair<pstring, pstring>("", plib::pfmt("/* solver doesn't support static compile */\n\n")); + return net == m_net; } - /* return number of floating point operations for solve */ - std::size_t ops() { return m_ops; } + void set_rail_start(std::size_t val) noexcept { m_rail_start = val; } - protected: + PALIGNAS_VECTOROPT() + + plib::arena_vector<arena_type, unsigned> m_nz; //!< all non zero for + //!< multiplication + plib::arena_vector<arena_type, unsigned> m_nzrd; //!< non zero right of + //!< the diagonal for + //!< elimination, may + //!< include RHS + //!< element + plib::arena_vector<arena_type, unsigned> m_nzbd; //!< non zero below of + //!< the diagonal for + //!< elimination - matrix_solver_t(netlist_state_t &anetlist, const pstring &name, - eSortType sort, const solver_parameters_t *params); + plib::arena_vector<arena_type, int> m_connected_net_idx; - void sort_terms(eSortType sort); + private: + plib::arena_vector<arena_type, terminal_t *> m_terms; + analog_net_t *m_net; + std::size_t m_rail_start; + }; - void setup_base(analog_net_t::list_t &nets); - void update_dynamic(); + class proxied_analog_output_t : public analog_output_t + { + public: + proxied_analog_output_t(core_device_t &dev, const pstring &aname, + analog_net_t *pnet) + : analog_output_t(dev, aname) + , m_proxied_net(pnet) + { + } - virtual void vsetup(analog_net_t::list_t &nets) = 0; - virtual unsigned vsolve_non_dynamic(const bool newton_raphson) = 0; + analog_net_t *proxied_net() const { return m_proxied_net; } - netlist_time compute_next_timestep(const double cur_ts); - /* virtual */ void add_term(std::size_t net_idx, terminal_t *term); + private: + analog_net_t *m_proxied_net; // only for proxy nets in analog input + // logic + }; - template <typename T> - void store(const T & V); + class matrix_solver_t : public device_t + { + public: + using list_t = std::vector<matrix_solver_t *>; + using fptype = nl_fptype; + using net_list_t = std::vector<analog_net_t *>; - template <typename T> - auto delta(const T & V) -> typename std::decay<decltype(V[0])>::type; + // after every call to solve, update inputs must be called. + // this can be done as well as a batch to ease parallel processing. - template <typename T> - void build_LE_A(T &child); - template <typename T> - void build_LE_RHS(T &child); + netlist_time solve(netlist_time_ext now, const char *source); + void update_inputs(); - void set_pointers() + std::size_t dynamic_device_count() const noexcept { - const std::size_t iN = this->m_nets.size(); + return m_dynamic_funcs.size(); + } + std::size_t time_step_device_count() const noexcept + { + return m_step_funcs.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]->m_railstart); - } + /// \brief reschedule solver execution + /// + /// Calls reschedule on main solver + /// + void reschedule(netlist_time ts); + + /// \brief Immediately solve system at current time + /// + /// This should only be called from update and update_param events. + /// It's purpose is to bring voltage values to the current time step. + /// This will be called BEFORE updating object properties. + void solve_now() + { + // this should only occur outside of execution and thus + // using time should be safe. - m_mat_ptr.resize(iN, max_rail+1); - 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); + [[maybe_unused]] const netlist_time new_time_step = solve( + exec().time(), "solve_now"); - for (std::size_t k = 0; k < iN; k++) - { - auto count = m_terms[k]->count(); + update_inputs(); - 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] = m_terms[k]->terms()[i]->connected_terminal()->net().Q_Analog_state_ptr(); - } + if (time_step_device_count() > 0) + { + this->reschedule(netlist_time::from_fp( + m_params.m_dynamic_ts ? m_params.m_min_time_step + : m_params.m_max_time_step)); } } - template <typename AP, typename FT> - void fill_matrix(std::size_t N, AP &tcr, FT &RHS) + template <typename F> + void change_state(F f) { - for (std::size_t k = 0; k < N; k++) + // We only need to update the net first if this is a time stepping + // net + if (time_step_device_count() > 0) { - auto *net = m_terms[k].get(); - auto **tcr_r = &(tcr[k][0]); - - const std::size_t term_count = net->count(); - const std::size_t railstart = net->m_railstart; - const auto &go = m_gonn[k]; - const auto > = m_gtn[k]; - const auto &Idr = m_Idrn[k]; - const auto &cnV = m_connected_net_Vn[k]; - - for (std::size_t i = 0; i < railstart; i++) - *tcr_r[i] += go[i]; - - typename FT::value_type gtot_t = 0.0; - typename FT::value_type RHS_t = 0.0; - - for (std::size_t i = 0; i < term_count; i++) - { - gtot_t += gt[i]; - RHS_t += Idr[i]; - } - // FIXME: Code above is faster than vec_sum - Check this - #if 0 - auto gtot_t = plib::vec_sum<FT>(term_count, m_gt); - auto RHS_t = plib::vec_sum<FT>(term_count, m_Idr); - #endif - - for (std::size_t i = railstart; i < term_count; i++) - { - RHS_t += (/*m_Idr[i]*/ (- go[i]) * *cnV[i]); - } - - RHS[k] = RHS_t; - // update diagonal element ... - *tcr_r[railstart] += gtot_t; //mat.A[mat.diag[k]] += gtot_t; + [[maybe_unused]] const netlist_time new_time_step = solve( + exec().time(), "change_state"); + update_inputs(); } - + f(); + if (time_step_device_count() > 0) + { + PFDEBUG(printf("here2\n");) + this->reschedule(netlist_time::from_fp(m_params.m_min_ts_ts())); + } + else + this->reschedule(netlist_time::quantum()); } - template <typename T> - using aligned_alloc = plib::aligned_allocator<T, PALIGN_VECTOROPT>; + NETLIB_RESETI(); - plib::pmatrix2d<nl_double, aligned_alloc<nl_double>> m_gonn; - plib::pmatrix2d<nl_double, aligned_alloc<nl_double>> m_gtn; - plib::pmatrix2d<nl_double, aligned_alloc<nl_double>> m_Idrn; - plib::pmatrix2d<nl_double *, aligned_alloc<nl_double *>> m_mat_ptr; - plib::pmatrix2d<nl_double *, aligned_alloc<nl_double *>> m_connected_net_Vn; + virtual void log_stats(); - plib::pmatrix2d<nl_double> m_test; + virtual std::pair<pstring, pstring> create_solver_code( + [[maybe_unused]] solver::static_compile_target target) + { + return {"", + plib::pfmt("// solver doesn't support static compile\n\n")}; + } - std::vector<plib::unique_ptr<terms_for_net_t>> m_terms; - std::vector<analog_net_t *> m_nets; - std::vector<pool_owned_ptr<proxied_analog_output_t>> m_inps; + // return number of floating point operations for solve + constexpr std::size_t ops() const { return m_ops; } - std::vector<plib::unique_ptr<terms_for_net_t>> m_rails_temp; + protected: + matrix_solver_t(devices::nld_solver &main_solver, const pstring &name, + const net_list_t &nets, + const solver_parameters_t *params); + + virtual void upstream_solve_non_dynamic() = 0; + virtual netlist_time + compute_next_time_step(fptype cur_ts, fptype min_ts, fptype max_ts) + = 0; + virtual bool check_err() const = 0; + virtual void store() = 0; + virtual void backup() = 0; + virtual void restore() = 0; + + std::size_t max_rail_start() const noexcept + { + std::size_t max_rail = 0; + for (const auto &term : m_terms) + max_rail = std::max(max_rail, term.rail_start()); + return max_rail; + } const solver_parameters_t &m_params; + arena_type m_arena; - state_var<int> m_stat_calculations; - state_var<int> m_stat_newton_raphson; - state_var<int> m_stat_vsolver_calls; - state_var<int> m_iterative_fail; - state_var<int> m_iterative_total; - - private: + plib::pmatrix2d_vrl<arena_type, fptype> m_gonn; + plib::pmatrix2d_vrl<arena_type, fptype> m_gtn; + plib::pmatrix2d_vrl<arena_type, fptype> m_Idrn; + plib::pmatrix2d_vrl<arena_type, fptype *> m_connected_net_Vn; - state_var<netlist_time> m_last_step; - std::vector<core_device_t *> m_step_devices; - std::vector<core_device_t *> m_dynamic_devices; + state_var<std::size_t> m_iterative_fail; + state_var<std::size_t> m_iterative_total; - logic_input_t m_fb_sync; - logic_output_t m_Q_sync; + std::vector<terms_for_net_t> m_terms; // setup only - /* calculate matrix */ - void setup_matrix(); - - void step(const netlist_time &delta); - - std::size_t m_ops; - const eSortType m_sort; - }; - - template <typename T> - auto matrix_solver_t::delta(const T & V) -> typename std::decay<decltype(V[0])>::type - { - /* NOTE: Ideally we should also include currents (RHS) here. This would - * need a reevaluation of the right hand side after voltages have been updated - * and thus belong into a different calculation. This applies to all solvers. - */ - - const std::size_t iN = this->m_terms.size(); - typename std::decay<decltype(V[0])>::type cerr = 0; - for (std::size_t i = 0; i < iN; i++) - cerr = std::max(cerr, std::abs(V[i] - this->m_nets[i]->Q_Analog())); - return cerr; - } - - template <typename T> - void matrix_solver_t::store(const T & V) - { - const std::size_t iN = this->m_terms.size(); - for (std::size_t i = 0; i < iN; i++) - this->m_nets[i]->set_Q_Analog(V[i]); - } + private: + // base setup - called from constructor + void setup_base(setup_t &setup, const net_list_t &nets) noexcept(false); - template <typename T> - void matrix_solver_t::build_LE_A(T &child) - { - using float_type = typename T::float_type; - static_assert(std::is_base_of<matrix_solver_t, T>::value, "T must derive from matrix_solver_t"); + bool solve_nr_base(); + netlist_time newton_loops_exceeded(netlist_time delta); - const std::size_t iN = child.size(); - for (std::size_t k = 0; k < iN; k++) - { - terms_for_net_t *terms = m_terms[k].get(); - float_type * Ak = &child.A(k, 0ul); + void sort_terms(matrix_sort_type_e sort); - for (std::size_t i=0; i < iN; i++) - Ak[i] = 0.0; + void update_dynamic() noexcept; + void step(detail::time_step_type ts_type, netlist_time delta) noexcept; - const std::size_t terms_count = terms->count(); - const std::size_t railstart = terms->m_railstart; - const float_type * const gt = m_gtn[k]; + int get_net_idx(const analog_net_t *net) const noexcept; + std::pair<int, int> + get_left_right_of_diagonal(std::size_t irow, std::size_t idiag); + fptype get_weight_around_diagonal(std::size_t row, std::size_t diag); - { - float_type akk = 0.0; - for (std::size_t i = 0; i < terms_count; i++) - akk += gt[i]; - - Ak[k] = akk; - } + void add_term(std::size_t net_idx, terminal_t *term) noexcept(false); - const float_type * const go = m_gonn[k]; - int * net_other = terms->m_connected_net_idx.data(); + // calculate matrix + void setup_matrix(); - for (std::size_t i = 0; i < railstart; i++) - Ak[net_other[i]] += go[i]; - } - } + void set_pointers(); - template <typename T> - void matrix_solver_t::build_LE_RHS(T &child) - { - static_assert(std::is_base_of<matrix_solver_t, T>::value, "T must derive from matrix_solver_t"); - using float_type = typename T::float_type; + analog_net_t *get_connected_net(terminal_t *term); - const std::size_t iN = child.size(); - for (std::size_t k = 0; k < iN; k++) - { - float_type rhsk_a = 0.0; - float_type rhsk_b = 0.0; + devices::nld_solver &m_main_solver; - const std::size_t terms_count = m_terms[k]->count(); - const float_type * const go = m_gonn[k]; - const float_type * const Idr = m_Idrn[k]; - const float_type * const * other_cur_analog = m_connected_net_Vn[k]; + state_var<std::size_t> m_stat_calculations; + state_var<std::size_t> m_stat_newton_raphson; + state_var<std::size_t> m_stat_newton_raphson_fail; + state_var<std::size_t> m_stat_vsolver_calls; - for (std::size_t i = 0; i < terms_count; i++) - rhsk_a = rhsk_a + Idr[i]; + state_var<netlist_time_ext> m_last_step; + plib::arena_vector<arena_type, nl_delegate_ts> m_step_funcs; + plib::arena_vector<arena_type, nl_delegate_dyn> m_dynamic_funcs; + plib::arena_vector<arena_type, + device_arena::unique_ptr<proxied_analog_output_t>> + m_inputs; - for (std::size_t i = m_terms[k]->m_railstart; i < terms_count; i++) - //rhsk = rhsk + go[i] * terms[i]->m_otherterm->net().as_analog().Q_Analog(); - rhsk_b = rhsk_b - go[i] * *other_cur_analog[i]; + std::size_t m_ops; - child.RHS(k) = rhsk_a + rhsk_b; - } - } + std::vector<terms_for_net_t> m_rails_temp; // setup only + }; -} //namespace devices -} // namespace netlist +} // namespace netlist::solver -#endif /* NLD_MS_DIRECT_H_ */ +#endif // NLD_MS_DIRECT_H_ diff --git a/src/lib/netlist/solver/nld_matrix_solver_ext.h b/src/lib/netlist/solver/nld_matrix_solver_ext.h new file mode 100644 index 00000000000..3627c0b8e91 --- /dev/null +++ b/src/lib/netlist/solver/nld_matrix_solver_ext.h @@ -0,0 +1,288 @@ +// license:BSD-3-Clause +// copyright-holders:Couriersud + +#ifndef NLD_MATRIX_SOLVER_EXT_H_ +#define NLD_MATRIX_SOLVER_EXT_H_ + +/// +/// \file nld_matrix_solver.h +/// + +#include "nld_matrix_solver.h" + +#include <numeric> + +namespace netlist::solver +{ + + template <typename FT, int SIZE> + class matrix_solver_ext_t: public matrix_solver_t + { + public: + + using float_type = FT; + + matrix_solver_ext_t(devices::nld_solver &main_solver, const pstring &name, + const net_list_t &nets, + const solver::solver_parameters_t *params, const std::size_t size) + : matrix_solver_t(main_solver, name, nets, params) + , m_new_V(size) + , m_RHS(size) + , m_mat_ptr(m_arena, size, this->max_rail_start() + 1) + , m_last_V(size, nlconst::zero()) + , m_DD_n_m_1(size, nlconst::zero()) + , m_h_n_m_1(size, nlconst::magic(1e-6)) // we need a non zero value here + , m_dim(size) + { + // + // save states + // + state().save(*this, m_last_V.as_base(), this->name(), "m_last_V"); + state().save(*this, m_DD_n_m_1.as_base(), this->name(), "m_DD_n_m_1"); + state().save(*this, m_h_n_m_1.as_base(), this->name(), "m_h_n_m_1"); + } + + protected: + static constexpr const std::size_t SIZEABS = plib::parray<FT, SIZE>::SIZEABS(); + static constexpr const std::size_t m_pitch_ABS = (((SIZEABS + 0) + 7) / 8) * 8; + + //PALIGNAS_VECTOROPT() `parray` defines alignment already + plib::parray<float_type, SIZE> m_new_V; + //PALIGNAS_VECTOROPT() `parray` defines alignment already + plib::parray<float_type, SIZE> m_RHS; + + //PALIGNAS_VECTOROPT() `parray` defines alignment already + plib::pmatrix2d<arena_type, float_type *> m_mat_ptr; + + template <typename T, typename M> + void log_fill(const T &fill, [[maybe_unused]] M &mat) + { + const std::size_t iN = fill.size(); + + // FIXME: Not yet working, mat_cr.h needs some more work +#if 0 + auto mat_GE = plib::dynamic_downcast<plib::pGEmatrix_cr_t<typename M::base> *>(&mat); +#endif + std::vector<unsigned> levL(iN, 0); + std::vector<unsigned> levU(iN, 0); + + // parallel scheme for L x = y + for (std::size_t k = 0; k < iN; k++) + { + unsigned lm=0; + for (std::size_t j = 0; j<k; j++) + if (fill[k][j] < M::FILL_INFINITY) + lm = std::max(lm, levL[j]); + levL[k] = 1+lm; + } + + // parallel scheme for U x = y + for (std::size_t k = iN; k-- > 0; ) + { + unsigned lm=0; + for (std::size_t j = iN; --j > k; ) + if (fill[k][j] < M::FILL_INFINITY) + lm = std::max(lm, levU[j]); + levU[k] = 1+lm; + } + for (std::size_t k = 0; k < iN; k++) + { + unsigned fm = 0; + pstring ml = ""; + for (std::size_t j = 0; j < iN; j++) + { + ml += fill[k][j] == 0 ? 'X' : fill[k][j] < M::FILL_INFINITY ? '+' : '.'; + if (fill[k][j] < M::FILL_INFINITY) + if (fill[k][j] > fm) + fm = fill[k][j]; + } +#if 0 + this->log().verbose("{1:4} {2} {3:4} {4:4} {5:4} {6:4}", k, ml, + levL[k], levU[k], mat_GE ? mat_GE->get_parallel_level(k) : 0, fm); +#else + this->log().verbose("{1:4} {2} {3:4} {4:4} {5:4} {6:4}", k, ml, + levL[k], levU[k], 0, fm); +#endif + } + } + + constexpr std::size_t size() const noexcept + { + return (SIZE > 0) ? static_cast<std::size_t>(SIZE) : m_dim; + } + +#if 1 + void store() override + { + const std::size_t iN = size(); + for (std::size_t i = 0; i < iN; i++) + this->m_terms[i].setV(static_cast<fptype>(m_new_V[i])); + } +#else + // global tanh damping (4.197) + // partially cures the symptoms but not the cause + void store() override + { + const std::size_t iN = size(); + for (std::size_t i = 0; i < iN; i++) + { + auto oldV = this->m_terms[i].template getV<fptype>(); + this->m_terms[i].setV(oldV + 0.02 * plib::tanh((m_new_V[i]-oldV)*50.0)); + } + } +#endif + bool check_err() const override + { + // NOTE: Ideally we should also include currents (RHS) here. This would + // need a reevaluation of the right hand side after voltages have been updated + // and thus belong into a different calculation. This applies to all solvers. + + const std::size_t iN = size(); + const float_type reltol(static_cast<float_type>(m_params.m_reltol)); + const float_type vntol(static_cast<float_type>(m_params.m_vntol)); + for (std::size_t i = 0; i < iN; i++) + { + const float_type vold(static_cast<float_type>(this->m_terms[i].getV())); + const float_type vnew(m_new_V[i]); + const float_type tol(vntol + reltol * std::max(plib::abs(vnew),plib::abs(vold))); + if (plib::abs(vnew - vold) > tol) + return true; + } + return false; + } + + void backup() override + { + const std::size_t iN = size(); + for (std::size_t i = 0; i < iN; i++) + m_last_V[i] = gsl::narrow_cast<fptype>(this->m_terms[i].getV()); + } + + void restore() override + { + const std::size_t iN = size(); + for (std::size_t i = 0; i < iN; i++) + this->m_terms[i].setV(static_cast<nl_fptype>(m_last_V[i])); + } + + netlist_time compute_next_time_step(fptype cur_ts, fptype min_ts, fptype max_ts) override + { + fptype new_solver_time_step_sq(max_ts * max_ts); + + for (std::size_t k = 0; k < size(); k++) + { + const auto &t = m_terms[k]; + const auto v(static_cast<fptype>(t.getV())); + // avoid floating point exceptions + const fptype DD_n = std::max(-fp_constants<fptype>::TIMESTEP_MAXDIFF(), + std::min(+fp_constants<fptype>::TIMESTEP_MAXDIFF(),(v - m_last_V[k]))); + + //m_last_V[k] = v; + const fptype hn = cur_ts; + + fptype DD2 = (DD_n / hn - m_DD_n_m_1[k] / m_h_n_m_1[k]) / (hn + m_h_n_m_1[k]); + + m_h_n_m_1[k] = hn; + m_DD_n_m_1[k] = DD_n; + if (plib::abs(DD2) > fp_constants<fptype>::TIMESTEP_MINDIV()) // avoid div-by-zero + { + // save the sqrt for the end + const fptype new_net_time_step_sq = m_params.m_dynamic_lte / plib::abs(nlconst::half()*DD2); + new_solver_time_step_sq = std::min(new_net_time_step_sq, new_solver_time_step_sq); + } + } + + new_solver_time_step_sq = std::max(plib::sqrt(new_solver_time_step_sq), min_ts); + + // FIXME: Factor 2 below is important. Without, we get timing issues. This must be a bug elsewhere. + return std::max(netlist_time::from_fp(new_solver_time_step_sq), netlist_time::quantum() * 2); + } + + template <typename M> + void build_mat_ptr(M &mat) + { + const std::size_t iN = size(); + + for (std::size_t k=0; k<iN; k++) + { + 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].rail_start();j++) + { + int other = this->m_terms[k].m_connected_net_idx[j]; + if (other >= 0) + { + m_mat_ptr[k][j] = &(mat[k][static_cast<std::size_t>(other)]); + cnt++; + } + } + nl_assert_always(cnt == this->m_terms[k].rail_start(), "Count and rail start mismatch"); + m_mat_ptr[k][this->m_terms[k].rail_start()] = &(mat[k][k]); + } + } + + template <typename M> + void clear_square_mat(M &m) + { + const std::size_t n = size(); + for (std::size_t k=0; k < n; k++) + { + auto *p = &(m[k][0]); + using mat_elem_type = typename std::decay<decltype(*p)>::type; + for (std::size_t i=0; i < n; i++) + p[i] = plib::constants<mat_elem_type>::zero(); + } + } + + void fill_matrix_and_rhs() + { + const std::size_t N = size(); + + for (std::size_t k = 0; k < N; k++) + { + auto &net = m_terms[k]; + auto **tcr_r = &(m_mat_ptr[k][0]); + + using source_type = typename decltype(m_gtn)::value_type; + const std::size_t term_count = net.count(); + const std::size_t rail_start = net.rail_start(); + const auto &go = m_gonn[k]; + const auto > = m_gtn[k]; + const auto &Idr = m_Idrn[k]; + const auto &cnV = m_connected_net_Vn[k]; + + //# FIXME: gonn, gtn and Idr - which float types should they have? + + auto gtot_t = std::accumulate(gt, gt + term_count, plib::constants<source_type>::zero()); + + // update diagonal element ... + *tcr_r[rail_start] = static_cast<FT>(gtot_t); //# mat.A[mat.diag[k]] += gtot_t; + + for (std::size_t i = 0; i < rail_start; i++) + *tcr_r[i] += static_cast<FT>(go[i]); + + auto RHS_t = std::accumulate(Idr, Idr + term_count, plib::constants<source_type>::zero()); + + for (std::size_t i = rail_start; i < term_count; i++) + RHS_t += (- go[i]) * *cnV[i]; + + m_RHS[k] = static_cast<FT>(RHS_t); + } + } + + private: + // state - variable time_stepping + //PALIGNAS_VECTOROPT() `parray` defines alignment already + plib::parray<fptype, SIZE> m_last_V; + //PALIGNAS_VECTOROPT() `parray` defines alignment already + plib::parray<fptype, SIZE> m_DD_n_m_1; + // PALIGNAS_VECTOROPT() parrays define alignment already + plib::parray<fptype, SIZE> m_h_n_m_1; + + const std::size_t m_dim; + + }; + +} // namespace netlist::solver + +#endif // NLD_MATRIX_SOLVER_EXT_H_ diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h index 2501742218d..a3e9b6345b5 100644 --- a/src/lib/netlist/solver/nld_ms_direct.h +++ b/src/lib/netlist/solver/nld_ms_direct.h @@ -1,63 +1,56 @@ -// license:GPL-2.0+ +// license:BSD-3-Clause // copyright-holders:Couriersud -/* - * nld_ms_direct.h - * - */ #ifndef NLD_MS_DIRECT_H_ #define NLD_MS_DIRECT_H_ +/// +/// \file nld_ms_direct.h +/// + #include "nld_matrix_solver.h" #include "nld_solver.h" -#include "plib/mat_cr.h" +#include "plib/parray.h" #include "plib/vector_ops.h" +#include "nld_matrix_solver_ext.h" + #include <algorithm> -#include <cmath> -namespace netlist -{ -namespace devices +namespace netlist::solver { template <typename FT, int SIZE> - class matrix_solver_direct_t: public matrix_solver_t + class matrix_solver_direct_t: public matrix_solver_ext_t<FT, SIZE> { - friend class matrix_solver_t; public: using float_type = FT; - matrix_solver_direct_t(netlist_state_t &anetlist, const pstring &name, const solver_parameters_t *params, const std::size_t size); - matrix_solver_direct_t(netlist_state_t &anetlist, const pstring &name, const eSortType sort, const solver_parameters_t *params, const std::size_t size); + matrix_solver_direct_t(devices::nld_solver &main_solver, const pstring &name, + const matrix_solver_t::net_list_t &nets, + const solver::solver_parameters_t *params, std::size_t size); - void vsetup(analog_net_t::list_t &nets) override; void reset() override { matrix_solver_t::reset(); } + private: + + const std::size_t m_pitch; + protected: - unsigned vsolve_non_dynamic(const bool newton_raphson) override; - unsigned solve_non_dynamic(const bool newton_raphson); + static constexpr const std::size_t SIZEABS = plib::parray<FT, SIZE>::SIZEABS(); + static constexpr const std::size_t m_pitch_ABS = (((SIZEABS + 0) + 7) / 8) * 8; - constexpr std::size_t size() const { return (SIZE > 0) ? static_cast<std::size_t>(SIZE) : m_dim; } + void upstream_solve_non_dynamic() override; + void solve_non_dynamic(); void LE_solve(); template <typename T> void LE_back_subst(T & x); - FT &A(std::size_t r, std::size_t c) { return m_A[r * m_pitch + c]; } - FT &RHS(std::size_t r) { return m_A[r * m_pitch + size()]; } - plib::parray<FT, SIZE> m_new_V; - - private: - static constexpr const std::size_t SIZEABS = plib::parray<FT, SIZE>::SIZEABS(); - static constexpr const std::size_t m_pitch_ABS = (((SIZEABS + 1) + 7) / 8) * 8; - - const std::size_t m_dim; - const std::size_t m_pitch; - plib::parray<FT, SIZE * int(m_pitch_ABS)> m_A; - + // PALIGNAS_VECTOROPT() `parray` defines alignment already + plib::parray2D<FT, SIZE, m_pitch_ABS> m_A; }; // ---------------------------------------------------------------------------------------- @@ -65,43 +58,26 @@ namespace devices // ---------------------------------------------------------------------------------------- template <typename FT, int SIZE> - void matrix_solver_direct_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets) - { - matrix_solver_t::setup_base(nets); - - /* add RHS element */ - for (std::size_t k = 0; k < size(); k++) - { - terms_for_net_t * t = m_terms[k].get(); - - if (!plib::container::contains(t->m_nzrd, static_cast<unsigned>(size()))) - t->m_nzrd.push_back(static_cast<unsigned>(size())); - } - - // FIXME: This shouldn't be necessary ... - for (std::size_t k = 0; k < size(); k++) - state().save(*this, RHS(k), this->name(), plib::pfmt("RHS.{1}")(k)); - } - - template <typename FT, int SIZE> void matrix_solver_direct_t<FT, SIZE>::LE_solve() { - const std::size_t kN = size(); - if (!m_params.m_pivot) + const std::size_t kN = this->size(); + if (!this->m_params.m_pivot) { for (std::size_t i = 0; i < kN; i++) { - /* FIXME: Singular matrix? */ - const FT f = 1.0 / A(i,i); - const auto &nzrd = m_terms[i]->m_nzrd; - const auto &nzbd = m_terms[i]->m_nzbd; + // FIXME: Singular matrix? + const auto &Ai = m_A[i]; + const FT f = plib::reciprocal(Ai[i]); + const auto &nzrd = this->m_terms[i].m_nzrd; + const auto &nzbd = this->m_terms[i].m_nzbd; - for (std::size_t j : nzbd) + for (auto &j : nzbd) { - const FT f1 = -f * A(j, i); - for (std::size_t k : nzrd) - A(j, k) += A(i, k) * f1; - //RHS(j) += RHS(i) * f1; + auto &Aj = m_A[j]; + const FT f1 = -f * Aj[i]; + for (auto &k : nzrd) + Aj[k] += Ai[k] * f1; + this->m_RHS[j] += this->m_RHS[i] * f1; } } } @@ -109,45 +85,47 @@ namespace devices { for (std::size_t i = 0; i < kN; i++) { - /* Find the row with the largest first value */ - std::size_t maxrow = i; + // Find the row with the largest first value + std::size_t max_row = i; for (std::size_t j = i + 1; j < kN; j++) { - //if (std::abs(m_A[j][i]) > std::abs(m_A[maxrow][i])) - if (A(j,i) * A(j,i) > A(maxrow,i) * A(maxrow,i)) - maxrow = j; + if (plib::abs(m_A[j][i]) > plib::abs(m_A[max_row][i])) + //#if (m_A[j][i] * m_A[j][i] > m_A[max_row][i] * m_A[max_row][i]) + max_row = j; } - if (maxrow != i) + if (max_row != i) { - /* Swap the maxrow and ith row */ - for (std::size_t k = 0; k < kN + 1; k++) { - std::swap(A(i,k), A(maxrow,k)); +#if 0 + // Swap the max_row and ith row + for (std::size_t k = 0; k < kN; k++) { + std::swap(m_A[i][k], m_A[max_row][k]); } - //std::swap(RHS(i), RHS(maxrow)); +#else + std::swap(m_A[i], m_A[max_row]); +#endif + std::swap(this->m_RHS[i], this->m_RHS[max_row]); } - /* FIXME: Singular matrix? */ - const FT f = 1.0 / A(i,i); + // FIXME: Singular matrix? + const auto &Ai = m_A[i]; + const FT f = plib::reciprocal(Ai[i]); - /* Eliminate column i from row j */ + // Eliminate column i from row j for (std::size_t j = i + 1; j < kN; j++) { - const FT f1 = - A(j,i) * f; + auto &Aj = m_A[j]; + const FT f1 = - m_A[j][i] * f; if (f1 != plib::constants<FT>::zero()) { - const FT * pi = &A(i,i+1); - FT * pj = &A(j,i+1); - #if 1 - plib::vec_add_mult_scalar_p(kN-i,pj, pi,f1); - #else - vec_add_mult_scalar_p1(kN-i-1,pj,pi,f1); + const FT * pi = &(Ai[i+1]); + FT * pj = &(Aj[i+1]); + plib::vec_add_mult_scalar_p(kN-i-1,pj,pi,f1); //for (unsigned k = i+1; k < kN; k++) // pj[k] = pj[k] + pi[k] * f1; //for (unsigned k = i+1; k < kN; k++) //A(j,k) += A(i,k) * f1; - RHS(j) += RHS(i) * f1; - #endif + this->m_RHS[j] += this->m_RHS[i] * f1; } } } @@ -159,77 +137,66 @@ namespace devices void matrix_solver_direct_t<FT, SIZE>::LE_back_subst( T & x) { - const std::size_t kN = size(); + const std::size_t kN = this->size(); - /* back substitution */ - if (m_params.m_pivot) + // back substitution + if (this->m_params.m_pivot) { for (std::size_t j = kN; j-- > 0; ) { - FT tmp = 0; + FT tmp(0); + const auto & Aj(m_A[j]); + for (std::size_t k = j+1; k < kN; k++) - tmp += A(j,k) * x[k]; - x[j] = (RHS(j) - tmp) / A(j,j); + tmp += Aj[k] * x[k]; + x[j] = (this->m_RHS[j] - tmp) / Aj[j]; } } else { for (std::size_t j = kN; j-- > 0; ) { - FT tmp = 0; - const auto &nzrd = m_terms[j]->m_nzrd; - const auto e = nzrd.size() - 1; /* exclude RHS element */ + FT tmp(0); + const auto &nzrd = this->m_terms[j].m_nzrd; + const auto & Aj(m_A[j]); + const auto e = nzrd.size(); + for ( std::size_t k = 0; k < e; k++) - tmp += A(j, nzrd[k]) * x[nzrd[k]]; - x[j] = (RHS(j) - tmp) / A(j,j); + tmp += Aj[nzrd[k]] * x[nzrd[k]]; + x[j] = (this->m_RHS[j] - tmp) / Aj[j]; } } } template <typename FT, int SIZE> - unsigned matrix_solver_direct_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphson) + void matrix_solver_direct_t<FT, SIZE>::solve_non_dynamic() { this->LE_solve(); - this->LE_back_subst(m_new_V); - - const FT err = (newton_raphson ? delta(m_new_V) : 0.0); - store(m_new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; + this->LE_back_subst(this->m_new_V); } template <typename FT, int SIZE> - unsigned matrix_solver_direct_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson) + void matrix_solver_direct_t<FT, SIZE>::upstream_solve_non_dynamic() { - this->build_LE_A(*this); - this->build_LE_RHS(*this); + // populate matrix + this->clear_square_mat(m_A); + this->fill_matrix_and_rhs(); - this->m_stat_calculations++; - return this->solve_non_dynamic(newton_raphson); - } - - template <typename FT, int SIZE> - matrix_solver_direct_t<FT, SIZE>::matrix_solver_direct_t(netlist_state_t &anetlist, const pstring &name, - const solver_parameters_t *params, const std::size_t size) - : matrix_solver_t(anetlist, name, ASCENDING, params) - , m_new_V(size) - , m_dim(size) - , m_pitch(m_pitch_ABS ? m_pitch_ABS : (((m_dim + 1) + 7) / 8) * 8) - , m_A(size * m_pitch) - { + this->solve_non_dynamic(); } template <typename FT, int SIZE> - matrix_solver_direct_t<FT, SIZE>::matrix_solver_direct_t(netlist_state_t &anetlist, const pstring &name, - const eSortType sort, const solver_parameters_t *params, const std::size_t size) - : matrix_solver_t(anetlist, name, sort, params) - , m_new_V(size) - , m_dim(size) - , m_pitch(m_pitch_ABS ? m_pitch_ABS : (((m_dim + 1) + 7) / 8) * 8) - , m_A(size * m_pitch) + matrix_solver_direct_t<FT, SIZE>::matrix_solver_direct_t(devices::nld_solver &main_solver, const pstring &name, + const matrix_solver_t::net_list_t &nets, + const solver::solver_parameters_t *params, + std::size_t size) + : matrix_solver_ext_t<FT, SIZE>(main_solver, name, nets, params, size) + , m_pitch(m_pitch_ABS ? m_pitch_ABS : (((size + 0) + 7) / 8) * 8) + , m_A(size, m_pitch) { + this->build_mat_ptr(m_A); } -} // namespace devices -} // namespace netlist +} // namespace netlist::solver -#endif /* NLD_MS_DIRECT_H_ */ +#endif // NLD_MS_DIRECT_H_ diff --git a/src/lib/netlist/solver/nld_ms_direct1.h b/src/lib/netlist/solver/nld_ms_direct1.h index fbbb8ecb098..899dcdf7b91 100644 --- a/src/lib/netlist/solver/nld_ms_direct1.h +++ b/src/lib/netlist/solver/nld_ms_direct1.h @@ -1,19 +1,18 @@ -// license:GPL-2.0+ +// license:BSD-3-Clause // copyright-holders:Couriersud -/* - * nld_ms_direct1.h - * - */ #ifndef NLD_MS_DIRECT1_H_ #define NLD_MS_DIRECT1_H_ +/// +/// \file nld_ms_direct1.h +/// + +#include "nld_matrix_solver_ext.h" #include "nld_ms_direct.h" #include "nld_solver.h" -namespace netlist -{ -namespace devices +namespace netlist::solver { template <typename FT> class matrix_solver_direct1_t: public matrix_solver_direct_t<FT, 1> @@ -23,32 +22,26 @@ namespace devices using float_type = FT; using base_type = matrix_solver_direct_t<FT, 1>; - matrix_solver_direct1_t(netlist_state_t &anetlist, const pstring &name, const solver_parameters_t *params) - : matrix_solver_direct_t<FT, 1>(anetlist, name, params, 1) + matrix_solver_direct1_t(devices::nld_solver &main_solver, const pstring &name, + const matrix_solver_t::net_list_t &nets, + const solver::solver_parameters_t *params) + : matrix_solver_direct_t<FT, 1>(main_solver, name, nets, params, 1) {} // ---------------------------------------------------------------------------------------- // matrix_solver - Direct1 // ---------------------------------------------------------------------------------------- - unsigned vsolve_non_dynamic(const bool newton_raphson) override + void upstream_solve_non_dynamic() override { - this->build_LE_A(*this); - this->build_LE_RHS(*this); - //NL_VERBOSE_OUT(("{1} {2}\n", new_val, m_RHS[0] / m_A[0][0]); - - std::array<FT, 1> new_V = { this->RHS(0) / this->A(0,0) }; + this->clear_square_mat(this->m_A); + this->fill_matrix_and_rhs(); - const FT err = (newton_raphson ? this->delta(new_V) : 0.0); - this->store(new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; + this->m_new_V[0] = this->m_RHS[0] / this->m_A[0][0]; } - }; - -} //namespace devices -} // namespace netlist +} // namespace netlist::solver -#endif /* NLD_MS_DIRECT1_H_ */ +#endif // NLD_MS_DIRECT1_H_ diff --git a/src/lib/netlist/solver/nld_ms_direct2.h b/src/lib/netlist/solver/nld_ms_direct2.h index 01f77c3bc3c..0f4702ce84b 100644 --- a/src/lib/netlist/solver/nld_ms_direct2.h +++ b/src/lib/netlist/solver/nld_ms_direct2.h @@ -1,19 +1,18 @@ -// license:GPL-2.0+ +// license:BSD-3-Clause // copyright-holders:Couriersud -/* - * nld_ms_direct1.h - * - */ #ifndef NLD_MS_DIRECT2_H_ #define NLD_MS_DIRECT2_H_ +/// +/// \file nld_ms_direct2.h +/// + +#include "nld_matrix_solver_ext.h" #include "nld_ms_direct.h" #include "nld_solver.h" -namespace netlist -{ -namespace devices +namespace netlist::solver { // ---------------------------------------------------------------------------------------- @@ -27,32 +26,29 @@ namespace devices using float_type = FT; - matrix_solver_direct2_t(netlist_state_t &anetlist, const pstring &name, const solver_parameters_t *params) - : matrix_solver_direct_t<double, 2>(anetlist, name, params, 2) - {} - unsigned vsolve_non_dynamic(const bool newton_raphson) override + matrix_solver_direct2_t(devices::nld_solver &main_solver, const pstring &name, + const matrix_solver_t::net_list_t &nets, + const solver::solver_parameters_t *params) + : matrix_solver_direct_t<FT, 2>(main_solver, name, nets, params, 2) + {} + void upstream_solve_non_dynamic() override { - this->build_LE_A(*this); - this->build_LE_RHS(*this); - - const float_type a = this->A(0,0); - const float_type b = this->A(0,1); - const float_type c = this->A(1,0); - const float_type d = this->A(1,1); - - const float_type v1 = (a * this->RHS(1) - c * this->RHS(0)) / (a * d - b * c); - const float_type v0 = (this->RHS(0) - b * v1) / a; - std::array<float_type, 2> new_V = {v0, v1}; - - this->m_stat_calculations++; - const float_type err = (newton_raphson ? this->delta(new_V) : 0.0); - this->store(new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; + this->clear_square_mat(this->m_A); + this->fill_matrix_and_rhs(); + + const float_type a = this->m_A[0][0]; + const float_type b = this->m_A[0][1]; + const float_type c = this->m_A[1][0]; + const float_type d = this->m_A[1][1]; + + const float_type v1 = (a * this->m_RHS[1] - c * this->m_RHS[0]) / (a * d - b * c); + const float_type v0 = (this->m_RHS[0] - b * v1) / a; + this->m_new_V[0] = v0; + this->m_new_V[1] = v1; } }; -} //namespace devices -} // namespace netlist +} // namespace netlist::solver -#endif /* NLD_MS_DIRECT2_H_ */ +#endif // NLD_MS_DIRECT2_H_ diff --git a/src/lib/netlist/solver/nld_ms_direct_lu.h b/src/lib/netlist/solver/nld_ms_direct_lu.h index e7cedc1dd29..c51ccac4da5 100644 --- a/src/lib/netlist/solver/nld_ms_direct_lu.h +++ b/src/lib/netlist/solver/nld_ms_direct_lu.h @@ -1,13 +1,18 @@ -// license:GPL-2.0+ +// license:BSD-3-Clause // copyright-holders:Couriersud -/* - * nld_ms_direct.h - * - */ + +/// +/// \file nld_ms_direct.h +/// +/// #if 0 #ifndef NLD_MS_DIRECT_H_ #define NLD_MS_DIRECT_H_ +// Names +// spell-checker: words Seidel,Crout + + #include "solver/nld_solver.h" #include "solver/nld_matrix_solver.h" @@ -38,12 +43,12 @@ public: unsigned N() const { if (m_N == 0) return m_dim; else return m_N; } - int vsolve_non_dynamic(const bool newton_raphson); + int upstream_solve_non_dynamic(bool newton_raphson); protected: virtual void add_term(int net_idx, terminal_t *term) override; - int solve_non_dynamic(const bool newton_raphson); + int solve_non_dynamic(bool newton_raphson); void build_LE_A(); void build_LE_RHS(nl_double * RESTRICT rhs); @@ -132,11 +137,11 @@ protected: nl_double delta(const nl_double * RESTRICT V); void store(const nl_double * RESTRICT V); - /* bring the whole system to the current time - * Don't schedule a new calculation time. The recalculation has to be - * triggered by the caller after the netlist element was changed. - */ - nl_double compute_next_timestep(); + // bring the whole system to the current time + // Don't schedule a new calculation time. The recalculation has to be + // triggered by the caller after the netlist element was changed. + + nl_double compute_next_time_step(); template <typename T1, typename T2> nl_ext_double &A(const T1 r, const T2 c) { return m_A[r][c]; } @@ -164,65 +169,64 @@ matrix_solver_direct_t<m_N, storage_N>::~matrix_solver_direct_t() } template <unsigned m_N, unsigned storage_N> -nl_double matrix_solver_direct_t<m_N, storage_N>::compute_next_timestep() +nl_double matrix_solver_direct_t<m_N, storage_N>::compute_next_time_step() { - nl_double new_solver_timestep = m_params.m_max_timestep; + nl_double new_solver_time_step = m_params.m_max_time_step; if (m_params.m_dynamic_ts) { - /* - * FIXME: We should extend the logic to use either all nets or - * only output nets. - */ + // + // FIXME: We should extend the logic to use either all nets or + // only output nets. for (unsigned k = 0, iN=N(); k < iN; k++) { analog_net_t *n = m_nets[k]; const nl_double DD_n = (n->Q_Analog() - m_last_V[k]); - const nl_double hn = current_timestep(); + const nl_double hn = current_time_step(); nl_double DD2 = (DD_n / hn - n->m_DD_n_m_1 / n->m_h_n_m_1) / (hn + n->m_h_n_m_1); - nl_double new_net_timestep; + nl_double new_net_time_step; n->m_h_n_m_1 = hn; n->m_DD_n_m_1 = DD_n; - if (std::abs(DD2) > NL_FCONST(1e-30)) // avoid div-by-zero - new_net_timestep = std::sqrt(m_params.m_dynamic_lte / std::abs(NL_FCONST(0.5)*DD2)); + if (plib::abs(DD2) > NL_FCONST(1e-30)) // avoid div-by-zero + new_net_time_step = std::sqrt(m_params.m_dynamic_lte / plib::abs(NL_FCONST(0.5)*DD2)); else - new_net_timestep = m_params.m_max_timestep; + new_net_time_step = m_params.m_max_time_step; - if (new_net_timestep < new_solver_timestep) - new_solver_timestep = new_net_timestep; + if (new_net_time_step < new_solver_time_step) + new_solver_time_step = new_net_time_step; } - if (new_solver_timestep < m_params.m_min_timestep) - new_solver_timestep = m_params.m_min_timestep; - if (new_solver_timestep > m_params.m_max_timestep) - new_solver_timestep = m_params.m_max_timestep; + if (new_solver_time_step < m_params.m_min_time_step) + new_solver_time_step = m_params.m_min_time_step; + if (new_solver_time_step > m_params.m_max_time_step) + new_solver_time_step = m_params.m_max_time_step; } - //if (new_solver_timestep > 10.0 * hn) - // new_solver_timestep = 10.0 * hn; - return new_solver_timestep; + //#if (new_solver_time_step > 10.0 * hn) + //# new_solver_time_step = 10.0 * hn; + return new_solver_time_step; } template <unsigned m_N, unsigned storage_N> void matrix_solver_direct_t<m_N, storage_N>::add_term(int k, terminal_t *term) { - if (term->m_otherterm->net().isRailNet()) + if (term->m_other_terminal->net().isRailNet()) { m_rails_temp[k].add(term, -1, false); } else { - int ot = get_net_idx(&term->m_otherterm->net()); + int ot = get_net_idx(&term->m_other_terminal->net()); if (ot>=0) { m_terms[k]->add(term, ot, true); } - /* Should this be allowed ? */ + // Should this be allowed ? else // if (ot<0) { m_rails_temp[k].add(term, ot, true); - netlist().error("found term with missing othernet {1}\n", term->name()); + netlist().error("found term with missing other net {1}\n", term->name()); } } } @@ -244,7 +248,7 @@ void matrix_solver_direct_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets) for (unsigned k = 0; k < N(); k++) { - m_terms[k]->m_railstart = m_terms[k]->count(); + m_terms[k]->m_rail_start = m_terms[k]->count(); for (unsigned i = 0; i < m_rails_temp[k].count(); i++) this->m_terms[k]->add(m_rails_temp[k].terms()[i], m_rails_temp[k].connected_net_idx()[i], false); @@ -254,32 +258,32 @@ void matrix_solver_direct_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets) #if 1 - /* 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. - * FIXME: Even better would be to sort on elements right of the matrix diagonal. - * - */ + // 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 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. + // FIXME: Even better would be to sort on elements right of the matrix diagonal. + // + // int sort_order = (type() == GAUSS_SEIDEL ? 1 : -1); for (unsigned k = 0; k < N() / 2; k++) for (unsigned i = 0; i < N() - 1; i++) { - if ((m_terms[i]->m_railstart - m_terms[i+1]->m_railstart) * sort_order < 0) + if ((m_terms[i]->m_rail_start - m_terms[i+1]->m_rail_start) * sort_order < 0) { std::swap(m_terms[i],m_terms[i+1]); std::swap(m_nets[i], m_nets[i+1]); @@ -291,19 +295,18 @@ void matrix_solver_direct_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets) int *other = m_terms[k]->connected_net_idx(); for (unsigned i = 0; i < m_terms[k]->count(); i++) if (other[i] != -1) - other[i] = get_net_idx(&m_terms[k]->terms()[i]->m_otherterm->net()); + other[i] = get_net_idx(&m_terms[k]->terms()[i]->m_other_terminal->net()); } #endif - /* 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 (unsigned k = 0; k < N(); k++) { terms_for_net_t * t = m_terms[k]; - /* pretty brutal */ int *other = t->connected_net_idx(); t->m_nz.clear(); @@ -325,7 +328,7 @@ void matrix_solver_direct_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets) for (unsigned j = 0; j < N(); j++) { - for (unsigned i = 0; i < t->m_railstart; i++) + for (unsigned i = 0; i < t->m_rail_start; i++) { if (!t->m_nzrd.contains(other[i]) && other[i] >= (int) (k + 1)) t->m_nzrd.add(other[i]); @@ -349,9 +352,10 @@ void matrix_solver_direct_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets) log("\n"); } - /* - * save states - */ + // + // save states + // + save(NLNAME(m_RHS)); save(NLNAME(m_last_V)); @@ -378,7 +382,7 @@ void matrix_solver_direct_t<m_N, storage_N>::build_LE_A() nl_double akk = 0.0; const unsigned terms_count = m_terms[k]->count(); - const unsigned railstart = m_terms[k]->m_railstart; + const unsigned rail_start = m_terms[k]->m_rail_start; const nl_double * RESTRICT gt = m_terms[k]->gt(); const nl_double * RESTRICT go = m_terms[k]->go(); const int * RESTRICT net_other = m_terms[k]->connected_net_idx(); @@ -388,7 +392,7 @@ void matrix_solver_direct_t<m_N, storage_N>::build_LE_A() A(k,k) += akk; - for (unsigned i = 0; i < railstart; i++) + for (unsigned i = 0; i < rail_start; i++) A(k, net_other[i]) -= go[i]; } } @@ -410,8 +414,8 @@ void matrix_solver_direct_t<m_N, storage_N>::build_LE_RHS(nl_double * RESTRICT r for (int i = 0; i < terms_count; i++) rhsk_a = rhsk_a + Idr[i]; - for (int i = m_terms[k]->m_railstart; i < terms_count; i++) - //rhsk = rhsk + go[i] * terms[i]->m_otherterm->net().as_analog().Q_Analog(); + for (int i = m_terms[k]->m_rail_start; i < terms_count; i++) + //#rhsk = rhsk + go[i] * terms[i]->m_other_terminal->net().as_analog().Q_Analog(); rhsk_b = rhsk_b + go[i] * *other_cur_analog[i]; rhs[k] = rhsk_a + rhsk_b; @@ -426,8 +430,8 @@ void matrix_solver_direct_t<m_N, storage_N>::LE_solve() { const unsigned kN = N(); - ATTR_UNUSED int imax; - ATTR_UNUSED double big,temp; + [[maybe_unused]] int imax; + [[maybe_unused]] double big,temp; #if 0 double vv[storage_N]; @@ -438,7 +442,7 @@ void matrix_solver_direct_t<m_N, storage_N>::LE_solve() for (j=0;j<kN;j++) if ((temp=fabs(m_A[i][j])) > big) big=temp; - //if (big == 0.0) nrerror("Singular matrix in routine LUDCMP"); + //#if (big == 0.0) nrerror("Singular matrix in routine LUDCMP"); vv[i]=1.0/big; } #endif @@ -489,7 +493,6 @@ void matrix_solver_direct_t<m_N, storage_N>::LE_solve() m_A[imax][k]=m_A[j][k]; m_A[j][k]=dum; } - //*d = -(*d); vv[imax]=vv[j]; } indx[j]=imax; @@ -509,7 +512,7 @@ void matrix_solver_direct_t<m_N, storage_N>::LE_back_subst( { const unsigned kN = N(); - /* back substitution */ + // back substitution // int ip; // ii=-1 @@ -541,15 +544,14 @@ template <unsigned m_N, unsigned storage_N> nl_double matrix_solver_direct_t<m_N, storage_N>::delta( const nl_double * RESTRICT V) { - /* FIXME: Ideally we should also include currents (RHS) here. This would - * need a revaluation of the right hand side after voltages have been updated - * and thus belong into a different calculation. This applies to all solvers. - */ + // FIXME: Ideally we should also include currents (RHS) here. This would + // need a revaluation of the right hand side after voltages have been updated + // and thus belong into a different calculation. This applies to all solvers. const unsigned iN = this->N(); nl_double cerr = 0; for (unsigned i = 0; i < iN; i++) - cerr = std::fmax(cerr, std::abs(V[i] - this->m_nets[i]->m_cur_Analog)); + cerr = std::fmax(cerr, plib::abs(V[i] - this->m_nets[i]->m_cur_Analog)); return cerr; } @@ -565,7 +567,7 @@ void matrix_solver_direct_t<m_N, storage_N>::store( template <unsigned m_N, unsigned storage_N> -unsigned matrix_solver_direct_t<m_N, storage_N>::solve_non_dynamic(const bool newton_raphson) +unsigned matrix_solver_direct_t<m_N, storage_N>::solve_non_dynamic(bool newton_raphson) { nl_double new_V[storage_N]; // = { 0.0 }; @@ -587,7 +589,7 @@ unsigned matrix_solver_direct_t<m_N, storage_N>::solve_non_dynamic(const bool ne } template <unsigned m_N, unsigned storage_N> -int matrix_solver_direct_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson) +int matrix_solver_direct_t<m_N, storage_N>::upstream_solve_non_dynamic(bool newton_raphson) { this->build_LE_A(); this->build_LE_RHS(m_RHS); @@ -621,5 +623,5 @@ matrix_solver_direct_t<m_N, storage_N>::matrix_solver_direct_t(const eSolverType } //namespace devices } // namespace netlist -#endif /* NLD_MS_DIRECT_H_ */ +#endif // NLD_MS_DIRECT_H_ #endif diff --git a/src/lib/netlist/solver/nld_ms_gcr.h b/src/lib/netlist/solver/nld_ms_gcr.h index e227802c1a2..a2c5d284c43 100644 --- a/src/lib/netlist/solver/nld_ms_gcr.h +++ b/src/lib/netlist/solver/nld_ms_gcr.h @@ -1,75 +1,134 @@ -// license:GPL-2.0+ +// license:BSD-3-Clause // copyright-holders:Couriersud -/* - * nld_ms_gcr.h - * - * Gaussian elimination using compressed row format. - * - */ #ifndef NLD_MS_GCR_H_ #define NLD_MS_GCR_H_ -#include "plib/mat_cr.h" +/// +/// \file nld_ms_gcr.h +/// +/// Gaussian elimination using compressed row format. +/// -#include "nld_ms_direct.h" +#include "nld_matrix_solver_ext.h" #include "nld_solver.h" #include "plib/pdynlib.h" +#include "plib/pmatrix_cr.h" #include "plib/pstream.h" #include "plib/vector_ops.h" #include <algorithm> -namespace netlist -{ -namespace devices +namespace netlist::solver { template <typename FT, int SIZE> - class matrix_solver_GCR_t: public matrix_solver_t + class matrix_solver_GCR_t: public matrix_solver_ext_t<FT, SIZE> { 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() + using mat_type = plib::pGEmatrix_cr<plib::pmatrix_cr<arena_type, FT, SIZE>>; + using base_type = matrix_solver_ext_t<FT, SIZE>; + 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<FT, SIZE>(main_solver, name, nets, params, size) + , mat(this->m_arena, static_cast<typename mat_type::index_type>(size)) + , m_proc() + { + const std::size_t iN = this->size(); + + // build the final matrix + + std::vector<std::vector<unsigned>> 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++; + } + } - constexpr std::size_t N() const { return m_dim; } + auto gr = mat.gaussian_extend_fill_mat(fill); + + this->log_fill(fill, mat); + + mat.build_from_fill_mat(fill); - void vsetup(analog_net_t::list_t &nets) override; - unsigned vsolve_non_dynamic(const bool newton_raphson) override; + for (mat_index_type k=0; k<iN; k++) + { + 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].rail_start();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])) + { + this->m_mat_ptr[k][j] = &mat.A[i]; + cnt++; + break; + } + } + nl_assert(cnt == this->m_terms[k].rail_start()); + this->m_mat_ptr[k][this->m_terms[k].rail_start()] = &mat.A[mat.diagonal[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<fptype>(mat.nz_num) / static_cast<fptype>(iN * iN)); + this->state().log().verbose(" Pre elimination occupancy ratio: {1}", + static_cast<fptype>(raw_elements) / static_cast<fptype>(iN * iN)); - std::pair<pstring, pstring> create_solver_code() override; + // FIXME: Move me + // + + if (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 upstream_solve_non_dynamic() override; + + std::pair<pstring, pstring> create_solver_code(static_compile_target target) override; private: - using mat_index_type = typename plib::matrix_compressed_rows_t<FT, SIZE>::index_type; + using mat_index_type = typename plib::pmatrix_cr<arena_type, FT, SIZE>::index_type; - void csc_private(plib::putf8_fmt_writer &strm); + template <typename T> + void stream_if_not_yet_done(plib::putf8_fmt_writer &strm, T &A, std::size_t i) + { + const pstring fptype(fp_constants<FT>::name()); + if (!A[i].empty()) + strm("\t{1} m_A{2} = {3};\n", fptype, i, A[i]); + A[i] = ""; + } - using extsolver = void (*)(double * m_A, double * RHS, double * V); + void generate_code(plib::putf8_fmt_writer &strm); pstring static_compile_name(); - const std::size_t m_dim; - plib::parray<FT, SIZE> RHS; - plib::parray<FT, SIZE> new_V; - mat_type mat; - - //extsolver m_proc; - plib::dynproc<void, double * , double * , double * > m_proc; + plib::dynamic_library::function<void, FT *, fptype *, fptype *, fptype *, fptype ** > m_proc; }; @@ -77,149 +136,116 @@ namespace devices // matrix_solver - GCR // ---------------------------------------------------------------------------------------- - // FIXME: namespace or static class member - template <typename V> - std::size_t inline get_level(const V &v, std::size_t k) - { - 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"); - } +#define COMPRESSED 0 template <typename FT, int SIZE> - void matrix_solver_GCR_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets) + void matrix_solver_GCR_t<FT, SIZE>::generate_code(plib::putf8_fmt_writer &strm) { - setup_base(nets); - - const std::size_t iN = this->N(); - - /* build the final matrix */ - - std::vector<std::vector<unsigned>> fill(iN); + const std::size_t iN = this->size(); + const pstring fptype(fp_constants<FT>::name()); + const pstring fp_suffix(fp_constants<FT>::suffix()); + std::vector<pstring> A(this->mat.nz_num); - std::size_t raw_elements = 0; + // avoid unused variable warnings + strm("\tplib::unused_var({1});\n", "cnV"); +#if !COMPRESSED + for (std::size_t i = 0; i < mat.nz_num; i++) + strm("\t{1} m_A{2}(0.0);\n", fptype, i, i); +#endif 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) + 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<FT>::zero()); + //# *tcr_r[railstart] = static_cast<FT>(gtot_t); //mat.A[mat.diag[k]] += gtot_t; + std::size_t pd = std::size_t(this->m_mat_ptr[k][net.rail_start()] - &this->mat.A[0]); + +#if COMPRESSED + //pstring terms = plib::pfmt("m_A{1} = gt[{2}]")(pd, this->m_gtn.didx(k,0)); + pstring terms = plib::pfmt("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)); + + A[pd] = terms; //strm("\t{1};\n", terms); + //auto RHS_t(std::accumulate(Idr, Idr + term_count, plib::constants<FT>::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 = rail_start; i < term_count; i++) + // RHS_t += (- go[i]) * *cnV[i]; + + for (std::size_t i = net.rail_start(); 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 + for (std::size_t i=0; i < net.count(); i++) + strm("\tm_A{1} += gt[{2}];\n", pd, this->m_gtn.didx(k,i)); + //for (std::size_t i = 0; i < rail_start; i++) + // *tcr_r[i] += static_cast<FT>(go[i]); + for (std::size_t i = 0; i < net.rail_start(); i++) { - fill[k][j] = 0; - raw_elements++; + 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)); } + //auto RHS_t(std::accumulate(Idr, Idr + term_count, plib::constants<FT>::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 = rail_start; i < term_count; i++) + // RHS_t += (- go[i]) * *cnV[i]; - } - - auto gr = mat.gaussian_extend_fill_mat(fill); + for (std::size_t i = net.rail_start(); 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)); - /* FIXME: move this to the cr matrix class and use computed - * parallel ordering once it makes sense. - */ - - std::vector<unsigned> levL(iN, 0); - std::vector<unsigned> levU(iN, 0); - - // parallel scheme for L x = y - for (std::size_t k = 0; k < iN; k++) - { - 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; +#endif } - - // parallel scheme for U x = y - for (std::size_t k = iN; k-- > 0; ) - { - 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; - } - - +#if COMPRESSED for (std::size_t k = 0; k < iN; k++) { - unsigned fm = 0; - pstring ml = ""; - for (std::size_t j = 0; j < iN; j++) - { - 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); - } - - - mat.build_from_fill_mat(fill); - - for (mat_index_type k=0; k<iN; k++) - { - 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++) + auto &net = this->m_terms[k]; + for (std::size_t i = 0; i < net.rail_start(); i++) { - 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; - } + std::size_t p = std::size_t(this->m_mat_ptr[k][i] - &this->mat.A[0]); + if (!A[p].empty()) + A[p] += " + "; + A[p] += plib::pfmt("go[{1}]")(this->m_gonn.didx(k,i)); } - nl_assert(cnt == this->m_terms[k]->m_railstart); - m_mat_ptr[k][this->m_terms[k]->m_railstart] = &mat.A[mat.diag[k]]; } - - 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)); - - // FIXME: Move me - - if (state().lib().isLoaded()) + for (std::size_t i = 0; i < mat.nz_num; i++) { - pstring symname = static_compile_name(); - m_proc.load(this->state().lib(), symname); - if (m_proc.resolved()) - this->log().info("External static solver {1} found ...", symname); - else - this->log().warning("External static solver {1} not found ...", symname); + if (A[i].empty()) + A[i] = plib::pfmt("0.0{1}")(fp_suffix); + //strm("\t{1} m_A{2} = {3};\n", fptype, i, A[i]); } - } - - 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); +#endif for (std::size_t i = 0; i < iN - 1; i++) { - const auto &nzbd = this->m_terms[i]->m_nzbd; + //#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.size() > 0) + if (nzbd_count > 0) { - std::size_t pi = mat.diag[i]; + std::size_t pi = mat.diagonal[i]; //const FT f = 1.0 / m_A[pi++]; - strm("const double f{1} = 1.0 / m_A{2};\n", i, pi); + if ((!COMPRESSED) || nzbd_count > 1) // keep code comparable to previous versions + { + stream_if_not_yet_done(strm, A, pi); + strm("\tconst {1} f{2} = 1.0{3} / m_A{4};\n", fptype, i, fp_suffix, pi); + } pi++; const std::size_t piie = mat.row_idx[i+1]; //for (auto & j : nzbd) - for (std::size_t 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]; @@ -227,100 +253,135 @@ namespace devices pj++; //const FT f1 = - m_A[pj++] * f; - strm("\tconst double f{1}_{2} = -f{3} * m_A{4};\n", i, j, i, pj); + stream_if_not_yet_done(strm, A, pi - 1); + stream_if_not_yet_done(strm, A, pj); + + if ((!COMPRESSED) || nzbd_count > 1) // keep code comparable to previous versions + strm("\tconst {1} f{2}_{3} = -f{4} * m_A{5};\n", fptype, i, j, i, pj); + else + strm("\tconst {1} f{2}_{3} = - m_A{4} / m_A{5};\n", fptype, i, j, pj, pi-1); pj++; - // subtract row i from j */ + // subtract row i from j for (std::size_t pii = pi; pii<piie; ) { while (mat.col_idx[pj] < mat.col_idx[pii]) pj++; //m_A[pj++] += m_A[pii++] * f1; + + stream_if_not_yet_done(strm, A, pj); + stream_if_not_yet_done(strm, A, pii); + 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); + strm("\tRHS{1} += f{2}_{3} * RHS{4};\n", j, i, j, i); } } } - //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]); + //#new_V[iN - 1] = RHS[iN - 1] / mat.A[mat.diag[iN - 1]]; + stream_if_not_yet_done(strm, A, mat.diagonal[iN - 1]); + strm("\tV[{1}] = RHS{2} / m_A{3};\n", iN - 1, iN - 1, mat.diagonal[iN - 1]); for (std::size_t j = iN - 1; j-- > 0;) { - strm("\tdouble tmp{1} = 0.0;\n", j); +#if COMPRESSED + pstring tmp; + const std::size_t e = mat.row_idx[j+1]; + for (std::size_t pk = mat.diagonal[j] + 1; pk < e; pk++) + { + stream_if_not_yet_done(strm, A, pk); + tmp = tmp + plib::pfmt(" + m_A{2} * V[{3}]")(j, pk, mat.col_idx[pk]); + } + + stream_if_not_yet_done(strm, A, mat.diagonal[j]); + if (tmp.empty()) + { + strm("\tV[{1}] = RHS{1} / m_A{2};\n", j, mat.diagonal[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.diagonal[j]); + } +#else + strm("\t{1} tmp{2} = 0.0{3};\n", fptype, j, fp_suffix); const std::size_t e = mat.row_idx[j+1]; - for (std::size_t pk = mat.diag[j] + 1; pk < e; pk++) + for (std::size_t pk = mat.diagonal[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]); + strm("\tV[{1}] = (RHS{1} - tmp{1}) / m_A{4};\n", j, j, j, mat.diagonal[j]); +#endif } } template <typename FT, int SIZE> pstring matrix_solver_GCR_t<FT, SIZE>::static_compile_name() { + pstring str_floattype(fp_constants<FT>::name()); + pstring str_fptype(fp_constants<fptype>::name()); std::stringstream t; t.imbue(std::locale::classic()); plib::putf8_fmt_writer w(&t); - csc_private(w); - std::hash<typename std::remove_const<std::remove_reference<decltype(t.str())>::type>::type> h; - - return plib::pfmt("nl_gcr_{1:x}_{2}")(h( t.str() ))(mat.nz_num); + generate_code(w); + //#std::hash<typename std::remove_const<std::remove_reference<decltype(t.str())>::type>::type> h; + return plib::pfmt("nl_gcr_{1}_{2}_{3}_{4:x}")(mat.nz_num)(str_fptype)(str_floattype)(plib::hash<uint64_t>( t.str().c_str(), t.str().size() )); } template <typename FT, int SIZE> - std::pair<pstring, pstring> matrix_solver_GCR_t<FT, SIZE>::create_solver_code() + std::pair<pstring, pstring> matrix_solver_GCR_t<FT, SIZE>::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(); - - 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, pstring(t.str())); + pstring str_float_type(fp_constants<FT>::name()); + pstring str_fptype(fp_constants<fptype>::name()); + + pstring external_qualifier; + if (target == CXX_EXTERNAL_C) + external_qualifier = "extern \"C\""; + else if (target == CXX_STATIC) + external_qualifier = "static"; + strm.write_line(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")(external_qualifier, name, str_float_type, str_fptype)); + strm.write_line("{\n"); + generate_code(strm); + strm.write_line("}\n"); + // some compilers (_WIN32, _WIN64, mac osx) need an explicit cast + return { name, putf8string(t.str()) }; } template <typename FT, int SIZE> - unsigned matrix_solver_GCR_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson) + void matrix_solver_GCR_t<FT, SIZE>::upstream_solve_non_dynamic() { - 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]); + 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 { - // mat.gaussian_elimination_parallel(RHS); - mat.gaussian_elimination(RHS); - /* backward substitution */ - mat.gaussian_back_substitution(new_V, RHS); - } + // clear matrix + mat.set_scalar(plib::constants<FT>::zero()); - this->m_stat_calculations++; + // populate matrix + this->fill_matrix_and_rhs(); - const FT err = (newton_raphson ? delta(new_V) : 0.0); - store(new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; + // 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 devices -} // namespace netlist +} // namespace netlist::solver -#endif /* NLD_MS_GCR_H_ */ +#endif // NLD_MS_GCR_H_ diff --git a/src/lib/netlist/solver/nld_ms_gmres.h b/src/lib/netlist/solver/nld_ms_gmres.h index 2ff515ebda7..836f6e61cbb 100644 --- a/src/lib/netlist/solver/nld_ms_gmres.h +++ b/src/lib/netlist/solver/nld_ms_gmres.h @@ -1,27 +1,24 @@ -// license:GPL-2.0+ +// license:BSD-3-Clause // copyright-holders:Couriersud -/* - * nld_ms_gmres.h - * - */ #ifndef NLD_MS_GMRES_H_ #define NLD_MS_GMRES_H_ +/// +/// \file nld_ms_gmres.h +/// + +#include "nld_matrix_solver_ext.h" #include "nld_ms_direct.h" #include "nld_solver.h" #include "plib/gmres.h" -#include "plib/mat_cr.h" #include "plib/parray.h" +#include "plib/pmatrix_cr.h" #include "plib/vector_ops.h" #include <algorithm> -#include <cmath> - -namespace netlist -{ -namespace devices +namespace netlist::solver { template <typename FT, int SIZE> @@ -31,26 +28,64 @@ namespace devices using float_type = FT; - /* Sort rows in ascending order. This should minimize fill-in and thus - * maximize the efficiency of the incomplete LUT. - * This is already preconditioning. - */ - matrix_solver_GMRES_t(netlist_state_t &anetlist, const pstring &name, const solver_parameters_t *params, const std::size_t size) - : matrix_solver_direct_t<FT, SIZE>(anetlist, name, matrix_solver_t::PREFER_BAND_MATRIX, params, size) - //, m_ops(size, 2) - , m_ops(size, 4) + // Sort rows in ascending order. This should minimize fill-in and thus + // maximize the efficiency of the incomplete LUT. + // This is already preconditioning. + + matrix_solver_GMRES_t(devices::nld_solver &main_solver, const pstring &name, + matrix_solver_t::net_list_t &nets, + const solver::solver_parameters_t *params, + const std::size_t size) + : matrix_solver_direct_t<FT, SIZE>(main_solver, name, nets, params, size) + , m_ops(this->m_arena, size, 0) , m_gmres(size) { + const std::size_t iN = this->size(); + + std::vector<std::vector<unsigned>> fill(iN); + + for (std::size_t k=0; k<iN; k++) + { + fill[k].resize(iN, decltype(m_ops.m_mat)::FILL_INFINITY); + terms_for_net_t & row = this->m_terms[k]; + for (const auto &nz_j : row.m_nz) + { + fill[k][static_cast<matrix_type>(nz_j)] = 0; + } } - void vsetup(analog_net_t::list_t &nets) override; - unsigned vsolve_non_dynamic(const bool newton_raphson) override; + m_ops.build(fill); + this->log_fill(fill, m_ops.m_mat); + + // build pointers into the compressed row format matrix for each terminal + + for (std::size_t k=0; k<iN; k++) + { + std::size_t cnt = 0; + for (std::size_t j=0; j< this->m_terms[k].rail_start();j++) + { + for (std::size_t i = m_ops.m_mat.row_idx[k]; i<m_ops.m_mat.row_idx[k+1]; i++) + if (this->m_terms[k].m_connected_net_idx[j] == static_cast<int>(m_ops.m_mat.col_idx[i])) + { + this->m_mat_ptr[k][j] = &m_ops.m_mat.A[i]; + cnt++; + break; + } + } + nl_assert(cnt == this->m_terms[k].rail_start()); + this->m_mat_ptr[k][this->m_terms[k].rail_start()] = &m_ops.m_mat.A[m_ops.m_mat.diagonal[k]]; + } + } + + void upstream_solve_non_dynamic() override; private: - using mattype = typename plib::matrix_compressed_rows_t<FT, SIZE>::index_type; + using matrix_type = typename plib::pmatrix_cr<arena_type, FT, SIZE>::index_type; - plib::mat_precondition_ILU<FT, SIZE> m_ops; + //plib::mat_precondition_none<FT, SIZE> m_ops; + plib::mat_precondition_ILU<arena_type, FT, SIZE> m_ops; + //plib::mat_precondition_diag<FT, SIZE> m_ops; plib::gmres_t<FT, SIZE> m_gmres; }; @@ -59,87 +94,38 @@ namespace devices // ---------------------------------------------------------------------------------------- template <typename FT, int SIZE> - void matrix_solver_GMRES_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets) + void matrix_solver_GMRES_t<FT, SIZE>::upstream_solve_non_dynamic() { - matrix_solver_direct_t<FT, SIZE>::vsetup(nets); - const std::size_t iN = this->size(); - std::vector<std::vector<unsigned>> fill(iN); - - for (std::size_t k=0; k<iN; k++) - { - fill[k].resize(iN, decltype(m_ops.m_mat)::FILL_INFINITY); - terms_for_net_t * row = this->m_terms[k].get(); - for (const auto &nz_j : row->m_nz) - { - fill[k][static_cast<mattype>(nz_j)] = 0; - } - } - - m_ops.build(fill); - - /* build pointers into the compressed row format matrix for each terminal */ - - for (std::size_t k=0; k<iN; k++) - { - std::size_t cnt = 0; - for (std::size_t j=0; j< this->m_terms[k]->m_railstart;j++) - { - for (std::size_t i = m_ops.m_mat.row_idx[k]; i<m_ops.m_mat.row_idx[k+1]; i++) - if (this->m_terms[k]->m_connected_net_idx[j] == static_cast<int>(m_ops.m_mat.col_idx[i])) - { - this->m_mat_ptr[k][j] = &m_ops.m_mat.A[i]; - cnt++; - break; - } - } - nl_assert(cnt == this->m_terms[k]->m_railstart); - this->m_mat_ptr[k][this->m_terms[k]->m_railstart] = &m_ops.m_mat.A[m_ops.m_mat.diag[k]]; - } - } - - template <typename FT, int SIZE> - unsigned matrix_solver_GMRES_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson) - { - const std::size_t iN = this->size(); - - plib::parray<FT, SIZE> RHS(iN); - //float_type new_V[storage_N]; - - m_ops.m_mat.set_scalar(0.0); + m_ops.m_mat.set_scalar(plib::constants<FT>::zero()); - /* populate matrix and V for first estimate */ - this->fill_matrix(iN, this->m_mat_ptr, RHS); + // populate matrix and V for first estimate + this->fill_matrix_and_rhs(); for (std::size_t k = 0; k < iN; k++) { - this->m_new_V[k] = this->m_nets[k]->Q_Analog(); + this->m_new_V[k] = static_cast<float_type>(this->m_terms[k].getV()); } - const float_type accuracy = this->m_params.m_accuracy; + const float_type accuracy(static_cast<float_type>(this->m_params.m_accuracy)); - auto iter = std::max(plib::constants<std::size_t>::one(), this->m_params.m_gs_loops); - auto gsl = m_gmres.solve(m_ops, this->m_new_V, RHS, iter, accuracy); + auto iter = std::max(plib::constants<std::size_t>::one(), this->m_params.m_gs_loops()); + auto gsl = m_gmres.solve(m_ops, this->m_new_V, this->m_RHS, iter, accuracy); this->m_iterative_total += gsl; - this->m_stat_calculations++; if (gsl > iter) { this->m_iterative_fail++; - return matrix_solver_direct_t<FT, SIZE>::vsolve_non_dynamic(newton_raphson); + matrix_solver_direct_t<FT, SIZE>::upstream_solve_non_dynamic(); } - const float_type err = (newton_raphson ? this->delta(this->m_new_V) : 0.0); - this->store(this->m_new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; } -} // namespace devices -} // namespace netlist +} // namespace netlist::solver -#endif /* NLD_MS_GMRES_H_ */ +#endif // NLD_MS_GMRES_H_ diff --git a/src/lib/netlist/solver/nld_ms_sm.h b/src/lib/netlist/solver/nld_ms_sm.h index 5b59effecf8..8a0b67f18ce 100644 --- a/src/lib/netlist/solver/nld_ms_sm.h +++ b/src/lib/netlist/solver/nld_ms_sm.h @@ -1,54 +1,50 @@ -// license:GPL-2.0+ +// license:BSD-3-Clause // copyright-holders:Couriersud -/* - * nld_ms_direct.h - * - * - * Sherman-Morrison Solver - * - * Computes the updated inverse of A given that the change in A is - * - * A <- A + (u x v) u,v vectors - * - * In this specific implementation, u is a unit vector specifying the row which - * changed. Thus v contains the changed column. - * - * Than z = A^-1 u , w = transpose(A^-1) v , lambda = v z - * - * A^-1 <- 1.0 / (1.0 + lambda) * (z x w) - * - * The approach is iterative and applied for each row changed. - * - * The performance for a typical circuit like kidniki compared to Gaussian - * elimination is poor: - * - * a) The code needs to be run for each row change. - * b) The inverse of A typically is fully occupied. - * - * It may have advantages for circuits with a high number of elements and only - * few dynamic/active components. - * - */ #ifndef NLD_MS_SM_H_ #define NLD_MS_SM_H_ +/// +/// \file nld_ms_sm.h +/// +/// Sherman-Morrison Solver +/// +/// Computes the updated inverse of A given that the change in A is +/// +/// A <- A + (u x v) u,v vectors +/// +/// In this specific implementation, u is a unit vector specifying the row which +/// changed. Thus v contains the changed column. +/// +/// Than z = A^-1 u , w = transpose(A^-1) v , lambda = v z +/// +/// A^-1 <- 1.0 / (1.0 + lambda) * (z x w) +/// +/// The approach is iterative and applied for each row changed. +/// +/// The performance for a typical circuit like `kidniki` compared to Gaussian +/// elimination is poor: +/// +/// a) The code needs to be run for each row change. +/// b) The inverse of A typically is fully occupied. +/// +/// It may have advantages for circuits with a high number of elements and only +/// few dynamic/active components. +/// + #include "nld_matrix_solver.h" +#include "nld_matrix_solver_ext.h" #include "nld_solver.h" #include "plib/vector_ops.h" #include <algorithm> -namespace netlist -{ -namespace devices +namespace netlist::solver { template <typename FT, int SIZE> - class matrix_solver_sm_t: public matrix_solver_t + class matrix_solver_sm_t: public matrix_solver_ext_t<FT, SIZE> { - friend class matrix_solver_t; - public: using float_ext_type = FT; @@ -56,17 +52,20 @@ namespace devices // FIXME: dirty hack to make this compile static constexpr const std::size_t storage_N = 100; - matrix_solver_sm_t(netlist_state_t &anetlist, const pstring &name, - const solver_parameters_t *params, const std::size_t size); + matrix_solver_sm_t(devices::nld_solver &main_solver, const pstring &name, + const matrix_solver_t::net_list_t &nets, + const solver_parameters_t *params, const std::size_t size) + : matrix_solver_ext_t<FT, SIZE>(main_solver, name, nets, params, size) + , m_cnt(0) + { + this->build_mat_ptr(m_A); + } - void vsetup(analog_net_t::list_t &nets) override; void reset() override { matrix_solver_t::reset(); } protected: - unsigned vsolve_non_dynamic(const bool newton_raphson) override; - unsigned solve_non_dynamic(const bool newton_raphson); - - constexpr std::size_t size() const { return m_dim; } + void upstream_solve_non_dynamic() override; + void solve_non_dynamic(); void LE_invert(); @@ -81,7 +80,7 @@ namespace devices template <typename T1, typename T2> float_ext_type &Ainv(const T1 &r, const T2 &c) { return m_Ainv[r][c]; } template <typename T1> - float_ext_type &RHS(const T1 &r) { return m_RHS[r]; } + float_ext_type &RHS(const T1 &r) { return this->m_RHS[r]; } template <typename T1, typename T2> @@ -90,18 +89,18 @@ namespace devices float_ext_type &lAinv(const T1 &r, const T2 &c) { return m_lAinv[r][c]; } private: + template <typename T, std::size_t N, std::size_t M> + using array2D = std::array<std::array<T, M>, N>; static constexpr std::size_t m_pitch = ((( storage_N) + 7) / 8) * 8; - float_ext_type m_A[storage_N][m_pitch]; - float_ext_type m_Ainv[storage_N][m_pitch]; - float_ext_type m_W[storage_N][m_pitch]; - std::array<float_ext_type, storage_N> m_RHS; // right hand side - contains currents + array2D<float_ext_type, storage_N, m_pitch> m_A; + array2D<float_ext_type, storage_N, m_pitch> m_Ainv; + array2D<float_ext_type, storage_N, m_pitch> m_W; - float_ext_type m_lA[storage_N][m_pitch]; - float_ext_type m_lAinv[storage_N][m_pitch]; + array2D<float_ext_type, storage_N, m_pitch> m_lA; + array2D<float_ext_type, storage_N, m_pitch> m_lAinv; //float_ext_type m_RHSx[storage_N]; - const std::size_t m_dim; std::size_t m_cnt; }; @@ -111,46 +110,37 @@ namespace devices // ---------------------------------------------------------------------------------------- template <typename FT, int SIZE> - void matrix_solver_sm_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets) - { - matrix_solver_t::setup_base(nets); - - /* FIXME: Shouldn't be necessary */ - for (std::size_t k = 0; k < size(); k++) - state().save(*this, RHS(k), this->name(), plib::pfmt("RHS.{1}")(k)); - } - - template <typename FT, int SIZE> void matrix_solver_sm_t<FT, SIZE>::LE_invert() { - const std::size_t kN = size(); + const std::size_t kN = this->size(); for (std::size_t i = 0; i < kN; i++) { for (std::size_t j = 0; j < kN; j++) { W(i,j) = lA(i,j) = A(i,j); - Ainv(i,j) = 0.0; + Ainv(i,j) = plib::constants<FT>::zero(); } - Ainv(i,i) = 1.0; + Ainv(i,i) = plib::constants<FT>::one(); } - /* down */ + // down for (std::size_t i = 0; i < kN; i++) { - /* FIXME: Singular matrix? */ - const float_type f = 1.0 / W(i,i); - const auto * const p = m_terms[i]->m_nzrd.data(); - const std::size_t e = m_terms[i]->m_nzrd.size(); + // FIXME: Singular matrix? + const float_type f = plib::reciprocal(W(i,i)); + const auto * const p = this->m_terms[i].m_nzrd.data(); + const std::size_t e = this->m_terms[i].m_nzrd.size(); - /* Eliminate column i from row j */ + // Eliminate column i from row j - const auto * const pb = m_terms[i]->m_nzbd.data(); - const std::size_t eb = m_terms[i]->m_nzbd.size(); + const auto * const pb = this->m_terms[i].m_nzbd.data(); + const std::size_t eb = this->m_terms[i].m_nzbd.size(); for (std::size_t jb = 0; jb < eb; jb++) { const unsigned j = pb[jb]; const float_type f1 = - W(j,i) * f; - if (f1 != 0.0) + // FIXME: comparison to zero + if (f1 != plib::constants<float_type>::zero()) { for (std::size_t k = 0; k < e; k++) W(j,p[k]) += W(i,p[k]) * f1; @@ -159,15 +149,16 @@ namespace devices } } } - /* up */ + // up for (std::size_t i = kN; i-- > 0; ) { - /* FIXME: Singular matrix? */ - const float_type f = 1.0 / W(i,i); + // FIXME: Singular matrix? + const float_type f = plib::reciprocal(W(i,i)); for (std::size_t j = i; j-- > 0; ) { const float_type f1 = - W(j,i) * f; - if (f1 != 0.0) + // FIXME: comparison to zero + if (f1 != plib::constants<float_type>::zero()) { for (std::size_t k = i; k < kN; k++) W(j,k) += W(i,k) * f1; @@ -188,10 +179,10 @@ namespace devices void matrix_solver_sm_t<FT, SIZE>::LE_compute_x( T & x) { - const std::size_t kN = size(); + const std::size_t kN = this->size(); for (std::size_t i=0; i<kN; i++) - x[i] = 0.0; + x[i] = plib::constants<FT>::zero(); for (std::size_t k=0; k<kN; k++) { @@ -203,16 +194,21 @@ namespace devices } template <typename FT, int SIZE> - unsigned matrix_solver_sm_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphson) + void matrix_solver_sm_t<FT, SIZE>::solve_non_dynamic() { static constexpr const bool incremental = true; - const std::size_t iN = size(); + const std::size_t iN = this->size(); - std::array<float_type, storage_N> new_V; // = { 0.0 }; + // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) + std::array<float_type, m_pitch> v; + // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) + std::array<std::size_t, m_pitch> cols; + // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) + std::array<float_type, m_pitch> z; if ((m_cnt % 50) == 0) { - /* complete calculation */ + // complete calculation this->LE_invert(); } else @@ -225,46 +221,45 @@ namespace devices } for (std::size_t row = 0; row < iN; row ++) { - std::array<float_type, m_pitch> v = {0}; - std::array<std::size_t, m_pitch> cols; - std::size_t colcount = 0; + std::size_t col_count = 0; - auto &nz = m_terms[row]->m_nz; + auto &nz = this->m_terms[row].m_nz; for (unsigned & col : nz) { v[col] = A(row,col) - lA(row,col); if (incremental) lA(row,col) = A(row,col); - if (v[col] != 0.0) - cols[colcount++] = col; + // FIXME: comparison to zero + if (v[col] != plib::constants<float_type>::zero()) + cols[col_count++] = col; } - if (colcount > 0) + if (col_count > 0) { - float_type lamba = 0.0; + auto lambda(plib::constants<FT>::zero()); std::array<float_type, m_pitch> w = {0}; - std::array<float_type, m_pitch> z; - /* compute w and lamba */ + // compute w and lambda for (std::size_t i = 0; i < iN; i++) - z[i] = Ainv(i, row); /* u is row'th column */ + z[i] = Ainv(i, row); // u is row'th column - for (std::size_t j = 0; j < colcount; j++) - lamba += v[cols[j]] * z[cols[j]]; + for (std::size_t j = 0; j < col_count; j++) + lambda += v[cols[j]] * z[cols[j]]; - for (std::size_t j=0; j<colcount; j++) + for (std::size_t j=0; j<col_count; j++) { std::size_t col = cols[j]; float_type f = v[col]; for (std::size_t k = 0; k < iN; k++) - w[k] += Ainv(col,k) * f; /* Transpose(Ainv) * v */ + w[k] += Ainv(col,k) * f; //# Transpose(Ainv) * v } - lamba = -1.0 / (1.0 + lamba); + lambda = -plib::reciprocal(plib::constants<float_type>::one() + lambda); for (std::size_t i=0; i<iN; i++) { - const float_type f = lamba * z[i]; - if (f != 0.0) + const float_type f = lambda * z[i]; + // FIXME: comparison to zero + if (f != plib::constants<float_type>::zero()) for (std::size_t k = 0; k < iN; k++) Ainv(i,k) += f * w[k]; } @@ -275,33 +270,20 @@ namespace devices m_cnt++; - this->LE_compute_x(new_V); - - const float_type err = (newton_raphson ? delta(new_V) : 0.0); - store(new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; + this->LE_compute_x(this->m_new_V); } template <typename FT, int SIZE> - unsigned matrix_solver_sm_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson) + void matrix_solver_sm_t<FT, SIZE>::upstream_solve_non_dynamic() { - this->build_LE_A(*this); - this->build_LE_RHS(*this); - this->m_stat_calculations++; - return this->solve_non_dynamic(newton_raphson); - } + this->clear_square_mat(this->m_A); + this->fill_matrix_and_rhs(); - template <typename FT, int SIZE> - matrix_solver_sm_t<FT, SIZE>::matrix_solver_sm_t(netlist_state_t &anetlist, const pstring &name, - const solver_parameters_t *params, const std::size_t size) - : matrix_solver_t(anetlist, name, NOSORT, params) - , m_dim(size) - , m_cnt(0) - { + this->solve_non_dynamic(); } -} // namespace devices -} // namespace netlist -#endif /* NLD_MS_DIRECT_H_ */ +} // namespace netlist::solver + +#endif // NLD_MS_SM_H_ diff --git a/src/lib/netlist/solver/nld_ms_sor.h b/src/lib/netlist/solver/nld_ms_sor.h index c31aaa6d46a..ca5a4839ce3 100644 --- a/src/lib/netlist/solver/nld_ms_sor.h +++ b/src/lib/netlist/solver/nld_ms_sor.h @@ -1,175 +1,164 @@ -// license:GPL-2.0+ +// license:BSD-3-Clause // copyright-holders:Couriersud -/* - * nld_ms_sor.h - * - * Generic successive over relaxation solver. - * - * Fow w==1 we will do the classic Gauss-Seidel approach - * - */ #ifndef NLD_MS_SOR_H_ #define NLD_MS_SOR_H_ +// Names +// spell-checker: words Seidel + +/// +/// \file nld_ms_sor.h +/// +/// Generic successive over relaxation solver. +/// +/// Fow w==1 we will do the classic Gauss-Seidel approach. +/// + +#include "nld_matrix_solver_ext.h" #include "nld_ms_direct.h" -#include "nld_solver.h" #include <algorithm> -namespace netlist -{ - namespace devices -{ - -template <typename FT, int SIZE> -class matrix_solver_SOR_t: public matrix_solver_direct_t<FT, SIZE> +namespace netlist::solver { -public: - using float_type = FT; - - matrix_solver_SOR_t(netlist_state_t &anetlist, const pstring &name, const solver_parameters_t *params, const std::size_t size) - : matrix_solver_direct_t<FT, SIZE>(anetlist, name, matrix_solver_t::ASCENDING, params, size) - , m_lp_fact(*this, "m_lp_fact", 0) - , w(size, 0.0) - , one_m_w(size, 0.0) - , RHS(size, 0.0) - //, new_V(size, 0.0) - { - } + template <typename FT, int SIZE> + class matrix_solver_SOR_t: public matrix_solver_direct_t<FT, SIZE> + { + public: - void vsetup(analog_net_t::list_t &nets) override; - unsigned vsolve_non_dynamic(const bool newton_raphson) override; + using float_type = FT; -private: - state_var<float_type> m_lp_fact; - std::vector<float_type> w; - std::vector<float_type> one_m_w; - std::vector<float_type> RHS; - //std::vector<float_type> new_V; -}; + matrix_solver_SOR_t(devices::nld_solver &main_solver, const pstring &name, + matrix_solver_t::net_list_t &nets, + const solver_parameters_t *params, const std::size_t size) + : matrix_solver_direct_t<FT, SIZE>(main_solver, name, nets, params, size) + , m_lp_fact(*this, "m_lp_fact", 0) + , w(size, plib::constants<FT>::zero()) + , one_m_w(size, plib::constants<FT>::zero()) + { + } -// ---------------------------------------------------------------------------------------- -// matrix_solver - Gauss - Seidel -// ---------------------------------------------------------------------------------------- + void upstream_solve_non_dynamic() override; + private: + state_var<float_type> m_lp_fact; + std::vector<float_type> w; + std::vector<float_type> one_m_w; + }; -template <typename FT, int SIZE> -void matrix_solver_SOR_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets) -{ - matrix_solver_direct_t<FT, SIZE>::vsetup(nets); -} + /// + /// \brief Gauss - Seidel matrix_solver + /// + /// + template <typename FT, int SIZE> + void matrix_solver_SOR_t<FT, SIZE>::upstream_solve_non_dynamic() + { + const std::size_t iN = this->size(); + bool resched = false; + unsigned resched_cnt = 0; -template <typename FT, int SIZE> -unsigned matrix_solver_SOR_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson) -{ - const std::size_t iN = this->size(); - bool resched = false; - unsigned resched_cnt = 0; + // ideally, we could get an estimate for the spectral radius of + // Inv(D - L) * U + // + // and estimate using + // + // omega = 2.0 / (1.0 + std::sqrt(1-rho)) + // - /* ideally, we could get an estimate for the spectral radius of - * Inv(D - L) * U - * - * and estimate using - * - * omega = 2.0 / (1.0 + std::sqrt(1-rho)) - */ + const auto ws(static_cast<float_type>(this->m_params.m_gs_sor)); - const float_type ws = this->m_params.m_gs_sor; + for (std::size_t k = 0; k < iN; k++) + { - for (std::size_t k = 0; k < iN; k++) - { - float_type gtot_t = 0.0; - float_type gabs_t = 0.0; - float_type RHS_t = 0.0; + const std::size_t term_count = this->m_terms[k].count(); + const auto * const gt = this->m_gtn[k]; + const auto * const go = this->m_gonn[k]; + const auto * const Idr = this->m_Idrn[k]; + auto other_cur_analog = this->m_connected_net_Vn[k]; - const std::size_t term_count = this->m_terms[k]->count(); - const float_type * const gt = this->m_gtn[k]; - const float_type * const go = this->m_gonn[k]; - const float_type * const Idr = this->m_Idrn[k]; - auto other_cur_analog = this->m_connected_net_Vn[k]; + using fpaggtype = std::remove_reference_t<std::remove_cv_t<decltype(this->m_gtn[0][0])>>; - this->m_new_V[k] = this->m_nets[k]->Q_Analog(); + fpaggtype gtot_t = nlconst_base<fpaggtype>::zero(); + fpaggtype gabs_t = nlconst_base<fpaggtype>::zero(); + fpaggtype RHS_t = nlconst_base<fpaggtype>::zero(); - for (std::size_t i = 0; i < term_count; i++) - { - gtot_t = gtot_t + gt[i]; - RHS_t = RHS_t + Idr[i]; - } + this->m_new_V[k] = static_cast<float_type>(this->m_terms[k].getV()); - for (std::size_t i = this->m_terms[k]->m_railstart; i < term_count; i++) - RHS_t = RHS_t - go[i] * *other_cur_analog[i]; + for (std::size_t i = 0; i < term_count; i++) + { + gtot_t = gtot_t + gt[i]; + RHS_t = RHS_t + Idr[i]; + } - RHS[k] = RHS_t; + for (std::size_t i = this->m_terms[k].rail_start(); i < term_count; i++) + RHS_t = RHS_t - go[i] * *other_cur_analog[i]; - if (this->m_params.m_use_gabs) - { - for (std::size_t i = 0; i < term_count; i++) - gabs_t = gabs_t + std::abs(go[i]); + this->m_RHS[k] = static_cast<float_type>(RHS_t); - gabs_t *= plib::constants<nl_double>::cast(0.5); // derived by try and error - if (gabs_t <= gtot_t) + if (this->m_params.m_use_gabs) { - w[k] = ws / gtot_t; - one_m_w[k] = plib::constants<FT>::one() - ws; + for (std::size_t i = 0; i < term_count; i++) + gabs_t = gabs_t + plib::abs(go[i]); + + gabs_t *= nlconst::half(); // derived by try and error + if (gabs_t <= gtot_t) + { + w[k] = ws / static_cast<float_type>(gtot_t); + one_m_w[k] = plib::constants<FT>::one() - ws; + } + else + { + w[k] = plib::reciprocal(static_cast<float_type>(gtot_t + gabs_t)); + one_m_w[k] = plib::constants<FT>::one() - plib::constants<FT>::one() * static_cast<FT>(gtot_t / (gtot_t + gabs_t)); + } } else { - w[k] = plib::constants<FT>::one() / (gtot_t + gabs_t); - one_m_w[k] = plib::constants<FT>::one() - plib::constants<FT>::one() * gtot_t / (gtot_t + gabs_t); + w[k] = ws / static_cast<float_type>(gtot_t); + one_m_w[k] = plib::constants<FT>::one() - ws; } } - else - { - w[k] = ws / gtot_t; - one_m_w[k] = plib::constants<FT>::one() - ws; - } - } - const float_type accuracy = this->m_params.m_accuracy; + const auto accuracy(static_cast<float_type>(this->m_params.m_accuracy)); - do { - resched = false; - float_type err = 0; - for (std::size_t k = 0; k < iN; k++) - { - const int * net_other = this->m_terms[k]->m_connected_net_idx.data(); - const std::size_t railstart = this->m_terms[k]->m_railstart; - const float_type * go = this->m_gonn[k]; + do { + resched = false; + float_type err = 0; + for (std::size_t k = 0; k < iN; k++) + { + const int * net_other = this->m_terms[k].m_connected_net_idx.data(); + const std::size_t rail_start = this->m_terms[k].rail_start(); + const auto * go = this->m_gonn[k]; - float_type Idrive = 0.0; - for (std::size_t i = 0; i < railstart; i++) - Idrive = Idrive - go[i] * this->m_new_V[static_cast<std::size_t>(net_other[i])]; + float_type Idrive = plib::constants<float_type>::zero(); + for (std::size_t i = 0; i < rail_start; i++) + Idrive = Idrive - static_cast<float_type>(go[i]) * this->m_new_V[static_cast<std::size_t>(net_other[i])]; - const float_type new_val = this->m_new_V[k] * one_m_w[k] + (Idrive + RHS[k]) * w[k]; + const float_type new_val = this->m_new_V[k] * one_m_w[k] + (Idrive + this->m_RHS[k]) * w[k]; - err = std::max(std::abs(new_val - this->m_new_V[k]), err); - this->m_new_V[k] = new_val; - } + err = std::max(plib::abs(new_val - this->m_new_V[k]), err); + this->m_new_V[k] = new_val; + } - if (err > accuracy) - resched = true; + if (err > accuracy) + resched = true; - resched_cnt++; - } while (resched && (resched_cnt < this->m_params.m_gs_loops)); + resched_cnt++; + } while (resched && (resched_cnt < this->m_params.m_gs_loops)); - this->m_iterative_total += resched_cnt; - this->m_stat_calculations++; + this->m_iterative_total += resched_cnt; - if (resched) - { - // Fallback to direct solver ... - this->m_iterative_fail++; - return matrix_solver_direct_t<FT, SIZE>::vsolve_non_dynamic(newton_raphson); - } + if (resched) + { + // Fallback to direct solver ... + this->m_iterative_fail++; + matrix_solver_direct_t<FT, SIZE>::upstream_solve_non_dynamic(); + } - const float_type err = (newton_raphson ? this->delta(this->m_new_V) : 0.0); - this->store(this->m_new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; -} + } - } //namespace devices -} // namespace netlist +} // namespace netlist::solver -#endif /* NLD_MS_SOR_H_ */ +#endif // NLD_MS_SOR_H_ diff --git a/src/lib/netlist/solver/nld_ms_sor_mat.h b/src/lib/netlist/solver/nld_ms_sor_mat.h index cdcd2e48847..dd370da4df9 100644 --- a/src/lib/netlist/solver/nld_ms_sor_mat.h +++ b/src/lib/netlist/solver/nld_ms_sor_mat.h @@ -1,56 +1,48 @@ -// license:GPL-2.0+ +// license:BSD-3-Clause // copyright-holders:Couriersud -/* - * nld_ms_sor.h - * - * Generic successive over relaxation solver. - * - * Fow w==1 we will do the classic Gauss-Seidel approach - * - */ #ifndef NLD_MS_SOR_MAT_H_ #define NLD_MS_SOR_MAT_H_ -#include "nld_matrix_solver.h" +// Names +// spell-checker: words Seidel, +// + +/// +/// \file nld_ms_sor.h +/// +/// Generic successive over relaxation solver. +/// +/// For w==1 we will do the classic Gauss-Seidel approach +/// + +#include "nld_matrix_solver_ext.h" #include "nld_ms_direct.h" -#include "nld_solver.h" #include <algorithm> -namespace netlist -{ -namespace devices +namespace netlist::solver { template <typename FT, int SIZE> class matrix_solver_SOR_mat_t: public matrix_solver_direct_t<FT, SIZE> { - friend class matrix_solver_t; - public: using float_type = FT; - matrix_solver_SOR_mat_t(netlist_state_t &anetlist, const pstring &name, const solver_parameters_t *params, std::size_t size) - : matrix_solver_direct_t<FT, SIZE>(anetlist, name, matrix_solver_t::ASCENDING, params, size) - , m_Vdelta(*this, "m_Vdelta", std::vector<float_type>(size)) - , m_omega(*this, "m_omega", params->m_gs_sor) - , m_lp_fact(*this, "m_lp_fact", 0) + matrix_solver_SOR_mat_t(devices::nld_solver &main_solver, const pstring &name, + const matrix_solver_t::net_list_t &nets, + const solver_parameters_t *params, std::size_t size) + : matrix_solver_direct_t<FT, SIZE>(main_solver, name, nets, params, size) + , m_omega(*this, "m_omega", static_cast<float_type>(params->m_gs_sor)) { } - void vsetup(analog_net_t::list_t &nets) override; - - unsigned vsolve_non_dynamic(const bool newton_raphson) override; + void upstream_solve_non_dynamic() override; private: - //state_var<float_type[storage_N]> m_Vdelta; - state_var<std::vector<float_type>> m_Vdelta; - state_var<float_type> m_omega; - state_var<float_type> m_lp_fact; - }; // ---------------------------------------------------------------------------------------- @@ -58,79 +50,21 @@ namespace devices // ---------------------------------------------------------------------------------------- template <typename FT, int SIZE> - void matrix_solver_SOR_mat_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets) - { - matrix_solver_direct_t<FT, SIZE>::vsetup(nets); - } - - #if 0 - //FIXME: move to solve_base - template <unsigned m_N, unsigned storage_N> - float_type matrix_solver_SOR_mat_t<m_N, storage_N>::vsolve() - { - /* - * enable linear prediction on first newton pass - */ - - if (this->m_params->use_linear_prediction) - for (unsigned k = 0; k < this->size(); k++) - { - this->m_last_V[k] = this->m_nets[k]->m_cur_Analog; - this->m_nets[k]->m_cur_Analog = this->m_nets[k]->m_cur_Analog + this->m_Vdelta[k] * this->current_timestep() * m_lp_fact; - } - else - for (unsigned k = 0; k < this->size(); k++) - { - this->m_last_V[k] = this->m_nets[k]->m_cur_Analog; - } - - this->solve_base(this); - - if (this->m_params->use_linear_prediction) - { - float_type sq = 0; - float_type sqo = 0; - const float_type rez_cts = 1.0 / this->current_timestep(); - for (unsigned k = 0; k < this->size(); k++) - { - const analog_net_t *n = this->m_nets[k]; - const float_type nv = (n->Q_Analog() - this->m_last_V[k]) * rez_cts ; - sq += nv * nv; - sqo += this->m_Vdelta[k] * this->m_Vdelta[k]; - this->m_Vdelta[k] = nv; - } - - // FIXME: used to be 1e90, but this would not be compatible with float - if (sqo > NL_FCONST(1e-20)) - m_lp_fact = std::min(std::sqrt(sq/sqo), (float_type) 2.0); - else - m_lp_fact = NL_FCONST(0.0); - } - - - return this->compute_next_timestep(); - } - #endif - - template <typename FT, int SIZE> - unsigned matrix_solver_SOR_mat_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson) + void matrix_solver_SOR_mat_t<FT, SIZE>::upstream_solve_non_dynamic() { - /* The matrix based code looks a lot nicer but actually is 30% slower than - * the optimized code which works directly on the data structures. - * Need something like that for gaussian elimination as well. - */ - + // The matrix based code looks a lot nicer but actually is 30% slower than + // the optimized code which works directly on the data structures. + // Need something like that for gaussian elimination as well. const std::size_t iN = this->size(); - this->build_LE_A(*this); - this->build_LE_RHS(*this); + this->clear_square_mat(this->m_A); + this->fill_matrix_and_rhs(); bool resched = false; unsigned resched_cnt = 0; - #if 0 static int ws_cnt = 0; ws_cnt++; @@ -142,16 +76,16 @@ namespace devices for (int k = 0; k < iN; k++) { #if 0 - float_type akk = std::abs(this->m_A[k][k]); + float_type akk = plib::abs(this->m_A[k][k]); if ( akk > lambdaN) lambdaN = akk; if (akk < lambda1) lambda1 = akk; #else - float_type akk = std::abs(this->m_A[k][k]); + float_type akk = plib::abs(this->m_A[k][k]); float_type s = 0.0; for (int i=0; i<iN; i++) - s = s + std::abs(this->m_A[k][i]); + s = s + plib::abs(this->m_A[k][i]); akk = s / akk - 1.0; if ( akk > lambdaN) lambdaN = akk; @@ -166,65 +100,58 @@ namespace devices #endif for (std::size_t k = 0; k < iN; k++) - this->m_new_V[k] = this->m_nets[k]->Q_Analog(); + this->m_new_V[k] = static_cast<float_type>(this->m_terms[k].getV()); do { resched = false; - float_type cerr = 0.0; + FT cerr = plib::constants<FT>::zero(); for (std::size_t k = 0; k < iN; k++) { float_type Idrive = 0; - const auto *p = this->m_terms[k]->m_nz.data(); - const std::size_t e = this->m_terms[k]->m_nz.size(); + const auto *p = this->m_terms[k].m_nz.data(); + const std::size_t e = this->m_terms[k].m_nz.size(); for (std::size_t i = 0; i < e; i++) - Idrive = Idrive + this->A(k,p[i]) * this->m_new_V[p[i]]; + Idrive = Idrive + this->m_A[k][p[i]] * this->m_new_V[p[i]]; - FT w = m_omega / this->A(k,k); + FT w = m_omega / this->m_A[k][k]; if (this->m_params.m_use_gabs) { - FT gabs_t = 0.0; + FT gabs_t = plib::constants<FT>::zero(); for (std::size_t i = 0; i < e; i++) if (p[i] != k) - gabs_t = gabs_t + std::abs(this->A(k,p[i])); + gabs_t = gabs_t + plib::abs(this->m_A[k][p[i]]); gabs_t *= plib::constants<FT>::one(); // derived by try and error - if (gabs_t > this->A(k,k)) + if (gabs_t > this->m_A[k][k]) { - w = plib::constants<FT>::one() / (this->A(k,k) + gabs_t); + w = plib::constants<FT>::one() / (this->m_A[k][k] + gabs_t); } } - const float_type delta = w * (this->RHS(k) - Idrive) ; - cerr = std::max(cerr, std::abs(delta)); + const float_type delta = w * (this->m_RHS[k] - Idrive) ; + cerr = std::max(cerr, plib::abs(delta)); this->m_new_V[k] += delta; } - if (cerr > this->m_params.m_accuracy) + if (cerr > static_cast<float_type>(this->m_params.m_accuracy)) { resched = true; } resched_cnt++; } while (resched && (resched_cnt < this->m_params.m_gs_loops)); - this->m_stat_calculations++; this->m_iterative_total += resched_cnt; if (resched) { this->m_iterative_fail++; - return matrix_solver_direct_t<FT, SIZE>::solve_non_dynamic(newton_raphson); + matrix_solver_direct_t<FT, SIZE>::solve_non_dynamic(); } - - const float_type err = (newton_raphson ? this->delta(this->m_new_V) : 0.0); - this->store(this->m_new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; - } -} // namespace devices -} // namespace netlist +} // namespace netlist::solver -#endif /* NLD_MS_GAUSS_SEIDEL_H_ */ +#endif // NLD_MS_SOR_MAT_ diff --git a/src/lib/netlist/solver/nld_ms_w.h b/src/lib/netlist/solver/nld_ms_w.h index a01fc239938..b5b58f4ebe7 100644 --- a/src/lib/netlist/solver/nld_ms_w.h +++ b/src/lib/netlist/solver/nld_ms_w.h @@ -1,379 +1,360 @@ -// license:GPL-2.0+ +// license:BSD-3-Clause // copyright-holders:Couriersud -/* - * nld_ms_direct.h - * - * - * Woodbury Solver - * - * Computes the updated solution of A given that the change in A is - * - * A <- A + (U x transpose(V)) U,V matrices - * - * The approach is describes in "Numerical Recipes in C", Second edition, Page 75ff - * - * Whilst the book proposes to invert the matrix R=(I+transpose(V)*Z) we define - * - * w = transpose(V)*y - * a = R^-1 * w - * - * and consequently - * - * R * a = w - * - * And solve for a using Gaussian elimination. This is a lot faster. - * - * One fact omitted in the book is the fact that actually the matrix Z which contains - * in it's columns the solutions of - * - * A * zk = uk - * - * for uk being unit vectors for full rank (max(k) == n) is identical to the - * inverse of A. - * - * The approach performs relatively well for matrices up to n ~ 40 (kidniki using frontiers). - * Kidniki without frontiers has n==88. Here, the average number of Newton-Raphson - * loops increase to 20. It looks like that the approach for larger matrices - * introduces numerical instability. - */ #ifndef NLD_MS_W_H_ #define NLD_MS_W_H_ -#include "nld_matrix_solver.h" -#include "nld_solver.h" +// Names +// spell-checker: words Woodbury, Raphson, +// +// Specific technical terms +// spell-checker: words Cgso, Cgdo, Cgbo, Cjsw, Mjsw, Ucrit, Uexp, Utra, Neff, Tnom, capval, Udsat, Utst + + +/// +/// \file nld_ms_direct.h +/// +/// Woodbury Solver +/// +/// Computes the updated solution of A given that the change in A is +/// +/// A <- A + (U x transpose(V)) U,V matrices +/// +/// The approach is describes in "Numerical Recipes in C", Second edition, Page 75ff +/// +/// Whilst the book proposes to invert the matrix R=(I+transpose(V)*Z) we define +/// +/// w = transpose(V)*y +/// a = R^-1 * w +/// +/// and consequently +/// +/// R * a = w +/// +/// And solve for a using Gaussian elimination. This is a lot faster. +/// +/// One fact omitted in the book is the fact that actually the matrix Z which contains +/// in it's columns the solutions of +/// +/// A * zk = uk +/// +/// for uk being unit vectors for full rank (max(k) == n) is identical to the +/// inverse of A. +/// +/// The approach performs relatively well for matrices up to n ~ 40 (`kidniki` using frontiers). +/// `Kidniki` without frontiers has n==88. Here, the average number of Newton-Raphson +/// loops increase to 20. It looks like that the approach for larger matrices +/// introduces numerical instability. +/// + +#include "nld_matrix_solver_ext.h" #include "plib/vector_ops.h" #include <algorithm> -namespace netlist +namespace netlist::solver { - namespace devices - { - -template <typename FT, int SIZE> -class matrix_solver_w_t: public matrix_solver_t -{ - friend class matrix_solver_t; - -public: - using float_ext_type = FT; - using float_type = FT; - - // FIXME: dirty hack to make this compile - static constexpr const std::size_t storage_N = 100; - - matrix_solver_w_t(netlist_state_t &anetlist, const pstring &name, const solver_parameters_t *params, const std::size_t size); - - void vsetup(analog_net_t::list_t &nets) override; - void reset() override { matrix_solver_t::reset(); } - -protected: - unsigned vsolve_non_dynamic(const bool newton_raphson) override; - unsigned solve_non_dynamic(const bool newton_raphson); - - constexpr std::size_t size() const { return m_dim; } - - void LE_invert(); - - template <typename T> - void LE_compute_x(T & x); - - - template <typename T1, typename T2> - float_ext_type &A(const T1 &r, const T2 &c) { return m_A[r][c]; } - template <typename T1, typename T2> - float_ext_type &W(const T1 &r, const T2 &c) { return m_W[r][c]; } - - /* access to Ainv for fixed columns over row, there store transposed */ - template <typename T1, typename T2> - float_ext_type &Ainv(const T1 &r, const T2 &c) { return m_Ainv[c][r]; } - template <typename T1> - float_ext_type &RHS(const T1 &r) { return m_RHS[r]; } + template <typename FT, int SIZE> + class matrix_solver_w_t: public matrix_solver_ext_t<FT, SIZE> + { + public: + using float_ext_type = FT; + using float_type = FT; + + // FIXME: dirty hack to make this compile + static constexpr const std::size_t storage_N = 100; + + matrix_solver_w_t(devices::nld_solver &main_solver, const pstring &name, + const matrix_solver_t::net_list_t &nets, + const solver_parameters_t *params, const std::size_t size) + : matrix_solver_ext_t<FT, SIZE>(main_solver, name, nets, params, size) + , m_cnt(0) + { + this->build_mat_ptr(m_A); + } - template <typename T1, typename T2> - float_ext_type &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; } + void reset() override { matrix_solver_t::reset(); } + protected: + void upstream_solve_non_dynamic() override; -private: - static constexpr std::size_t m_pitch = ((( storage_N) + 7) / 8) * 8; - float_ext_type m_A[storage_N][m_pitch]; - float_ext_type m_Ainv[storage_N][m_pitch]; - float_ext_type m_W[storage_N][m_pitch]; - std::array<float_ext_type, storage_N> m_RHS; // right hand side - contains currents + void LE_invert(); - float_ext_type m_lA[storage_N][m_pitch]; + template <typename T> + void LE_compute_x(T & x); - /* temporary */ - float_type H[storage_N][m_pitch] ; - std::array<unsigned, storage_N> rows; - unsigned cols[storage_N][m_pitch]; - std::array<unsigned, storage_N> colcount; + template <typename T1, typename T2> + float_ext_type &A(const T1 &r, const T2 &c) { return m_A[r][c]; } + template <typename T1, typename T2> + float_ext_type &W(const T1 &r, const T2 &c) { return m_W[r][c]; } - unsigned m_cnt; + // access to the inverted matrix for fixed columns over row, values stored transposed + template <typename T1, typename T2> + float_ext_type &Ainv(const T1 &r, const T2 &c) { return m_Ainv[c][r]; } + template <typename T1> + float_ext_type &RHS(const T1 &r) { return this->m_RHS[r]; } - //float_ext_type m_RHSx[storage_N]; - const std::size_t m_dim; + template <typename T1, typename T2> + float_ext_type &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; } -}; -// ---------------------------------------------------------------------------------------- -// matrix_solver_direct -// ---------------------------------------------------------------------------------------- + private: + void solve_non_dynamic(); -template <typename FT, int SIZE> -void matrix_solver_w_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets) -{ - matrix_solver_t::setup_base(nets); + template <typename T, std::size_t N, std::size_t M> + using array2D = std::array<std::array<T, M>, N>; + static constexpr std::size_t m_pitch = ((( storage_N) + 7) / 8) * 8; + array2D<float_ext_type, storage_N, m_pitch> m_A; + array2D<float_ext_type, storage_N, m_pitch> m_Ainv; + array2D<float_ext_type, storage_N, m_pitch> m_W; - // FIXME: This shouldn't be necessary, recalculate on each entry ... - for (std::size_t k = 0; k < size(); k++) - state().save(*this, RHS(k), this->name(), plib::pfmt("RHS.{1}")(k)); -} + array2D<float_ext_type, storage_N, m_pitch> m_lA; + // temporary + array2D<float_ext_type, storage_N, m_pitch> H; + std::array<unsigned, storage_N> rows; + array2D<unsigned, storage_N, m_pitch> cols; + std::array<unsigned, storage_N> col_count; + unsigned m_cnt; + }; -template <typename FT, int SIZE> -void matrix_solver_w_t<FT, SIZE>::LE_invert() -{ - const std::size_t kN = size(); + // ---------------------------------------------------------------------------------------- + // matrix_solver_direct + // ---------------------------------------------------------------------------------------- - for (std::size_t i = 0; i < kN; i++) + template <typename FT, int SIZE> + void matrix_solver_w_t<FT, SIZE>::LE_invert() { - for (std::size_t j = 0; j < kN; j++) - { - W(i,j) = lA(i,j) = A(i,j); - Ainv(i,j) = 0.0; - } - Ainv(i,i) = 1.0; - } - /* down */ - for (std::size_t i = 0; i < kN; i++) - { - /* FIXME: Singular matrix? */ - const float_type f = 1.0 / W(i,i); - const auto * const p = m_terms[i]->m_nzrd.data(); - const size_t e = m_terms[i]->m_nzrd.size(); - - /* Eliminate column i from row j */ + const std::size_t kN = this->size(); - const auto * const pb = m_terms[i]->m_nzbd.data(); - const size_t eb = m_terms[i]->m_nzbd.size(); - for (std::size_t jb = 0; jb < eb; jb++) + for (std::size_t i = 0; i < kN; i++) { - const auto j = pb[jb]; - const float_type f1 = - W(j,i) * f; - if (f1 != 0.0) + for (std::size_t j = 0; j < kN; j++) { - for (std::size_t k = 0; k < e; k++) - W(j,p[k]) += W(i,p[k]) * f1; - for (std::size_t k = 0; k <= i; k ++) - Ainv(j,k) += Ainv(i,k) * f1; + W(i,j) = lA(i,j) = A(i,j); + Ainv(i,j) = plib::constants<FT>::zero(); } + Ainv(i,i) = plib::constants<FT>::one(); } - } - /* up */ - for (std::size_t i = kN; i-- > 0; ) - { - /* FIXME: Singular matrix? */ - const float_type f = 1.0 / W(i,i); - for (std::size_t j = i; j-- > 0; ) + // down + for (std::size_t i = 0; i < kN; i++) { - const float_type f1 = - W(j,i) * f; - if (f1 != 0.0) + // FIXME: Singular matrix? + const float_type f = plib::reciprocal(W(i,i)); + const auto * const p = this->m_terms[i].m_nzrd.data(); + const size_t e = this->m_terms[i].m_nzrd.size(); + + // Eliminate column i from row j + + const auto * const pb = this->m_terms[i].m_nzbd.data(); + const size_t eb = this->m_terms[i].m_nzbd.size(); + for (std::size_t jb = 0; jb < eb; jb++) { - for (std::size_t k = i; k < kN; k++) - W(j,k) += W(i,k) * f1; - for (std::size_t k = 0; k < kN; k++) - Ainv(j,k) += Ainv(i,k) * f1; + const auto j = pb[jb]; + const float_type f1 = - W(j,i) * f; + // FIXME: comparison to zero + if (f1 != plib::constants<float_type>::zero()) + { + for (std::size_t k = 0; k < e; k++) + W(j,p[k]) += W(i,p[k]) * f1; + for (std::size_t k = 0; k <= i; k ++) + Ainv(j,k) += Ainv(i,k) * f1; + } } } - for (std::size_t k = 0; k < kN; k++) + // up + for (std::size_t i = kN; i-- > 0; ) { - Ainv(i,k) *= f; + // FIXME: Singular matrix? + const float_type f = plib::reciprocal(W(i,i)); + for (std::size_t j = i; j-- > 0; ) + { + const float_type f1 = - W(j,i) * f; + // FIXME: comparison to zero + if (f1 != plib::constants<float_type>::zero()) + { + for (std::size_t k = i; k < kN; k++) + W(j,k) += W(i,k) * f1; + for (std::size_t k = 0; k < kN; k++) + Ainv(j,k) += Ainv(i,k) * f1; + } + } + for (std::size_t k = 0; k < kN; k++) + { + Ainv(i,k) *= f; + } } } -} - -template <typename FT, int SIZE> -template <typename T> -void matrix_solver_w_t<FT, SIZE>::LE_compute_x( - T & x) -{ - const std::size_t kN = size(); - - for (std::size_t i=0; i<kN; i++) - x[i] = 0.0; - for (std::size_t k=0; k<kN; k++) + template <typename FT, int SIZE> + template <typename T> + void matrix_solver_w_t<FT, SIZE>::LE_compute_x( + T & x) { - const float_type f = RHS(k); + const std::size_t kN = this->size(); for (std::size_t i=0; i<kN; i++) - x[i] += Ainv(i,k) * f; - } -} + x[i] = plib::constants<FT>::zero(); + for (std::size_t k=0; k<kN; k++) + { + const float_type f = RHS(k); -template <typename FT, int SIZE> -unsigned matrix_solver_w_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphson) -{ - const auto iN = size(); + for (std::size_t i=0; i<kN; i++) + x[i] += Ainv(i,k) * f; + } + } - std::array<float_type, storage_N> new_V; // = { 0.0 }; - if ((m_cnt % 50) == 0) + template <typename FT, int SIZE> + void matrix_solver_w_t<FT, SIZE>::solve_non_dynamic() { - /* complete calculation */ - this->LE_invert(); - this->LE_compute_x(new_V); - } - else - { - /* Solve Ay = b for y */ - this->LE_compute_x(new_V); - - /* determine changed rows */ + const auto iN = this->size(); - unsigned rowcount=0; - #define VT(r,c) (A(r,c) - lA(r,c)) + // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) + std::array<float_type, storage_N> t; // FIXME: convert to member + // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) + std::array<float_type, storage_N> w; - for (unsigned row = 0; row < iN; row ++) + if ((m_cnt % 50) == 0) { - unsigned cc=0; - auto &nz = m_terms[row]->m_nz; - for (auto & col : nz) - { - if (A(row,col) != lA(row,col)) - cols[rowcount][cc++] = col; - } - if (cc > 0) - { - colcount[rowcount] = cc; - rows[rowcount++] = row; - } + // complete calculation + this->LE_invert(); + this->LE_compute_x(this->m_new_V); } - if (rowcount > 0) + else { - /* construct w = transform(V) * y - * dim: rowcount x iN - * */ - std::array<float_type, storage_N> w; - for (unsigned i = 0; i < rowcount; i++) + // Solve Ay = b for y + this->LE_compute_x(this->m_new_V); + + // determine changed rows + + unsigned row_count=0; + #define VT(r,c) (A(r,c) - lA(r,c)) + + for (unsigned row = 0; row < iN; row ++) { - const unsigned r = rows[i]; - double tmp = 0.0; - for (unsigned k = 0; k < iN; k++) - tmp += VT(r,k) * new_V[k]; - w[i] = tmp; + unsigned cc=0; + auto &nz = this->m_terms[row].m_nz; + for (auto & col : nz) + { + if (A(row,col) != lA(row,col)) + cols[row_count][cc++] = col; + } + if (cc > 0) + { + col_count[row_count] = cc; + rows[row_count++] = row; + } } + if (row_count > 0) + { + // construct w = transform(V) * y + // dim: row_count x iN + // + for (unsigned i = 0; i < row_count; i++) + { + const unsigned r = rows[i]; + FT tmp = plib::constants<FT>::zero(); + for (unsigned k = 0; k < iN; k++) + tmp += VT(r,k) * this->m_new_V[k]; + w[i] = tmp; + } + + for (unsigned i = 0; i < row_count; i++) + for (unsigned k=0; k< row_count; k++) + H[i][k] = plib::constants<FT>::zero(); - for (unsigned i = 0; i < rowcount; i++) - for (unsigned k=0; k< rowcount; k++) - H[i][k] = 0.0; + for (unsigned i = 0; i < row_count; i++) + H[i][i] = plib::constants<FT>::one(); + // Construct H = (I + VT*Z) + for (unsigned i = 0; i < row_count; i++) + for (unsigned k=0; k< col_count[i]; k++) + { + const unsigned col = cols[i][k]; + float_type f = VT(rows[i],col); + // FIXME: comparison to zero + if (f != plib::constants<float_type>::zero()) + for (unsigned j= 0; j < row_count; j++) + H[i][j] += f * Ainv(col,rows[j]); + } - for (unsigned i = 0; i < rowcount; i++) - H[i][i] = 1.0; - /* Construct H = (I + VT*Z) */ - for (unsigned i = 0; i < rowcount; i++) - for (unsigned k=0; k< colcount[i]; k++) + // Gaussian elimination of H + for (unsigned i = 0; i < row_count; i++) { - const unsigned col = cols[i][k]; - float_type f = VT(rows[i],col); - if (f!=0.0) - for (unsigned j= 0; j < rowcount; j++) - H[i][j] += f * Ainv(col,rows[j]); + // FIXME: comparison to zero + if (H[i][i] == plib::constants<float_type>::zero()) + plib::perrlogger("{} H singular\n", this->name()); + const float_type f = plib::reciprocal(H[i][i]); + for (unsigned j = i+1; j < row_count; j++) + { + const float_type f1 = - f * H[j][i]; + + // FIXME: comparison to zero + if (f1 != plib::constants<float_type>::zero()) + { + float_type *pj = &H[j][i+1]; + const float_type *pi = &H[i][i+1]; + for (unsigned k = 0; k < row_count-i-1; k++) + pj[k] += f1 * pi[k]; + //H[j][k] += f1 * H[i][k]; + w[j] += f1 * w[i]; + } + } } - - /* Gaussian elimination of H */ - for (unsigned i = 0; i < rowcount; i++) - { - if (H[i][i] == 0.0) - plib::perrlogger("{} H singular\n", this->name()); - const float_type f = 1.0 / H[i][i]; - for (unsigned j = i+1; j < rowcount; j++) + // Back substitution + //inv(H) w = t w = H t + for (unsigned j = row_count; j-- > 0; ) { - const float_type f1 = - f * H[j][i]; + float_type tmp = 0; + const float_type *pj = &H[j][j+1]; + const float_type *tj = &t[j+1]; + for (unsigned k = 0; k < row_count-j-1; k++) + tmp += pj[k] * tj[k]; + //tmp += H[j][k] * t[k]; + t[j] = (w[j] - tmp) / H[j][j]; + } - if (f1!=0.0) + // x = y - Zt + for (unsigned i=0; i<iN; i++) + { + float_type tmp = plib::constants<FT>::zero(); + for (unsigned j=0; j<row_count;j++) { - float_type *pj = &H[j][i+1]; - const float_type *pi = &H[i][i+1]; - for (unsigned k = 0; k < rowcount-i-1; k++) - pj[k] += f1 * pi[k]; - //H[j][k] += f1 * H[i][k]; - w[j] += f1 * w[i]; + const unsigned row = rows[j]; + tmp += Ainv(i,row) * t[j]; } + this->m_new_V[i] -= tmp; } } - /* Back substitution */ - //inv(H) w = t w = H t - std::array<float_type, storage_N> t; // FIXME: convert to member - for (unsigned j = rowcount; j-- > 0; ) - { - float_type tmp = 0; - const float_type *pj = &H[j][j+1]; - const float_type *tj = &t[j+1]; - for (unsigned k = 0; k < rowcount-j-1; k++) - tmp += pj[k] * tj[k]; - //tmp += H[j][k] * t[k]; - t[j] = (w[j] - tmp) / H[j][j]; - } + } + m_cnt++; - /* x = y - Zt */ + if (false) // NOLINT for (unsigned i=0; i<iN; i++) { - float_type tmp = 0.0; - for (unsigned j=0; j<rowcount;j++) + float_type tmp = plib::constants<FT>::zero(); + for (unsigned j=0; j<iN; j++) { - const unsigned row = rows[j]; - tmp += Ainv(i,row) * t[j]; + tmp += A(i,j) * this->m_new_V[j]; } - new_V[i] -= tmp; + if (plib::abs(tmp-RHS(i)) > static_cast<float_type>(1e-6)) + plib::perrlogger("{} failed on row {}: {} RHS: {}\n", this->name(), i, plib::abs(tmp-RHS(i)), RHS(i)); } - } } - m_cnt++; - - if (false) - for (unsigned i=0; i<iN; i++) - { - float_type tmp = 0.0; - for (unsigned j=0; j<iN; j++) - { - tmp += A(i,j) * new_V[j]; - } - if (std::abs(tmp-RHS(i)) > 1e-6) - plib::perrlogger("{} failed on row {}: {} RHS: {}\n", this->name(), i, std::abs(tmp-RHS(i)), RHS(i)); - } - const float_type 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> + void matrix_solver_w_t<FT, SIZE>::upstream_solve_non_dynamic() + { + this->clear_square_mat(this->m_A); + this->fill_matrix_and_rhs(); -template <typename FT, int SIZE> -unsigned matrix_solver_w_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson) -{ - this->build_LE_A(*this); - this->build_LE_RHS(*this); - - this->m_stat_calculations++; - return this->solve_non_dynamic(newton_raphson); -} - -template <typename FT, int SIZE> -matrix_solver_w_t<FT, SIZE>::matrix_solver_w_t(netlist_state_t &anetlist, const pstring &name, - const solver_parameters_t *params, const std::size_t size) - : matrix_solver_t(anetlist, name, NOSORT, params) - , m_cnt(0) - , m_dim(size) -{ -} + this->solve_non_dynamic(); + } - } //namespace devices -} // namespace netlist +} // namespace netlist::solver -#endif /* NLD_MS_DIRECT_H_ */ +#endif // NLD_MS_DIRECT_H_ diff --git a/src/lib/netlist/solver/nld_solver.cpp b/src/lib/netlist/solver/nld_solver.cpp index 0aa7e6f584c..d7381cc813a 100644 --- a/src/lib/netlist/solver/nld_solver.cpp +++ b/src/lib/netlist/solver/nld_solver.cpp @@ -1,38 +1,14 @@ -// license:GPL-2.0+ +// license:BSD-3-Clause // copyright-holders:Couriersud -/* - * nld_solver.c - * - */ - -/* Commented out for now. Relatively low number of terminals / nets make - * the vectorizations fast-math enables pretty expensive - */ - -#if 0 -#pragma GCC optimize "-ftree-vectorize" -#pragma GCC optimize "-ffast-math" -#pragma GCC optimize "-funsafe-math-optimizations" -#pragma GCC optimize "-funroll-loops" -#pragma GCC optimize "-funswitch-loops" -#pragma GCC optimize "-fstrict-aliasing" -#pragma GCC optimize "tree-vectorizer-verbose=7" -#pragma GCC optimize "opt-info-vec" -#pragma GCC optimize "opt-info-vec-missed" -//#pragma GCC optimize "tree-parallelize-loops=4" -#pragma GCC optimize "variable-expansion-in-unroller" -#pragma GCC optimize "unsafe-loop-optimizations" -#pragma GCC optimize "vect-cost-model" -#pragma GCC optimize "variable-expansion-in-unroller" -#pragma GCC optimize "tree-loop-if-convert-stores" -#pragma GCC optimize "tree-loop-distribution" -#pragma GCC optimize "tree-loop-im" -#pragma GCC optimize "tree-loop-ivcanon" -#pragma GCC optimize "ivopts" -#endif -#include "netlist/nl_lists.h" -#include "netlist/nl_factory.h" +// 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" @@ -43,25 +19,48 @@ #include "nld_ms_sor.h" #include "nld_ms_sor_mat.h" #include "nld_ms_w.h" -#include "nld_solver.h" + +#include "core/setup.h" + #include "plib/pomp.h" +#include "plib/ptimed_queue.h" #include <algorithm> -#include <cmath> +#include <type_traits> -namespace netlist -{ -namespace devices +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<false>)) + , 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<plib::state_manager_t::callback_t &>(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<true>)); for (auto &s : m_mat_solvers) s->reset(); + for (auto &s : m_mat_solvers) + m_queue.push<false>({netlist_time_ext::zero(), s.get()}); } void NETLIB_NAME(solver)::stop() @@ -70,352 +69,541 @@ namespace devices s->log_stats(); } - NETLIB_UPDATE(solver) +#if 1 + + template <bool KEEP_STATS> + 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<std::size_t>( + 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<solver::matrix_solver_t *, + config::max_solver_queue_size::value> + tmp; // NOLINT + plib::uninitialised_array<netlist_time, + config::max_solver_queue_size::value> + 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<false>({now + nt[i], tmp[i]}); + tmp[i]->update_inputs(); + } + } + else + { + plib::omp::set_num_threads(nthreads); + plib::omp::for_static(static_cast<std::size_t>(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<false>({now + nt[i], tmp[i]}); + tmp[i]->update_inputs(); + } + } + if (!m_queue.empty()) + m_Q_step.net().toggle_and_push_to_queue( + static_cast<netlist_time>(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<false>(solv); + m_queue.push<false>({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 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::from_double(2 * m_params.m_max_timestep)); + 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<std::size_t>(m_parallel()), plib::omp::get_max_threads()); + std::size_t nthreads = std::min( + static_cast<std::size_t>(m_params.m_parallel()), + plib::omp::get_max_threads()); - std::vector<matrix_solver_t *> &solvers = (force_solve ? m_mat_solvers_all : m_mat_solvers_timestepping); + std::vector<solver_entry *> + &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<std::size_t>(0), solvers.size(), [&solvers, now](std::size_t i) - { - const netlist_time ts = solvers[i]->solve(now); - plib::unused_var(ts); + plib::omp::for_static( + static_cast<std::size_t>(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) + for (auto &solver : solvers) { - const netlist_time ts = solver->solve(now); - plib::unused_var(ts); + [[maybe_unused]] const netlist_time ts = solver->ptr->solve( + now); } - for (auto & solver : solvers) - solver->update_inputs(); + for (auto &solver : solvers) + solver->ptr->update_inputs(); - /* step circuit */ + // step circuit if (!m_Q_step.net().is_queued()) { - m_Q_step.net().toggle_and_push_to_queue(netlist_time::from_double(m_params.m_max_timestep)); + m_Q_step.net().toggle_and_push_to_queue( + netlist_time::from_fp(m_params.m_max_time_step)); } } +#endif - template <class C> - pool_owned_ptr<matrix_solver_t> create_it(netlist_state_t &nl, pstring name, solver_parameters_t ¶ms, std::size_t size) + // FIXME: should be created in device space + template <class C, class A> + 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 pool().make_poolptr<C>(nl, name, ¶ms, size); + return plib::make_unique<C>(arena, main_solver, name, nets, params, + size); } template <typename FT, int SIZE> - pool_owned_ptr<matrix_solver_t> NETLIB_NAME(solver)::create_solver(std::size_t size, const pstring &solvername) + 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) { - if (m_method() == "SOR_MAT") - { - return create_it<matrix_solver_SOR_mat_t<FT, SIZE>>(state(), solvername, m_params, size); - //typedef matrix_solver_SOR_mat_t<m_N,storage_N> solver_sor_mat; - //return plib::make_unique<solver_sor_mat>(state(), solvername, &m_params, size); - } - else if (m_method() == "MAT_CR") - { - if (size > 0) // GCR always outperforms MAT solver - { - return create_it<matrix_solver_GCR_t<FT, SIZE>>(state(), solvername, m_params, size); - } - else - { - return create_it<matrix_solver_direct_t<FT, SIZE>>(state(), solvername, m_params, size); - } - } - else if (m_method() == "MAT") - { - return create_it<matrix_solver_direct_t<FT, SIZE>>(state(), solvername, m_params, size); - } - else if (m_method() == "SM") - { - /* Sherman-Morrison Formula */ - return create_it<matrix_solver_sm_t<FT, SIZE>>(state(), solvername, m_params, size); - } - else if (m_method() == "W") - { - /* Woodbury Formula */ - return create_it<matrix_solver_w_t<FT, SIZE>>(state(), solvername, m_params, size); - } - else if (m_method() == "SOR") - { - return create_it<matrix_solver_SOR_t<FT, SIZE>>(state(), solvername, m_params, size); - } - else if (m_method() == "GMRES") - { - return create_it<matrix_solver_GMRES_t<FT, SIZE>>(state(), solvername, m_params, size); - } - else + switch (params->m_method()) { - log().fatal(MF_UNKNOWN_SOLVER_TYPE(m_method())); - return pool_owned_ptr<matrix_solver_t>(); + case solver::matrix_type_e::MAT_CR: + return create_it<solver::matrix_solver_GCR_t<FT, SIZE>>( + state().pool(), *this, solver_name, nets, params, size); + case solver::matrix_type_e::MAT: + return create_it<solver::matrix_solver_direct_t<FT, SIZE>>( + state().pool(), *this, solver_name, nets, params, size); + case solver::matrix_type_e::GMRES: + return create_it<solver::matrix_solver_GMRES_t<FT, SIZE>>( + state().pool(), *this, solver_name, nets, params, size); +#if (NL_USE_ACADEMIC_SOLVERS) + case solver::matrix_type_e::SOR: + return create_it<solver::matrix_solver_SOR_t<FT, SIZE>>( + state().pool(), *this, solver_name, nets, params, size); + case solver::matrix_type_e::SOR_MAT: + return create_it<solver::matrix_solver_SOR_mat_t<FT, SIZE>>( + state().pool(), *this, solver_name, nets, params, size); + case solver::matrix_type_e::SM: + // Sherman-Morrison Formula + return create_it<solver::matrix_solver_sm_t<FT, SIZE>>( + state().pool(), *this, solver_name, nets, params, size); + case solver::matrix_type_e::W: + // Woodbury Formula + return create_it<solver::matrix_solver_w_t<FT, SIZE>>( + 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<solver::matrix_solver_GCR_t<FT, SIZE>>( + state().pool(), *this, solver_name, nets, params, size); +#endif } + return solver_ptr(); } - template <typename FT, int SIZE> - pool_owned_ptr<matrix_solver_t> NETLIB_NAME(solver)::create_solver_x(std::size_t size, const pstring &solvername) + template <typename FT> + NETLIB_NAME(solver)::solver_ptr NETLIB_NAME(solver)::create_solvers( + const pstring &sname, const solver::solver_parameters_t *params, + net_list_t &nets) { - if (SIZE > 0) - { - if (size == SIZE) - return create_solver<FT, SIZE>(size, solvername); - else - return this->create_solver_x<FT, SIZE-1>(size, solvername); - } - else + std::size_t net_count = nets.size(); + switch (net_count) { - if (size * 2 > -SIZE ) - return create_solver<FT, SIZE>(size, solvername); - else - return this->create_solver_x<FT, SIZE / 2>(size, solvername); +#if !defined(__EMSCRIPTEN__) + case 1: + return plib::make_unique<solver::matrix_solver_direct1_t<FT>>( + state().pool(), *this, sname, nets, params); + case 2: + return plib::make_unique<solver::matrix_solver_direct2_t<FT>>( + state().pool(), *this, sname, nets, params); + case 3: return create_solver<FT, 3>(3, sname, params, nets); + case 4: return create_solver<FT, 4>(4, sname, params, nets); + case 5: return create_solver<FT, 5>(5, sname, params, nets); + case 6: return create_solver<FT, 6>(6, sname, params, nets); + case 7: return create_solver<FT, 7>(7, sname, params, nets); + case 8: return create_solver<FT, 8>(8, sname, params, nets); +#endif + default: + log().info(MI_NO_SPECIFIC_SOLVER(net_count)); + if (net_count <= 16) + { + return create_solver<FT, -16>(net_count, sname, params, + nets); + } + if (net_count <= 32) + { + return create_solver<FT, -32>(net_count, sname, params, + nets); + } + if (net_count <= 64) + { + return create_solver<FT, -64>(net_count, sname, params, + nets); + } + if (net_count <= 128) + { + return create_solver<FT, -128>(net_count, sname, params, + nets); + } + if (net_count <= 256) + { + return create_solver<FT, -256>(net_count, sname, params, + nets); + } + if (net_count <= 512) + { + return create_solver<FT, -512>(net_count, sname, params, + nets); + } + return create_solver<FT, 0>(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<analog_net_t *>(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<NETLIB_NAME(solver)::net_list_t> groups; + private: bool already_processed(const analog_net_t &n) const { - /* no need to process rail nets - these are known variables */ - if (n.isRailNet()) + // 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 (auto & grp : groups) + // if it's already processed - no need to continue + for (const auto &grp : groups) if (plib::container::contains(grp, &n)) return true; return false; } - void process_net(analog_net_t &n) + bool check_if_processed_and_join(const analog_net_t &n) { - /* ignore empty nets. FIXME: print a warning message */ - if (n.num_cons() == 0) - return; - /* add the net */ - groups.back().push_back(&n); - /* process all terminals connected to this net */ - for (auto &term : n.core_terms()) + // 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) { - /* only process analog terminals */ - if (term->is_type(detail::terminal_type::TERMINAL)) - { - auto *pt = static_cast<terminal_t *>(term); - /* check the connected terminal */ - analog_net_t &connected_net = pt->connected_terminal()->net(); - if (!already_processed(connected_net)) - process_net(connected_net); - } + 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; } - void run(netlist_state_t &netlist) + // NOLINTNEXTLINE(misc-no-recursion) + void process_net(netlist_state_t &nl_state, analog_net_t &n) { - for (auto & net : netlist.nets()) + // ignore empty nets. FIXME: print a warning message + nl_state.log().verbose("Net {}", n.name()); + auto terminals(n.core_terms_copy()); + + if (!terminals.empty()) { - netlist.log().debug("processing {1}\n", net->name()); - if (!net->isRailNet() && net->num_cons() > 0) + // add the net + groupspre.back().push_back(&n); + // process all terminals connected to this net + for (detail::core_terminal_t *term : terminals) { - netlist.log().debug(" ==> not a rail net\n"); - /* Must be an analog net */ - auto &n = *static_cast<analog_net_t *>(net.get()); - if (!already_processed(n)) + nl_state.log().verbose("Term {} {}", term->name(), + static_cast<int>(term->type())); + // only process analog terminals + if (term->is_type(detail::terminal_type::TERMINAL)) { - groups.emplace_back(analog_net_t::list_t()); - process_net(n); + auto pt = plib::dynamic_downcast<terminal_t *>(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<analog_net_t::list_t> groups; + std::vector<NETLIB_NAME(solver)::net_list_t> groupspre; }; void NETLIB_NAME(solver)::post_start() { - m_params.m_pivot = m_pivot(); - m_params.m_accuracy = m_accuracy(); - /* FIXME: Throw when negative */ - m_params.m_gs_loops = static_cast<unsigned>(m_gs_loops()); - m_params.m_nr_loops = static_cast<unsigned>(m_nr_loops()); - m_params.m_nr_recalc_delay = netlist_time::from_double(m_nr_recalc_delay()); - m_params.m_dynamic_lte = m_dynamic_lte(); - m_params.m_gs_sor = m_gs_sor(); - - m_params.m_min_timestep = m_dynamic_min_ts(); - m_params.m_dynamic_ts = (m_dynamic_ts() == 1 ? true : false); - m_params.m_max_timestep = netlist_time::from_double(1.0 / m_freq()).as_double(); - - m_params.m_use_gabs = m_use_gabs(); - m_params.m_use_linear_prediction = m_use_linear_prediction(); - - - if (m_params.m_dynamic_ts) - { - m_params.m_max_timestep *= 1;//NL_FCONST(1000.0); - } - else - { - m_params.m_min_timestep = m_params.m_max_timestep; - } - - //m_params.m_max_timestep = std::max(m_params.m_max_timestep, m_params.m_max_timestep::) - 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()); - // setup the solvers - log().verbose("Found {1} net groups in {2} nets\n", splitter.groups.size(), state().nets().size()); - for (auto & grp : splitter.groups) - { - pool_owned_ptr<matrix_solver_t> ms; - std::size_t net_count = grp.size(); - pstring sname = plib::pfmt("Solver_{1}")(m_mat_solvers.size()); + int num_errors = 0; - switch (net_count) + log().verbose("checking net consistency ..."); + for (const auto &grp : splitter.groups) + { + int rail_terminals = 0; + pstring nets_in_grp; + for (const auto &n : grp) { - #if 1 - case 1: - ms = pool().make_poolptr<matrix_solver_direct1_t<double>>(state(), sname, &m_params); - break; - case 2: - ms = pool().make_poolptr<matrix_solver_direct2_t<double>>(state(), sname, &m_params); - break; - case 3: - ms = create_solver<double, 3>(3, sname); - break; - case 4: - ms = create_solver<double, 4>(4, sname); - break; - case 5: - ms = create_solver<double, 5>(5, sname); - break; - case 6: - ms = create_solver<double, 6>(6, sname); - break; - case 7: - ms = create_solver<double, 7>(7, sname); - break; - case 8: - ms = create_solver<double, 8>(8, sname); - break; - case 9: - ms = create_solver<double, 9>(9, sname); - break; - case 10: - ms = create_solver<double, 10>(10, sname); - break; - #if 0 - case 11: - ms = create_solver<double, 11>(11, sname); - break; - case 12: - ms = create_solver<double, 12>(12, sname); - break; - case 15: - ms = create_solver<double, 15>(15, sname); - break; - case 31: - ms = create_solver<double, 31>(31, sname); - break; - case 35: - ms = create_solver<double, 35>(35, sname); - break; - case 43: - ms = create_solver<double, 43>(43, sname); - break; - case 49: - ms = create_solver<double, 49>(49, sname); - break; - #endif - #if 1 - case 86: - ms = create_solver<double,86>(86, sname); - break; - #endif - #endif - default: - log().info(MI_NO_SPECIFIC_SOLVER(net_count)); - if (net_count <= 8) - { - ms = create_solver<double, -8>(net_count, sname); - } - else if (net_count <= 16) - { - ms = create_solver<double, -16>(net_count, sname); - } - else if (net_count <= 32) - { - ms = create_solver<double, -32>(net_count, sname); - } - else - if (net_count <= 64) - { - ms = create_solver<double, -64>(net_count, sname); - } - else - if (net_count <= 128) + 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()) { - ms = create_solver<double, -128>(net_count, sname); + state().log().error( + ME_SOLVER_TERMINAL_NO_NET(t->name())); + num_errors++; } else { - log().fatal(MF_NETGROUP_SIZE_EXCEEDED_1(128)); - return; /* tease compilers */ + 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<solver::solver_parameters_t>( + 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<std::conditional_t< + config::use_float_matrix::value, float, double>>( + sname, params.get(), grp); + break; + case solver::matrix_fp_type_e::DOUBLE: + ms = create_solvers<double>(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<std::conditional_t< + config::use_long_double_matrix::value, long double, + double>>(sname, params.get(), grp); + break; + case solver::matrix_fp_type_e::FLOATQ128: +#if (NL_USE_FLOAT128) + ms = create_solvers<FLOAT128>(sname, params.get(), grp); +#else + log().info("FPTYPE {1} not supported. Using DOUBLE", + params->m_fp_type().name()); + ms = create_solvers<double>(sname, params.get(), grp); +#endif break; } - // FIXME ... - ms->setup(grp); + state().register_device( + ms->name(), + device_arena::owned_ptr<core_device_t>(ms.get(), false)); log().verbose("Solver {1}", ms->name()); log().verbose(" ==> {1} nets", grp.size()); - log().verbose(" has {1} elements", ms->has_dynamic_devices() ? "dynamic" : "no dynamic"); - log().verbose(" has {1} elements", ms->has_timestep_devices() ? "timestep" : "no timestep"); + 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 auto &pcore : n->core_terms()) + for (const detail::core_terminal_t *t : n->core_terms_copy()) { - log().verbose(" {1}", pcore->name()); + log().verbose(" {1}", t->name()); } } - m_mat_solvers_all.push_back(ms.get()); - if (ms->has_timestep_devices()) - m_mat_solvers_timestepping.push_back(ms.get()); - - m_mat_solvers.emplace_back(std::move(ms)); + m_mat_params.push_back(std::move(params)); + m_mat_solvers.push_back(std::move(ms)); } } - void NETLIB_NAME(solver)::create_solver_code(std::map<pstring, pstring> &mp) + solver::static_compile_container NETLIB_NAME(solver)::create_solver_code( + solver::static_compile_target target) { - for (auto & s : m_mat_solvers) + solver::static_compile_container mp; + for (auto &s : m_mat_solvers) { - auto r = s->create_solver_code(); - mp[r.first] = r.second; // automatically overwrites identical names + 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<std::size_t>::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 devices -} // namespace netlist +} // namespace netlist::devices diff --git a/src/lib/netlist/solver/nld_solver.h b/src/lib/netlist/solver/nld_solver.h index 99937f72812..c4f6c2e5bf1 100644 --- a/src/lib/netlist/solver/nld_solver.h +++ b/src/lib/netlist/solver/nld_solver.h @@ -1,118 +1,91 @@ -// license:GPL-2.0+ +// license:BSD-3-Clause // copyright-holders:Couriersud -/* - * nld_solver.h - * - */ #ifndef NLD_SOLVER_H_ #define NLD_SOLVER_H_ -#include "netlist/nl_base.h" +/// +/// \file nld_solver.h +/// + #include "nld_matrix_solver.h" -#include "plib/pstream.h" + +#include "core/core_device.h" +#include "core/logic.h" +#include "core/state_var.h" + +#include "../plib/pstream.h" #include <map> #include <memory> #include <vector> -//#define ATTR_ALIGNED(N) __attribute__((aligned(N))) -#define ATTR_ALIGNED(N) ATTR_ALIGN - -// ---------------------------------------------------------------------------------------- -// solver -// ---------------------------------------------------------------------------------------- - -namespace netlist +namespace netlist::devices { -namespace devices -{ - class NETLIB_NAME(solver); - - class matrix_solver_t; + // ------------------------------------------------------------------------- + // solver + // ------------------------------------------------------------------------- - NETLIB_OBJECT(solver) + class nld_solver : public device_t { - NETLIB_CONSTRUCTOR(solver) - , m_fb_step(*this, "FB_step") - , m_Q_step(*this, "Q_step") - , m_freq(*this, "FREQ", 48000.0) - - /* iteration parameters */ - , m_gs_sor(*this, "SOR_FACTOR", 1.059) - , m_method(*this, "METHOD", "MAT_CR") - , m_accuracy(*this, "ACCURACY", 1e-7) - , m_gs_loops(*this, "GS_LOOPS", 9) // Gauss-Seidel loops - - /* general parameters */ - , m_gmin(*this, "GMIN", 1e-9) - , m_pivot(*this, "PIVOT", false) // use pivoting - on supported solvers - , m_nr_loops(*this, "NR_LOOPS", 250) // Newton-Raphson loops - , m_nr_recalc_delay(*this, "NR_RECALC_DELAY", netlist_time::quantum().as_double()) // Delay to next solve attempt if nr loops exceeded - , m_parallel(*this, "PARALLEL", 0) - - /* automatic time step */ - , m_dynamic_ts(*this, "DYNAMIC_TS", false) - , m_dynamic_lte(*this, "DYNAMIC_LTE", 1e-5) // diff/timestep - , m_dynamic_min_ts(*this, "DYNAMIC_MIN_TIMESTEP", 1e-6) // nl_double timestep resolution - - /* special */ - , m_use_gabs(*this, "USE_GABS", true) - , m_use_linear_prediction(*this, "USE_LINEAR_PREDICTION", false) // // savings are eaten up by effort - - , m_params() - { - // internal staff + public: + using solver_arena = device_arena; + using queue_type = detail::queue_base<solver_arena, + solver::matrix_solver_t>; - connect(m_fb_step, m_Q_step); - } + nld_solver(constructor_param_t data); void post_start(); void stop(); - nl_double gmin() const { return m_gmin(); } + auto gmin() const -> decltype(solver::solver_parameters_t::m_gmin()) + { + return m_params.m_gmin(); + } - void create_solver_code(std::map<pstring, pstring> &mp); + solver::static_compile_container + create_solver_code(solver::static_compile_target target); - NETLIB_UPDATEI(); NETLIB_RESETI(); // NETLIB_UPDATE_PARAMI(); + using solver_ptr = solver_arena::unique_ptr<solver::matrix_solver_t>; + + using net_list_t = solver::matrix_solver_t::net_list_t; + + void reschedule(solver::matrix_solver_t *solv, netlist_time ts); + private: - logic_input_t m_fb_step; + using params_uptr = solver_arena::unique_ptr< + solver::solver_parameters_t>; + + template <bool KEEP_STATS> + NETLIB_HANDLERI(fb_step); + + logic_input_t m_fb_step; logic_output_t m_Q_step; - param_double_t m_freq; - param_double_t m_gs_sor; - param_str_t m_method; - param_double_t m_accuracy; - param_int_t m_gs_loops; - param_double_t m_gmin; - param_logic_t m_pivot; - param_int_t m_nr_loops; - param_double_t m_nr_recalc_delay; - param_int_t m_parallel; - param_logic_t m_dynamic_ts; - param_double_t m_dynamic_lte; - param_double_t m_dynamic_min_ts; - - param_logic_t m_use_gabs; - param_logic_t m_use_linear_prediction; - - std::vector<pool_owned_ptr<matrix_solver_t>> m_mat_solvers; - std::vector<matrix_solver_t *> m_mat_solvers_all; - std::vector<matrix_solver_t *> m_mat_solvers_timestepping; - - solver_parameters_t m_params; + // FIXME: these should be created in device space + std::vector<params_uptr> m_mat_params; + std::vector<solver_ptr> m_mat_solvers; - template <typename FT, int SIZE> - pool_owned_ptr<matrix_solver_t> create_solver(std::size_t size, const pstring &solvername); + solver::solver_parameters_t m_params; + queue_type m_queue; template <typename FT, int SIZE> - pool_owned_ptr<matrix_solver_t> create_solver_x(std::size_t size, const pstring &solvername); + solver_ptr create_solver(std::size_t size, const pstring &solver_name, + const solver::solver_parameters_t *params, + net_list_t & nets); + + template <typename FT> + solver_ptr create_solvers(const pstring & sname, + const solver::solver_parameters_t *params, + net_list_t & nets); + + std::size_t get_solver_id(const solver::matrix_solver_t *net) const; + solver::matrix_solver_t *solver_by_id(std::size_t id) const; }; -} //namespace devices -} // namespace netlist +} // namespace netlist::devices -#endif /* NLD_SOLVER_H_ */ +#endif // NLD_SOLVER_H_ |