summaryrefslogtreecommitdiffstatshomepage
path: root/src/lib/netlist/solver
diff options
context:
space:
mode:
author couriersud <couriersud@gmx.org>2020-06-13 15:49:35 +0200
committer couriersud <couriersud@gmx.org>2020-06-13 15:49:35 +0200
commite949e9c29de82ee7c32692e7b65da05dd22bdc9d (patch)
tree8770696a48827e3b63d94063f8162eaadeff6539 /src/lib/netlist/solver
parent75681d760c478f772fdb1603222498f96a9b61e7 (diff)
netlist: Performance improvement and refactoring. [Couriersud]
Kidniki now achieves up to 910% when run with static solvers and with nltool. That is significant better than the 860% we have seen previously. This increase is driven by using a global memory pool in the solver code. In addition the following refactoring and code maintenance work is included. Please excuse the large commit, some of this took interfered with other work and the detail development steps were ugly. - gsl support: This commit adds pgsl.h which implements a very limited number of the functionality of the gsl header described in the c++ core guidelines. - clang-tidy fixes - A significant refactoring of palloc.h. Aligned hints were removed, they added complexity without a significant performance gain. Vector operations should better be done on special spans/views. The code has been tested on linux with g++-7, g++-9, clang-11. On Windows mingw-10 and VS2019, OSX clang-11.
Diffstat (limited to 'src/lib/netlist/solver')
-rw-r--r--src/lib/netlist/solver/nld_matrix_solver.cpp22
-rw-r--r--src/lib/netlist/solver/nld_matrix_solver.h322
-rw-r--r--src/lib/netlist/solver/nld_matrix_solver_ext.h288
-rw-r--r--src/lib/netlist/solver/nld_ms_direct.h4
-rw-r--r--src/lib/netlist/solver/nld_ms_direct1.h1
-rw-r--r--src/lib/netlist/solver/nld_ms_direct2.h1
-rw-r--r--src/lib/netlist/solver/nld_ms_gcr.h15
-rw-r--r--src/lib/netlist/solver/nld_ms_gmres.h1
-rw-r--r--src/lib/netlist/solver/nld_ms_sm.h1
-rw-r--r--src/lib/netlist/solver/nld_ms_sor.h19
-rw-r--r--src/lib/netlist/solver/nld_ms_sor_mat.h3
-rw-r--r--src/lib/netlist/solver/nld_ms_w.h3
-rw-r--r--src/lib/netlist/solver/nld_solver.cpp31
-rw-r--r--src/lib/netlist/solver/nld_solver.h2
14 files changed, 364 insertions, 349 deletions
diff --git a/src/lib/netlist/solver/nld_matrix_solver.cpp b/src/lib/netlist/solver/nld_matrix_solver.cpp
index aeb5d11055f..87ae2de71d7 100644
--- a/src/lib/netlist/solver/nld_matrix_solver.cpp
+++ b/src/lib/netlist/solver/nld_matrix_solver.cpp
@@ -121,7 +121,7 @@ namespace solver
{
pstring nname(this->name() + "." + pstring(plib::pfmt("m{1}")(m_inps.size())));
nl_assert(p->net().is_analog());
- auto net_proxy_output_u = state().make_pool_object<proxied_analog_output_t>(*this, nname, static_cast<analog_net_t *>(&p->net()));
+ auto net_proxy_output_u = state().make_pool_object<proxied_analog_output_t>(*this, nname, &dynamic_cast<analog_net_t &>(p->net()));
net_proxy_output = net_proxy_output_u.get();
m_inps.emplace_back(std::move(net_proxy_output_u));
}
@@ -432,7 +432,7 @@ namespace solver
void matrix_solver_t::step(netlist_time delta) noexcept
{
- const auto dd(delta.as_fp<nl_fptype>());
+ const auto dd(delta.as_fp<fptype>());
for (auto &d : m_step_funcs)
d(dd);
}
@@ -489,7 +489,7 @@ namespace solver
if (m_params.m_dynamic_ts)
- return compute_next_timestep(delta.as_fp<nl_fptype>(), m_params.m_max_timestep);
+ return compute_next_timestep(delta.as_fp<fptype>(), m_params.m_max_timestep);
return netlist_time::from_fp(m_params.m_max_timestep);
}
@@ -534,7 +534,7 @@ namespace solver
return {colmax, colmin};
}
- nl_fptype matrix_solver_t::get_weight_around_diag(std::size_t row, std::size_t diag)
+ matrix_solver_t::fptype matrix_solver_t::get_weight_around_diag(std::size_t row, std::size_t diag)
{
{
//
@@ -543,7 +543,7 @@ namespace solver
std::vector<bool> touched(1024, false); // FIXME!
- nl_fptype weight = nlconst::zero();
+ fptype weight = nlconst::zero();
auto &term = m_terms[row];
for (std::size_t i = 0; i < term.count(); i++)
{
@@ -556,7 +556,7 @@ namespace solver
if (colu==row) colu = static_cast<unsigned>(diag);
else if (colu==diag) colu = static_cast<unsigned>(row);
- weight = weight + plib::abs(static_cast<nl_fptype>(colu) - static_cast<nl_fptype>(diag));
+ weight = weight + plib::abs(static_cast<fptype>(colu) - static_cast<fptype>(diag));
touched[colu] = true;
}
}
@@ -596,14 +596,14 @@ namespace solver
log().verbose(" has {1} dynamic elements", this->dynamic_device_count());
log().verbose(" has {1} timestep elements", this->timestep_device_count());
log().verbose(" {1:6.3} average newton raphson loops",
- static_cast<nl_fptype>(this->m_stat_newton_raphson) / static_cast<nl_fptype>(this->m_stat_vsolver_calls));
+ static_cast<fptype>(this->m_stat_newton_raphson) / static_cast<fptype>(this->m_stat_vsolver_calls));
log().verbose(" {1:10} invocations ({2:6.0} Hz) {3:10} gs fails ({4:6.2} %) {5:6.3} average",
this->m_stat_calculations,
- static_cast<nl_fptype>(this->m_stat_calculations) / this->exec().time().as_fp<nl_fptype>(),
+ static_cast<fptype>(this->m_stat_calculations) / this->exec().time().as_fp<fptype>(),
this->m_iterative_fail,
- nlconst::hundred() * static_cast<nl_fptype>(this->m_iterative_fail)
- / static_cast<nl_fptype>(this->m_stat_calculations),
- static_cast<nl_fptype>(this->m_iterative_total) / static_cast<nl_fptype>(this->m_stat_calculations));
+ nlconst::hundred() * static_cast<fptype>(this->m_iterative_fail)
+ / static_cast<fptype>(this->m_stat_calculations),
+ static_cast<fptype>(this->m_iterative_total) / static_cast<fptype>(this->m_stat_calculations));
}
}
diff --git a/src/lib/netlist/solver/nld_matrix_solver.h b/src/lib/netlist/solver/nld_matrix_solver.h
index 606ccbb9ee2..6a769775248 100644
--- a/src/lib/netlist/solver/nld_matrix_solver.h
+++ b/src/lib/netlist/solver/nld_matrix_solver.h
@@ -14,6 +14,7 @@
#include "netlist/plib/palloc.h"
#include "netlist/plib/penum.h"
#include "netlist/plib/pmatrix2d.h"
+#include "netlist/plib/pmempool.h"
#include "netlist/plib/putil.h"
#include "netlist/plib/vector_ops.h"
@@ -183,6 +184,8 @@ namespace solver
{
public:
using list_t = std::vector<matrix_solver_t *>;
+ using fptype = nl_fptype;
+ using arena_type = plib::mempool_arena<plib::aligned_arena>;
// after every call to solve, update inputs must be called.
// this can be done as well as a batch to ease parallel processing.
@@ -259,57 +262,31 @@ namespace solver
}
// return number of floating point operations for solve
- std::size_t ops() const { return m_ops; }
+ constexpr std::size_t ops() const { return m_ops; }
protected:
- template <typename T>
- using aligned_alloc = plib::aligned_allocator<T, PALIGN_VECTOROPT>;
-
matrix_solver_t(netlist_state_t &anetlist, const pstring &name,
const analog_net_t::list_t &nets,
const solver_parameters_t *params);
virtual void vsolve_non_dynamic() = 0;
- virtual netlist_time compute_next_timestep(nl_fptype cur_ts, nl_fptype max_ts) = 0;
- virtual bool check_err() = 0;
+ virtual netlist_time compute_next_timestep(fptype cur_ts, fptype max_ts) = 0;
+ virtual bool check_err() const = 0;
virtual void store() = 0;
-#if 0
- plib::pmatrix2d<nl_fptype, aligned_alloc<nl_fptype>> m_gonn;
- plib::pmatrix2d<nl_fptype, aligned_alloc<nl_fptype>> m_gtn;
- plib::pmatrix2d<nl_fptype, aligned_alloc<nl_fptype>> m_Idrn;
- plib::pmatrix2d<nl_fptype *, aligned_alloc<nl_fptype *>> m_connected_net_Vn;
-#else
- plib::pmatrix2d_vrl<nl_fptype, aligned_alloc<nl_fptype>> m_gonn;
- plib::pmatrix2d_vrl<nl_fptype, aligned_alloc<nl_fptype>> m_gtn;
- plib::pmatrix2d_vrl<nl_fptype, aligned_alloc<nl_fptype>> m_Idrn;
- plib::pmatrix2d_vrl<nl_fptype *, aligned_alloc<nl_fptype *>> m_connected_net_Vn;
-#endif
-
- plib::aligned_vector<terms_for_net_t> m_terms;
+ plib::pmatrix2d_vrl<fptype, arena_type> m_gonn;
+ plib::pmatrix2d_vrl<fptype, arena_type> m_gtn;
+ plib::pmatrix2d_vrl<fptype, arena_type> m_Idrn;
+ plib::pmatrix2d_vrl<fptype *, arena_type> m_connected_net_Vn;
const solver_parameters_t &m_params;
state_var<std::size_t> m_iterative_fail;
state_var<std::size_t> m_iterative_total;
- private:
-
- plib::aligned_vector<terms_for_net_t> m_rails_temp;
- std::vector<device_arena::unique_ptr<proxied_analog_output_t>> m_inps;
-
- state_var<std::size_t> m_stat_calculations;
- state_var<std::size_t> m_stat_newton_raphson;
- state_var<std::size_t> m_stat_vsolver_calls;
-
- state_var<netlist_time_ext> m_last_step;
- plib::aligned_vector<nldelegate_ts> m_step_funcs;
- plib::aligned_vector<nldelegate_dyn> m_dynamic_funcs;
+ plib::aligned_vector<terms_for_net_t> m_terms; // setup only
- logic_input_t m_fb_sync;
- logic_output_t m_Q_sync;
-
- std::size_t m_ops;
+ private:
// base setup - called from constructor
void setup_base(setup_t &setup, const analog_net_t::list_t &nets) noexcept(false);
@@ -321,7 +298,7 @@ namespace solver
int get_net_idx(const analog_net_t *net) const noexcept;
std::pair<int, int> get_left_right_of_diag(std::size_t irow, std::size_t idiag);
- nl_fptype get_weight_around_diag(std::size_t row, std::size_t diag);
+ fptype get_weight_around_diag(std::size_t row, std::size_t diag);
void add_term(std::size_t net_idx, terminal_t *term) noexcept(false);
@@ -332,272 +309,21 @@ namespace solver
analog_net_t *get_connected_net(terminal_t *term);
- };
-
- template <typename FT, int SIZE>
- class matrix_solver_ext_t: public matrix_solver_t
- {
- public:
-
- using float_type = FT;
-
- matrix_solver_ext_t(netlist_state_t &anetlist, const pstring &name,
- const analog_net_t::list_t &nets,
- const solver_parameters_t *params, const std::size_t size)
- : matrix_solver_t(anetlist, name, nets, params)
- , m_dim(size)
- , m_new_V(size)
- , m_RHS(size)
- , m_mat_ptr(size, this->max_railstart() + 1)
- , m_last_V(size, nlconst::zero())
- , m_DD_n_m_1(size, nlconst::zero())
- , m_h_n_m_1(size, nlconst::magic(1e-6)) // we need a non zero value here
- {
- //
- // save states
- //
- state().save(*this, m_last_V.as_base(), this->name(), "m_last_V");
- state().save(*this, m_DD_n_m_1.as_base(), this->name(), "m_DD_n_m_1");
- state().save(*this, m_h_n_m_1.as_base(), this->name(), "m_h_n_m_1");
- }
-
-
- private:
- const std::size_t m_dim;
-
- protected:
- static constexpr const std::size_t SIZEABS = plib::parray<FT, SIZE>::SIZEABS();
- static constexpr const std::size_t m_pitch_ABS = (((SIZEABS + 0) + 7) / 8) * 8;
-
- PALIGNAS_VECTOROPT()
- plib::parray<float_type, SIZE> m_new_V;
- PALIGNAS_VECTOROPT()
- plib::parray<float_type, SIZE> m_RHS;
-
- PALIGNAS_VECTOROPT()
- plib::pmatrix2d<float_type *> m_mat_ptr;
-
- // FIXME: below should be private
- // state - variable time_stepping
- PALIGNAS_VECTOROPT()
- plib::parray<nl_fptype, SIZE> m_last_V;
- PALIGNAS_VECTOROPT()
- plib::parray<nl_fptype, SIZE> m_DD_n_m_1;
- PALIGNAS_VECTOROPT()
- plib::parray<nl_fptype, SIZE> m_h_n_m_1;
-
- std::size_t max_railstart() const noexcept
- {
- std::size_t max_rail = 0;
- for (std::size_t k = 0; k < m_terms.size(); k++)
- max_rail = std::max(max_rail, m_terms[k].railstart());
- return max_rail;
- }
-
-
- template <typename T, typename M>
- void log_fill(const T &fill, M &mat)
- {
- const std::size_t iN = fill.size();
-
- // FIXME: Not yet working, mat_cr.h needs some more work
-#if 0
- auto mat_GE = dynamic_cast<plib::pGEmatrix_cr_t<typename M::base> *>(&mat);
-#else
- plib::unused_var(mat);
-#endif
- std::vector<unsigned> levL(iN, 0);
- std::vector<unsigned> levU(iN, 0);
-
- // parallel scheme for L x = y
- for (std::size_t k = 0; k < iN; k++)
- {
- unsigned lm=0;
- for (std::size_t j = 0; j<k; j++)
- if (fill[k][j] < M::FILL_INFINITY)
- lm = std::max(lm, levL[j]);
- levL[k] = 1+lm;
- }
-
- // parallel scheme for U x = y
- for (std::size_t k = iN; k-- > 0; )
- {
- unsigned lm=0;
- for (std::size_t j = iN; --j > k; )
- if (fill[k][j] < M::FILL_INFINITY)
- lm = std::max(lm, levU[j]);
- levU[k] = 1+lm;
- }
- for (std::size_t k = 0; k < iN; k++)
- {
- unsigned fm = 0;
- pstring ml = "";
- for (std::size_t j = 0; j < iN; j++)
- {
- ml += fill[k][j] == 0 ? 'X' : fill[k][j] < M::FILL_INFINITY ? '+' : '.';
- if (fill[k][j] < M::FILL_INFINITY)
- if (fill[k][j] > fm)
- fm = fill[k][j];
- }
-#if 0
- this->log().verbose("{1:4} {2} {3:4} {4:4} {5:4} {6:4}", k, ml,
- levL[k], levU[k], mat_GE ? mat_GE->get_parallel_level(k) : 0, fm);
-#else
- this->log().verbose("{1:4} {2} {3:4} {4:4} {5:4} {6:4}", k, ml,
- levL[k], levU[k], 0, fm);
-#endif
- }
- }
-
- constexpr std::size_t size() const noexcept
- {
- return (SIZE > 0) ? static_cast<std::size_t>(SIZE) : m_dim;
- }
-
-#if 1
- void store() override
- {
- const std::size_t iN = size();
- for (std::size_t i = 0; i < iN; i++)
- this->m_terms[i].setV(static_cast<nl_fptype>(m_new_V[i]));
- }
-#else
- // global tanh damping (4.197)
- // partially cures the symptoms but not the cause
- void store() override
- {
- const std::size_t iN = size();
- for (std::size_t i = 0; i < iN; i++)
- {
- auto oldV = this->m_terms[i].template getV<nl_fptype>();
- this->m_terms[i].setV(oldV + 0.02 * plib::tanh((m_new_V[i]-oldV)*50.0));
- }
- }
-#endif
- bool check_err() override
- {
- // NOTE: Ideally we should also include currents (RHS) here. This would
- // need a reevaluation of the right hand side after voltages have been updated
- // and thus belong into a different calculation. This applies to all solvers.
-
- const std::size_t iN = size();
- const auto reltol(static_cast<float_type>(m_params.m_reltol));
- const auto vntol(static_cast<float_type>(m_params.m_vntol));
- for (std::size_t i = 0; i < iN; i++)
- {
- const auto vold(static_cast<float_type>(this->m_terms[i].getV()));
- const auto vnew(m_new_V[i]);
- const auto tol(vntol + reltol * std::max(plib::abs(vnew),plib::abs(vold)));
- if (plib::abs(vnew - vold) > tol)
- return true;
- }
- return false;
- }
-
- netlist_time compute_next_timestep(nl_fptype cur_ts, nl_fptype max_ts) override
- {
- nl_fptype new_solver_timestep(max_ts);
-
- for (std::size_t k = 0; k < size(); k++)
- {
- const auto &t = m_terms[k];
- const auto v(static_cast<nl_fptype>(t.getV()));
- // avoid floating point exceptions
- const nl_fptype DD_n = std::max(-fp_constants<nl_fptype>::TIMESTEP_MAXDIFF(),
- std::min(+fp_constants<nl_fptype>::TIMESTEP_MAXDIFF(),(v - m_last_V[k])));
-
- m_last_V[k] = v;
- const nl_fptype hn = cur_ts;
-
- nl_fptype DD2 = (DD_n / hn - m_DD_n_m_1[k] / m_h_n_m_1[k]) / (hn + m_h_n_m_1[k]);
- nl_fptype new_net_timestep(0);
-
- m_h_n_m_1[k] = hn;
- m_DD_n_m_1[k] = DD_n;
- if (plib::abs(DD2) > fp_constants<nl_fptype>::TIMESTEP_MINDIV()) // avoid div-by-zero
- new_net_timestep = plib::sqrt(m_params.m_dynamic_lte / plib::abs(nlconst::half()*DD2));
- else
- new_net_timestep = m_params.m_max_timestep;
-
- new_solver_timestep = std::min(new_net_timestep, new_solver_timestep);
- }
- new_solver_timestep = std::max(new_solver_timestep, m_params.m_min_timestep);
-
- // FIXME: Factor 2 below is important. Without, we get timing issues. This must be a bug elsewhere.
- return std::max(netlist_time::from_fp(new_solver_timestep), netlist_time::quantum() * 2);
- }
-
- template <typename M>
- void build_mat_ptr(M &mat)
- {
- const std::size_t iN = size();
-
- for (std::size_t k=0; k<iN; k++)
- {
- std::size_t cnt(0);
- // build pointers into the compressed row format matrix for each terminal
- for (std::size_t j=0; j< this->m_terms[k].railstart();j++)
- {
- int other = this->m_terms[k].m_connected_net_idx[j];
- if (other >= 0)
- {
- m_mat_ptr[k][j] = &(mat[k][static_cast<std::size_t>(other)]);
- cnt++;
- }
- }
- nl_assert_always(cnt == this->m_terms[k].railstart(), "Count and railstart mismatch");
- m_mat_ptr[k][this->m_terms[k].railstart()] = &(mat[k][k]);
- }
- }
-
- template <typename M>
- void clear_square_mat(M &m)
- {
- const std::size_t n = size();
- for (std::size_t k=0; k < n; k++)
- {
- auto *p = &(m[k][0]);
- using mat_elem_type = typename std::decay<decltype(*p)>::type;
- for (std::size_t i=0; i < n; i++)
- p[i] = plib::constants<mat_elem_type>::zero();
- }
- }
-
- void fill_matrix_and_rhs()
- {
- const std::size_t N = size();
-
- for (std::size_t k = 0; k < N; k++)
- {
- auto &net = m_terms[k];
- auto **tcr_r = &(m_mat_ptr[k][0]);
-
- using source_type = typename decltype(m_gtn)::value_type;
- const std::size_t term_count = net.count();
- const std::size_t railstart = net.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];
-
- // FIXME: gonn, gtn and Idr - which float types should they have?
-
- auto gtot_t = std::accumulate(gt, gt + term_count, plib::constants<source_type>::zero());
-
- // update diagonal element ...
- *tcr_r[railstart] = static_cast<FT>(gtot_t); //mat.A[mat.diag[k]] += gtot_t;
+ state_var<std::size_t> m_stat_calculations;
+ state_var<std::size_t> m_stat_newton_raphson;
+ state_var<std::size_t> m_stat_vsolver_calls;
- for (std::size_t i = 0; i < railstart; i++)
- *tcr_r[i] += static_cast<FT>(go[i]);
+ state_var<netlist_time_ext> m_last_step;
+ plib::aligned_vector<nldelegate_ts> m_step_funcs;
+ plib::aligned_vector<nldelegate_dyn> m_dynamic_funcs;
+ plib::aligned_vector<device_arena::unique_ptr<proxied_analog_output_t>> m_inps;
- auto RHS_t(std::accumulate(Idr, Idr + term_count, plib::constants<source_type>::zero()));
+ logic_input_t m_fb_sync;
+ logic_output_t m_Q_sync;
- for (std::size_t i = railstart; i < term_count; i++)
- RHS_t += (- go[i]) * *cnV[i];
+ std::size_t m_ops;
- m_RHS[k] = static_cast<FT>(RHS_t);
- }
- }
+ plib::aligned_vector<terms_for_net_t> m_rails_temp; // setup only
};
diff --git a/src/lib/netlist/solver/nld_matrix_solver_ext.h b/src/lib/netlist/solver/nld_matrix_solver_ext.h
new file mode 100644
index 00000000000..275ef135c69
--- /dev/null
+++ b/src/lib/netlist/solver/nld_matrix_solver_ext.h
@@ -0,0 +1,288 @@
+// license:GPL-2.0+
+// copyright-holders:Couriersud
+
+#ifndef NLD_MATRIX_SOLVER_EXT_H_
+#define NLD_MATRIX_SOLVER_EXT_H_
+
+///
+/// \file nld_matrix_solver.h
+///
+
+#include "nld_matrix_solver.h"
+
+#include <numeric>
+
+namespace netlist
+{
+namespace solver
+{
+
+ template <typename FT, int SIZE>
+ class matrix_solver_ext_t: public matrix_solver_t
+ {
+ public:
+
+ using float_type = FT;
+
+ matrix_solver_ext_t(netlist_state_t &anetlist, const pstring &name,
+ const analog_net_t::list_t &nets,
+ const solver_parameters_t *params, const std::size_t size)
+ : matrix_solver_t(anetlist, name, nets, params)
+ , m_new_V(size)
+ , m_RHS(size)
+ , m_mat_ptr(size, this->max_railstart() + 1)
+ , m_last_V(size, nlconst::zero())
+ , m_DD_n_m_1(size, nlconst::zero())
+ , m_h_n_m_1(size, nlconst::magic(1e-6)) // we need a non zero value here
+ , m_dim(size)
+ {
+ //
+ // save states
+ //
+ state().save(*this, m_last_V.as_base(), this->name(), "m_last_V");
+ state().save(*this, m_DD_n_m_1.as_base(), this->name(), "m_DD_n_m_1");
+ state().save(*this, m_h_n_m_1.as_base(), this->name(), "m_h_n_m_1");
+ }
+
+ protected:
+ static constexpr const std::size_t SIZEABS = plib::parray<FT, SIZE>::SIZEABS();
+ static constexpr const std::size_t m_pitch_ABS = (((SIZEABS + 0) + 7) / 8) * 8;
+
+ //PALIGNAS_VECTOROPT() parrays define alignment already
+ plib::parray<float_type, SIZE> m_new_V;
+ //PALIGNAS_VECTOROPT() parrays define alignment already
+ plib::parray<float_type, SIZE> m_RHS;
+
+ //PALIGNAS_VECTOROPT() parrays define alignment already
+ plib::pmatrix2d<float_type *> m_mat_ptr;
+
+ std::size_t max_railstart() const noexcept
+ {
+ std::size_t max_rail = 0;
+ for (std::size_t k = 0; k < m_terms.size(); k++)
+ max_rail = std::max(max_rail, m_terms[k].railstart());
+ return max_rail;
+ }
+
+
+ template <typename T, typename M>
+ void log_fill(const T &fill, M &mat)
+ {
+ const std::size_t iN = fill.size();
+
+ // FIXME: Not yet working, mat_cr.h needs some more work
+#if 0
+ auto mat_GE = dynamic_cast<plib::pGEmatrix_cr_t<typename M::base> *>(&mat);
+#else
+ plib::unused_var(mat);
+#endif
+ std::vector<unsigned> levL(iN, 0);
+ std::vector<unsigned> levU(iN, 0);
+
+ // parallel scheme for L x = y
+ for (std::size_t k = 0; k < iN; k++)
+ {
+ unsigned lm=0;
+ for (std::size_t j = 0; j<k; j++)
+ if (fill[k][j] < M::FILL_INFINITY)
+ lm = std::max(lm, levL[j]);
+ levL[k] = 1+lm;
+ }
+
+ // parallel scheme for U x = y
+ for (std::size_t k = iN; k-- > 0; )
+ {
+ unsigned lm=0;
+ for (std::size_t j = iN; --j > k; )
+ if (fill[k][j] < M::FILL_INFINITY)
+ lm = std::max(lm, levU[j]);
+ levU[k] = 1+lm;
+ }
+ for (std::size_t k = 0; k < iN; k++)
+ {
+ unsigned fm = 0;
+ pstring ml = "";
+ for (std::size_t j = 0; j < iN; j++)
+ {
+ ml += fill[k][j] == 0 ? 'X' : fill[k][j] < M::FILL_INFINITY ? '+' : '.';
+ if (fill[k][j] < M::FILL_INFINITY)
+ if (fill[k][j] > fm)
+ fm = fill[k][j];
+ }
+#if 0
+ this->log().verbose("{1:4} {2} {3:4} {4:4} {5:4} {6:4}", k, ml,
+ levL[k], levU[k], mat_GE ? mat_GE->get_parallel_level(k) : 0, fm);
+#else
+ this->log().verbose("{1:4} {2} {3:4} {4:4} {5:4} {6:4}", k, ml,
+ levL[k], levU[k], 0, fm);
+#endif
+ }
+ }
+
+ constexpr std::size_t size() const noexcept
+ {
+ return (SIZE > 0) ? static_cast<std::size_t>(SIZE) : m_dim;
+ }
+
+#if 1
+ void store() override
+ {
+ const std::size_t iN = size();
+ for (std::size_t i = 0; i < iN; i++)
+ this->m_terms[i].setV(static_cast<fptype>(m_new_V[i]));
+ }
+#else
+ // global tanh damping (4.197)
+ // partially cures the symptoms but not the cause
+ void store() override
+ {
+ const std::size_t iN = size();
+ for (std::size_t i = 0; i < iN; i++)
+ {
+ auto oldV = this->m_terms[i].template getV<fptype>();
+ this->m_terms[i].setV(oldV + 0.02 * plib::tanh((m_new_V[i]-oldV)*50.0));
+ }
+ }
+#endif
+ bool check_err() const override
+ {
+ // NOTE: Ideally we should also include currents (RHS) here. This would
+ // need a reevaluation of the right hand side after voltages have been updated
+ // and thus belong into a different calculation. This applies to all solvers.
+
+ const std::size_t iN = size();
+ const auto reltol(static_cast<float_type>(m_params.m_reltol));
+ const auto vntol(static_cast<float_type>(m_params.m_vntol));
+ for (std::size_t i = 0; i < iN; i++)
+ {
+ const auto vold(static_cast<float_type>(this->m_terms[i].getV()));
+ const auto vnew(m_new_V[i]);
+ const auto tol(vntol + reltol * std::max(plib::abs(vnew),plib::abs(vold)));
+ if (plib::abs(vnew - vold) > tol)
+ return true;
+ }
+ return false;
+ }
+
+ netlist_time compute_next_timestep(fptype cur_ts, fptype max_ts) override
+ {
+ fptype new_solver_timestep(max_ts);
+
+ for (std::size_t k = 0; k < size(); k++)
+ {
+ const auto &t = m_terms[k];
+ const auto v(static_cast<fptype>(t.getV()));
+ // avoid floating point exceptions
+ const fptype DD_n = std::max(-fp_constants<fptype>::TIMESTEP_MAXDIFF(),
+ std::min(+fp_constants<fptype>::TIMESTEP_MAXDIFF(),(v - m_last_V[k])));
+
+ m_last_V[k] = v;
+ const fptype hn = cur_ts;
+
+ fptype DD2 = (DD_n / hn - m_DD_n_m_1[k] / m_h_n_m_1[k]) / (hn + m_h_n_m_1[k]);
+ fptype new_net_timestep(0);
+
+ m_h_n_m_1[k] = hn;
+ m_DD_n_m_1[k] = DD_n;
+ if (plib::abs(DD2) > fp_constants<fptype>::TIMESTEP_MINDIV()) // avoid div-by-zero
+ new_net_timestep = plib::sqrt(m_params.m_dynamic_lte / plib::abs(nlconst::half()*DD2));
+ else
+ new_net_timestep = m_params.m_max_timestep;
+
+ new_solver_timestep = std::min(new_net_timestep, new_solver_timestep);
+ }
+ new_solver_timestep = std::max(new_solver_timestep, m_params.m_min_timestep);
+
+ // FIXME: Factor 2 below is important. Without, we get timing issues. This must be a bug elsewhere.
+ return std::max(netlist_time::from_fp(new_solver_timestep), netlist_time::quantum() * 2);
+ }
+
+ template <typename M>
+ void build_mat_ptr(M &mat)
+ {
+ const std::size_t iN = size();
+
+ for (std::size_t k=0; k<iN; k++)
+ {
+ std::size_t cnt(0);
+ // build pointers into the compressed row format matrix for each terminal
+ for (std::size_t j=0; j< this->m_terms[k].railstart();j++)
+ {
+ int other = this->m_terms[k].m_connected_net_idx[j];
+ if (other >= 0)
+ {
+ m_mat_ptr[k][j] = &(mat[k][static_cast<std::size_t>(other)]);
+ cnt++;
+ }
+ }
+ nl_assert_always(cnt == this->m_terms[k].railstart(), "Count and railstart mismatch");
+ m_mat_ptr[k][this->m_terms[k].railstart()] = &(mat[k][k]);
+ }
+ }
+
+ template <typename M>
+ void clear_square_mat(M &m)
+ {
+ const std::size_t n = size();
+ for (std::size_t k=0; k < n; k++)
+ {
+ auto *p = &(m[k][0]);
+ using mat_elem_type = typename std::decay<decltype(*p)>::type;
+ for (std::size_t i=0; i < n; i++)
+ p[i] = plib::constants<mat_elem_type>::zero();
+ }
+ }
+
+ void fill_matrix_and_rhs()
+ {
+ const std::size_t N = size();
+
+ for (std::size_t k = 0; k < N; k++)
+ {
+ auto &net = m_terms[k];
+ auto **tcr_r = &(m_mat_ptr[k][0]);
+
+ using source_type = typename decltype(m_gtn)::value_type;
+ const std::size_t term_count = net.count();
+ const std::size_t railstart = net.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];
+
+ // FIXME: gonn, gtn and Idr - which float types should they have?
+
+ auto gtot_t = std::accumulate(gt, gt + term_count, plib::constants<source_type>::zero());
+
+ // update diagonal element ...
+ *tcr_r[railstart] = static_cast<FT>(gtot_t); //mat.A[mat.diag[k]] += gtot_t;
+
+ for (std::size_t i = 0; i < railstart; i++)
+ *tcr_r[i] += static_cast<FT>(go[i]);
+
+ auto RHS_t(std::accumulate(Idr, Idr + term_count, plib::constants<source_type>::zero()));
+
+ for (std::size_t i = railstart; i < term_count; i++)
+ RHS_t += (- go[i]) * *cnV[i];
+
+ m_RHS[k] = static_cast<FT>(RHS_t);
+ }
+ }
+
+ private:
+ // state - variable time_stepping
+ //PALIGNAS_VECTOROPT() parrays define alignment already
+ plib::parray<fptype, SIZE> m_last_V;
+ // PALIGNAS_VECTOROPT() parrays define alignment already
+ plib::parray<fptype, SIZE> m_DD_n_m_1;
+ // PALIGNAS_VECTOROPT() parrays define alignment already
+ plib::parray<fptype, SIZE> m_h_n_m_1;
+
+ const std::size_t m_dim;
+
+ };
+
+} // namespace solver
+} // namespace netlist
+
+#endif // NLD_MATRIX_SOLVER_EXT_H_
diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h
index ab4df0758e1..b32e4d2172b 100644
--- a/src/lib/netlist/solver/nld_ms_direct.h
+++ b/src/lib/netlist/solver/nld_ms_direct.h
@@ -13,6 +13,8 @@
#include "plib/parray.h"
#include "plib/vector_ops.h"
+#include "nld_matrix_solver_ext.h"
+
#include <algorithm>
namespace netlist
@@ -49,7 +51,7 @@ namespace solver
template <typename T>
void LE_back_subst(T & x);
- PALIGNAS_VECTOROPT()
+ // PALIGNAS_VECTOROPT() parrays define alignment already
plib::parray2D<FT, SIZE, m_pitch_ABS> m_A;
};
diff --git a/src/lib/netlist/solver/nld_ms_direct1.h b/src/lib/netlist/solver/nld_ms_direct1.h
index 095d35a30c0..e33846ad403 100644
--- a/src/lib/netlist/solver/nld_ms_direct1.h
+++ b/src/lib/netlist/solver/nld_ms_direct1.h
@@ -8,6 +8,7 @@
/// \file nld_ms_direct1.h
///
+#include "nld_matrix_solver_ext.h"
#include "nld_ms_direct.h"
#include "nld_solver.h"
diff --git a/src/lib/netlist/solver/nld_ms_direct2.h b/src/lib/netlist/solver/nld_ms_direct2.h
index 744737b3892..365c1cac33e 100644
--- a/src/lib/netlist/solver/nld_ms_direct2.h
+++ b/src/lib/netlist/solver/nld_ms_direct2.h
@@ -8,6 +8,7 @@
/// \file nld_ms_direct2.h
///
+#include "nld_matrix_solver_ext.h"
#include "nld_ms_direct.h"
#include "nld_solver.h"
diff --git a/src/lib/netlist/solver/nld_ms_gcr.h b/src/lib/netlist/solver/nld_ms_gcr.h
index 6c7f5d7e502..3de955627ee 100644
--- a/src/lib/netlist/solver/nld_ms_gcr.h
+++ b/src/lib/netlist/solver/nld_ms_gcr.h
@@ -12,7 +12,8 @@
#include "plib/mat_cr.h"
-#include "nld_ms_direct.h"
+//#include "nld_ms_direct.h"
+#include "nld_matrix_solver_ext.h"
#include "nld_solver.h"
#include "plib/pdynlib.h"
#include "plib/pstream.h"
@@ -31,6 +32,8 @@ namespace solver
public:
using mat_type = plib::pGEmatrix_cr_t<plib::pmatrix_cr_t<FT, SIZE>>;
+ using base_type = matrix_solver_ext_t<FT, SIZE>;
+ using fptype = typename base_type::fptype;
matrix_solver_GCR_t(netlist_state_t &anetlist, const pstring &name,
const analog_net_t::list_t &nets,
@@ -85,9 +88,9 @@ namespace solver
anetlist.log().verbose("maximum fill: {1}", gr.first);
anetlist.log().verbose("Post elimination occupancy ratio: {2} Ops: {1}", gr.second,
- static_cast<nl_fptype>(mat.nz_num) / static_cast<nl_fptype>(iN * iN));
+ static_cast<fptype>(mat.nz_num) / static_cast<fptype>(iN * iN));
anetlist.log().verbose(" Pre elimination occupancy ratio: {2}",
- static_cast<nl_fptype>(raw_elements) / static_cast<nl_fptype>(iN * iN));
+ static_cast<fptype>(raw_elements) / static_cast<fptype>(iN * iN));
// FIXME: Move me
//
@@ -123,7 +126,7 @@ namespace solver
pstring static_compile_name();
mat_type mat;
- plib::dynproc<void, FT *, nl_fptype *, nl_fptype *, nl_fptype *, nl_fptype ** > m_proc;
+ plib::dynproc<void, FT *, fptype *, fptype *, fptype *, fptype ** > m_proc;
};
@@ -273,7 +276,7 @@ namespace solver
pstring matrix_solver_GCR_t<FT, SIZE>::static_compile_name()
{
pstring str_floattype(fp_constants<FT>::name());
- pstring str_fptype(fp_constants<nl_fptype>::name());
+ pstring str_fptype(fp_constants<fptype>::name());
std::stringstream t;
t.imbue(std::locale::classic());
plib::putf8_fmt_writer w(&t);
@@ -290,7 +293,7 @@ namespace solver
plib::putf8_fmt_writer strm(&t);
pstring name = static_compile_name();
pstring str_floattype(fp_constants<FT>::name());
- pstring str_fptype(fp_constants<nl_fptype>::name());
+ pstring str_fptype(fp_constants<fptype>::name());
pstring extqual;
if (target == CXX_EXTERNAL_C)
diff --git a/src/lib/netlist/solver/nld_ms_gmres.h b/src/lib/netlist/solver/nld_ms_gmres.h
index add885c2bfb..837c0bc780e 100644
--- a/src/lib/netlist/solver/nld_ms_gmres.h
+++ b/src/lib/netlist/solver/nld_ms_gmres.h
@@ -8,6 +8,7 @@
/// \file nld_ms_gmres.h
///
+#include "nld_matrix_solver_ext.h"
#include "nld_ms_direct.h"
#include "nld_solver.h"
#include "plib/gmres.h"
diff --git a/src/lib/netlist/solver/nld_ms_sm.h b/src/lib/netlist/solver/nld_ms_sm.h
index 599252adf6a..fd11a32a102 100644
--- a/src/lib/netlist/solver/nld_ms_sm.h
+++ b/src/lib/netlist/solver/nld_ms_sm.h
@@ -33,6 +33,7 @@
///
#include "nld_matrix_solver.h"
+#include "nld_matrix_solver_ext.h"
#include "nld_solver.h"
#include "plib/vector_ops.h"
diff --git a/src/lib/netlist/solver/nld_ms_sor.h b/src/lib/netlist/solver/nld_ms_sor.h
index a253f160c4a..6b1a21afd39 100644
--- a/src/lib/netlist/solver/nld_ms_sor.h
+++ b/src/lib/netlist/solver/nld_ms_sor.h
@@ -12,8 +12,8 @@
/// Fow w==1 we will do the classic Gauss-Seidel approach.
///
+#include "nld_matrix_solver_ext.h"
#include "nld_ms_direct.h"
-#include "nld_solver.h"
#include <algorithm>
@@ -70,16 +70,19 @@ namespace solver
for (std::size_t k = 0; k < iN; k++)
{
- nl_fptype gtot_t = nlconst::zero();
- nl_fptype gabs_t = nlconst::zero();
- nl_fptype RHS_t = nlconst::zero();
const std::size_t term_count = this->m_terms[k].count();
- const nl_fptype * const gt = this->m_gtn[k];
- const nl_fptype * const go = this->m_gonn[k];
- const nl_fptype * const Idr = this->m_Idrn[k];
+ const auto * const gt = this->m_gtn[k];
+ const auto * const go = this->m_gonn[k];
+ const auto * const Idr = this->m_Idrn[k];
auto other_cur_analog = this->m_connected_net_Vn[k];
+ using fpaggtype = std::remove_reference_t<std::remove_cv_t<decltype(this->m_gtn[0][0])>>;
+
+ fpaggtype gtot_t = nlconst_base<fpaggtype>::zero();
+ fpaggtype gabs_t = nlconst_base<fpaggtype>::zero();
+ fpaggtype RHS_t = nlconst_base<fpaggtype>::zero();
+
this->m_new_V[k] = static_cast<float_type>(this->m_terms[k].getV());
for (std::size_t i = 0; i < term_count; i++)
@@ -126,7 +129,7 @@ namespace solver
{
const int * net_other = this->m_terms[k].m_connected_net_idx.data();
const std::size_t railstart = this->m_terms[k].railstart();
- const nl_fptype * go = this->m_gonn[k];
+ const auto * go = this->m_gonn[k];
float_type Idrive = plib::constants<float_type>::zero();
for (std::size_t i = 0; i < railstart; i++)
diff --git a/src/lib/netlist/solver/nld_ms_sor_mat.h b/src/lib/netlist/solver/nld_ms_sor_mat.h
index d575cbb1f37..2d29d0f01cd 100644
--- a/src/lib/netlist/solver/nld_ms_sor_mat.h
+++ b/src/lib/netlist/solver/nld_ms_sor_mat.h
@@ -12,9 +12,8 @@
/// For w==1 we will do the classic Gauss-Seidel approach
///
-#include "nld_matrix_solver.h"
+#include "nld_matrix_solver_ext.h"
#include "nld_ms_direct.h"
-#include "nld_solver.h"
#include <algorithm>
diff --git a/src/lib/netlist/solver/nld_ms_w.h b/src/lib/netlist/solver/nld_ms_w.h
index 8571d459ad8..598a4b9e516 100644
--- a/src/lib/netlist/solver/nld_ms_w.h
+++ b/src/lib/netlist/solver/nld_ms_w.h
@@ -40,8 +40,7 @@
/// introduces numerical instability.
///
-#include "nld_matrix_solver.h"
-#include "nld_solver.h"
+#include "nld_matrix_solver_ext.h"
#include "plib/vector_ops.h"
#include <algorithm>
diff --git a/src/lib/netlist/solver/nld_solver.cpp b/src/lib/netlist/solver/nld_solver.cpp
index 8787247e505..4e0d8f4c4b3 100644
--- a/src/lib/netlist/solver/nld_solver.cpp
+++ b/src/lib/netlist/solver/nld_solver.cpp
@@ -87,7 +87,7 @@ namespace devices
analog_net_t::list_t &nets,
solver::solver_parameters_t &params, std::size_t size)
{
- return host_arena::make_unique<C>(nl, name, nets, &params, size);
+ return plib::make_unique<C, host_arena>(nl, name, nets, &params, size);
}
template <typename FT, int SIZE>
@@ -136,29 +136,21 @@ namespace devices
switch (net_count)
{
case 1:
- return host_arena::make_unique<solver::matrix_solver_direct1_t<FT>>(state(), sname, nets, &m_params);
- break;
+ return plib::make_unique<solver::matrix_solver_direct1_t<FT>, host_arena>(state(), sname, nets, &m_params);
case 2:
- return host_arena::make_unique<solver::matrix_solver_direct2_t<FT>>(state(), sname, nets, &m_params);
- break;
+ return plib::make_unique<solver::matrix_solver_direct2_t<FT>, host_arena>(state(), sname, nets, &m_params);
case 3:
return create_solver<FT, 3>(3, sname, nets);
- break;
case 4:
return create_solver<FT, 4>(4, sname, nets);
- break;
case 5:
return create_solver<FT, 5>(5, sname, nets);
- break;
case 6:
return create_solver<FT, 6>(6, sname, nets);
- break;
case 7:
return create_solver<FT, 7>(7, sname, nets);
- break;
case 8:
return create_solver<FT, 8>(8, sname, nets);
- break;
default:
log().info(MI_NO_SPECIFIC_SOLVER(net_count));
if (net_count <= 16)
@@ -186,7 +178,6 @@ namespace devices
return create_solver<FT, -512>(net_count, sname, nets);
}
return create_solver<FT, 0>(net_count, sname, nets);
- break;
}
}
@@ -201,7 +192,7 @@ namespace devices
{
netlist.log().verbose(" ==> not a rail net");
// Must be an analog net
- auto &n = *static_cast<analog_net_t *>(net.get());
+ auto &n = dynamic_cast<analog_net_t &>(*net);
if (!already_processed(n))
{
groupspre.emplace_back(analog_net_t::list_t());
@@ -273,9 +264,9 @@ namespace devices
// only process analog terminals
if (term->is_type(detail::terminal_type::TERMINAL))
{
- auto *pt = static_cast<terminal_t *>(term);
+ auto &pt = dynamic_cast<terminal_t &>(*term);
// check the connected terminal
- analog_net_t &connected_net = netlist.setup().get_connected_terminal(*pt)->net();
+ analog_net_t &connected_net = netlist.setup().get_connected_terminal(pt)->net();
netlist.log().verbose(" Connected net {}", connected_net.name());
if (!check_if_processed_and_join(connected_net))
process_net(netlist, connected_net);
@@ -306,17 +297,17 @@ namespace devices
switch (m_params.m_fp_type())
{
case solver::matrix_fp_type_e::FLOAT:
- if (!NL_USE_FLOAT_MATRIX)
+ if (!config::use_float_matrix())
log().info("FPTYPE {1} not supported. Using DOUBLE", m_params.m_fp_type().name());
- ms = create_solvers<std::conditional_t<NL_USE_FLOAT_MATRIX,float, double>>(sname, grp);
+ ms = create_solvers<std::conditional_t<config::use_float_matrix::value, float, double>>(sname, grp);
break;
case solver::matrix_fp_type_e::DOUBLE:
ms = create_solvers<double>(sname, grp);
break;
case solver::matrix_fp_type_e::LONGDOUBLE:
- if (!NL_USE_LONG_DOUBLE_MATRIX)
+ if (!config::use_long_double_matrix())
log().info("FPTYPE {1} not supported. Using DOUBLE", m_params.m_fp_type().name());
- ms = create_solvers<std::conditional_t<NL_USE_LONG_DOUBLE_MATRIX, long double, double>>(sname, grp);
+ ms = create_solvers<std::conditional_t<config::use_long_double_matrix::value, long double, double>>(sname, grp);
break;
case solver::matrix_fp_type_e::FLOATQ128:
#if (NL_USE_FLOAT128)
@@ -355,7 +346,7 @@ namespace devices
for (auto & s : m_mat_solvers)
{
auto r = s->create_solver_code(target);
- if (r.first != "") // ignore solvers not supporting static compile
+ if (!r.first.empty()) // ignore solvers not supporting static compile
mp.push_back(r);
}
return mp;
diff --git a/src/lib/netlist/solver/nld_solver.h b/src/lib/netlist/solver/nld_solver.h
index a0c28e6b302..e9734af1520 100644
--- a/src/lib/netlist/solver/nld_solver.h
+++ b/src/lib/netlist/solver/nld_solver.h
@@ -39,7 +39,7 @@ namespace devices
void post_start();
void stop();
- nl_fptype gmin() const { return m_params.m_gmin(); }
+ auto gmin() const -> decltype(solver::solver_parameters_t::m_gmin()) { return m_params.m_gmin(); }
solver::static_compile_container create_solver_code(solver::static_compile_target target);