summaryrefslogtreecommitdiffstatshomepage
path: root/src/lib/netlist/solver
diff options
context:
space:
mode:
Diffstat (limited to 'src/lib/netlist/solver')
-rw-r--r--src/lib/netlist/solver/nld_matrix_solver.cpp839
-rw-r--r--src/lib/netlist/solver/nld_matrix_solver.h659
-rw-r--r--src/lib/netlist/solver/nld_matrix_solver_ext.h288
-rw-r--r--src/lib/netlist/solver/nld_ms_direct.h217
-rw-r--r--src/lib/netlist/solver/nld_ms_direct1.h41
-rw-r--r--src/lib/netlist/solver/nld_ms_direct2.h58
-rw-r--r--src/lib/netlist/solver/nld_ms_direct_lu.h170
-rw-r--r--src/lib/netlist/solver/nld_ms_gcr.h453
-rw-r--r--src/lib/netlist/solver/nld_ms_gmres.h152
-rw-r--r--src/lib/netlist/solver/nld_ms_sm.h226
-rw-r--r--src/lib/netlist/solver/nld_ms_sor.h249
-rw-r--r--src/lib/netlist/solver/nld_ms_sor_mat.h165
-rw-r--r--src/lib/netlist/solver/nld_ms_w.h585
-rw-r--r--src/lib/netlist/solver/nld_solver.cpp762
-rw-r--r--src/lib/netlist/solver/nld_solver.h143
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 &gt = 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 &gt = 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 &params, 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, &params, 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_