summaryrefslogtreecommitdiffstatshomepage
path: root/src/lib/netlist/solver/nld_matrix_solver.h
diff options
context:
space:
mode:
Diffstat (limited to 'src/lib/netlist/solver/nld_matrix_solver.h')
-rw-r--r--src/lib/netlist/solver/nld_matrix_solver.h322
1 files changed, 24 insertions, 298 deletions
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
};