summaryrefslogtreecommitdiffstats
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/mat_cr.h175
-rw-r--r--src/lib/netlist/solver/nld_matrix_solver.cpp937
-rw-r--r--src/lib/netlist/solver/nld_matrix_solver.h471
-rw-r--r--src/lib/netlist/solver/nld_ms_direct.h386
-rw-r--r--src/lib/netlist/solver/nld_ms_direct1.h50
-rw-r--r--src/lib/netlist/solver/nld_ms_direct2.h59
-rw-r--r--src/lib/netlist/solver/nld_ms_direct_lu.h26
-rw-r--r--src/lib/netlist/solver/nld_ms_gcr.h562
-rw-r--r--src/lib/netlist/solver/nld_ms_gmres.h412
-rw-r--r--src/lib/netlist/solver/nld_ms_sm.h416
-rw-r--r--src/lib/netlist/solver/nld_ms_sor.h112
-rw-r--r--src/lib/netlist/solver/nld_ms_sor_mat.h338
-rw-r--r--src/lib/netlist/solver/nld_ms_w.h165
-rw-r--r--src/lib/netlist/solver/nld_solver.cpp641
-rw-r--r--src/lib/netlist/solver/nld_solver.h178
-rw-r--r--src/lib/netlist/solver/vector_base.h144
16 files changed, 2841 insertions, 2231 deletions
diff --git a/src/lib/netlist/solver/mat_cr.h b/src/lib/netlist/solver/mat_cr.h
new file mode 100644
index 00000000000..8693c49c3e1
--- /dev/null
+++ b/src/lib/netlist/solver/mat_cr.h
@@ -0,0 +1,175 @@
+// license:GPL-2.0+
+// copyright-holders:Couriersud
+/*
+ * mat_cr.h
+ *
+ * Compressed row format matrices
+ *
+ */
+
+#ifndef MAT_CR_H_
+#define MAT_CR_H_
+
+#include <algorithm>
+#include "../plib/pconfig.h"
+#include "../plib/palloc.h"
+
+template<std::size_t N, typename C = uint16_t, typename T = double>
+struct mat_cr_t
+{
+ typedef C index_type;
+ typedef T value_type;
+
+ C diag[N]; // diagonal index pointer n
+ C ia[N+1]; // row index pointer n + 1
+ C ja[N*N]; // column index array nz_num, initially (n * n)
+ T A[N*N]; // Matrix elements nz_num, initially (n * n)
+
+ std::size_t size;
+ std::size_t nz_num;
+
+ explicit mat_cr_t(const std::size_t n)
+ : size(n)
+ , nz_num(0)
+ {
+#if 0
+#if 0
+ ia = plib::palloc_array<C>(n + 1);
+ ja = plib::palloc_array<C>(n * n);
+ diag = plib::palloc_array<C>(n);
+#else
+ diag = plib::palloc_array<C>(n + (n + 1) + n * n);
+ ia = diag + n;
+ ja = ia + (n+1);
+ A = plib::palloc_array<T>(n * n);
+#endif
+#endif
+ }
+
+ ~mat_cr_t()
+ {
+#if 0
+ plib::pfree_array(diag);
+#if 0
+ plib::pfree_array(ia);
+ plib::pfree_array(ja);
+#endif
+ plib::pfree_array(A);
+#endif
+ }
+
+ void set_scalar(const T scalar)
+ {
+ for (std::size_t i=0, e=nz_num; i<e; i++)
+ A[i] = scalar;
+ }
+
+ void mult_vec(const T * RESTRICT x, T * RESTRICT res)
+ {
+ /*
+ * res = A * x
+ */
+
+ std::size_t i = 0;
+ std::size_t k = 0;
+ const std::size_t oe = nz_num;
+
+ while (k < oe)
+ {
+ T tmp = 0.0;
+ const std::size_t e = ia[i+1];
+ for (; k < e; k++)
+ tmp += A[k] * x[ja[k]];
+ res[i++] = tmp;
+ }
+ }
+
+ void incomplete_LU_factorization(T * RESTRICT LU)
+ {
+ /*
+ * incomplete LU Factorization according to http://de.wikipedia.org/wiki/ILU-Zerlegung
+ *
+ * Result is stored in matrix LU
+ *
+ */
+
+ const std::size_t lnz = nz_num;
+
+ for (std::size_t k = 0; k < lnz; k++)
+ LU[k] = A[k];
+
+ for (std::size_t i = 1; ia[i] < lnz; i++) // row i
+ {
+ const std::size_t iai1 = ia[i + 1];
+ const std::size_t pke = diag[i];
+ for (std::size_t pk = ia[i]; pk < pke; pk++) // all columns left of diag in row i
+ {
+ // pk == (i, k)
+ const std::size_t k = ja[pk];
+ const std::size_t iak1 = ia[k + 1];
+ const T LUpk = LU[pk] = LU[pk] / LU[diag[k]];
+
+ std::size_t pt = ia[k];
+
+ for (std::size_t pj = pk + 1; pj < iai1; pj++) // pj = (i, j)
+ {
+ // we can assume that within a row ja increases continuously */
+ const std::size_t ej = ja[pj];
+ while (ja[pt] < ej && pt < iak1)
+ pt++;
+ if (pt < iak1 && ja[pt] == ej)
+ LU[pj] = LU[pj] - LUpk * LU[pt];
+ }
+ }
+ }
+ }
+
+ void solveLUx (const T * RESTRICT LU, T * RESTRICT r)
+ {
+ /*
+ * Solve a linear equation Ax = r
+ * where
+ * A = L*U
+ *
+ * L unit lower triangular
+ * U upper triangular
+ *
+ * ==> LUx = r
+ *
+ * ==> Ux = L⁻¹ r = w
+ *
+ * ==> r = Lw
+ *
+ * This can be solved for w using backwards elimination in L.
+ *
+ * Now Ux = w
+ *
+ * This can be solved for x using backwards elimination in U.
+ *
+ */
+
+ for (std::size_t i = 1; ia[i] < nz_num; ++i )
+ {
+ T tmp = 0.0;
+ const std::size_t j1 = ia[i];
+ const std::size_t j2 = diag[i];
+
+ for (std::size_t j = j1; j < j2; ++j )
+ tmp += LU[j] * r[ja[j]];
+
+ r[i] -= tmp;
+ }
+ // i now is equal to n;
+ for (std::size_t i = size; i-- > 0; )
+ {
+ T tmp = 0.0;
+ const std::size_t di = diag[i];
+ const std::size_t j2 = ia[i+1];
+ for (std::size_t j = di + 1; j < j2; j++ )
+ tmp += LU[j] * r[ja[j]];
+ r[i] = (r[i] - tmp) / LU[di];
+ }
+ }
+};
+
+#endif /* MAT_CR_H_ */
diff --git a/src/lib/netlist/solver/nld_matrix_solver.cpp b/src/lib/netlist/solver/nld_matrix_solver.cpp
index 1fe14b0f5e1..87c9b46ba21 100644
--- a/src/lib/netlist/solver/nld_matrix_solver.cpp
+++ b/src/lib/netlist/solver/nld_matrix_solver.cpp
@@ -6,623 +6,552 @@
*/
#include "nld_matrix_solver.h"
-#include "plib/putil.h"
+#include "../plib/putil.h"
#include <cmath> // <<= needed by windows build
namespace netlist
{
-namespace devices
+ namespace devices
+ {
+
+proxied_analog_output_t::~proxied_analog_output_t()
{
+}
- 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-9)
- {
- }
+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-9)
+{
+}
- void terms_for_net_t::add(terminal_t *term, int net_other, bool sorted)
- {
- if (sorted)
- for (std::size_t i=0; i < m_connected_net_idx.size(); i++)
+void terms_for_net_t::clear()
+{
+ m_terms.clear();
+ m_connected_net_idx.clear();
+ m_gt.clear();
+ m_go.clear();
+ m_Idr.clear();
+ m_connected_net_V.clear();
+}
+
+void terms_for_net_t::add(terminal_t *term, int net_other, bool sorted)
+{
+ if (sorted)
+ for (unsigned i=0; i < m_connected_net_idx.size(); i++)
+ {
+ if (m_connected_net_idx[i] > net_other)
{
- if (m_connected_net_idx[i] > net_other)
- {
- plib::container::insert_at(m_terms, i, term);
- plib::container::insert_at(m_connected_net_idx, i, net_other);
- return;
- }
+ plib::container::insert_at(m_terms, i, term);
+ plib::container::insert_at(m_connected_net_idx, i, net_other);
+ plib::container::insert_at(m_gt, i, 0.0);
+ plib::container::insert_at(m_go, i, 0.0);
+ plib::container::insert_at(m_Idr, i, 0.0);
+ plib::container::insert_at(m_connected_net_V, i, nullptr);
+ return;
}
- m_terms.push_back(term);
- m_connected_net_idx.push_back(net_other);
- }
-
- // ----------------------------------------------------------------------------------------
- // matrix_solver
- // ----------------------------------------------------------------------------------------
-
- matrix_solver_t::matrix_solver_t(netlist_state_t &anetlist, const pstring &name,
- const 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)
+ }
+ m_terms.push_back(term);
+ m_connected_net_idx.push_back(net_other);
+ m_gt.push_back(0.0);
+ m_go.push_back(0.0);
+ m_Idr.push_back(0.0);
+ m_connected_net_V.push_back(nullptr);
+}
+
+void terms_for_net_t::set_pointers()
+{
+ for (unsigned i = 0; i < count(); i++)
{
- connect_post_start(m_fb_sync, m_Q_sync);
+ m_terms[i]->set_ptrs(&m_gt[i], &m_go[i], &m_Idr[i]);
+ m_connected_net_V[i] = m_terms[i]->m_otherterm->net().Q_Analog_state_ptr();
}
+}
+
+// ----------------------------------------------------------------------------------------
+// matrix_solver
+// ----------------------------------------------------------------------------------------
+
+matrix_solver_t::matrix_solver_t(netlist_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)
+{
+ connect_post_start(m_fb_sync, m_Q_sync);
+}
- void matrix_solver_t::setup_base(analog_net_t::list_t &nets)
- {
+matrix_solver_t::~matrix_solver_t()
+{
+}
- log().debug("New solver setup\n");
+void matrix_solver_t::setup_base(analog_net_t::list_t &nets)
+{
- m_nets.clear();
- m_terms.clear();
+ log().debug("New solver setup\n");
- for (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_nets.clear();
+ m_terms.clear();
- for (std::size_t k = 0; k < nets.size(); k++)
- {
- analog_net_t *net = nets[k];
+ for (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::palloc<terms_for_net_t>());
+ }
- log().debug("setting up net\n");
+ for (std::size_t k = 0; k < nets.size(); k++)
+ {
+ analog_net_t *net = nets[k];
- net->set_solver(this);
+ log().debug("setting up net\n");
- for (auto &p : net->core_terms())
- {
- log().debug("{1} {2} {3}\n", p->name(), net->name(), net->isRailNet());
- 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_dynamic())
- if (!plib::container::contains(m_dynamic_devices, &p->device()))
- m_dynamic_devices.push_back(&p->device());
- {
- auto *pterm = dynamic_cast<terminal_t *>(p);
- add_term(k, pterm);
- }
- log().debug("Added terminal {1}\n", p->name());
- break;
- case detail::terminal_type::INPUT:
- {
- proxied_analog_output_t *net_proxy_output = nullptr;
- for (auto & input : m_inps)
- if (input->proxied_net() == &p->net())
- {
- net_proxy_output = input.get();
- break;
- }
-
- if (net_proxy_output == nullptr)
- {
- pstring nname = this->name() + "." + pstring(plib::pfmt("m{1}")(m_inps.size()));
- nl_assert(p->net().is_analog());
- auto net_proxy_output_u = 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->net().add_terminal(*p);
- // FIXME: repeated calling - kind of brute force
- net_proxy_output->net().rebuild_list();
- log().debug("Added input\n");
- }
- break;
- case detail::terminal_type::OUTPUT:
- log().fatal(MF_1_UNHANDLED_ELEMENT_1_FOUND,
- p->name());
- break;
- }
- }
- log().debug("added net with {1} populated connections\n", net->core_terms().size());
- }
+ net->set_solver(this);
- /* now setup the matrix */
- setup_matrix();
- }
-
- void matrix_solver_t::sort_terms(eSortType 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();
-
- switch (sort)
+ for (auto &p : net->m_core_terms)
{
- case PREFER_BAND_MATRIX:
- {
- for (std::size_t k = 0; k < iN - 1; k++)
+ log().debug("{1} {2} {3}\n", p->name(), net->name(), net->isRailNet());
+ 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_dynamic())
+ if (!plib::container::contains(m_dynamic_devices, &p->device()))
+ m_dynamic_devices.push_back(&p->device());
{
- auto pk = get_weight_around_diag(k,k);
- for (std::size_t i = k+1; i < iN; i++)
- {
- auto pi = get_weight_around_diag(i,k);
- if (pi < pk)
- {
- std::swap(m_terms[i], m_terms[k]);
- std::swap(m_nets[i], m_nets[k]);
- pk = get_weight_around_diag(k,k);
- }
- }
+ terminal_t *pterm = dynamic_cast<terminal_t *>(p);
+ add_term(k, pterm);
}
- }
- break;
- case PREFER_IDENTITY_TOP_LEFT:
- {
- for (std::size_t k = 0; k < iN - 1; k++)
+ log().debug("Added terminal {1}\n", p->name());
+ break;
+ case detail::terminal_type::INPUT:
{
- auto pk = get_left_right_of_diag(k,k);
- for (std::size_t i = k+1; i < iN; i++)
- {
- auto pi = get_left_right_of_diag(i,k);
- if (pi.first <= pk.first && pi.second >= pk.second)
+ proxied_analog_output_t *net_proxy_output = nullptr;
+ for (auto & input : m_inps)
+ if (input->m_proxied_net == &p->net())
{
- std::swap(m_terms[i], m_terms[k]);
- std::swap(m_nets[i], m_nets[k]);
- pk = get_left_right_of_diag(k,k);
+ net_proxy_output = input.get();
+ break;
}
- }
- }
- }
- break;
- case ASCENDING:
- case DESCENDING:
- {
- int sort_order = (m_sort == DESCENDING ? 1 : -1);
- for (std::size_t k = 0; k < iN - 1; k++)
- for (std::size_t i = k+1; i < iN; i++)
+ if (net_proxy_output == nullptr)
{
- 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]);
- }
+ pstring nname = this->name() + "." + pstring(plib::pfmt("m{1}")(m_inps.size()));
+ auto net_proxy_output_u = plib::make_unique<proxied_analog_output_t>(*this, nname);
+ net_proxy_output = net_proxy_output_u.get();
+ m_inps.push_back(std::move(net_proxy_output_u));
+ nl_assert(p->net().is_analog());
+ net_proxy_output->m_proxied_net = static_cast<analog_net_t *>(&p->net());
}
- }
- break;
- case NOSORT:
- break;
- }
- /* rebuild */
- for (auto &term : m_terms)
- {
- int *other = term->m_connected_net_idx.data();
- for (std::size_t i = 0; i < term->count(); i++)
- //FIXME: this is weird
- if (other[i] != -1)
- other[i] = get_net_idx(&term->terms()[i]->connected_terminal()->net());
+ 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");
+ }
+ break;
+ case detail::terminal_type::OUTPUT:
+ log().fatal(MF_1_UNHANDLED_ELEMENT_1_FOUND,
+ p->name());
+ break;
+ }
}
+ log().debug("added net with {1} populated connections\n", net->m_core_terms.size());
}
- void matrix_solver_t::setup_matrix()
- {
- const std::size_t iN = m_nets.size();
+ /* now setup the matrix */
+ setup_matrix();
+}
- 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);
- }
-
- // free all - no longer needed
- m_rails_temp.clear();
-
- sort_terms(m_sort);
-
- this->set_pointers();
+void matrix_solver_t::setup_matrix()
+{
+ const std::size_t iN = m_nets.size();
- /* 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();
+ 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]->connected_net_idx()[i], false);
- t->m_nz.clear();
+ m_terms[k]->set_pointers();
+ }
- 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 (terms_for_net_t *rt : m_rails_temp)
+ {
+ rt->clear(); // no longer needed
+ plib::pfree(rt); // no longer needed
+ }
- t->m_nz.push_back(k); // add diagonal
+ m_rails_temp.clear();
- /* and sort */
- std::sort(t->m_nz.begin(), t->m_nz.end());
- }
+ /* 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.
+ *
+ */
- /* 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();
+ if (m_sort != NOSORT)
+ {
+ int sort_order = (m_sort == DESCENDING ? 1 : -1);
- if (k==0)
- t->m_nzrd.clear();
- else
+ for (unsigned k = 0; k < iN - 1; k++)
+ for (unsigned i = k+1; i < iN; i++)
{
- t->m_nzrd = m_terms[k-1]->m_nzrd;
- for (auto j = t->m_nzrd.begin(); j != t->m_nzrd.end(); )
+ if ((static_cast<int>(m_terms[k]->m_railstart) - static_cast<int>(m_terms[i]->m_railstart)) * sort_order < 0)
{
- if (*j < k + 1)
- j = t->m_nzrd.erase(j);
- else
- ++j;
+ std::swap(m_terms[i], m_terms[k]);
+ std::swap(m_nets[i], m_nets[k]);
}
}
- 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]));
-
- /* and sort */
- std::sort(t->m_nzrd.begin(), t->m_nzrd.end());
- }
-
- /* create a list of non zero elements below diagonal k
- * This should reduce cache misses ...
- */
-
- std::vector<std::vector<bool>> touched(iN, std::vector<bool>(iN));
-
- for (std::size_t k = 0; k < iN; k++)
- {
- for (std::size_t j = 0; j < iN; j++)
- touched[k][j] = false;
- for (std::size_t j = 0; j < m_terms[k]->m_nz.size(); j++)
- touched[k][m_terms[k]->m_nz[j]] = true;
- }
-
- m_ops = 0;
- for (unsigned k = 0; k < iN; k++)
- {
- m_ops++; // 1/A(k,k)
- for (unsigned row = k + 1; row < iN; row++)
- {
- if (touched[row][k])
- {
- m_ops++;
- if (!plib::container::contains(m_terms[k]->m_nzbd, row))
- m_terms[k]->m_nzbd.push_back(row);
- for (std::size_t col = k + 1; col < iN; col++)
- if (touched[k][col])
- {
- touched[row][col] = true;
- m_ops += 2;
- }
- }
- }
- }
- log().verbose("Number of mults/adds for {1}: {2}", name(), m_ops);
-
- if ((false))
- for (std::size_t k = 0; k < iN; k++)
- {
- pstring line = plib::pfmt("{1:3}")(k);
- for (const auto & nzrd : m_terms[k]->m_nzrd)
- line += plib::pfmt(" {1:3}")(nzrd);
- log().verbose("{1}", line);
- }
-
- /*
- * save states
- */
- for (std::size_t k = 0; k < iN; k++)
+ for (auto &term : m_terms)
{
- 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);
-
- // 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());
+ int *other = term->connected_net_idx();
+ for (unsigned i = 0; i < term->count(); i++)
+ if (other[i] != -1)
+ other[i] = get_net_idx(&term->terms()[i]->m_otherterm->net());
}
}
- void matrix_solver_t::update_inputs()
+ /* create a list of non zero elements. */
+ for (unsigned k = 0; k < iN; k++)
{
- // avoid recursive calls. Inputs are updated outside this call
- for (auto &inp : m_inps)
- inp->push(inp->proxied_net()->Q_Analog());
- }
+ terms_for_net_t * t = m_terms[k].get();
+ /* pretty brutal */
+ int *other = t->connected_net_idx();
- void matrix_solver_t::update_dynamic()
- {
- /* update all non-linear devices */
- for (auto &dyn : m_dynamic_devices)
- dyn->update_terminals();
- }
+ t->m_nz.clear();
- void matrix_solver_t::reset()
- {
- m_last_step = netlist_time::zero();
- }
+ for (unsigned 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]));
- void matrix_solver_t::update() NL_NOEXCEPT
- {
- const netlist_time new_timestep = solve(exec().time());
- update_inputs();
+ t->m_nz.push_back(k); // add diagonal
- if (m_params.m_dynamic_ts && has_timestep_devices() && new_timestep > netlist_time::zero())
- {
- m_Q_sync.net().toggle_and_push_to_queue(new_timestep);
- }
+ /* and sort */
+ std::sort(t->m_nz.begin(), t->m_nz.end());
}
- /* update_forced is called from within param_update
- *
- * this should only occur outside of execution and thus
- * using time should be safe.
- *
+ /* create a list of non zero elements right of the diagonal
+ * These list anticipate the population of array elements by
+ * Gaussian elimination.
*/
- void matrix_solver_t::update_forced()
+ for (unsigned k = 0; k < iN; k++)
{
- const netlist_time new_timestep = solve(exec().time());
- plib::unused_var(new_timestep);
-
- update_inputs();
+ terms_for_net_t * t = m_terms[k].get();
+ /* pretty brutal */
+ int *other = t->connected_net_idx();
- if (m_params.m_dynamic_ts && has_timestep_devices())
+ if (k==0)
+ t->m_nzrd.clear();
+ else
{
- m_Q_sync.net().toggle_and_push_to_queue(netlist_time::from_double(m_params.m_min_timestep));
+ t->m_nzrd = m_terms[k-1]->m_nzrd;
+ for (auto j = t->m_nzrd.begin(); j != t->m_nzrd.end(); )
+ {
+ if (*j < k + 1)
+ j = t->m_nzrd.erase(j);
+ else
+ ++j;
+ }
}
+
+ for (unsigned 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]));
+
+ /* and sort */
+ std::sort(t->m_nzrd.begin(), t->m_nzrd.end());
}
- void matrix_solver_t::step(const netlist_time &delta)
+ /* create a list of non zero elements below diagonal k
+ * This should reduce cache misses ...
+ */
+
+ bool **touched = plib::palloc_array<bool *>(iN);
+ for (unsigned k=0; k<iN; k++)
+ touched[k] = plib::palloc_array<bool>(iN);
+
+ for (unsigned k = 0; k < iN; k++)
{
- const nl_double dd = delta.as_double();
- for (auto &d : m_step_devices)
- d->timestep(dd);
+ for (unsigned j = 0; j < iN; j++)
+ touched[k][j] = false;
+ for (unsigned j = 0; j < m_terms[k]->m_nz.size(); j++)
+ touched[k][m_terms[k]->m_nz[j]] = true;
}
- void matrix_solver_t::solve_base()
+ m_ops = 0;
+ for (unsigned k = 0; k < iN; k++)
{
- ++m_stat_vsolver_calls;
- if (has_dynamic_devices())
+ m_ops++; // 1/A(k,k)
+ for (unsigned row = k + 1; row < iN; row++)
{
- 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())
+ if (touched[row][k])
{
- log().warning(MW_1_NEWTON_LOOPS_EXCEEDED_ON_NET_1, this->name());
- m_Q_sync.net().toggle_and_push_to_queue(m_params.m_nr_recalc_delay);
+ m_ops++;
+ if (!plib::container::contains(m_terms[k]->m_nzbd, row))
+ m_terms[k]->m_nzbd.push_back(row);
+ for (unsigned col = k + 1; col < iN; col++)
+ if (touched[k][col])
+ {
+ touched[row][col] = true;
+ m_ops += 2;
+ }
}
}
- else
+ }
+ log().verbose("Number of mults/adds for {1}: {2}", name(), m_ops);
+
+ if ((0))
+ for (unsigned k = 0; k < iN; k++)
{
- this->vsolve_non_dynamic(false);
+ pstring line = plib::pfmt("{1:3}")(k);
+ for (unsigned j = 0; j < m_terms[k]->m_nzrd.size(); j++)
+ line += plib::pfmt(" {1:3}")(m_terms[k]->m_nzrd[j]);
+ log().verbose("{1}", line);
}
- }
- const netlist_time matrix_solver_t::solve(netlist_time now)
+ /*
+ * save states
+ */
+ for (unsigned k = 0; k < iN; k++)
{
- const netlist_time delta = now - m_last_step;
-
- // We are already up to date. Avoid oscillations.
- // FIXME: Make this a parameter!
- if (delta < netlist_time::quantum())
- return netlist_time::zero();
+ pstring num = plib::pfmt("{1}")(k);
- /* 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());
+ netlist().save(*this, m_terms[k]->m_last_V, "lastV." + num);
+ netlist().save(*this, m_terms[k]->m_DD_n_m_1, "m_DD_n_m_1." + num);
+ netlist().save(*this, m_terms[k]->m_h_n_m_1, "m_h_n_m_1." + num);
- return next_time_step;
+ netlist().save(*this, m_terms[k]->go(),"GO" + num, m_terms[k]->count());
+ netlist().save(*this, m_terms[k]->gt(),"GT" + num, m_terms[k]->count());
+ netlist().save(*this, m_terms[k]->Idr(),"IDR" + num , m_terms[k]->count());
}
- int matrix_solver_t::get_net_idx(detail::net_t *net)
- {
- for (std::size_t k = 0; k < m_nets.size(); k++)
- if (m_nets[k] == net)
- return static_cast<int>(k);
- return -1;
- }
+ for (unsigned k=0; k<iN; k++)
+ plib::pfree_array(touched[k]);
+ plib::pfree_array(touched);
+}
- std::pair<int, int> matrix_solver_t::get_left_right_of_diag(std::size_t irow, std::size_t idiag)
- {
- /*
- * return the maximum column left of the diagonal (-1 if no cols found)
- * return the minimum column right of the diagonal (999999 if no cols found)
- */
+void matrix_solver_t::update_inputs()
+{
+ // avoid recursive calls. Inputs are updated outside this call
+ for (auto &inp : m_inps)
+ inp->push(inp->m_proxied_net->Q_Analog());
+}
- const auto row = static_cast<int>(irow);
- const auto diag = static_cast<int>(idiag);
+void matrix_solver_t::update_dynamic()
+{
+ /* update all non-linear devices */
+ for (auto &dyn : m_dynamic_devices)
+ dyn->update_terminals();
+}
- int colmax = -1;
- int colmin = 999999;
+void matrix_solver_t::reset()
+{
+ m_last_step = netlist_time::zero();
+}
- auto &term = m_terms[irow];
+void matrix_solver_t::update() NL_NOEXCEPT
+{
+ const netlist_time new_timestep = solve();
+ update_inputs();
- for (std::size_t i = 0; i < term->count(); i++)
- {
- auto col = get_net_idx(&term->terms()[i]->connected_terminal()->net());
- if (col != -1)
- {
- if (col==row) col = diag;
- else if (col==diag) col = row;
+ if (m_params.m_dynamic_ts && has_timestep_devices() && new_timestep > netlist_time::zero())
+ {
+ m_Q_sync.net().toggle_and_push_to_queue(new_timestep);
+ }
+}
- if (col > diag && col < colmin)
- colmin = col;
- else if (col < diag && col > colmax)
- colmax = col;
- }
- }
- return {colmax, colmin};
+void matrix_solver_t::update_forced()
+{
+ ATTR_UNUSED const netlist_time new_timestep = solve();
+ update_inputs();
+
+ if (m_params.m_dynamic_ts && has_timestep_devices())
+ {
+ m_Q_sync.net().toggle_and_push_to_queue(netlist_time::from_double(m_params.m_min_timestep));
}
+}
- double matrix_solver_t::get_weight_around_diag(std::size_t row, std::size_t diag)
+void matrix_solver_t::step(const netlist_time &delta)
+{
+ const nl_double dd = delta.as_double();
+ for (std::size_t k=0; k < m_step_devices.size(); k++)
+ m_step_devices[k]->timestep(dd);
+}
+
+void matrix_solver_t::solve_base()
+{
+ ++m_stat_vsolver_calls;
+ if (has_dynamic_devices())
{
+ unsigned this_resched;
+ unsigned newton_loops = 0;
+ do
{
- /*
- * return average absolute distance
- */
+ 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_1_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);
+ }
+}
- std::vector<bool> touched(1024, false); // FIXME!
+const netlist_time matrix_solver_t::solve()
+{
+ const netlist_time now = netlist().time();
+ const netlist_time delta = now - m_last_step;
- double weight = 0.0;
- 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());
- 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);
+ // We are already up to date. Avoid oscillations.
+ // FIXME: Make this a parameter!
+ if (delta < netlist_time::quantum())
+ return netlist_time::zero();
- weight = weight + std::abs(static_cast<double>(colu) - static_cast<double>(diag));
- touched[colu] = true;
- }
- }
- }
- return weight; // / static_cast<double>(term->m_railstart);
- }
- }
+ /* 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;
+}
- void matrix_solver_t::add_term(std::size_t k, terminal_t *term)
+int matrix_solver_t::get_net_idx(detail::net_t *net)
+{
+ for (std::size_t k = 0; k < m_nets.size(); k++)
+ if (m_nets[k] == net)
+ return static_cast<int>(k);
+ return -1;
+}
+
+void matrix_solver_t::add_term(std::size_t k, terminal_t *term)
+{
+ if (term->m_otherterm->net().isRailNet())
{
- if (term->connected_terminal()->net().isRailNet())
+ m_rails_temp[k]->add(term, -1, false);
+ }
+ else
+ {
+ int ot = get_net_idx(&term->m_otherterm->net());
+ if (ot>=0)
{
- m_rails_temp[k]->add(term, -1, false);
+ m_terms[k]->add(term, ot, true);
}
- else
+ /* Should this be allowed ? */
+ else // if (ot<0)
{
- int ot = get_net_idx(&term->connected_terminal()->net());
- if (ot>=0)
- {
- m_terms[k]->add(term, ot, true);
- }
- /* Should this be allowed ? */
- else // if (ot<0)
- {
- m_rails_temp[k]->add(term, ot, true);
- log().fatal(MF_1_FOUND_TERM_WITH_MISSING_OTHERNET, term->name());
- }
+ m_rails_temp[k]->add(term, ot, true);
+ log().fatal(MF_1_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;
+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)
+ if (m_params.m_dynamic_ts)
+ {
+ for (std::size_t k = 0, iN=m_terms.size(); k < iN; k++)
{
- 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();
+ 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);
- const nl_double hn = cur_ts;
+ const nl_double DD_n = (n->Q_Analog() - t->m_last_V);
+ const nl_double hn = cur_ts;
- 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;
+ //printf("%f %f %f %f\n", DD_n, t->m_DD_n_m_1, hn, 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;
+ t->m_h_n_m_1 = hn;
+ t->m_DD_n_m_1 = DD_n;
+ if (std::fabs(DD2) > NL_FCONST(1e-60)) // avoid div-by-zero
+ new_net_timestep = std::sqrt(m_params.m_dynamic_lte / std::fabs(NL_FCONST(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;
+ 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)
- {
- //log().warning("Dynamic timestep below min timestep. Consider decreasing MIN_TIMESTEP: {1} us", new_solver_timestep*1.0e6);
- new_solver_timestep = m_params.m_min_timestep;
- }
+ t->m_last_V = n->Q_Analog();
+ }
+ if (new_solver_timestep < m_params.m_min_timestep)
+ {
+ //log().warning("Dynamic timestep below min timestep. Consider decreasing MIN_TIMESTEP: {1} us", new_solver_timestep*1.0e6);
+ 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);
}
+ //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()
+
+
+void matrix_solver_t::log_stats()
+{
+ if (this->m_stat_calculations != 0 && this->m_stat_vsolver_calls && this->m_params.m_log_stats)
{
- if (this->m_stat_calculations != 0 && this->m_stat_vsolver_calls && this->m_params.m_log_stats)
- {
- 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("==============================================");
+ 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->netlist().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));
}
+}
+
-} // namespace devices
+ } //namespace devices
} // namespace netlist
diff --git a/src/lib/netlist/solver/nld_matrix_solver.h b/src/lib/netlist/solver/nld_matrix_solver.h
index f76660e1cb9..da44370d0d9 100644
--- a/src/lib/netlist/solver/nld_matrix_solver.h
+++ b/src/lib/netlist/solver/nld_matrix_solver.h
@@ -10,367 +10,282 @@
#include "netlist/nl_base.h"
#include "netlist/nl_errstr.h"
-#include "plib/palloc.h"
-#include "plib/pmatrix2d.h"
-#include "plib/putil.h"
-#include "plib/vector_ops.h"
-
-#include <cmath>
+#include "netlist/plib/putil.h"
namespace netlist
{
-namespace devices
-{
+ namespace devices
+ {
/* FIXME: these should become proper devices */
struct solver_parameters_t
{
- bool m_pivot;
+ int 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;
+ unsigned m_gs_loops;
+ unsigned m_nr_loops;
netlist_time m_nr_recalc_delay;
- bool m_use_gabs;
- bool m_use_linear_prediction;
bool m_log_stats;
};
- class terms_for_net_t : plib::nocopyassignmove
- {
- public:
- terms_for_net_t();
-
- void clear();
+class terms_for_net_t : plib::nocopyassignmove
+{
+public:
+ terms_for_net_t();
- void add(terminal_t *term, int net_other, bool sorted);
+ void clear();
- std::size_t count() const { return m_terms.size(); }
+ void add(terminal_t *term, int net_other, bool sorted);
- terminal_t **terms() { return m_terms.data(); }
+ inline std::size_t count() const { return m_terms.size(); }
- std::size_t m_railstart;
+ inline terminal_t **terms() { return m_terms.data(); }
+ inline int *connected_net_idx() { return m_connected_net_idx.data(); }
+ inline nl_double *gt() { return m_gt.data(); }
+ inline nl_double *go() { return m_go.data(); }
+ inline nl_double *Idr() { return m_Idr.data(); }
+ inline nl_double * const *connected_net_V() const { return m_connected_net_V.data(); }
- 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 */
+ void set_pointers();
- /* state */
- nl_double m_last_V;
- nl_double m_DD_n_m_1;
- nl_double m_h_n_m_1;
+ std::size_t m_railstart;
- std::vector<int> m_connected_net_idx;
- private:
- std::vector<terminal_t *> m_terms;
+ 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 */
- };
+ /* state */
+ nl_double m_last_V;
+ nl_double m_DD_n_m_1;
+ nl_double m_h_n_m_1;
- class proxied_analog_output_t : public analog_output_t
- {
- public:
+private:
+ std::vector<int> m_connected_net_idx;
+ std::vector<nl_double> m_go;
+ std::vector<nl_double> m_gt;
+ std::vector<nl_double> m_Idr;
+ std::vector<nl_double *> m_connected_net_V;
+ std::vector<terminal_t *> m_terms;
- proxied_analog_output_t(core_device_t &dev, const pstring &aname, analog_net_t *pnet)
- : analog_output_t(dev, aname)
- , m_proxied_net(pnet)
- { }
+};
- 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
- };
+class proxied_analog_output_t : public analog_output_t
+{
+public:
- class matrix_solver_t : public device_t
- {
- public:
- using list_t = std::vector<matrix_solver_t *>;
+ proxied_analog_output_t(core_device_t &dev, const pstring &aname)
+ : analog_output_t(dev, aname)
+ , m_proxied_net(nullptr)
+ { }
+ virtual ~proxied_analog_output_t();
- enum eSortType
- {
- NOSORT,
- ASCENDING,
- DESCENDING,
- PREFER_IDENTITY_TOP_LEFT,
- PREFER_BAND_MATRIX
- };
-
- void setup(analog_net_t::list_t &nets)
- {
- vsetup(nets);
- }
+ analog_net_t *m_proxied_net; // only for proxy nets in analog input logic
+};
- void solve_base();
- /* 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();
+class matrix_solver_t : public device_t
+{
+public:
+ using list_t = std::vector<matrix_solver_t *>;
- bool has_dynamic_devices() const { return m_dynamic_devices.size() > 0; }
- bool has_timestep_devices() const { return m_step_devices.size() > 0; }
+ enum eSortType
+ {
+ NOSORT,
+ ASCENDING,
+ DESCENDING
+ };
- void update_forced();
- void update_after(const netlist_time after)
- {
- m_Q_sync.net().toggle_and_push_to_queue(after);
- }
+ virtual ~matrix_solver_t() override;
- /* netdevice functions */
- NETLIB_UPDATEI();
- NETLIB_RESETI();
+ void setup(analog_net_t::list_t &nets)
+ {
+ vsetup(nets);
+ }
- 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);
+ void solve_base();
- virtual void log_stats();
+ /* 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();
+ void update_inputs();
- virtual std::pair<pstring, pstring> create_solver_code()
- {
- return std::pair<pstring, pstring>("", plib::pfmt("/* solver doesn't support static compile */\n\n"));
- }
+ inline bool has_dynamic_devices() const { return m_dynamic_devices.size() > 0; }
+ inline bool has_timestep_devices() const { return m_step_devices.size() > 0; }
- /* return number of floating point operations for solve */
- std::size_t ops() { return m_ops; }
+ void update_forced();
+ void update_after(const netlist_time &after)
+ {
+ m_Q_sync.net().toggle_and_push_to_queue(after);
+ }
- protected:
+ /* netdevice functions */
+ NETLIB_UPDATEI();
+ NETLIB_RESETI();
- matrix_solver_t(netlist_state_t &anetlist, const pstring &name,
- eSortType sort, const solver_parameters_t *params);
+public:
+ int get_net_idx(detail::net_t *net);
- void sort_terms(eSortType sort);
+ virtual void log_stats();
- void setup_base(analog_net_t::list_t &nets);
- void update_dynamic();
+ virtual std::pair<pstring, pstring> create_solver_code()
+ {
+ return std::pair<pstring, pstring>("", plib::pfmt("/* solver doesn't support static compile */\n\n"));
+ }
- virtual void vsetup(analog_net_t::list_t &nets) = 0;
- virtual unsigned vsolve_non_dynamic(const bool newton_raphson) = 0;
+ /* return number of floating point operations for solve */
+ std::size_t ops() { return m_ops; }
- netlist_time compute_next_timestep(const double cur_ts);
- /* virtual */ void add_term(std::size_t net_idx, terminal_t *term);
+protected:
- template <typename T>
- void store(const T & V);
+ matrix_solver_t(netlist_t &anetlist, const pstring &name,
+ const eSortType sort, const solver_parameters_t *params);
- template <typename T>
- auto delta(const T & V) -> typename std::decay<decltype(V[0])>::type;
+ void setup_base(analog_net_t::list_t &nets);
+ void update_dynamic();
- template <typename T>
- void build_LE_A(T &child);
- template <typename T>
- void build_LE_RHS(T &child);
+ virtual void vsetup(analog_net_t::list_t &nets) = 0;
+ virtual unsigned vsolve_non_dynamic(const bool newton_raphson) = 0;
- void set_pointers()
- {
- const std::size_t iN = this->m_nets.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);
- }
-
- 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);
-
- 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] = m_terms[k]->terms()[i]->connected_terminal()->net().Q_Analog_state_ptr();
- }
- }
- }
+ netlist_time compute_next_timestep(const double cur_ts);
+ /* virtual */ void add_term(std::size_t net_idx, terminal_t *term);
- template <typename AP, typename FT>
- void fill_matrix(std::size_t N, AP &tcr, FT &RHS)
- {
- for (std::size_t k = 0; k < N; k++)
- {
- 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;
- }
+ template <typename T>
+ void store(const T * RESTRICT V);
+ template <typename T>
+ T delta(const T * RESTRICT V);
- }
+ template <typename T>
+ void build_LE_A();
+ template <typename T>
+ void build_LE_RHS();
- template <typename T>
- using aligned_alloc = plib::aligned_allocator<T, PALIGN_VECTOROPT>;
+ std::vector<std::unique_ptr<terms_for_net_t>> m_terms;
+ std::vector<analog_net_t *> m_nets;
+ std::vector<std::unique_ptr<proxied_analog_output_t>> m_inps;
- 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;
+ std::vector<terms_for_net_t *> m_rails_temp;
- plib::pmatrix2d<nl_double> m_test;
+ const solver_parameters_t &m_params;
- 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;
+ 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;
- std::vector<plib::unique_ptr<terms_for_net_t>> m_rails_temp;
+private:
- const solver_parameters_t &m_params;
+ 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<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;
+ logic_input_t m_fb_sync;
+ logic_output_t m_Q_sync;
- private:
+ /* calculate matrix */
+ void setup_matrix();
- state_var<netlist_time> m_last_step;
- std::vector<core_device_t *> m_step_devices;
- std::vector<core_device_t *> m_dynamic_devices;
+ void step(const netlist_time &delta);
- logic_input_t m_fb_sync;
- logic_output_t m_Q_sync;
+ std::size_t m_ops;
+ const eSortType m_sort;
+};
- /* calculate matrix */
- void setup_matrix();
+template <typename T>
+T matrix_solver_t::delta(const T * RESTRICT V)
+{
+ /* 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();
+ T cerr = 0;
+ for (std::size_t i = 0; i < iN; i++)
+ cerr = std::max(cerr, std::abs(V[i] - static_cast<T>(this->m_nets[i]->Q_Analog())));
+ return cerr;
+}
+
+template <typename T>
+void matrix_solver_t::store(const T * RESTRICT 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]);
+}
- void step(const netlist_time &delta);
+template <typename T>
+void matrix_solver_t::build_LE_A()
+{
+ static_assert(std::is_base_of<matrix_solver_t, T>::value, "T must derive from matrix_solver_t");
- std::size_t m_ops;
- const eSortType m_sort;
- };
+ T &child = static_cast<T &>(*this);
- template <typename T>
- auto matrix_solver_t::delta(const T & V) -> typename std::decay<decltype(V[0])>::type
+ const std::size_t iN = child.N();
+ for (std::size_t k = 0; k < iN; k++)
{
- /* 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;
- }
+ terms_for_net_t *terms = m_terms[k].get();
+ nl_double * Ak = &child.A(k, 0);
- 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]);
- }
+ for (std::size_t i=0; i < iN; i++)
+ Ak[i] = 0.0;
- 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");
+ const std::size_t terms_count = terms->count();
+ const std::size_t railstart = terms->m_railstart;
+ const nl_double * const RESTRICT gt = terms->gt();
- 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);
-
- for (std::size_t i=0; i < iN; i++)
- Ak[i] = 0.0;
-
- const std::size_t terms_count = terms->count();
- const std::size_t railstart = terms->m_railstart;
- const float_type * const gt = m_gtn[k];
-
- {
- float_type akk = 0.0;
- for (std::size_t i = 0; i < terms_count; i++)
- akk += gt[i];
+ nl_double akk = 0.0;
+ for (std::size_t i = 0; i < terms_count; i++)
+ akk += gt[i];
- Ak[k] = akk;
- }
+ Ak[k] = akk;
+ }
- const float_type * const go = m_gonn[k];
- int * net_other = terms->m_connected_net_idx.data();
+ const nl_double * const RESTRICT go = terms->go();
+ int * RESTRICT net_other = terms->connected_net_idx();
- for (std::size_t i = 0; i < railstart; i++)
- Ak[net_other[i]] += go[i];
- }
+ for (std::size_t i = 0; i < railstart; i++)
+ Ak[net_other[i]] -= go[i];
}
+}
- 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;
+template <typename T>
+void matrix_solver_t::build_LE_RHS()
+{
+ static_assert(std::is_base_of<matrix_solver_t, T>::value, "T must derive from matrix_solver_t");
+ T &child = static_cast<T &>(*this);
- 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;
+ const std::size_t iN = child.N();
+ for (std::size_t k = 0; k < iN; k++)
+ {
+ nl_double rhsk_a = 0.0;
+ nl_double rhsk_b = 0.0;
- 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];
+ const std::size_t terms_count = m_terms[k]->count();
+ const nl_double * const RESTRICT go = m_terms[k]->go();
+ const nl_double * const RESTRICT Idr = m_terms[k]->Idr();
+ const nl_double * const * RESTRICT other_cur_analog = m_terms[k]->connected_net_V();
- for (std::size_t i = 0; i < terms_count; i++)
- rhsk_a = rhsk_a + Idr[i];
+ for (std::size_t i = 0; i < terms_count; i++)
+ rhsk_a = rhsk_a + Idr[i];
- 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];
+ 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];
- child.RHS(k) = rhsk_a + rhsk_b;
- }
+ child.RHS(k) = rhsk_a + rhsk_b;
}
+}
-} //namespace devices
+ } //namespace devices
} // namespace netlist
#endif /* NLD_MS_DIRECT_H_ */
diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h
index 2501742218d..f03cb738028 100644
--- a/src/lib/netlist/solver/nld_ms_direct.h
+++ b/src/lib/netlist/solver/nld_ms_direct.h
@@ -8,228 +8,286 @@
#ifndef NLD_MS_DIRECT_H_
#define NLD_MS_DIRECT_H_
-#include "nld_matrix_solver.h"
+#include <algorithm>
+
#include "nld_solver.h"
-#include "plib/mat_cr.h"
-#include "plib/vector_ops.h"
+#include "nld_matrix_solver.h"
+#include "vector_base.h"
-#include <algorithm>
-#include <cmath>
+/* Disabling dynamic allocation gives a ~10% boost in performance
+ * This flag has been added to support continuous storage for arrays
+ * going forward in case we implement cuda solvers in the future.
+ */
+#define NL_USE_DYNAMIC_ALLOCATION (1)
namespace netlist
{
-namespace devices
-{
-
- template <typename FT, int SIZE>
- class matrix_solver_direct_t: public matrix_solver_t
+ namespace devices
{
- friend class matrix_solver_t;
- public:
-
- using float_type = FT;
+//#define nl_ext_double _float128 // slow, very slow
+//#define nl_ext_double long double // slightly slower
+#define nl_ext_double nl_double
- 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);
- 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);
+template <std::size_t m_N, std::size_t storage_N>
+class matrix_solver_direct_t: public matrix_solver_t
+{
+ friend class matrix_solver_t;
+public:
- constexpr std::size_t size() const { return (SIZE > 0) ? static_cast<std::size_t>(SIZE) : m_dim; }
+ matrix_solver_direct_t(netlist_t &anetlist, const pstring &name, const solver_parameters_t *params, const std::size_t size);
+ matrix_solver_direct_t(netlist_t &anetlist, const pstring &name, const eSortType sort, const solver_parameters_t *params, const std::size_t size);
- void LE_solve();
+ virtual ~matrix_solver_direct_t() override;
- template <typename T>
- void LE_back_subst(T & x);
+ virtual void vsetup(analog_net_t::list_t &nets) override;
+ virtual void reset() override { matrix_solver_t::reset(); }
- 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;
+protected:
+ virtual unsigned vsolve_non_dynamic(const bool newton_raphson) override;
+ unsigned solve_non_dynamic(const bool newton_raphson);
- 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;
+ constexpr std::size_t N() const { return (m_N == 0) ? m_dim : m_N; }
- const std::size_t m_dim;
- const std::size_t m_pitch;
- plib::parray<FT, SIZE * int(m_pitch_ABS)> m_A;
+ void LE_solve();
- };
+ template <typename T>
+ void LE_back_subst(T * RESTRICT x);
+
+#if (NL_USE_DYNAMIC_ALLOCATION)
+ template <typename T1, typename T2>
+ nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r * m_pitch + c]; }
+ template <typename T1>
+ nl_ext_double &RHS(const T1 &r) { return m_A[r * m_pitch + N()]; }
+#else
+ template <typename T1, typename T2>
+ nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r][c]; }
+ template <typename T1>
+ nl_ext_double &RHS(const T1 &r) { return m_A[r][N()]; }
+#endif
+ nl_double m_last_RHS[storage_N]; // right hand side - contains currents
+
+private:
+ //static const std::size_t m_pitch = (((storage_N + 1) + 0) / 1) * 1;
+ static constexpr std::size_t m_pitch = (((storage_N + 1) + 7) / 8) * 8;
+ //static const std::size_t m_pitch = (((storage_N + 1) + 15) / 16) * 16;
+ //static const std::size_t m_pitch = (((storage_N + 1) + 31) / 32) * 32;
+#if (NL_USE_DYNAMIC_ALLOCATION)
+ //nl_ext_double * RESTRICT m_A;
+ std::vector<nl_ext_double> m_A;
+#else
+ nl_ext_double m_A[storage_N][m_pitch];
+#endif
+ //nl_ext_double m_RHSx[storage_N];
+
+ const std::size_t m_dim;
+
+};
+
+// ----------------------------------------------------------------------------------------
+// matrix_solver_direct
+// ----------------------------------------------------------------------------------------
+
+template <std::size_t m_N, std::size_t storage_N>
+matrix_solver_direct_t<m_N, storage_N>::~matrix_solver_direct_t()
+{
+#if (NL_USE_DYNAMIC_ALLOCATION)
+ //plib::pfree_array(m_A);
+#endif
+}
- // ----------------------------------------------------------------------------------------
- // matrix_solver_direct
- // ----------------------------------------------------------------------------------------
+template <std::size_t m_N, std::size_t storage_N>
+void matrix_solver_direct_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets)
+{
+ matrix_solver_t::setup_base(nets);
- template <typename FT, int SIZE>
- void matrix_solver_direct_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets)
+ /* add RHS element */
+ for (std::size_t k = 0; k < N(); k++)
{
- matrix_solver_t::setup_base(nets);
+ terms_for_net_t * t = m_terms[k].get();
- /* 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>(N())))
+ t->m_nzrd.push_back(static_cast<unsigned>(N()));
+ }
- if (!plib::container::contains(t->m_nzrd, static_cast<unsigned>(size())))
- t->m_nzrd.push_back(static_cast<unsigned>(size()));
- }
+ netlist().save(*this, m_last_RHS, "m_last_RHS");
- // 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));
- }
+ for (std::size_t k = 0; k < N(); k++)
+ netlist().save(*this, RHS(k), plib::pfmt("RHS.{1}")(k));
+}
- template <typename FT, int SIZE>
- void matrix_solver_direct_t<FT, SIZE>::LE_solve()
+
+template <std::size_t m_N, std::size_t storage_N>
+void matrix_solver_direct_t<m_N, storage_N>::LE_solve()
+{
+ const std::size_t kN = N();
+ if (!m_params.m_pivot)
{
- const std::size_t kN = size();
- if (!m_params.m_pivot)
+ for (std::size_t i = 0; i < kN; i++)
{
- 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;
- for (std::size_t 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;
- }
+ /* FIXME: Singular matrix? */
+ nl_double *Ai = &A(i, 0);
+ const nl_double f = 1.0 / A(i,i);
+ const auto &nzrd = m_terms[i]->m_nzrd;
+ const auto &nzbd = m_terms[i]->m_nzbd;
+
+ for (std::size_t j : nzbd)
+ {
+ nl_double *Aj = &A(j, 0);
+ const nl_double f1 = -f * Aj[i];
+ for (std::size_t k : nzrd)
+ Aj[k] += Ai[k] * f1;
+ //RHS(j) += RHS(i) * f1;
}
}
- else
+ }
+ else
+ {
+ for (std::size_t i = 0; i < kN; i++)
{
- for (std::size_t i = 0; i < kN; i++)
+ /* Find the row with the largest first value */
+ std::size_t maxrow = i;
+ for (std::size_t j = i + 1; j < kN; j++)
{
- /* Find the row with the largest first value */
- std::size_t maxrow = 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 (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 (maxrow != 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));
- }
- //std::swap(RHS(i), RHS(maxrow));
+ if (maxrow != 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));
}
- /* FIXME: Singular matrix? */
- const FT f = 1.0 / A(i,i);
+ //std::swap(RHS(i), RHS(maxrow));
+ }
+ /* FIXME: Singular matrix? */
+ const nl_double f = 1.0 / A(i,i);
- /* Eliminate column i from row j */
+ /* Eliminate column i from row j */
- for (std::size_t j = i + 1; j < kN; j++)
+ for (std::size_t j = i + 1; j < kN; j++)
+ {
+ const nl_double f1 = - A(j,i) * f;
+ if (f1 != NL_FCONST(0.0))
{
- const FT f1 = - 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);
- //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
- }
+ const nl_double * RESTRICT pi = &A(i,i+1);
+ nl_double * RESTRICT pj = &A(j,i+1);
+#if 1
+ vec_add_mult_scalar_p(kN-i,pi,f1,pj);
+#else
+ vec_add_mult_scalar_p(kN-i-1,pj,f1,pi);
+ //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
}
}
}
}
+}
- template <typename FT, int SIZE>
- template <typename T>
- void matrix_solver_direct_t<FT, SIZE>::LE_back_subst(
- T & x)
- {
- const std::size_t kN = size();
+template <std::size_t m_N, std::size_t storage_N>
+template <typename T>
+void matrix_solver_direct_t<m_N, storage_N>::LE_back_subst(
+ T * RESTRICT x)
+{
+ const std::size_t kN = N();
- /* back substitution */
- if (m_params.m_pivot)
+ /* back substitution */
+ if (m_params.m_pivot)
+ {
+ for (std::size_t j = kN; j-- > 0; )
{
- for (std::size_t j = kN; j-- > 0; )
- {
- FT tmp = 0;
- for (std::size_t k = j+1; k < kN; k++)
- tmp += A(j,k) * x[k];
- x[j] = (RHS(j) - tmp) / A(j,j);
- }
+ T tmp = 0;
+ for (std::size_t k = j+1; k < kN; k++)
+ tmp += A(j,k) * x[k];
+ x[j] = (RHS(j) - tmp) / A(j,j);
}
- else
+ }
+ else
+ {
+ for (std::size_t j = kN; j-- > 0; )
{
- for (std::size_t j = kN; j-- > 0; )
+ T tmp = 0;
+
+ const auto *p = m_terms[j]->m_nzrd.data();
+ const auto e = m_terms[j]->m_nzrd.size() - 1; /* exclude RHS element */
+ T * Aj = &A(j,0);
+ for (std::size_t k = 0; k < e; k++)
{
- FT tmp = 0;
- const auto &nzrd = m_terms[j]->m_nzrd;
- const auto e = nzrd.size() - 1; /* exclude RHS element */
- 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);
+ const auto pk = p[k];
+ tmp += Aj[pk] * x[pk];
}
+ x[j] = (RHS(j) - tmp) / A(j,j);
}
}
+}
- template <typename FT, int SIZE>
- unsigned matrix_solver_direct_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphson)
- {
- 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;
- }
+template <std::size_t m_N, std::size_t storage_N>
+unsigned matrix_solver_direct_t<m_N, storage_N>::solve_non_dynamic(const bool newton_raphson)
+{
+ nl_double new_V[storage_N]; // = { 0.0 };
- template <typename FT, int SIZE>
- unsigned matrix_solver_direct_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson)
- {
- this->build_LE_A(*this);
- this->build_LE_RHS(*this);
+ this->LE_solve();
+ this->LE_back_subst(new_V);
- this->m_stat_calculations++;
- return this->solve_non_dynamic(newton_raphson);
- }
+ const nl_double 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>
- 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)
+template <std::size_t m_N, std::size_t storage_N>
+inline unsigned matrix_solver_direct_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson)
+{
+ build_LE_A<matrix_solver_direct_t>();
+ build_LE_RHS<matrix_solver_direct_t>();
+
+ for (std::size_t i=0, iN=N(); i < iN; i++)
+ m_last_RHS[i] = RHS(i);
+
+ this->m_stat_calculations++;
+ return this->solve_non_dynamic(newton_raphson);
+}
+
+template <std::size_t m_N, std::size_t storage_N>
+matrix_solver_direct_t<m_N, storage_N>::matrix_solver_direct_t(netlist_t &anetlist, const pstring &name,
+ const solver_parameters_t *params, const std::size_t size)
+: matrix_solver_t(anetlist, name, ASCENDING, params)
+, m_dim(size)
+{
+#if (NL_USE_DYNAMIC_ALLOCATION)
+ m_A.resize(N() * m_pitch);
+ //m_A = plib::palloc_array<nl_ext_double>(N() * m_pitch);
+#endif
+ for (unsigned k = 0; k < N(); k++)
{
+ m_last_RHS[k] = 0.0;
}
+}
- 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)
+template <std::size_t m_N, std::size_t storage_N>
+matrix_solver_direct_t<m_N, storage_N>::matrix_solver_direct_t(netlist_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_dim(size)
+{
+#if (NL_USE_DYNAMIC_ALLOCATION)
+ m_A.resize(N() * m_pitch);
+ //m_A = plib::palloc_array<nl_ext_double>(N() * m_pitch);
+#endif
+ for (unsigned k = 0; k < N(); k++)
{
+ m_last_RHS[k] = 0.0;
}
+}
-} // namespace devices
+ } //namespace devices
} // namespace netlist
#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..6e1f99bad7d 100644
--- a/src/lib/netlist/solver/nld_ms_direct1.h
+++ b/src/lib/netlist/solver/nld_ms_direct1.h
@@ -13,41 +13,37 @@
namespace netlist
{
-namespace devices
-{
- template <typename FT>
- class matrix_solver_direct1_t: public matrix_solver_direct_t<FT, 1>
+ namespace devices
{
- public:
-
- 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)
- {}
+class matrix_solver_direct1_t: public matrix_solver_direct_t<1,1>
+{
+public:
- // ----------------------------------------------------------------------------------------
- // matrix_solver - Direct1
- // ----------------------------------------------------------------------------------------
- unsigned vsolve_non_dynamic(const bool newton_raphson) 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]);
+ matrix_solver_direct1_t(netlist_t &anetlist, const pstring &name, const solver_parameters_t *params)
+ : matrix_solver_direct_t<1, 1>(anetlist, name, params, 1)
+ {}
+ virtual unsigned vsolve_non_dynamic(const bool newton_raphson) override;
- std::array<FT, 1> new_V = { this->RHS(0) / this->A(0,0) };
+};
- const FT err = (newton_raphson ? this->delta(new_V) : 0.0);
- this->store(new_V);
- return (err > this->m_params.m_accuracy) ? 2 : 1;
- }
+// ----------------------------------------------------------------------------------------
+// matrix_solver - Direct1
+// ----------------------------------------------------------------------------------------
- };
+inline unsigned matrix_solver_direct1_t::vsolve_non_dynamic(ATTR_UNUSED const bool newton_raphson)
+{
+ build_LE_A<matrix_solver_direct1_t>();
+ build_LE_RHS<matrix_solver_direct1_t>();
+ //NL_VERBOSE_OUT(("{1} {2}\n", new_val, m_RHS[0] / m_A[0][0]);
+ nl_double new_V[1] = { RHS(0) / A(0,0) };
+ const nl_double err = (newton_raphson ? delta(new_V) : 0.0);
+ store(new_V);
+ return (err > this->m_params.m_accuracy) ? 2 : 1;
+}
-} //namespace devices
+ } //namespace devices
} // namespace netlist
diff --git a/src/lib/netlist/solver/nld_ms_direct2.h b/src/lib/netlist/solver/nld_ms_direct2.h
index 01f77c3bc3c..4004bce9cc4 100644
--- a/src/lib/netlist/solver/nld_ms_direct2.h
+++ b/src/lib/netlist/solver/nld_ms_direct2.h
@@ -13,46 +13,43 @@
namespace netlist
{
-namespace devices
+ namespace devices
+ {
+class matrix_solver_direct2_t: public matrix_solver_direct_t<2,2>
{
+public:
- // ----------------------------------------------------------------------------------------
- // matrix_solver - Direct2
- // ----------------------------------------------------------------------------------------
-
- template <typename FT>
- class matrix_solver_direct2_t: public matrix_solver_direct_t<FT, 2>
- {
- public:
+ matrix_solver_direct2_t(netlist_t &anetlist, const pstring &name, const solver_parameters_t *params)
+ : matrix_solver_direct_t<2, 2>(anetlist, name, params, 2)
+ {}
+ virtual unsigned vsolve_non_dynamic(const bool newton_raphson) override;
- 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
- {
- this->build_LE_A(*this);
- this->build_LE_RHS(*this);
+// ----------------------------------------------------------------------------------------
+// matrix_solver - Direct2
+// ----------------------------------------------------------------------------------------
- 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);
+inline unsigned matrix_solver_direct2_t::vsolve_non_dynamic(ATTR_UNUSED const bool newton_raphson)
+{
+ build_LE_A<matrix_solver_direct2_t>();
+ build_LE_RHS<matrix_solver_direct2_t>();
- 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};
+ const nl_double a = A(0,0);
+ const nl_double b = A(0,1);
+ const nl_double c = A(1,0);
+ const nl_double d = A(1,1);
- 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;
- }
+ nl_double new_V[2];
+ new_V[1] = (a * RHS(1) - c * RHS(0)) / (a * d - b * c);
+ new_V[0] = (RHS(0) - b * new_V[1]) / a;
- };
+ const nl_double err = (newton_raphson ? delta(new_V) : 0.0);
+ store(new_V);
+ return (err > this->m_params.m_accuracy) ? 2 : 1;
+}
-} //namespace devices
+ } //namespace devices
} // namespace netlist
#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..c379cf8d055 100644
--- a/src/lib/netlist/solver/nld_ms_direct_lu.h
+++ b/src/lib/netlist/solver/nld_ms_direct_lu.h
@@ -8,11 +8,11 @@
#ifndef NLD_MS_DIRECT_H_
#define NLD_MS_DIRECT_H_
+#include <algorithm>
+
#include "solver/nld_solver.h"
#include "solver/nld_matrix_solver.h"
-#include <algorithm>
-
//#define A(r, c) m_A[_r][_c]
namespace netlist
@@ -36,9 +36,9 @@ public:
virtual void vsetup(analog_net_t::list_t &nets) override;
virtual void reset() override { matrix_solver_t::reset(); }
- unsigned N() const { if (m_N == 0) return m_dim; else return m_N; }
+ inline unsigned N() const { if (m_N == 0) return m_dim; else return m_N; }
- int vsolve_non_dynamic(const bool newton_raphson);
+ inline int vsolve_non_dynamic(const bool newton_raphson);
protected:
virtual void add_term(int net_idx, terminal_t *term) override;
@@ -139,10 +139,12 @@ protected:
nl_double compute_next_timestep();
template <typename T1, typename T2>
- nl_ext_double &A(const T1 r, const T2 c) { return m_A[r][c]; }
+ inline nl_ext_double &A(const T1 r, const T2 c) { return m_A[r][c]; }
//nl_double m_A[storage_N][((storage_N + 7) / 8) * 8];
nl_double m_RHS[storage_N];
+ nl_double m_last_RHS[storage_N]; // right hand side - contains currents
+ nl_double m_last_V[storage_N];
terms_for_net_t *m_rails_temp;
@@ -353,6 +355,7 @@ void matrix_solver_direct_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets)
* save states
*/
save(NLNAME(m_RHS));
+ save(NLNAME(m_last_RHS));
save(NLNAME(m_last_V));
for (unsigned k = 0; k < N(); k++)
@@ -587,10 +590,13 @@ 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)
+inline int matrix_solver_direct_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson)
{
this->build_LE_A();
- this->build_LE_RHS(m_RHS);
+ this->build_LE_RHS(m_last_RHS);
+
+ for (unsigned i=0, iN=N(); i < iN; i++)
+ m_RHS[i] = m_last_RHS[i];
return this->solve_non_dynamic(newton_raphson);
}
@@ -602,6 +608,11 @@ matrix_solver_direct_t<m_N, storage_N>::matrix_solver_direct_t(const solver_para
, m_lp_fact(0)
{
m_rails_temp = palloc_array(terms_for_net_t, N());
+
+ for (unsigned k = 0; k < N(); k++)
+ {
+ m_last_RHS[k] = 0.0;
+ }
}
template <unsigned m_N, unsigned storage_N>
@@ -615,6 +626,7 @@ matrix_solver_direct_t<m_N, storage_N>::matrix_solver_direct_t(const eSolverType
for (unsigned k = 0; k < N(); k++)
{
m_terms[k] = palloc(terms_for_net_t);
+ m_last_RHS[k] = 0.0;
}
}
diff --git a/src/lib/netlist/solver/nld_ms_gcr.h b/src/lib/netlist/solver/nld_ms_gcr.h
index f3e56a342d6..a17ca1ec282 100644
--- a/src/lib/netlist/solver/nld_ms_gcr.h
+++ b/src/lib/netlist/solver/nld_ms_gcr.h
@@ -10,201 +10,390 @@
#ifndef NLD_MS_GCR_H_
#define NLD_MS_GCR_H_
-#include "plib/mat_cr.h"
+#include <algorithm>
+#include "../plib/pdynlib.h"
+#include "mat_cr.h"
#include "nld_ms_direct.h"
#include "nld_solver.h"
-#include "plib/pdynlib.h"
-#include "plib/pstream.h"
-#include "plib/vector_ops.h"
-
-#include <algorithm>
+#include "vector_base.h"
+#include "../plib/pstream.h"
namespace netlist
{
-namespace devices
+ namespace devices
+ {
+template <std::size_t m_N, std::size_t storage_N>
+class matrix_solver_GCR_t: public matrix_solver_t
{
+public:
+
+ matrix_solver_GCR_t(netlist_t &anetlist, const pstring &name,
+ const solver_parameters_t *params, const std::size_t size)
+ : matrix_solver_t(anetlist, name, matrix_solver_t::ASCENDING, params)
+ , m_dim(size)
+ , mat(size)
+ , m_proc()
+ {
+ }
- template <typename FT, int SIZE>
- class matrix_solver_GCR_t: public matrix_solver_t
+ virtual ~matrix_solver_GCR_t() override
{
- 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()
- {
- }
+ }
+
+ constexpr std::size_t N() const { return (m_N == 0) ? m_dim : m_N; }
- constexpr std::size_t N() const { return m_dim; }
+ virtual void vsetup(analog_net_t::list_t &nets) override;
+ virtual unsigned vsolve_non_dynamic(const bool newton_raphson) override;
- void vsetup(analog_net_t::list_t &nets) override;
- unsigned vsolve_non_dynamic(const bool newton_raphson) override;
+ virtual std::pair<pstring, pstring> create_solver_code() override;
- std::pair<pstring, pstring> create_solver_code() override;
+private:
- private:
+ //typedef typename mat_cr_t<storage_N>::type mattype;
+ typedef typename mat_cr_t<storage_N>::index_type mattype;
- using mat_index_type = typename plib::matrix_compressed_rows_t<FT, SIZE>::index_type;
+ void csc_private(plib::putf8_fmt_writer &strm);
- void csc_private(plib::putf8_fmt_writer &strm);
+ using extsolver = void (*)(double * RESTRICT m_A, double * RESTRICT RHS, double * RESTRICT V);
- using extsolver = void (*)(double * m_A, double * RHS, double * V);
+ pstring static_compile_name();
- pstring static_compile_name();
+ const std::size_t m_dim;
+ std::vector<unsigned> m_term_cr[storage_N];
+ mat_cr_t<storage_N> mat;
- const std::size_t m_dim;
- plib::parray<FT, SIZE> RHS;
- plib::parray<FT, SIZE> new_V;
+ //extsolver m_proc;
+ plib::dynproc<void, double * RESTRICT, double * RESTRICT, double * RESTRICT> m_proc;
- mat_type mat;
+};
- //extsolver m_proc;
- plib::dynproc<void, double * , double * , double * > m_proc;
+// ----------------------------------------------------------------------------------------
+// matrix_solver - GCR
+// ----------------------------------------------------------------------------------------
+
+template <std::size_t m_N, std::size_t storage_N>
+void matrix_solver_GCR_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets)
+{
+ setup_base(nets);
- };
+ mattype nz = 0;
+ const std::size_t iN = this->N();
- // ----------------------------------------------------------------------------------------
- // matrix_solver - GCR
- // ----------------------------------------------------------------------------------------
+ /* build the final matrix */
- // FIXME: namespace or static class member
- template <typename V>
- std::size_t inline get_level(const V &v, std::size_t k)
+ bool touched[storage_N][storage_N] = { { false } };
+ for (std::size_t k = 0; k < iN; 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");
+ for (auto &j : this->m_terms[k]->m_nz)
+ touched[k][j] = true;
}
- template <typename FT, int SIZE>
- void matrix_solver_GCR_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets)
- {
- setup_base(nets);
+ unsigned fc = 0;
- const std::size_t iN = this->N();
+ unsigned ops = 0;
- /* build the final matrix */
+ for (std::size_t k = 0; k < iN; k++)
+ {
+ ops++; // 1/A(k,k)
+ for (std::size_t row = k + 1; row < iN; row++)
+ {
+ if (touched[row][k])
+ {
+ ops++;
+ fc++;
+ for (std::size_t col = k + 1; col < iN; col++)
+ if (touched[k][col])
+ {
+ touched[row][col] = true;
+ ops += 2;
+ }
+ }
+ }
+ }
- std::vector<std::vector<unsigned>> fill(iN);
- std::size_t raw_elements = 0;
+ for (mattype k=0; k<iN; k++)
+ {
+ mat.ia[k] = nz;
- for (std::size_t k = 0; k < iN; k++)
+ for (mattype j=0; j<iN; j++)
{
- fill[k].resize(iN, decltype(mat)::FILL_INFINITY);
- for (auto &j : this->m_terms[k]->m_nz)
+ if (touched[k][j])
{
- fill[k][j] = 0;
- raw_elements++;
+ mat.ja[nz] = j;
+ if (j == k)
+ mat.diag[k] = nz;
+ nz++;
}
+ }
+ m_term_cr[k].clear();
+ /* 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++)
+ {
+ int other = this->m_terms[k]->connected_net_idx()[j];
+ for (auto i = mat.ia[k]; i < nz; i++)
+ if (other == static_cast<int>(mat.ja[i]))
+ {
+ m_term_cr[k].push_back(i);
+ break;
+ }
}
+ nl_assert(m_term_cr[k].size() == this->m_terms[k]->m_railstart);
+ }
- auto gr = mat.gaussian_extend_fill_mat(fill);
+ mat.ia[iN] = nz;
+ mat.nz_num = nz;
- /* FIXME: move this to the cr matrix class and use computed
- * parallel ordering once it makes sense.
- */
+ this->log().verbose("Ops: {1} Occupancy ratio: {2}\n", ops,
+ static_cast<double>(nz) / static_cast<double>(iN * iN));
+
+ // FIXME: Move me
+
+ if (netlist().lib().isLoaded())
+ {
+ pstring symname = static_compile_name();
+#if 0
+ m_proc = this->netlist().lib().template getsym<extsolver>(symname);
+ if (m_proc != nullptr)
+ this->log().verbose("External static solver {1} found ...", symname);
+ else
+ this->log().warning("External static solver {1} not found ...", symname);
+#else
+ m_proc.load(this->netlist().lib(), symname);
+ if (m_proc.resolved())
+ this->log().warning("External static solver {1} found ...", symname);
+ else
+ this->log().warning("External static solver {1} not found ...", symname);
+#endif
+ }
- std::vector<unsigned> levL(iN, 0);
- std::vector<unsigned> levU(iN, 0);
+}
+#if 0
+template <std::size_t m_N, std::size_t storage_N>
+void matrix_solver_GCR_t<m_N, storage_N>::csc_private(plib::putf8_fmt_writer &strm)
+{
+ const std::size_t iN = N();
+ for (std::size_t i = 0; i < iN - 1; i++)
+ {
+ const auto &nzbd = this->m_terms[i]->m_nzbd;
- // parallel scheme for L x = y
- for (std::size_t k = 0; k < iN; k++)
+ if (nzbd.size() > 0)
{
- 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;
+ std::size_t pi = mat.diag[i];
+
+ //const nl_double f = 1.0 / m_A[pi++];
+ strm("const double f{1} = 1.0 / m_A[{2}];\n", i, pi);
+ pi++;
+ const std::size_t piie = mat.ia[i+1];
+
+ //for (auto & j : nzbd)
+ for (std::size_t j : nzbd)
+ {
+ // proceed to column i
+ std::size_t pj = mat.ia[j];
+
+ while (mat.ja[pj] < i)
+ pj++;
+
+ //const nl_double f1 = - m_A[pj++] * f;
+ strm("\tconst double f{1}_{2} = -f{3} * m_A[{4}];\n", i, j, i, pj);
+ pj++;
+
+ // subtract row i from j */
+ for (std::size_t pii = pi; pii<piie; )
+ {
+ while (mat.ja[pj] < mat.ja[pii])
+ pj++;
+ //m_A[pj++] += m_A[pii++] * f1;
+ 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);
+ }
}
+ }
- // parallel scheme for U x = y
- for (std::size_t k = iN; k-- > 0; )
+ //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]);
+ for (std::size_t j = iN - 1; j-- > 0;)
+ {
+ strm("\tdouble tmp{1} = 0.0;\n", j);
+ const std::size_t e = mat.ia[j+1];
+ for (std::size_t pk = mat.diag[j] + 1; pk < e; pk++)
{
- 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;
+ strm("\ttmp{1} += m_A[{2}] * V[{3}];\n", j, pk, mat.ja[pk]);
}
+ strm("\tV[{1}] = (RHS[{1}] - tmp{1}) / m_A[{4}];\n", j, j, j, mat.diag[j]);
+ }
+}
+#else
+template <std::size_t m_N, std::size_t storage_N>
+void matrix_solver_GCR_t<m_N, storage_N>::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);
- for (std::size_t k = 0; k < iN; k++)
+ for (std::size_t i = 0; i < iN - 1; i++)
+ {
+ const auto &nzbd = this->m_terms[i]->m_nzbd;
+
+ if (nzbd.size() > 0)
{
- unsigned fm = 0;
- pstring ml = "";
- for (std::size_t j = 0; j < iN; j++)
+ std::size_t pi = mat.diag[i];
+
+ //const nl_double f = 1.0 / m_A[pi++];
+ strm("const double f{1} = 1.0 / m_A{2};\n", i, pi);
+ pi++;
+ const std::size_t piie = mat.ia[i+1];
+
+ //for (auto & j : nzbd)
+ for (std::size_t j : nzbd)
{
- 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];
+ // proceed to column i
+ std::size_t pj = mat.ia[j];
+
+ while (mat.ja[pj] < i)
+ pj++;
+
+ //const nl_double f1 = - m_A[pj++] * f;
+ strm("\tconst double f{1}_{2} = -f{3} * m_A{4};\n", i, j, i, pj);
+ pj++;
+
+ // subtract row i from j */
+ for (std::size_t pii = pi; pii<piie; )
+ {
+ while (mat.ja[pj] < mat.ja[pii])
+ pj++;
+ //m_A[pj++] += m_A[pii++] * f1;
+ 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);
}
- 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);
}
+ }
+
+ //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]);
+ for (std::size_t j = iN - 1; j-- > 0;)
+ {
+ strm("\tdouble tmp{1} = 0.0;\n", j);
+ const std::size_t e = mat.ia[j+1];
+ for (std::size_t pk = mat.diag[j] + 1; pk < e; pk++)
+ {
+ strm("\ttmp{1} += m_A{2} * V[{3}];\n", j, pk, mat.ja[pk]);
+ }
+ strm("\tV[{1}] = (RHS[{1}] - tmp{1}) / m_A{4};\n", j, j, j, mat.diag[j]);
+ }
+}
+#endif
+
+template <std::size_t m_N, std::size_t storage_N>
+pstring matrix_solver_GCR_t<m_N, storage_N>::static_compile_name()
+{
+ plib::postringstream t;
+ plib::putf8_fmt_writer w(t);
+ csc_private(w);
+ std::hash<pstring> h;
+
+ return plib::pfmt("nl_gcr_{1:x}_{2}")(h( t.str() ))(mat.nz_num);
+}
+
+template <std::size_t m_N, std::size_t storage_N>
+std::pair<pstring, pstring> matrix_solver_GCR_t<m_N, storage_N>::create_solver_code()
+{
+ plib::postringstream t;
+ 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, t.str());
+}
- mat.build_from_fill_mat(fill);
+template <std::size_t m_N, std::size_t storage_N>
+unsigned matrix_solver_GCR_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson)
+{
+ const std::size_t iN = this->N();
+
+ nl_double RHS[storage_N];
+ nl_double new_V[storage_N];
- for (mat_index_type k=0; k<iN; k++)
+ mat.set_scalar(0.0);
+
+ for (std::size_t k = 0; k < iN; k++)
+ {
+ terms_for_net_t *t = this->m_terms[k].get();
+ nl_double gtot_t = 0.0;
+ nl_double RHS_t = 0.0;
+
+ const std::size_t term_count = t->count();
+ const std::size_t railstart = t->m_railstart;
+ const nl_double * const RESTRICT gt = t->gt();
+ const nl_double * const RESTRICT go = t->go();
+ const nl_double * const RESTRICT Idr = t->Idr();
+ const nl_double * const * RESTRICT other_cur_analog = t->connected_net_V();
+ const unsigned * const RESTRICT tcr = m_term_cr[k].data();
+
+#if 0
+ for (std::size_t i = 0; i < term_count; i++)
{
- 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++)
- {
- 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;
- }
- }
- nl_assert(cnt == this->m_terms[k]->m_railstart);
- m_mat_ptr[k][this->m_terms[k]->m_railstart] = &mat.A[mat.diag[k]];
+ gtot_t += gt[i];
+ RHS_t += Idr[i];
}
- 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));
+ for (std::size_t i = railstart; i < term_count; i++)
+ RHS_t += go[i] * *other_cur_analog[i];
- // FIXME: Move me
+ RHS[k] = RHS_t;
- if (state().lib().isLoaded())
+ // add diagonal element
+ mat.A[mat.diag[k]] = gtot_t;
+
+ for (std::size_t i = 0; i < railstart; i++)
+ mat.A[tcr[i]] -= go[i];
+ }
+#else
+ for (std::size_t i = 0; i < railstart; i++)
+ mat.A[tcr[i]] -= go[i];
+
+ for (std::size_t i = 0; i < railstart; i++)
{
- pstring symname = static_compile_name();
- m_proc.load(this->state().lib(), symname);
- if (m_proc.resolved())
- this->log().warning("External static solver {1} found ...", symname);
- else
- this->log().warning("External static solver {1} not found ...", symname);
+ gtot_t += gt[i];
+ RHS_t += Idr[i];
}
- }
+ for (std::size_t i = railstart; i < term_count; i++)
+ {
+ RHS_t += (Idr[i] + go[i] * *other_cur_analog[i]);
+ gtot_t += gt[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();
+ RHS[k] = RHS_t;
+ mat.A[mat.diag[k]] += gtot_t;
+ }
+#endif
+ mat.ia[iN] = static_cast<mattype>(mat.nz_num);
- for (std::size_t i = 0; i < mat.nz_num; i++)
- strm("double m_A{1} = m_A[{2}];\n", i, i);
+ /* 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]);
+ }
+ else
+ {
for (std::size_t i = 0; i < iN - 1; i++)
{
const auto &nzbd = this->m_terms[i]->m_nzbd;
@@ -212,114 +401,61 @@ namespace devices
if (nzbd.size() > 0)
{
std::size_t pi = mat.diag[i];
+ const nl_double f = 1.0 / mat.A[pi++];
+ const std::size_t piie = mat.ia[i+1];
- //const FT f = 1.0 / m_A[pi++];
- strm("const double f{1} = 1.0 / m_A{2};\n", i, 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 j : nzbd) // for (std::size_t j = i + 1; j < iN; j++)
{
// proceed to column i
- std::size_t pj = mat.row_idx[j];
+ //__builtin_prefetch(&m_A[mat.diag[j+1]], 1);
+ std::size_t pj = mat.ia[j];
- while (mat.col_idx[pj] < i)
+ while (mat.ja[pj] < i)
pj++;
- //const FT f1 = - m_A[pj++] * f;
- strm("\tconst double f{1}_{2} = -f{3} * m_A{4};\n", i, j, i, pj);
- pj++;
+ const nl_double f1 = - mat.A[pj++] * f;
// subtract row i from j */
for (std::size_t pii = pi; pii<piie; )
{
- while (mat.col_idx[pj] < mat.col_idx[pii])
+ while (mat.ja[pj] < mat.ja[pii])
pj++;
- //m_A[pj++] += m_A[pii++] * f1;
- strm("\tm_A{1} += m_A{2} * f{3}_{4};\n", pj, pii, i, j);
- pj++; pii++;
+ mat.A[pj++] += mat.A[pii++] * f1;
}
- //RHS[j] += f1 * RHS[i];
- strm("\tRHS[{1}] += f{2}_{3} * RHS[{4}];\n", j, i, j, i);
+ RHS[j] += f1 * RHS[i];
}
}
}
+ /* backward substitution
+ *
+ */
+
+ /* row n-1 */
+ new_V[iN - 1] = RHS[iN - 1] / mat.A[mat.diag[iN - 1]];
- //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]);
for (std::size_t j = iN - 1; j-- > 0;)
{
- strm("\tdouble tmp{1} = 0.0;\n", j);
- const std::size_t e = mat.row_idx[j+1];
- for (std::size_t pk = mat.diag[j] + 1; pk < e; pk++)
+ //__builtin_prefetch(&new_V[j-1], 1);
+ //if (j>0)__builtin_prefetch(&m_A[mat.diag[j-1]], 0);
+ double tmp = 0;
+ auto jdiag = mat.diag[j];
+ const std::size_t e = mat.ia[j+1];
+ for (std::size_t pk = jdiag + 1; pk < e; pk++)
{
- strm("\ttmp{1} += m_A{2} * V[{3}];\n", j, pk, mat.col_idx[pk]);
+ tmp += mat.A[pk] * new_V[mat.ja[pk]];
}
- strm("\tV[{1}] = (RHS[{1}] - tmp{1}) / m_A{4};\n", j, j, j, mat.diag[j]);
+ new_V[j] = (RHS[j] - tmp) / mat.A[jdiag];
}
}
- template <typename FT, int SIZE>
- pstring matrix_solver_GCR_t<FT, SIZE>::static_compile_name()
- {
- plib::postringstream t;
- plib::putf8_fmt_writer w(&t);
- csc_private(w);
- std::hash<pstring> h;
-
- return plib::pfmt("nl_gcr_{1:x}_{2}")(h( t.str() ))(mat.nz_num);
- }
+ this->m_stat_calculations++;
- template <typename FT, int SIZE>
- std::pair<pstring, pstring> matrix_solver_GCR_t<FT, SIZE>::create_solver_code()
- {
- plib::postringstream t;
- 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, t.str());
- }
-
- template <typename FT, int SIZE>
- unsigned matrix_solver_GCR_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson)
- {
- 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]);
- }
- else
- {
- // mat.gaussian_elimination_parallel(RHS);
- mat.gaussian_elimination(RHS);
- /* backward substitution */
- mat.gaussian_back_substitution(new_V, RHS);
- }
-
- this->m_stat_calculations++;
-
- const FT err = (newton_raphson ? delta(new_V) : 0.0);
- store(new_V);
- return (err > this->m_params.m_accuracy) ? 2 : 1;
- }
+ const nl_double err = (newton_raphson ? delta(new_V) : 0.0);
+ store(new_V);
+ return (err > this->m_params.m_accuracy) ? 2 : 1;
+}
-} // namespace devices
+ } //namespace devices
} // namespace netlist
#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..2e4e447d14f 100644
--- a/src/lib/netlist/solver/nld_ms_gmres.h
+++ b/src/lib/netlist/solver/nld_ms_gmres.h
@@ -1,145 +1,387 @@
// license:GPL-2.0+
// copyright-holders:Couriersud
/*
- * nld_ms_gmres.h
+ * nld_ms_sor.h
+ *
+ * Generic successive over relaxation solver.
+ *
+ * Fow w==1 we will do the classic Gauss-Seidel approach
*
*/
#ifndef NLD_MS_GMRES_H_
#define NLD_MS_GMRES_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/vector_ops.h"
-
#include <algorithm>
-#include <cmath>
+#include "mat_cr.h"
+#include "nld_ms_direct.h"
+#include "nld_solver.h"
+#include "vector_base.h"
namespace netlist
{
-namespace devices
+ namespace devices
+ {
+template <std::size_t m_N, std::size_t storage_N>
+class matrix_solver_GMRES_t: public matrix_solver_direct_t<m_N, storage_N>
{
+public:
+
+ matrix_solver_GMRES_t(netlist_t &anetlist, const pstring &name, const solver_parameters_t *params, const std::size_t size)
+ : matrix_solver_direct_t<m_N, storage_N>(anetlist, name, matrix_solver_t::ASCENDING, params, size)
+ , m_use_iLU_preconditioning(true)
+ , m_use_more_precise_stop_condition(false)
+ , m_accuracy_mult(1.0)
+ , mat(size)
+ {
+ }
- template <typename FT, int SIZE>
- class matrix_solver_GMRES_t: public matrix_solver_direct_t<FT, SIZE>
+ virtual ~matrix_solver_GMRES_t() override
{
- public:
+ }
- using float_type = FT;
+ virtual void vsetup(analog_net_t::list_t &nets) override;
+ virtual unsigned vsolve_non_dynamic(const bool newton_raphson) override;
- /* 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)
- , m_gmres(size)
- {
- }
+private:
+
+ //typedef typename mat_cr_t<storage_N>::type mattype;
+ typedef typename mat_cr_t<storage_N>::index_type mattype;
+
+ unsigned solve_ilu_gmres(nl_double (& RESTRICT x)[storage_N], const nl_double (& RESTRICT rhs)[storage_N], const unsigned restart_max, std::size_t mr, nl_double accuracy);
+
+ std::vector<unsigned> m_term_cr[storage_N];
+
+ bool m_use_iLU_preconditioning;
+ bool m_use_more_precise_stop_condition;
+ nl_double m_accuracy_mult; // FXIME: Save state
- void vsetup(analog_net_t::list_t &nets) override;
- unsigned vsolve_non_dynamic(const bool newton_raphson) override;
+ mat_cr_t<storage_N> mat;
- private:
+ nl_double m_LU[storage_N * storage_N];
- using mattype = typename plib::matrix_compressed_rows_t<FT, SIZE>::index_type;
+ nl_double m_c[storage_N + 1]; /* mr + 1 */
+ nl_double m_g[storage_N + 1]; /* mr + 1 */
+ nl_double m_ht[storage_N + 1][storage_N]; /* (mr + 1), mr */
+ nl_double m_s[storage_N + 1]; /* mr + 1 */
+ nl_double m_v[storage_N + 1][storage_N]; /*(mr + 1), n */
+ nl_double m_y[storage_N + 1]; /* mr + 1 */
- plib::mat_precondition_ILU<FT, SIZE> m_ops;
- plib::gmres_t<FT, SIZE> m_gmres;
- };
+};
- // ----------------------------------------------------------------------------------------
- // matrix_solver - GMRES
- // ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+// matrix_solver - GMRES
+// ----------------------------------------------------------------------------------------
- template <typename FT, int SIZE>
- void matrix_solver_GMRES_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets)
+template <std::size_t m_N, std::size_t storage_N>
+void matrix_solver_GMRES_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets)
+{
+ matrix_solver_direct_t<m_N, storage_N>::vsetup(nets);
+
+ mattype nz = 0;
+ const std::size_t iN = this->N();
+
+ for (std::size_t k=0; k<iN; k++)
{
- matrix_solver_direct_t<FT, SIZE>::vsetup(nets);
+ terms_for_net_t * RESTRICT row = this->m_terms[k].get();
+ mat.ia[k] = nz;
- const std::size_t iN = this->size();
+ for (std::size_t j=0; j<row->m_nz.size(); j++)
+ {
+ mat.ja[nz] = static_cast<mattype>(row->m_nz[j]);
+ if (row->m_nz[j] == k)
+ mat.diag[k] = nz;
+ nz++;
+ }
- std::vector<std::vector<unsigned>> fill(iN);
+ /* build pointers into the compressed row format matrix for each terminal */
- for (std::size_t k=0; k<iN; k++)
+ for (unsigned j=0; j< this->m_terms[k]->m_railstart;j++)
{
- 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;
- }
+ for (unsigned i = mat.ia[k]; i<nz; i++)
+ if (this->m_terms[k]->connected_net_idx()[j] == static_cast<int>(mat.ja[i]))
+ {
+ m_term_cr[k].push_back(i);
+ break;
+ }
+ nl_assert(m_term_cr[k].size() == this->m_terms[k]->m_railstart);
}
+ }
- m_ops.build(fill);
+ mat.ia[iN] = nz;
+ mat.nz_num = nz;
+}
- /* build pointers into the compressed row format matrix for each terminal */
+template <std::size_t m_N, std::size_t storage_N>
+unsigned matrix_solver_GMRES_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson)
+{
+ const std::size_t iN = this->N();
+
+ /* 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))
+ */
- for (std::size_t k=0; k<iN; k++)
+ //nz_num = 0;
+ nl_double RHS[storage_N];
+ nl_double new_V[storage_N];
+
+ mat.set_scalar(0.0);
+
+ for (std::size_t k = 0; k < iN; k++)
+ {
+ nl_double gtot_t = 0.0;
+ nl_double RHS_t = 0.0;
+
+ const std::size_t term_count = this->m_terms[k]->count();
+ const std::size_t railstart = this->m_terms[k]->m_railstart;
+ const nl_double * const RESTRICT gt = this->m_terms[k]->gt();
+ const nl_double * const RESTRICT go = this->m_terms[k]->go();
+ const nl_double * const RESTRICT Idr = this->m_terms[k]->Idr();
+ const nl_double * const * RESTRICT other_cur_analog = this->m_terms[k]->connected_net_V();
+
+ for (std::size_t i = 0; i < term_count; i++)
{
- 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]];
+ gtot_t = gtot_t + gt[i];
+ RHS_t = RHS_t + Idr[i];
}
+
+ for (std::size_t i = railstart; i < term_count; i++)
+ RHS_t = RHS_t + go[i] * *other_cur_analog[i];
+
+ RHS[k] = RHS_t;
+
+ // add diagonal element
+ mat.A[mat.diag[k]] = gtot_t;
+
+ for (std::size_t i = 0; i < railstart; i++)
+ {
+ const std::size_t pi = m_term_cr[k][i];
+ mat.A[pi] -= go[i];
+ }
+
+ new_V[k] = this->m_nets[k]->Q_Analog();
+
}
+ mat.ia[iN] = static_cast<mattype>(mat.nz_num);
+
+ const nl_double accuracy = this->m_params.m_accuracy;
+
+ unsigned mr = iN;
+ if (iN > 3 )
+ mr = static_cast<unsigned>(std::sqrt(iN) * 2.0);
+ unsigned iter = std::max(1u, this->m_params.m_gs_loops);
+ unsigned gsl = solve_ilu_gmres(new_V, RHS, iter, mr, accuracy);
+ unsigned failed = mr * iter;
- template <typename FT, int SIZE>
- unsigned matrix_solver_GMRES_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson)
+ this->m_iterative_total += gsl;
+ this->m_stat_calculations++;
+
+ if (gsl>=failed)
+ {
+ this->m_iterative_fail++;
+ return matrix_solver_direct_t<m_N, storage_N>::vsolve_non_dynamic(newton_raphson);
+ }
+
+ const nl_double err = (newton_raphson ? this->delta(new_V) : 0.0);
+ this->store(new_V);
+ return (err > this->m_params.m_accuracy) ? 2 : 1;
+}
+
+template <typename T>
+inline static void givens_mult( const T c, const T s, T & g0, T & g1 )
+{
+ const T tg0 = c * g0 - s * g1;
+ const T tg1 = s * g0 + c * g1;
+
+ g0 = tg0;
+ g1 = tg1;
+}
+
+template <std::size_t m_N, std::size_t storage_N>
+unsigned matrix_solver_GMRES_t<m_N, storage_N>::solve_ilu_gmres (nl_double (& RESTRICT x)[storage_N], const nl_double (& RESTRICT rhs)[storage_N], const unsigned restart_max, std::size_t mr, nl_double accuracy)
+{
+ /*-------------------------------------------------------------------------
+ * The code below was inspired by code published by John Burkardt under
+ * the LPGL here:
+ *
+ * http://people.sc.fsu.edu/~jburkardt/cpp_src/mgmres/mgmres.html
+ *
+ * The code below was completely written from scratch based on the pseudo code
+ * found here:
+ *
+ * http://de.wikipedia.org/wiki/GMRES-Verfahren
+ *
+ * The Algorithm itself is described in
+ *
+ * Yousef Saad,
+ * Iterative Methods for Sparse Linear Systems,
+ * Second Edition,
+ * SIAM, 20003,
+ * ISBN: 0898715342,
+ * LC: QA188.S17.
+ *
+ *------------------------------------------------------------------------*/
+
+ unsigned itr_used = 0;
+ double rho_delta = 0.0;
+
+ const std::size_t n = this->N();
+
+ if (mr > n) mr = n;
+
+ if (m_use_iLU_preconditioning)
+ mat.incomplete_LU_factorization(m_LU);
+
+ if (m_use_more_precise_stop_condition)
{
- const std::size_t iN = this->size();
+ /* derive residual for a given delta x
+ *
+ * LU y = A dx
+ *
+ * ==> rho / accuracy = sqrt(y * y)
+ *
+ * This approach will approximate the iterative stop condition
+ * based |xnew - xold| pretty precisely. But it is slow, or expressed
+ * differently: The invest doesn't pay off.
+ * Therefore we use the approach in the else part.
+ */
+ nl_double t[storage_N];
+ nl_double Ax[storage_N];
+ vec_set(n, accuracy, t);
+ mat.mult_vec(t, Ax);
+
+ mat.solveLUx(m_LU, Ax);
+
+ const nl_double rho_to_accuracy = std::sqrt(vec_mult2(n, Ax)) / accuracy;
+
+ rho_delta = accuracy * rho_to_accuracy;
+ }
+ else
+ rho_delta = accuracy * std::sqrt(n) * m_accuracy_mult;
+
+ for (unsigned itr = 0; itr < restart_max; itr++)
+ {
+ std::size_t last_k = mr;
+ nl_double rho;
+
+ nl_double Ax[storage_N];
+ nl_double residual[storage_N];
+
+ mat.mult_vec(x, Ax);
+
+ vec_sub(n, rhs, Ax, residual);
+
+ if (m_use_iLU_preconditioning)
+ {
+ mat.solveLUx(m_LU, residual);
+ }
+
+ rho = std::sqrt(vec_mult2(n, residual));
- plib::parray<FT, SIZE> RHS(iN);
- //float_type new_V[storage_N];
+ if (rho < rho_delta)
+ return itr_used + 1;
- m_ops.m_mat.set_scalar(0.0);
+ vec_set(mr+1, NL_FCONST(0.0), m_g);
+ m_g[0] = rho;
- /* populate matrix and V for first estimate */
- this->fill_matrix(iN, this->m_mat_ptr, RHS);
+ for (std::size_t i = 0; i < mr; i++)
+ vec_set(mr + 1, NL_FCONST(0.0), m_ht[i]);
- for (std::size_t k = 0; k < iN; k++)
+ vec_mult_scalar(n, residual, NL_FCONST(1.0) / rho, m_v[0]);
+
+ for (std::size_t k = 0; k < mr; k++)
{
- this->m_new_V[k] = this->m_nets[k]->Q_Analog();
+ const std::size_t k1 = k + 1;
+
+ mat.mult_vec(m_v[k], m_v[k1]);
+
+ if (m_use_iLU_preconditioning)
+ mat.solveLUx(m_LU, m_v[k1]);
+
+ for (std::size_t j = 0; j <= k; j++)
+ {
+ m_ht[j][k] = vec_mult(n, m_v[k1], m_v[j]);
+ vec_add_mult_scalar(n, m_v[j], -m_ht[j][k], m_v[k1]);
+ }
+ m_ht[k1][k] = std::sqrt(vec_mult2(n, m_v[k1]));
+
+ if (m_ht[k1][k] != 0.0)
+ vec_scale(n, m_v[k1], NL_FCONST(1.0) / m_ht[k1][k]);
+
+ for (std::size_t j = 0; j < k; j++)
+ givens_mult(m_c[j], m_s[j], m_ht[j][k], m_ht[j+1][k]);
+
+ const nl_double mu = 1.0 / std::hypot(m_ht[k][k], m_ht[k1][k]);
+
+ m_c[k] = m_ht[k][k] * mu;
+ m_s[k] = -m_ht[k1][k] * mu;
+ m_ht[k][k] = m_c[k] * m_ht[k][k] - m_s[k] * m_ht[k1][k];
+ m_ht[k1][k] = 0.0;
+
+ givens_mult(m_c[k], m_s[k], m_g[k], m_g[k1]);
+
+ rho = std::abs(m_g[k1]);
+
+ itr_used = itr_used + 1;
+
+ if (rho <= rho_delta)
+ {
+ last_k = k;
+ break;
+ }
}
- const float_type accuracy = this->m_params.m_accuracy;
+ if (last_k >= mr)
+ /* didn't converge within accuracy */
+ last_k = mr - 1;
- 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);
+ /* Solve the system H * y = g */
+ /* x += m_v[j] * m_y[j] */
+ for (std::size_t i = last_k + 1; i-- > 0;)
+ {
+ double tmp = m_g[i];
+ for (std::size_t j = i + 1; j <= last_k; j++)
+ {
+ tmp -= m_ht[i][j] * m_y[j];
+ }
+ m_y[i] = tmp / m_ht[i][i];
+ }
+
+ for (std::size_t i = 0; i <= last_k; i++)
+ vec_add_mult_scalar(n, m_v[i], m_y[i], x);
- this->m_iterative_total += gsl;
- this->m_stat_calculations++;
+#if 1
+ if (rho <= rho_delta)
+ {
+ break;
+ }
+#else
+ /* we try to approximate the x difference between to steps using m_v[last_k] */
- if (gsl > iter)
+ double xdelta = m_y[last_k] * vec_maxabs(n, m_v[last_k]);
+ if (xdelta < accuracy)
{
- this->m_iterative_fail++;
- return matrix_solver_direct_t<FT, SIZE>::vsolve_non_dynamic(newton_raphson);
+ if (m_accuracy_mult < 16384.0)
+ m_accuracy_mult = m_accuracy_mult * 2.0;
+ break;
}
+ else
+ m_accuracy_mult = m_accuracy_mult / 2.0;
- 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;
+#endif
}
+ return itr_used;
+}
-} // namespace devices
+ } //namespace devices
} // namespace netlist
#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 d85ab0044f4..ed95cd2dc03 100644
--- a/src/lib/netlist/solver/nld_ms_sm.h
+++ b/src/lib/netlist/solver/nld_ms_sm.h
@@ -33,275 +33,293 @@
#ifndef NLD_MS_SM_H_
#define NLD_MS_SM_H_
-#include "nld_matrix_solver.h"
-#include "nld_solver.h"
-#include "plib/vector_ops.h"
-
#include <algorithm>
+#include "nld_solver.h"
+#include "nld_matrix_solver.h"
+#include "vector_base.h"
+
namespace netlist
{
-namespace devices
+ namespace devices
+ {
+//#define nl_ext_double _float128 // slow, very slow
+//#define nl_ext_double long double // slightly slower
+#define nl_ext_double nl_double
+
+template <std::size_t m_N, std::size_t storage_N>
+class matrix_solver_sm_t: public matrix_solver_t
{
+ friend class matrix_solver_t;
- template <typename FT, int SIZE>
- class matrix_solver_sm_t: public matrix_solver_t
- {
- friend class matrix_solver_t;
+public:
- public:
+ matrix_solver_sm_t(netlist_t &anetlist, const pstring &name,
+ const solver_parameters_t *params, const std::size_t size);
- 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;
+ virtual ~matrix_solver_sm_t() override;
- matrix_solver_sm_t(netlist_state_t &anetlist, const pstring &name,
- const solver_parameters_t *params, const std::size_t size);
+ virtual void vsetup(analog_net_t::list_t &nets) override;
+ virtual void reset() override { matrix_solver_t::reset(); }
- void vsetup(analog_net_t::list_t &nets) override;
- void reset() override { matrix_solver_t::reset(); }
+protected:
+ virtual unsigned vsolve_non_dynamic(const bool newton_raphson) override;
+ unsigned solve_non_dynamic(const bool newton_raphson);
- protected:
- unsigned vsolve_non_dynamic(const bool newton_raphson) override;
- unsigned solve_non_dynamic(const bool newton_raphson);
+ constexpr std::size_t N() const { return (m_N == 0) ? m_dim : m_N; }
- constexpr std::size_t size() const { return m_dim; }
+ void LE_invert();
- void LE_invert();
+ template <typename T>
+ void LE_compute_x(T * RESTRICT x);
- template <typename T>
- void LE_compute_x(T * x);
+ template <typename T1, typename T2>
+ nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r][c]; }
+ template <typename T1, typename T2>
+ nl_ext_double &W(const T1 &r, const T2 &c) { return m_W[r][c]; }
+ template <typename T1, typename T2>
+ nl_ext_double &Ainv(const T1 &r, const T2 &c) { return m_Ainv[r][c]; }
+ template <typename T1>
+ nl_ext_double &RHS(const T1 &r) { return m_RHS[r]; }
- 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]; }
- 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]; }
+ template <typename T1, typename T2>
+ nl_ext_double &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; }
+ template <typename T1, typename T2>
+ nl_ext_double &lAinv(const T1 &r, const T2 &c) { return m_lAinv[r][c]; }
- template <typename T1, typename T2>
- float_ext_type &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; }
- template <typename T1, typename T2>
- float_ext_type &lAinv(const T1 &r, const T2 &c) { return m_lAinv[r][c]; }
+ nl_double m_last_RHS[storage_N]; // right hand side - contains currents
- 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];
- float_ext_type m_RHS[storage_N]; // right hand side - contains currents
+private:
+ static constexpr std::size_t m_pitch = ((( storage_N) + 7) / 8) * 8;
+ nl_ext_double m_A[storage_N][m_pitch];
+ nl_ext_double m_Ainv[storage_N][m_pitch];
+ nl_ext_double m_W[storage_N][m_pitch];
+ nl_ext_double m_RHS[storage_N]; // right hand side - contains currents
- float_ext_type m_lA[storage_N][m_pitch];
- float_ext_type m_lAinv[storage_N][m_pitch];
+ nl_ext_double m_lA[storage_N][m_pitch];
+ nl_ext_double m_lAinv[storage_N][m_pitch];
- //float_ext_type m_RHSx[storage_N];
+ //nl_ext_double m_RHSx[storage_N];
- const std::size_t m_dim;
- std::size_t m_cnt;
+ const std::size_t m_dim;
+ std::size_t m_cnt;
- };
+};
- // ----------------------------------------------------------------------------------------
- // matrix_solver_direct
- // ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+// matrix_solver_direct
+// ----------------------------------------------------------------------------------------
- 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);
+template <std::size_t m_N, std::size_t storage_N>
+matrix_solver_sm_t<m_N, storage_N>::~matrix_solver_sm_t()
+{
+}
+
+template <std::size_t m_N, std::size_t storage_N>
+void matrix_solver_sm_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets)
+{
+ matrix_solver_t::setup_base(nets);
+
+ netlist().save(*this, m_last_RHS, "m_last_RHS");
+
+ for (unsigned k = 0; k < N(); k++)
+ netlist().save(*this, RHS(k), plib::pfmt("RHS.{1}")(k));
+}
- /* 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();
- for (std::size_t i = 0; i < kN; i++)
+template <std::size_t m_N, std::size_t storage_N>
+void matrix_solver_sm_t<m_N, storage_N>::LE_invert()
+{
+ const std::size_t kN = N();
+
+ for (std::size_t i = 0; i < kN; i++)
+ {
+ for (std::size_t j = 0; j < kN; j++)
{
- 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;
+ W(i,j) = lA(i,j) = A(i,j);
+ Ainv(i,j) = 0.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 std::size_t e = m_terms[i]->m_nzrd.size();
+ Ainv(i,i) = 1.0;
+ }
+ /* down */
+ for (std::size_t i = 0; i < kN; i++)
+ {
+ /* FIXME: Singular matrix? */
+ const nl_double f = 1.0 / W(i,i);
+ const auto * RESTRICT const p = m_terms[i]->m_nzrd.data();
+ const std::size_t e = 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();
- for (std::size_t jb = 0; jb < eb; jb++)
+ const auto * RESTRICT const pb = m_terms[i]->m_nzbd.data();
+ const std::size_t eb = m_terms[i]->m_nzbd.size();
+ for (std::size_t jb = 0; jb < eb; jb++)
+ {
+ const unsigned j = pb[jb];
+ const nl_double f1 = - W(j,i) * f;
+ if (f1 != 0.0)
{
- const unsigned j = pb[jb];
- const float_type f1 = - W(j,i) * f;
- if (f1 != 0.0)
- {
- 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 < 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;
}
}
- /* up */
- for (std::size_t i = kN; i-- > 0; )
+ }
+ /* up */
+ for (std::size_t i = kN; i-- > 0; )
+ {
+ /* FIXME: Singular matrix? */
+ const nl_double f = 1.0 / W(i,i);
+ for (std::size_t j = i; j-- > 0; )
{
- /* FIXME: Singular matrix? */
- const float_type f = 1.0 / W(i,i);
- for (std::size_t j = i; j-- > 0; )
+ const nl_double f1 = - W(j,i) * f;
+ if (f1 != 0.0)
{
- const float_type f1 = - W(j,i) * f;
- if (f1 != 0.0)
- {
- 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;
- lAinv(i,k) = Ainv(i,k);
+ 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;
+ lAinv(i,k) = Ainv(i,k);
+ }
}
+}
- template <typename FT, int SIZE>
- template <typename T>
- void matrix_solver_sm_t<FT, SIZE>::LE_compute_x(
- T * x)
- {
- const std::size_t kN = size();
+template <std::size_t m_N, std::size_t storage_N>
+template <typename T>
+void matrix_solver_sm_t<m_N, storage_N>::LE_compute_x(
+ T * RESTRICT x)
+{
+ const std::size_t kN = N();
- for (std::size_t i=0; i<kN; i++)
- x[i] = 0.0;
+ for (std::size_t i=0; i<kN; i++)
+ x[i] = 0.0;
- for (std::size_t k=0; k<kN; k++)
- {
- const float_type f = RHS(k);
+ for (std::size_t k=0; k<kN; k++)
+ {
+ const nl_double f = RHS(k);
- for (std::size_t i=0; i<kN; i++)
- x[i] += Ainv(i,k) * f;
- }
+ for (std::size_t i=0; i<kN; i++)
+ x[i] += Ainv(i,k) * f;
}
+}
- template <typename FT, int SIZE>
- unsigned matrix_solver_sm_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphson)
- {
- static constexpr const bool incremental = true;
- const std::size_t iN = size();
- float_type new_V[storage_N]; // = { 0.0 };
+template <std::size_t m_N, std::size_t storage_N>
+unsigned matrix_solver_sm_t<m_N, storage_N>::solve_non_dynamic(const bool newton_raphson)
+{
+ static constexpr const bool incremental = true;
+ const std::size_t iN = N();
+
+ nl_double new_V[storage_N]; // = { 0.0 };
- if ((m_cnt % 50) == 0)
+ if ((m_cnt % 50) == 0)
+ {
+ /* complete calculation */
+ this->LE_invert();
+ }
+ else
+ {
+ if (!incremental)
{
- /* complete calculation */
- this->LE_invert();
+ for (std::size_t row = 0; row < iN; row ++)
+ for (std::size_t k = 0; k < iN; k++)
+ Ainv(row,k) = lAinv(row, k);
}
- else
+ for (std::size_t row = 0; row < iN; row ++)
{
- if (!incremental)
+ nl_double v[m_pitch] = {0};
+ std::size_t cols[m_pitch];
+ std::size_t colcount = 0;
+
+ auto &nz = m_terms[row]->m_nz;
+ for (unsigned & col : nz)
{
- for (std::size_t row = 0; row < iN; row ++)
- for (std::size_t k = 0; k < iN; k++)
- Ainv(row,k) = lAinv(row, k);
+ v[col] = A(row,col) - lA(row,col);
+ if (incremental)
+ lA(row,col) = A(row,col);
+ if (v[col] != 0.0)
+ cols[colcount++] = col;
}
- for (std::size_t row = 0; row < iN; row ++)
+
+ if (colcount > 0)
{
- float_type v[m_pitch] = {0};
- std::size_t cols[m_pitch];
- std::size_t colcount = 0;
+ nl_double lamba = 0.0;
+ nl_double w[m_pitch] = {0};
+
+ nl_double z[m_pitch];
+ /* compute w and lamba */
+ for (std::size_t i = 0; i < iN; i++)
+ z[i] = Ainv(i, row); /* u is row'th column */
- auto &nz = m_terms[row]->m_nz;
- for (unsigned & col : nz)
+ for (std::size_t j = 0; j < colcount; j++)
+ lamba += v[cols[j]] * z[cols[j]];
+
+ for (std::size_t j=0; j<colcount; j++)
{
- v[col] = A(row,col) - lA(row,col);
- if (incremental)
- lA(row,col) = A(row,col);
- if (v[col] != 0.0)
- cols[colcount++] = col;
+ std::size_t col = cols[j];
+ nl_double f = v[col];
+ for (std::size_t k = 0; k < iN; k++)
+ w[k] += Ainv(col,k) * f; /* Transpose(Ainv) * v */
}
- if (colcount > 0)
+ lamba = -1.0 / (1.0 + lamba);
+ for (std::size_t i=0; i<iN; i++)
{
- float_type lamba = 0.0;
- float_type w[m_pitch] = {0};
-
- float_type z[m_pitch];
- /* compute w and lamba */
- for (std::size_t i = 0; i < iN; i++)
- 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<colcount; j++)
- {
- std::size_t col = cols[j];
- float_type f = v[col];
+ const nl_double f = lamba * z[i];
+ if (f != 0.0)
for (std::size_t k = 0; k < iN; k++)
- w[k] += Ainv(col,k) * f; /* Transpose(Ainv) * v */
- }
-
- lamba = -1.0 / (1.0 + lamba);
- for (std::size_t i=0; i<iN; i++)
- {
- const float_type f = lamba * z[i];
- if (f != 0.0)
- for (std::size_t k = 0; k < iN; k++)
- Ainv(i,k) += f * w[k];
- }
+ Ainv(i,k) += f * w[k];
}
-
}
- }
-
- 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;
+ }
}
- template <typename FT, int SIZE>
- unsigned matrix_solver_sm_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson)
- {
- this->build_LE_A(*this);
- this->build_LE_RHS(*this);
+ m_cnt++;
- this->m_stat_calculations++;
- return this->solve_non_dynamic(newton_raphson);
- }
+ this->LE_compute_x(new_V);
+
+ const nl_double 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>
- 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)
+template <std::size_t m_N, std::size_t storage_N>
+inline unsigned matrix_solver_sm_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson)
+{
+ build_LE_A<matrix_solver_sm_t>();
+ build_LE_RHS<matrix_solver_sm_t>();
+
+ for (std::size_t i=0, iN=N(); i < iN; i++)
+ m_last_RHS[i] = RHS(i);
+
+ this->m_stat_calculations++;
+ return this->solve_non_dynamic(newton_raphson);
+}
+
+template <std::size_t m_N, std::size_t storage_N>
+matrix_solver_sm_t<m_N, storage_N>::matrix_solver_sm_t(netlist_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)
+{
+ for (std::size_t k = 0; k < N(); k++)
{
+ m_last_RHS[k] = 0.0;
}
+}
-} // namespace devices
+ } //namespace devices
} // namespace netlist
#endif /* NLD_MS_DIRECT_H_ */
diff --git a/src/lib/netlist/solver/nld_ms_sor.h b/src/lib/netlist/solver/nld_ms_sor.h
index c31aaa6d46a..eea692d6c47 100644
--- a/src/lib/netlist/solver/nld_ms_sor.h
+++ b/src/lib/netlist/solver/nld_ms_sor.h
@@ -12,42 +12,33 @@
#ifndef NLD_MS_SOR_H_
#define NLD_MS_SOR_H_
+#include <algorithm>
+
#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>
+ {
+template <std::size_t m_N, std::size_t storage_N>
+class matrix_solver_SOR_t: public matrix_solver_direct_t<m_N, storage_N>
{
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)
+ matrix_solver_SOR_t(netlist_t &anetlist, const pstring &name, const solver_parameters_t *params, const std::size_t size)
+ : matrix_solver_direct_t<m_N, storage_N>(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)
{
}
- void vsetup(analog_net_t::list_t &nets) override;
- unsigned vsolve_non_dynamic(const bool newton_raphson) override;
+ virtual ~matrix_solver_SOR_t() override {}
+
+ virtual void vsetup(analog_net_t::list_t &nets) override;
+ virtual unsigned vsolve_non_dynamic(const bool newton_raphson) override;
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;
+ state_var<nl_double> m_lp_fact;
};
// ----------------------------------------------------------------------------------------
@@ -55,16 +46,16 @@ private:
// ----------------------------------------------------------------------------------------
-template <typename FT, int SIZE>
-void matrix_solver_SOR_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets)
+template <std::size_t m_N, std::size_t storage_N>
+void matrix_solver_SOR_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets)
{
- matrix_solver_direct_t<FT, SIZE>::vsetup(nets);
+ matrix_solver_direct_t<m_N, storage_N>::vsetup(nets);
}
-template <typename FT, int SIZE>
-unsigned matrix_solver_SOR_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson)
+template <std::size_t m_N, std::size_t storage_N>
+unsigned matrix_solver_SOR_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson)
{
- const std::size_t iN = this->size();
+ const std::size_t iN = this->N();
bool resched = false;
unsigned resched_cnt = 0;
@@ -76,21 +67,26 @@ unsigned matrix_solver_SOR_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_rap
* omega = 2.0 / (1.0 + std::sqrt(1-rho))
*/
- const float_type ws = this->m_params.m_gs_sor;
+ const nl_double ws = this->m_params.m_gs_sor;
+
+ nl_double w[storage_N];
+ nl_double one_m_w[storage_N];
+ nl_double RHS[storage_N];
+ nl_double new_V[storage_N];
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;
+ nl_double gtot_t = 0.0;
+ nl_double gabs_t = 0.0;
+ nl_double RHS_t = 0.0;
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];
+ const nl_double * const RESTRICT gt = this->m_terms[k]->gt();
+ const nl_double * const RESTRICT go = this->m_terms[k]->go();
+ const nl_double * const RESTRICT Idr = this->m_terms[k]->Idr();
+ const nl_double * const *other_cur_analog = this->m_terms[k]->connected_net_V();
- this->m_new_V[k] = this->m_nets[k]->Q_Analog();
+ new_V[k] = this->m_nets[k]->Q_Analog();
for (std::size_t i = 0; i < term_count; i++)
{
@@ -99,60 +95,61 @@ unsigned matrix_solver_SOR_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_rap
}
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];
+ RHS_t = RHS_t + go[i] * *other_cur_analog[i];
RHS[k] = RHS_t;
- if (this->m_params.m_use_gabs)
+ if (USE_GABS)
{
for (std::size_t i = 0; i < term_count; i++)
gabs_t = gabs_t + std::abs(go[i]);
- gabs_t *= plib::constants<nl_double>::cast(0.5); // derived by try and error
+ gabs_t *= NL_FCONST(0.5); // derived by try and error
if (gabs_t <= gtot_t)
{
w[k] = ws / gtot_t;
- one_m_w[k] = plib::constants<FT>::one() - ws;
+ one_m_w[k] = NL_FCONST(1.0) - ws;
}
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] = NL_FCONST(1.0) / (gtot_t + gabs_t);
+ one_m_w[k] = NL_FCONST(1.0) - NL_FCONST(1.0) * gtot_t / (gtot_t + gabs_t);
}
}
else
{
w[k] = ws / gtot_t;
- one_m_w[k] = plib::constants<FT>::one() - ws;
+ one_m_w[k] = NL_FCONST(1.0) - ws;
}
}
- const float_type accuracy = this->m_params.m_accuracy;
+ const nl_double accuracy = this->m_params.m_accuracy;
do {
resched = false;
- float_type err = 0;
+ nl_double 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 int * RESTRICT net_other = this->m_terms[k]->connected_net_idx();
const std::size_t railstart = this->m_terms[k]->m_railstart;
- const float_type * go = this->m_gonn[k];
+ const nl_double * RESTRICT go = this->m_terms[k]->go();
- float_type Idrive = 0.0;
+ nl_double 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])];
+ Idrive = Idrive + go[i] * new_V[net_other[i]];
- const float_type new_val = this->m_new_V[k] * one_m_w[k] + (Idrive + RHS[k]) * w[k];
+ const nl_double new_val = new_V[k] * one_m_w[k] + (Idrive + 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(std::abs(new_val - new_V[k]), err);
+ new_V[k] = new_val;
}
if (err > accuracy)
resched = true;
resched_cnt++;
- } while (resched && (resched_cnt < this->m_params.m_gs_loops));
+ //} while (resched && (resched_cnt < this->m_params.m_gs_loops));
+ } while (resched && ((resched_cnt < this->m_params.m_gs_loops)));
this->m_iterative_total += resched_cnt;
this->m_stat_calculations++;
@@ -161,12 +158,13 @@ unsigned matrix_solver_SOR_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_rap
{
// Fallback to direct solver ...
this->m_iterative_fail++;
- return matrix_solver_direct_t<FT, SIZE>::vsolve_non_dynamic(newton_raphson);
+ return matrix_solver_direct_t<m_N, storage_N>::vsolve_non_dynamic(newton_raphson);
}
- 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;
+ for (std::size_t k = 0; k < iN; k++)
+ this->m_nets[k]->set_Q_Analog(new_V[k]);
+
+ return resched_cnt;
}
} //namespace devices
diff --git a/src/lib/netlist/solver/nld_ms_sor_mat.h b/src/lib/netlist/solver/nld_ms_sor_mat.h
index 83e4870cf28..50bcac1a52d 100644
--- a/src/lib/netlist/solver/nld_ms_sor_mat.h
+++ b/src/lib/netlist/solver/nld_ms_sor_mat.h
@@ -12,220 +12,212 @@
#ifndef NLD_MS_SOR_MAT_H_
#define NLD_MS_SOR_MAT_H_
-#include "nld_matrix_solver.h"
+#include <algorithm>
+
#include "nld_ms_direct.h"
+#include "nld_matrix_solver.h"
#include "nld_solver.h"
-#include <algorithm>
-
namespace netlist
{
-namespace devices
-{
-
- template <typename FT, int SIZE>
- class matrix_solver_SOR_mat_t: public matrix_solver_direct_t<FT, SIZE>
+ namespace devices
{
- friend class matrix_solver_t;
+template <std::size_t m_N, std::size_t storage_N>
+class matrix_solver_SOR_mat_t: public matrix_solver_direct_t<m_N, storage_N>
+{
+ friend class matrix_solver_t;
- public:
+public:
- using float_type = FT;
+ matrix_solver_SOR_mat_t(netlist_t &anetlist, const pstring &name, const solver_parameters_t *params, std::size_t size)
+ : matrix_solver_direct_t<m_N, storage_N>(anetlist, name, matrix_solver_t::DESCENDING, params, size)
+ , m_Vdelta(*this, "m_Vdelta", 0.0)
+ , m_omega(*this, "m_omega", params->m_gs_sor)
+ , m_lp_fact(*this, "m_lp_fact", 0)
+ , m_gs_fail(*this, "m_gs_fail", 0)
+ , m_gs_total(*this, "m_gs_total", 0)
+ {
+ }
- 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)
- {
- }
+ virtual ~matrix_solver_SOR_mat_t() override {}
- void vsetup(analog_net_t::list_t &nets) override;
+ virtual void vsetup(analog_net_t::list_t &nets) override;
- unsigned vsolve_non_dynamic(const bool newton_raphson) override;
+ virtual unsigned vsolve_non_dynamic(const bool newton_raphson) override;
- private:
- //state_var<float_type[storage_N]> m_Vdelta;
- state_var<std::vector<float_type>> m_Vdelta;
+private:
+ state_var<nl_double[storage_N]> m_Vdelta;
- state_var<float_type> m_omega;
- state_var<float_type> m_lp_fact;
+ state_var<nl_double> m_omega;
+ state_var<nl_double> m_lp_fact;
+ state_var<int> m_gs_fail;
+ state_var<int> m_gs_total;
+};
- };
+// ----------------------------------------------------------------------------------------
+// matrix_solver - Gauss - Seidel
+// ----------------------------------------------------------------------------------------
- // ----------------------------------------------------------------------------------------
- // matrix_solver - Gauss - Seidel
- // ----------------------------------------------------------------------------------------
+template <std::size_t m_N, std::size_t storage_N>
+void matrix_solver_SOR_mat_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets)
+{
+ matrix_solver_direct_t<m_N, storage_N>::vsetup(nets);
+}
- 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>
+nl_double matrix_solver_SOR_mat_t<m_N, storage_N>::vsolve()
+{
+ /*
+ * enable linear prediction on first newton pass
+ */
- #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;
- }
+ if (USE_LINEAR_PREDICTION)
+ for (unsigned k = 0; k < this->N(); 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->N(); k++)
+ {
+ this->m_last_V[k] = this->m_nets[k]->m_cur_Analog;
+ }
- this->solve_base(this);
+ this->solve_base(this);
- if (this->m_params->use_linear_prediction)
+ if (USE_LINEAR_PREDICTION)
+ {
+ nl_double sq = 0;
+ nl_double sqo = 0;
+ const nl_double rez_cts = 1.0 / this->current_timestep();
+ for (unsigned k = 0; k < this->N(); k++)
{
- 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);
+ const analog_net_t *n = this->m_nets[k];
+ const nl_double 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;
}
-
- return this->compute_next_timestep();
+ // 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), (nl_double) 2.0);
+ else
+ m_lp_fact = NL_FCONST(0.0);
}
- #endif
- template <typename FT, int SIZE>
- unsigned matrix_solver_SOR_mat_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson)
- {
- /* 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.
- */
+ return this->compute_next_timestep();
+}
+#endif
+
+template <std::size_t m_N, std::size_t storage_N>
+unsigned matrix_solver_SOR_mat_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson)
+{
+ /* 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);
+ nl_double new_v[storage_N] = { 0.0 };
+ const std::size_t iN = this->N();
- bool resched = false;
+ matrix_solver_t::build_LE_A<matrix_solver_SOR_mat_t>();
+ matrix_solver_t::build_LE_RHS<matrix_solver_SOR_mat_t>();
- unsigned resched_cnt = 0;
+ bool resched = false;
+ unsigned resched_cnt = 0;
- #if 0
- static int ws_cnt = 0;
- ws_cnt++;
- if (1 && ws_cnt % 200 == 0)
+
+#if 0
+ static int ws_cnt = 0;
+ ws_cnt++;
+ if (1 && ws_cnt % 200 == 0)
+ {
+ // update omega
+ nl_double lambdaN = 0;
+ nl_double lambda1 = 1e9;
+ for (int k = 0; k < iN; k++)
{
- // update omega
- float_type lambdaN = 0;
- float_type lambda1 = 1e9;
- for (int k = 0; k < iN; k++)
- {
- #if 0
- float_type akk = std::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 s = 0.0;
- for (int i=0; i<iN; i++)
- s = s + std::abs(this->m_A[k][i]);
- akk = s / akk - 1.0;
- if ( akk > lambdaN)
- lambdaN = akk;
- if (akk < lambda1)
- lambda1 = akk;
- #endif
- }
-
- //ws = 2.0 / (2.0 - lambdaN - lambda1);
- m_omega = 2.0 / (2.0 - lambda1);
- }
+ #if 0
+ nl_double akk = std::abs(this->m_A[k][k]);
+ if ( akk > lambdaN)
+ lambdaN = akk;
+ if (akk < lambda1)
+ lambda1 = akk;
+ #else
+ nl_double akk = std::abs(this->m_A[k][k]);
+ nl_double s = 0.0;
+ for (int i=0; i<iN; i++)
+ s = s + std::abs(this->m_A[k][i]);
+ akk = s / akk - 1.0;
+ if ( akk > lambdaN)
+ lambdaN = akk;
+ if (akk < lambda1)
+ lambda1 = akk;
#endif
+ }
+ //printf("lambda: %f %f\n", lambda, 2.0 / (1.0 + 2 * sqrt(lambda)) );
+
+ //ws = 2.0 / (2.0 - lambdaN - lambda1);
+ m_omega = 2.0 / (2.0 - lambda1);
+ //printf("%f %f %f\n", m_omega, lambda1, lambdaN);
+ }
+#endif
+
+ for (std::size_t k = 0; k < iN; k++)
+ new_v[k] = this->m_nets[k]->Q_Analog();
+
+ do {
+ resched = false;
+ nl_double cerr = 0.0;
for (std::size_t k = 0; k < iN; k++)
- this->m_new_V[k] = this->m_nets[k]->Q_Analog();
-
- do {
- resched = false;
- float_type cerr = 0.0;
-
- 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();
-
- for (std::size_t i = 0; i < e; i++)
- Idrive = Idrive + this->A(k,p[i]) * this->m_new_V[p[i]];
-
- FT w = m_omega / this->A(k,k);
- if (this->m_params.m_use_gabs)
- {
- FT gabs_t = 0.0;
- 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 *= plib::constants<FT>::one(); // derived by try and error
- if (gabs_t > this->A(k,k))
- {
- w = plib::constants<FT>::one() / (this->A(k,k) + gabs_t);
- }
- }
-
- const float_type delta = w * (this->RHS(k) - Idrive) ;
- cerr = std::max(cerr, std::abs(delta));
- this->m_new_V[k] += delta;
- }
-
- if (cerr > 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++;
- //this->netlist().warning("Falling back to direct solver .. Consider increasing RESCHED_LOOPS");
- return matrix_solver_direct_t<FT, SIZE>::solve_non_dynamic(newton_raphson);
+ nl_double Idrive = 0;
+
+ 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]) * new_v[p[i]];
+
+ const nl_double delta = m_omega * (this->RHS(k) - Idrive) / this->A(k,k);
+ cerr = std::max(cerr, std::abs(delta));
+ new_v[k] += delta;
+ }
+
+ if (cerr > 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;
+ this->m_gs_total += resched_cnt;
- 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;
+ if (resched)
+ {
+ this->m_iterative_fail++;
+ //this->netlist().warning("Falling back to direct solver .. Consider increasing RESCHED_LOOPS");
+ this->m_gs_fail++;
+ return matrix_solver_direct_t<m_N, storage_N>::solve_non_dynamic(newton_raphson);
+ }
+ else {
+ this->store(new_v);
+ return resched_cnt;
}
-} // namespace devices
+}
+
+ } //namespace devices
} // namespace netlist
#endif /* NLD_MS_GAUSS_SEIDEL_H_ */
diff --git a/src/lib/netlist/solver/nld_ms_w.h b/src/lib/netlist/solver/nld_ms_w.h
index 3372b50c7c5..0d4e7781eed 100644
--- a/src/lib/netlist/solver/nld_ms_w.h
+++ b/src/lib/netlist/solver/nld_ms_w.h
@@ -40,80 +40,80 @@
#ifndef NLD_MS_W_H_
#define NLD_MS_W_H_
-#include "nld_matrix_solver.h"
-#include "nld_solver.h"
-#include "plib/vector_ops.h"
-
#include <algorithm>
+#include "nld_solver.h"
+#include "nld_matrix_solver.h"
+#include "vector_base.h"
+
namespace netlist
{
namespace devices
{
+//#define nl_ext_double _float128 // slow, very slow
+//#define nl_ext_double long double // slightly slower
+#define nl_ext_double nl_double
-template <typename FT, int SIZE>
+template <std::size_t m_N, std::size_t storage_N>
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_t &anetlist, const pstring &name, const solver_parameters_t *params, const std::size_t size);
- matrix_solver_w_t(netlist_state_t &anetlist, const pstring &name, const solver_parameters_t *params, const std::size_t size);
+ virtual ~matrix_solver_w_t() override;
- void vsetup(analog_net_t::list_t &nets) override;
- void reset() override { matrix_solver_t::reset(); }
+ virtual void vsetup(analog_net_t::list_t &nets) override;
+ virtual void reset() override { matrix_solver_t::reset(); }
protected:
- unsigned vsolve_non_dynamic(const bool newton_raphson) override;
+ virtual 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; }
+ constexpr std::size_t N() const { return (m_N == 0) ? m_dim : m_N; }
void LE_invert();
template <typename T>
- void LE_compute_x(T * x);
+ void LE_compute_x(T * RESTRICT x);
template <typename T1, typename T2>
- float_ext_type &A(const T1 &r, const T2 &c) { return m_A[r][c]; }
+ inline nl_ext_double &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]; }
+ inline nl_ext_double &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]; }
+ inline nl_ext_double &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]; }
+ inline nl_ext_double &RHS(const T1 &r) { return m_RHS[r]; }
template <typename T1, typename T2>
- float_ext_type &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; }
+ inline nl_ext_double &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; }
+ nl_double m_last_RHS[storage_N]; // right hand side - contains currents
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];
- float_ext_type m_RHS[storage_N]; // right hand side - contains currents
+ nl_ext_double m_A[storage_N][m_pitch];
+ nl_ext_double m_Ainv[storage_N][m_pitch];
+ nl_ext_double m_W[storage_N][m_pitch];
+ nl_ext_double m_RHS[storage_N]; // right hand side - contains currents
- float_ext_type m_lA[storage_N][m_pitch];
+ nl_ext_double m_lA[storage_N][m_pitch];
/* temporary */
- float_type H[storage_N][m_pitch] ;
+ nl_double H[storage_N][m_pitch] ;
unsigned rows[storage_N];
unsigned cols[storage_N][m_pitch];
unsigned colcount[storage_N];
unsigned m_cnt;
- //float_ext_type m_RHSx[storage_N];
+ //nl_ext_double m_RHSx[storage_N];
const std::size_t m_dim;
@@ -123,22 +123,28 @@ private:
// matrix_solver_direct
// ----------------------------------------------------------------------------------------
-template <typename FT, int SIZE>
-void matrix_solver_w_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets)
+template <std::size_t m_N, std::size_t storage_N>
+matrix_solver_w_t<m_N, storage_N>::~matrix_solver_w_t()
+{
+}
+
+template <std::size_t m_N, std::size_t storage_N>
+void matrix_solver_w_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets)
{
matrix_solver_t::setup_base(nets);
- // 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));
+ netlist().save(*this, m_last_RHS, "m_last_RHS");
+
+ for (unsigned k = 0; k < N(); k++)
+ netlist().save(*this, RHS(k), plib::pfmt("RHS.{1}")(k));
}
-template <typename FT, int SIZE>
-void matrix_solver_w_t<FT, SIZE>::LE_invert()
+template <std::size_t m_N, std::size_t storage_N>
+void matrix_solver_w_t<m_N, storage_N>::LE_invert()
{
- const std::size_t kN = size();
+ const std::size_t kN = N();
for (std::size_t i = 0; i < kN; i++)
{
@@ -153,18 +159,18 @@ void matrix_solver_w_t<FT, SIZE>::LE_invert()
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 nl_double f = 1.0 / W(i,i);
+ const auto * RESTRICT 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 auto * const pb = m_terms[i]->m_nzbd.data();
+ const auto * RESTRICT 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++)
{
const auto j = pb[jb];
- const float_type f1 = - W(j,i) * f;
+ const nl_double f1 = - W(j,i) * f;
if (f1 != 0.0)
{
for (std::size_t k = 0; k < e; k++)
@@ -178,10 +184,10 @@ void matrix_solver_w_t<FT, SIZE>::LE_invert()
for (std::size_t i = kN; i-- > 0; )
{
/* FIXME: Singular matrix? */
- const float_type f = 1.0 / W(i,i);
+ const nl_double f = 1.0 / W(i,i);
for (std::size_t j = i; j-- > 0; )
{
- const float_type f1 = - W(j,i) * f;
+ const nl_double f1 = - W(j,i) * f;
if (f1 != 0.0)
{
for (std::size_t k = i; k < kN; k++)
@@ -197,19 +203,19 @@ void matrix_solver_w_t<FT, SIZE>::LE_invert()
}
}
-template <typename FT, int SIZE>
+template <std::size_t m_N, std::size_t storage_N>
template <typename T>
-void matrix_solver_w_t<FT, SIZE>::LE_compute_x(
- T * x)
+void matrix_solver_w_t<m_N, storage_N>::LE_compute_x(
+ T * RESTRICT x)
{
- const std::size_t kN = size();
+ const std::size_t kN = N();
for (std::size_t i=0; i<kN; i++)
x[i] = 0.0;
for (std::size_t k=0; k<kN; k++)
{
- const float_type f = RHS(k);
+ const nl_double f = RHS(k);
for (std::size_t i=0; i<kN; i++)
x[i] += Ainv(i,k) * f;
@@ -217,14 +223,14 @@ void matrix_solver_w_t<FT, SIZE>::LE_compute_x(
}
-template <typename FT, int SIZE>
-unsigned matrix_solver_w_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphson)
+template <std::size_t m_N, std::size_t storage_N>
+unsigned matrix_solver_w_t<m_N, storage_N>::solve_non_dynamic(const bool newton_raphson)
{
- const auto iN = size();
+ const auto iN = N();
- float_type new_V[storage_N]; // = { 0.0 };
+ nl_double new_V[storage_N]; // = { 0.0 };
- if ((m_cnt % 50) == 0)
+ if ((m_cnt % 100) == 0)
{
/* complete calculation */
this->LE_invert();
@@ -260,7 +266,7 @@ unsigned matrix_solver_w_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphso
/* construct w = transform(V) * y
* dim: rowcount x iN
* */
- float_type w[storage_N];
+ nl_double w[storage_N];
for (unsigned i = 0; i < rowcount; i++)
{
const unsigned r = rows[i];
@@ -281,7 +287,7 @@ unsigned matrix_solver_w_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphso
for (unsigned k=0; k< colcount[i]; k++)
{
const unsigned col = cols[i][k];
- float_type f = VT(rows[i],col);
+ nl_double f = VT(rows[i],col);
if (f!=0.0)
for (unsigned j= 0; j < rowcount; j++)
H[i][j] += f * Ainv(col,rows[j]);
@@ -291,16 +297,16 @@ unsigned matrix_solver_w_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphso
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];
+ printf("%s H singular\n", this->name().c_str());
+ const nl_double f = 1.0 / H[i][i];
for (unsigned j = i+1; j < rowcount; j++)
{
- const float_type f1 = - f * H[j][i];
+ const nl_double f1 = - f * H[j][i];
if (f1!=0.0)
{
- float_type *pj = &H[j][i+1];
- const float_type *pi = &H[i][i+1];
+ nl_double *pj = &H[j][i+1];
+ const nl_double *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];
@@ -310,12 +316,12 @@ unsigned matrix_solver_w_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphso
}
/* Back substitution */
//inv(H) w = t w = H t
- float_type t[storage_N]; // FIXME: convert to member
+ nl_double t[storage_N]; // 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];
+ nl_double tmp = 0;
+ const nl_double *pj = &H[j][j+1];
+ const nl_double *tj = &t[j+1];
for (unsigned k = 0; k < rowcount-j-1; k++)
tmp += pj[k] * tj[k];
//tmp += H[j][k] * t[k];
@@ -325,7 +331,7 @@ unsigned matrix_solver_w_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphso
/* x = y - Zt */
for (unsigned i=0; i<iN; i++)
{
- float_type tmp = 0.0;
+ nl_double tmp = 0.0;
for (unsigned j=0; j<rowcount;j++)
{
const unsigned row = rows[j];
@@ -337,40 +343,47 @@ unsigned matrix_solver_w_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphso
}
m_cnt++;
- if (false)
+ if (0)
for (unsigned i=0; i<iN; i++)
{
- float_type tmp = 0.0;
+ nl_double 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));
+ printf("%s failed on row %d: %f RHS: %f\n", this->name().c_str(), i, std::abs(tmp-RHS(i)), RHS(i));
}
- const float_type err = (newton_raphson ? delta(new_V) : 0.0);
+ const nl_double 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>
-unsigned matrix_solver_w_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson)
+template <std::size_t m_N, std::size_t storage_N>
+inline unsigned matrix_solver_w_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson)
{
- this->build_LE_A(*this);
- this->build_LE_RHS(*this);
+ build_LE_A<matrix_solver_w_t>();
+ build_LE_RHS<matrix_solver_w_t>();
+
+ for (std::size_t i=0, iN=N(); i < iN; i++)
+ m_last_RHS[i] = RHS(i);
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,
+template <std::size_t m_N, std::size_t storage_N>
+matrix_solver_w_t<m_N, storage_N>::matrix_solver_w_t(netlist_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)
+: matrix_solver_t(anetlist, name, NOSORT, params)
+ ,m_cnt(0)
, m_dim(size)
{
+ for (std::size_t k = 0; k < N(); k++)
+ {
+ m_last_RHS[k] = 0.0;
+ }
}
} //namespace devices
diff --git a/src/lib/netlist/solver/nld_solver.cpp b/src/lib/netlist/solver/nld_solver.cpp
index 4734cc3624d..cae76caab1a 100644
--- a/src/lib/netlist/solver/nld_solver.cpp
+++ b/src/lib/netlist/solver/nld_solver.cpp
@@ -31,398 +31,393 @@
#pragma GCC optimize "ivopts"
#endif
-#include "netlist/nl_lists.h"
-#include "netlist/nl_factory.h"
+#include <algorithm>
+#include <cmath> // <<= needed by windows build
+
+#include "../nl_lists.h"
+
+#include "../plib/pomp.h"
+
+#include "../nl_factory.h"
+
+#include "nld_solver.h"
#include "nld_matrix_solver.h"
+
+#if 1
#include "nld_ms_direct.h"
-#include "nld_ms_direct1.h"
-#include "nld_ms_direct2.h"
#include "nld_ms_gcr.h"
-#include "nld_ms_gmres.h"
+#else
+#include "nld_ms_direct_lu.h"
+#endif
+#include "nld_ms_w.h"
#include "nld_ms_sm.h"
+#include "nld_ms_direct1.h"
+#include "nld_ms_direct2.h"
#include "nld_ms_sor.h"
#include "nld_ms_sor_mat.h"
-#include "nld_ms_w.h"
-#include "nld_solver.h"
-#include "plib/pomp.h"
-
-#include <algorithm>
-#include <cmath>
+#include "nld_ms_gmres.h"
namespace netlist
{
-namespace devices
-{
+ namespace devices
+ {
- // ----------------------------------------------------------------------------------------
- // solver
- // ----------------------------------------------------------------------------------------
- NETLIB_RESET(solver)
- {
- for (auto &s : m_mat_solvers)
- s->reset();
- }
- void NETLIB_NAME(solver)::stop()
- {
- for (auto &s : m_mat_solvers)
- s->log_stats();
- }
+// ----------------------------------------------------------------------------------------
+// solver
+// ----------------------------------------------------------------------------------------
- NETLIB_UPDATE(solver)
- {
- if (m_params.m_dynamic_ts)
- return;
+NETLIB_RESET(solver)
+{
+ for (std::size_t i = 0; i < m_mat_solvers.size(); i++)
+ m_mat_solvers[i]->do_reset();
+}
- 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));
+void NETLIB_NAME(solver)::stop()
+{
+ for (std::size_t i = 0; i < m_mat_solvers.size(); i++)
+ m_mat_solvers[i]->log_stats();
+}
- std::size_t nthreads = std::min(static_cast<std::size_t>(m_parallel()), plib::omp::get_max_threads());
+NETLIB_NAME(solver)::~NETLIB_NAME(solver)()
+{
+}
- std::vector<matrix_solver_t *> &solvers = (force_solve ? m_mat_solvers_all : m_mat_solvers_timestepping);
+NETLIB_UPDATE(solver)
+{
+ if (m_params.m_dynamic_ts)
+ return;
- 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);
- });
- }
- else
- for (auto & solver : solvers)
- {
- const netlist_time ts = solver->solve(now);
- plib::unused_var(ts);
- }
+ /* force solving during start up if there are no time-step devices */
+ /* FIXME: Needs a more elegant solution */
+ bool force_solve = (netlist().time() < netlist_time::from_double(2 * m_params.m_max_timestep));
- for (auto & solver : solvers)
- solver->update_inputs();
+ std::size_t nthreads = std::min(m_parallel(), plib::omp::get_max_threads());
+ std::size_t t_cnt = 0;
+ int solv[128];
+ for (int i = 0; i < m_mat_solvers.size(); i++)
+ if (m_mat_solvers[i]->has_timestep_devices() || force_solve)
+ solv[t_cnt++] = i;
- /* 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));
- }
+ if (nthreads > 1 && t_cnt > 1)
+ {
+ plib::omp::set_num_threads(nthreads);
+ plib::omp::for_static(0, t_cnt, [this, &solv](int i) { ATTR_UNUSED const netlist_time ts = this->m_mat_solvers[solv[i]]->solve(); });
}
+ else
+ for (auto & solver : m_mat_solvers)
+ if (solver->has_timestep_devices() || force_solve)
+ ATTR_UNUSED const netlist_time ts = solver->solve();
+
+ for (auto & solver : m_mat_solvers)
+ if (solver->has_timestep_devices() || force_solve)
+ solver->update_inputs();
- 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)
+ /* step circuit */
+ if (!m_Q_step.net().is_queued())
{
- return pool().make_poolptr<C>(nl, name, &params, size);
+ m_Q_step.net().toggle_and_push_to_queue(netlist_time::from_double(m_params.m_max_timestep));
}
+}
- template <typename FT, int SIZE>
- pool_owned_ptr<matrix_solver_t> NETLIB_NAME(solver)::create_solver(std::size_t size, const pstring &solvername)
+template <class C>
+std::unique_ptr<matrix_solver_t> create_it(netlist_t &nl, pstring name, solver_parameters_t &params, std::size_t size)
+{
+ typedef C solver;
+ return plib::make_unique<solver>(nl, name, &params, size);
+}
+
+template <std::size_t m_N, std::size_t storage_N>
+std::unique_ptr<matrix_solver_t> NETLIB_NAME(solver)::create_solver(std::size_t size, const pstring &solvername)
+{
+ if (pstring("SOR_MAT").equals(m_method()))
{
- 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
- {
- log().fatal(MF_1_UNKNOWN_SOLVER_TYPE, m_method());
- return pool_owned_ptr<matrix_solver_t>();
- }
+ return create_it<matrix_solver_SOR_mat_t<m_N, storage_N>>(netlist(), solvername, m_params, size);
+ //typedef matrix_solver_SOR_mat_t<m_N,storage_N> solver_sor_mat;
+ //return plib::make_unique<solver_sor_mat>(netlist(), solvername, &m_params, size);
}
-
- template <typename FT, int SIZE>
- pool_owned_ptr<matrix_solver_t> NETLIB_NAME(solver)::create_solver_x(std::size_t size, const pstring &solvername)
+ else if (pstring("MAT_CR").equals(m_method()))
{
- if (SIZE > 0)
+ if (size > 0) // GCR always outperforms MAT solver
{
- if (size == SIZE)
- return create_solver<FT, SIZE>(size, solvername);
- else
- return this->create_solver_x<FT, SIZE-1>(size, solvername);
+ typedef matrix_solver_GCR_t<m_N,storage_N> solver_mat;
+ return plib::make_unique<solver_mat>(netlist(), solvername, &m_params, size);
}
else
{
- if (size * 2 > -SIZE )
- return create_solver<FT, SIZE>(size, solvername);
- else
- return this->create_solver_x<FT, SIZE / 2>(size, solvername);
+ typedef matrix_solver_direct_t<m_N,storage_N> solver_mat;
+ return plib::make_unique<solver_mat>(netlist(), solvername, &m_params, size);
}
}
-
- struct net_splitter
+ else if (pstring("MAT").equals(m_method()))
+ {
+ typedef matrix_solver_direct_t<m_N,storage_N> solver_mat;
+ return plib::make_unique<solver_mat>(netlist(), solvername, &m_params, size);
+ }
+ else if (pstring("SM").equals(m_method()))
+ {
+ /* Sherman-Morrison Formula */
+ typedef matrix_solver_sm_t<m_N,storage_N> solver_mat;
+ return plib::make_unique<solver_mat>(netlist(), solvername, &m_params, size);
+ }
+ else if (pstring("W").equals(m_method()))
+ {
+ /* Woodbury Formula */
+ typedef matrix_solver_w_t<m_N,storage_N> solver_mat;
+ return plib::make_unique<solver_mat>(netlist(), solvername, &m_params, size);
+ }
+ else if (pstring("SOR").equals(m_method()))
+ {
+ typedef matrix_solver_SOR_t<m_N,storage_N> solver_GS;
+ return plib::make_unique<solver_GS>(netlist(), solvername, &m_params, size);
+ }
+ else if (pstring("GMRES").equals(m_method()))
{
+ typedef matrix_solver_GMRES_t<m_N,storage_N> solver_GMRES;
+ return plib::make_unique<solver_GMRES>(netlist(), solvername, &m_params, size);
+ }
+ else
+ {
+ log().fatal(MF_1_UNKNOWN_SOLVER_TYPE, m_method());
+ return nullptr;
+ }
+}
- bool already_processed(const analog_net_t &n) const
- {
- /* no need to process rail nets - these are known variables */
- if (n.isRailNet())
+struct net_splitter
+{
+
+ bool already_processed(analog_net_t *n)
+ {
+ if (n->isRailNet())
+ return true;
+ for (auto & grp : groups)
+ if (plib::container::contains(grp, n))
return true;
- /* if it's already processed - no need to continue */
- for (auto & grp : groups)
- if (plib::container::contains(grp, &n))
- return true;
- return false;
- }
+ return false;
+ }
- void process_net(analog_net_t &n)
+ void process_net(analog_net_t *n)
+ {
+ if (n->num_cons() == 0)
+ return;
+ /* add the net */
+ groups.back().push_back(n);
+ for (auto &p : n->m_core_terms)
{
- /* 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())
+ if (p->is_type(detail::terminal_type::TERMINAL))
{
- /* 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);
- }
+ terminal_t *pt = static_cast<terminal_t *>(p);
+ analog_net_t *other_net = &pt->m_otherterm->net();
+ if (!already_processed(other_net))
+ process_net(other_net);
}
}
+ }
- void run(netlist_state_t &netlist)
+ void run(netlist_t &netlist)
+ {
+ for (auto & net : netlist.m_nets)
{
- for (auto & net : netlist.nets())
+ netlist.log().debug("processing {1}\n", net->name());
+ if (!net->isRailNet() && net->num_cons() > 0)
{
- netlist.log().debug("processing {1}\n", net->name());
- if (!net->isRailNet() && net->num_cons() > 0)
+ netlist.log().debug(" ==> not a rail net\n");
+ /* Must be an analog net */
+ analog_net_t *n = static_cast<analog_net_t *>(net.get());
+ if (!already_processed(n))
{
- 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))
- {
- groups.emplace_back(analog_net_t::list_t());
- process_net(n);
- }
+ groups.push_back(analog_net_t::list_t());
+ process_net(n);
}
}
}
+ }
- std::vector<analog_net_t::list_t> groups;
- };
+ std::vector<analog_net_t::list_t> groups;
+};
- void NETLIB_NAME(solver)::post_start()
+void NETLIB_NAME(solver)::post_start()
+{
+ const bool use_specific = true;
+
+ 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();
+
+ if (m_params.m_dynamic_ts)
{
- 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();
-
+ m_params.m_max_timestep *= 1;//NL_FCONST(1000.0);
+ }
+ else
+ {
+ m_params.m_min_timestep = m_params.m_max_timestep;
+ }
- 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::)
- //m_params.m_max_timestep = std::max(m_params.m_max_timestep, m_params.m_max_timestep::)
+ // Override log statistics
+ pstring p = plib::util::environment("NL_STATS", "");
+ if (p != "")
+ m_params.m_log_stats = p.as_long();
+ else
+ m_params.m_log_stats = m_log_stats();
- // Override log statistics
- pstring p = plib::util::environment("NL_STATS", "");
- if (p != "")
- m_params.m_log_stats = plib::pstonum<decltype(m_params.m_log_stats)>(p);
- else
- m_params.m_log_stats = m_log_stats();
+ log().verbose("Scanning net groups ...");
+ // determine net groups
- log().verbose("Scanning net groups ...");
- // determine net groups
+ net_splitter splitter;
- net_splitter splitter;
+ splitter.run(netlist());
- splitter.run(state());
+ // setup the solvers
+ log().verbose("Found {1} net groups in {2} nets\n", splitter.groups.size(), netlist().m_nets.size());
+ for (auto & grp : splitter.groups)
+ {
+ std::unique_ptr<matrix_solver_t> ms;
+ std::size_t net_count = grp.size();
+ pstring sname = plib::pfmt("Solver_{1}")(m_mat_solvers.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)
+ switch (net_count)
{
- pool_owned_ptr<matrix_solver_t> ms;
- std::size_t net_count = grp.size();
- pstring sname = plib::pfmt("Solver_{1}")(m_mat_solvers.size());
+#if 1
+ case 1:
+ if (use_specific)
+ ms = plib::make_unique<matrix_solver_direct1_t>(netlist(), sname, &m_params);
+ else
+ ms = create_solver<1,1>(1, sname);
+ break;
+ case 2:
+ if (use_specific)
+ ms = plib::make_unique<matrix_solver_direct2_t>(netlist(), sname, &m_params);
+ else
+ ms = create_solver<2,2>(2, sname);
+ break;
+#if 0
+ case 3:
+ ms = create_solver<3,3>(3, sname);
+ break;
+ case 4:
+ ms = create_solver<4,4>(4, sname);
+ break;
+ case 5:
+ ms = create_solver<5,5>(5, sname);
+ break;
+ case 6:
+ ms = create_solver<6,6>(6, sname);
+ break;
+ case 7:
+ ms = create_solver<7,7>(7, sname);
+ break;
+ case 8:
+ ms = create_solver<8,8>(8, sname);
+ break;
+ case 9:
+ ms = create_solver<9,9>(9, sname);
+ break;
+ case 10:
+ ms = create_solver<10,10>(10, sname);
+ break;
+ case 11:
+ ms = create_solver<11,11>(11, sname);
+ break;
+ case 12:
+ ms = create_solver<12,12>(12, sname);
+ break;
+ case 15:
+ ms = create_solver<15,15>(15, sname);
+ break;
+ case 31:
+ ms = create_solver<31,31>(31, sname);
+ break;
+ case 35:
+ ms = create_solver<35,35>(35, sname);
+ break;
+ case 43:
+ ms = create_solver<43,43>(43, sname);
+ break;
+ case 49:
+ ms = create_solver<49,49>(49, sname);
+ break;
+#endif
+#if 0
+ case 87:
+ ms = create_solver<87,87>(87, sname);
+ break;
+#endif
+#endif
+ default:
+ log().warning(MW_1_NO_SPECIFIC_SOLVER, net_count);
+ if (net_count <= 8)
+ {
+ ms = create_solver<0, 8>(net_count, sname);
+ }
+ else if (net_count <= 16)
+ {
+ ms = create_solver<0,16>(net_count, sname);
+ }
+ else if (net_count <= 32)
+ {
+ ms = create_solver<0,32>(net_count, sname);
+ }
+ else
+ if (net_count <= 64)
+ {
+ ms = create_solver<0,64>(net_count, sname);
+ }
+ else
+ if (net_count <= 128)
+ {
+ ms = create_solver<0,128>(net_count, sname);
+ }
+ else
+ {
+ log().fatal(MF_1_NETGROUP_SIZE_EXCEEDED_1, 128);
+ ms = nullptr; /* tease compilers */
+ }
- switch (net_count)
- {
- #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().warning(MW_1_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)
- {
- ms = create_solver<double, -128>(net_count, sname);
- }
- else
- {
- log().fatal(MF_1_NETGROUP_SIZE_EXCEEDED_1, 128);
- return; /* tease compilers */
- }
- break;
- }
+ break;
+ }
- // FIXME ...
- ms->setup(grp);
+ // FIXME ...
+ ms->setup(grp);
- log().verbose("Solver {1}", ms->name());
- log().verbose(" ==> {2} 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");
- for (auto &n : grp)
+ log().verbose("Solver {1}", ms->name());
+ log().verbose(" ==> {2} 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");
+ for (auto &n : grp)
+ {
+ log().verbose("Net {1}", n->name());
+ for (const auto &pcore : n->m_core_terms)
{
- log().verbose("Net {1}", n->name());
- for (const auto &pcore : n->core_terms())
- {
- log().verbose(" {1}", pcore->name());
- }
+ log().verbose(" {1}", pcore->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_solvers.push_back(std::move(ms));
}
+}
- void NETLIB_NAME(solver)::create_solver_code(std::map<pstring, pstring> &mp)
+void NETLIB_NAME(solver)::create_solver_code(std::map<pstring, pstring> &mp)
+{
+ for (auto & s : m_mat_solvers)
{
- 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();
+ mp[r.first] = r.second; // automatically overwrites identical names
}
+}
- NETLIB_DEVICE_IMPL(solver, "SOLVER", "FREQ")
+ NETLIB_DEVICE_IMPL(solver)
-} // namespace devices
+ } //namespace devices
} // namespace netlist
diff --git a/src/lib/netlist/solver/nld_solver.h b/src/lib/netlist/solver/nld_solver.h
index c9ec967a72a..986d14f401a 100644
--- a/src/lib/netlist/solver/nld_solver.h
+++ b/src/lib/netlist/solver/nld_solver.h
@@ -8,13 +8,11 @@
#ifndef NLD_SOLVER_H_
#define NLD_SOLVER_H_
-#include "netlist/nl_base.h"
-#include "nld_matrix_solver.h"
-#include "plib/pstream.h"
-
#include <map>
-#include <memory>
-#include <vector>
+
+#include "../nl_base.h"
+#include "../plib/pstream.h"
+#include "nld_matrix_solver.h"
//#define ATTR_ALIGNED(N) __attribute__((aligned(N)))
#define ATTR_ALIGNED(N) ATTR_ALIGN
@@ -25,97 +23,89 @@
namespace netlist
{
-namespace devices
-{
- class NETLIB_NAME(solver);
+ namespace devices
+ {
+class NETLIB_NAME(solver);
- class matrix_solver_t;
- NETLIB_OBJECT(solver)
+class matrix_solver_t;
+
+NETLIB_OBJECT(solver)
+{
+ 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", NETLIST_GMIN_DEFAULT)
+ , m_pivot(*this, "PIVOT", 0) // use pivoting - on supported solvers
+ , m_nr_loops(*this, "NR_LOOPS", 250) // Newton-Raphson loops
+ , m_nr_recalc_delay(*this, "NR_RECALC_DELAY", NLTIME_FROM_NS(10).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", 0)
+ , m_dynamic_lte(*this, "DYNAMIC_LTE", 1e-5) // diff/timestep
+ , m_dynamic_min_ts(*this, "DYNAMIC_MIN_TIMESTEP", 1e-6) // nl_double timestep resolution
+
+ , m_log_stats(*this, "LOG_STATS", 0) // log statistics on shutdown
+ , m_params()
{
- 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", NLTIME_FROM_NS(10).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_log_stats(*this, "LOG_STATS", true) // log statistics on shutdown
- , m_params()
- {
- // internal staff
-
- connect(m_fb_step, m_Q_step);
- }
-
- void post_start();
- void stop();
-
- nl_double gmin() const { return m_gmin(); }
-
- void create_solver_code(std::map<pstring, pstring> &mp);
-
- NETLIB_UPDATEI();
- NETLIB_RESETI();
- // NETLIB_UPDATE_PARAMI();
-
- private:
- 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;
-
- param_logic_t m_log_stats;
-
- 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;
-
- template <typename FT, int SIZE>
- pool_owned_ptr<matrix_solver_t> create_solver(std::size_t size, const pstring &solvername);
-
- template <typename FT, int SIZE>
- pool_owned_ptr<matrix_solver_t> create_solver_x(std::size_t size, const pstring &solvername);
- };
-
-} //namespace devices
+ // internal staff
+
+ connect(m_fb_step, m_Q_step);
+ }
+
+ virtual ~NETLIB_NAME(solver)() override;
+
+ void post_start();
+ void stop();
+
+ inline nl_double gmin() { return m_gmin(); }
+
+ void create_solver_code(std::map<pstring, pstring> &mp);
+
+ NETLIB_UPDATEI();
+ NETLIB_RESETI();
+ // NETLIB_UPDATE_PARAMI();
+
+protected:
+ 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_log_stats;
+
+ std::vector<std::unique_ptr<matrix_solver_t>> m_mat_solvers;
+private:
+
+ solver_parameters_t m_params;
+
+ template <std::size_t m_N, std::size_t storage_N>
+ std::unique_ptr<matrix_solver_t> create_solver(std::size_t size, const pstring &solvername);
+};
+
+ } //namespace devices
} // namespace netlist
#endif /* NLD_SOLVER_H_ */
diff --git a/src/lib/netlist/solver/vector_base.h b/src/lib/netlist/solver/vector_base.h
new file mode 100644
index 00000000000..28f5fa0dd44
--- /dev/null
+++ b/src/lib/netlist/solver/vector_base.h
@@ -0,0 +1,144 @@
+// license:GPL-2.0+
+// copyright-holders:Couriersud
+/*
+ * vector_base.h
+ *
+ * Base vector operations
+ *
+ */
+
+#ifndef VECTOR_BASE_H_
+#define VECTOR_BASE_H_
+
+#include <algorithm>
+#include "../plib/pconfig.h"
+
+#if 0
+template <unsigned storage_N>
+struct pvector
+{
+ pvector(unsigned size)
+ : m_N(size) { }
+
+ unsigned size() {
+ if (storage_N)
+ }
+
+ double m_V[storage_N];
+private:
+ unsigned m_N;
+};
+#endif
+
+#if !defined(__clang__) && !defined(_MSC_VER) && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ > 6))
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
+#endif
+
+template<typename T, std::size_t N>
+inline static void vec_set (const std::size_t n, const T scalar, T (& RESTRICT v)[N])
+{
+ if (n != N)
+ for ( std::size_t i = 0; i < n; i++ )
+ v[i] = scalar;
+ else
+ for ( std::size_t i = 0; i < N; i++ )
+ v[i] = scalar;
+}
+
+template<typename T, std::size_t N>
+inline static T vec_mult (const std::size_t n, const T (& RESTRICT v1)[N], const T (& RESTRICT v2)[N] )
+{
+ T value = 0.0;
+ if (n != N)
+ for ( std::size_t i = 0; i < n; i++ )
+ value += v1[i] * v2[i];
+ else
+ for ( std::size_t i = 0; i < N; i++ )
+ value += v1[i] * v2[i];
+ return value;
+}
+
+template<typename T, std::size_t N>
+inline static T vec_mult2 (const std::size_t n, const T (& RESTRICT v)[N])
+{
+ T value = 0.0;
+ if (n != N)
+ for ( std::size_t i = 0; i < n; i++ )
+ value += v[i] * v[i];
+ else
+ for ( std::size_t i = 0; i < N; i++ )
+ value += v[i] * v[i];
+ return value;
+}
+
+template<typename T, std::size_t N>
+inline static void vec_mult_scalar (const std::size_t n, const T (& RESTRICT v)[N], const T & scalar, T (& RESTRICT result)[N])
+{
+ if (n != N)
+ for ( std::size_t i = 0; i < n; i++ )
+ result[i] = scalar * v[i];
+ else
+ for ( std::size_t i = 0; i < N; i++ )
+ result[i] = scalar * v[i];
+}
+
+template<typename T, std::size_t N>
+inline static void vec_add_mult_scalar (const std::size_t n, const T (& RESTRICT v)[N], const T scalar, T (& RESTRICT result)[N])
+{
+ if (n != N)
+ for ( std::size_t i = 0; i < n; i++ )
+ result[i] = result[i] + scalar * v[i];
+ else
+ for ( std::size_t i = 0; i < N; i++ )
+ result[i] = result[i] + scalar * v[i];
+}
+
+template<typename T>
+inline static void vec_add_mult_scalar_p(const std::size_t & n, const T * RESTRICT v, const T scalar, T * RESTRICT result)
+{
+ for ( std::size_t i = 0; i < n; i++ )
+ result[i] += scalar * v[i];
+}
+
+inline static void vec_add_ip(const std::size_t n, const double * RESTRICT v, double * RESTRICT result)
+{
+ for ( std::size_t i = 0; i < n; i++ )
+ result[i] += v[i];
+}
+
+template<typename T, std::size_t N>
+inline void vec_sub(const std::size_t n, const T (& RESTRICT v1)[N], const T (& RESTRICT v2)[N], T (& RESTRICT result)[N])
+{
+ if (n != N)
+ for ( std::size_t i = 0; i < n; i++ )
+ result[i] = v1[i] - v2[i];
+ else
+ for ( std::size_t i = 0; i < N; i++ )
+ result[i] = v1[i] - v2[i];
+}
+
+template<typename T, std::size_t N>
+inline void vec_scale(const std::size_t n, T (& RESTRICT v)[N], const T scalar)
+{
+ if (n != N)
+ for ( std::size_t i = 0; i < n; i++ )
+ v[i] = scalar * v[i];
+ else
+ for ( std::size_t i = 0; i < N; i++ )
+ v[i] = scalar * v[i];
+}
+
+inline double vec_maxabs(const std::size_t n, const double * RESTRICT v)
+{
+ double ret = 0.0;
+ for ( std::size_t i = 0; i < n; i++ )
+ ret = std::max(ret, std::abs(v[i]));
+
+ return ret;
+}
+#if !defined(__clang__) && !defined(_MSC_VER) && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ > 6))
+#pragma GCC diagnostic pop
+#endif
+
+#endif /* MAT_CR_H_ */