mirror of
https://github.com/holub/mame
synced 2025-07-05 18:08:04 +03:00
netlist: further solver refactoring. (nw)
This commit is contained in:
parent
e37b9b7b88
commit
47d938b149
@ -105,6 +105,8 @@ namespace plib {
|
||||
throw plib::pexception("parray: size error " + plib::to_string(size) + ">" + plib::to_string(SIZE));
|
||||
}
|
||||
|
||||
base_type &as_base() noexcept { return m_a; }
|
||||
|
||||
inline size_type size() const noexcept { return SIZE <= 0 ? m_size : SIZEABS(); }
|
||||
|
||||
constexpr size_type max_size() const noexcept { return base_type::max_size(); }
|
||||
|
@ -341,14 +341,6 @@ namespace solver
|
||||
* save states
|
||||
*/
|
||||
|
||||
m_last_V.resize(iN, plib::constants<nl_fptype>::zero());
|
||||
m_DD_n_m_1.resize(iN, plib::constants<nl_fptype>::zero());
|
||||
m_h_n_m_1.resize(iN, plib::constants<nl_fptype>::zero());
|
||||
|
||||
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");
|
||||
|
||||
for (std::size_t k = 0; k < iN; k++)
|
||||
{
|
||||
pstring num = plib::pfmt("{1}")(k);
|
||||
@ -553,50 +545,6 @@ namespace solver
|
||||
}
|
||||
}
|
||||
|
||||
netlist_time matrix_solver_t::compute_next_timestep(const nl_fptype cur_ts)
|
||||
{
|
||||
nl_fptype new_solver_timestep = m_params.m_max_timestep;
|
||||
|
||||
if (m_params.m_dynamic_ts)
|
||||
{
|
||||
for (std::size_t k = 0; k < m_terms.size(); k++)
|
||||
{
|
||||
auto &t = m_terms[k];
|
||||
//const nl_fptype DD_n = (n->Q_Analog() - t->m_last_V);
|
||||
// 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,(t.getV() - m_last_V[k])));
|
||||
const nl_fptype hn = cur_ts;
|
||||
|
||||
//printf("%g %g %g %g\n", DD_n, hn, t.m_DD_n_m_1, t.m_h_n_m_1);
|
||||
nl_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 (std::fabs(DD2) > fp_constants<nl_fptype>::TIMESTEP_MINDIV) // avoid div-by-zero
|
||||
new_net_timestep = std::sqrt(m_params.m_dynamic_lte / std::fabs(plib::constants<nl_fptype>::cast(0.5)*DD2));
|
||||
else
|
||||
new_net_timestep = m_params.m_max_timestep;
|
||||
|
||||
if (new_net_timestep < new_solver_timestep)
|
||||
new_solver_timestep = new_net_timestep;
|
||||
|
||||
m_last_V[k] = t.getV();
|
||||
}
|
||||
if (new_solver_timestep < m_params.m_min_timestep)
|
||||
{
|
||||
new_solver_timestep = m_params.m_min_timestep;
|
||||
}
|
||||
}
|
||||
//if (new_solver_timestep > 10.0 * hn)
|
||||
// new_solver_timestep = 10.0 * hn;
|
||||
/*
|
||||
* FIXME: Factor 2 below is important. Without, we get timing issues. This must be a bug elsewhere.
|
||||
*/
|
||||
return std::max(netlist_time::from_double(new_solver_timestep), netlist_time::quantum() * 2);
|
||||
}
|
||||
|
||||
void matrix_solver_t::log_stats()
|
||||
{
|
||||
if (this->m_stat_calculations != 0 && this->m_stat_vsolver_calls && log().verbose.is_enabled())
|
||||
|
@ -207,15 +207,42 @@ namespace solver
|
||||
void update_dynamic();
|
||||
|
||||
virtual unsigned vsolve_non_dynamic(const bool newton_raphson) = 0;
|
||||
virtual netlist_time compute_next_timestep(const nl_fptype cur_ts) = 0;
|
||||
|
||||
netlist_time compute_next_timestep(const nl_fptype cur_ts);
|
||||
/* virtual */ void add_term(std::size_t net_idx, terminal_t *term);
|
||||
void add_term(std::size_t net_idx, terminal_t *term);
|
||||
|
||||
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>
|
||||
void store(const T & V);
|
||||
void store(const T & V)
|
||||
{
|
||||
const std::size_t iN = this->m_terms.size();
|
||||
for (std::size_t i = 0; i < iN; i++)
|
||||
this->m_terms[i].setV(V[i]);
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
auto delta(const T & V) -> typename std::decay<decltype(V[0])>::type;
|
||||
auto delta(const T & V) -> typename std::decay<decltype(V[0])>::type
|
||||
{
|
||||
/* NOTE: Ideally we should also include currents (RHS) here. This would
|
||||
* need a reevaluation of the right hand side after voltages have been updated
|
||||
* and thus belong into a different calculation. This applies to all solvers.
|
||||
*/
|
||||
|
||||
const std::size_t iN = this->m_terms.size();
|
||||
using vtype = typename std::decay<decltype(V[0])>::type;
|
||||
vtype cerr = 0;
|
||||
for (std::size_t i = 0; i < iN; i++)
|
||||
cerr = std::max(cerr, std::abs(V[i] - static_cast<vtype>(this->m_terms[i].getV())));
|
||||
return cerr;
|
||||
}
|
||||
|
||||
void set_pointers()
|
||||
{
|
||||
@ -229,7 +256,6 @@ namespace solver
|
||||
max_rail = std::max(max_rail, m_terms[k].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);
|
||||
@ -247,49 +273,6 @@ namespace solver
|
||||
}
|
||||
}
|
||||
|
||||
template <typename FT>
|
||||
void fill_matrix(std::size_t N, FT &RHS)
|
||||
{
|
||||
for (std::size_t k = 0; k < N; k++)
|
||||
{
|
||||
auto &net = m_terms[k];
|
||||
auto **tcr_r = &(m_mat_ptr[k][0]);
|
||||
|
||||
const std::size_t term_count = net.count();
|
||||
const std::size_t railstart = net.railstart();
|
||||
const auto &go = m_gonn[k];
|
||||
const auto > = m_gtn[k];
|
||||
const auto &Idr = m_Idrn[k];
|
||||
const auto &cnV = m_connected_net_Vn[k];
|
||||
|
||||
for (std::size_t i = 0; i < railstart; i++)
|
||||
*tcr_r[i] += go[i];
|
||||
|
||||
typename FT::value_type gtot_t = 0.0;
|
||||
typename FT::value_type RHS_t = 0.0;
|
||||
|
||||
for (std::size_t i = 0; i < term_count; i++)
|
||||
{
|
||||
gtot_t += gt[i];
|
||||
RHS_t += Idr[i];
|
||||
}
|
||||
// FIXME: Code above is faster than vec_sum - Check this
|
||||
#if 0
|
||||
auto gtot_t = plib::vec_sum<FT>(term_count, m_gt);
|
||||
auto RHS_t = plib::vec_sum<FT>(term_count, m_Idr);
|
||||
#endif
|
||||
|
||||
for (std::size_t i = railstart; i < term_count; i++)
|
||||
{
|
||||
RHS_t += (/*m_Idr[i]*/ (- go[i]) * *cnV[i]);
|
||||
}
|
||||
|
||||
RHS[k] = RHS_t;
|
||||
// update diagonal element ...
|
||||
*tcr_r[railstart] += gtot_t; //mat.A[mat.diag[k]] += gtot_t;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
template <typename T, typename M>
|
||||
void log_fill(const T &fill, M &mat)
|
||||
@ -345,61 +328,17 @@ namespace solver
|
||||
}
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void clear_square_mat(std::size_t n, M &m)
|
||||
{
|
||||
for (std::size_t k=0; k < n; k++)
|
||||
{
|
||||
auto *p = &(m[k][0]);
|
||||
for (std::size_t i=0; i < n; i++)
|
||||
p[i] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void build_mat_ptr(std::size_t iN, M &mat)
|
||||
{
|
||||
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 T>
|
||||
using aligned_alloc = plib::aligned_allocator<T, PALIGN_VECTOROPT>;
|
||||
|
||||
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_mat_fptype *, aligned_alloc<nl_mat_fptype *>> m_mat_ptr;
|
||||
plib::pmatrix2d<nl_fptype *, aligned_alloc<nl_fptype *>> m_connected_net_Vn;
|
||||
|
||||
plib::aligned_vector<terms_for_net_t> m_terms;
|
||||
plib::aligned_vector<terms_for_net_t> m_rails_temp;
|
||||
|
||||
/* state - variable time_stepping */
|
||||
plib::aligned_vector<nl_fptype> m_last_V;
|
||||
plib::aligned_vector<nl_fptype> m_DD_n_m_1;
|
||||
plib::aligned_vector<nl_fptype> m_h_n_m_1;
|
||||
|
||||
// FIXME: it should be like this, however dimensions are determined
|
||||
// in vsetup.
|
||||
//state_container<std::vector<nl_fptype>> m_last_V;
|
||||
//state_container<std::vector<nl_fptype>> m_DD_n_m_1;
|
||||
//state_container<std::vector<nl_fptype>> m_h_n_m_1;
|
||||
|
||||
std::vector<unique_pool_ptr<proxied_analog_output_t>> m_inps;
|
||||
|
||||
const solver_parameters_t &m_params;
|
||||
@ -430,29 +369,182 @@ namespace solver
|
||||
std::size_t m_ops;
|
||||
};
|
||||
|
||||
template <typename T>
|
||||
auto matrix_solver_t::delta(const T & V) -> typename std::decay<decltype(V[0])>::type
|
||||
template <typename FT, int SIZE>
|
||||
class matrix_solver_ext_t: public matrix_solver_t
|
||||
{
|
||||
/* 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.
|
||||
*/
|
||||
friend class matrix_solver_t;
|
||||
public:
|
||||
|
||||
const std::size_t iN = this->m_terms.size();
|
||||
using vtype = typename std::decay<decltype(V[0])>::type;
|
||||
vtype cerr = 0;
|
||||
for (std::size_t i = 0; i < iN; i++)
|
||||
cerr = std::max(cerr, std::abs(V[i] - static_cast<vtype>(this->m_terms[i].getV())));
|
||||
return cerr;
|
||||
}
|
||||
using float_type = FT;
|
||||
|
||||
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_terms[i].setV(V[i]);
|
||||
}
|
||||
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_mat_ptr(size, this->max_railstart() + 1)
|
||||
, m_last_V(size, plib::constants<float_type>::zero())
|
||||
, m_DD_n_m_1(size, plib::constants<float_type>::zero())
|
||||
, m_h_n_m_1(size, plib::constants<float_type>::zero())
|
||||
{
|
||||
/*
|
||||
* 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;
|
||||
|
||||
plib::parray2D<float_type *, SIZE, 0> m_mat_ptr;
|
||||
/* state - variable time_stepping */
|
||||
plib::parray<float_type, SIZE> m_last_V;
|
||||
plib::parray<float_type, SIZE> m_DD_n_m_1;
|
||||
plib::parray<float_type, SIZE> m_h_n_m_1;
|
||||
|
||||
// FIXME: it should be like this, however dimensions are determined
|
||||
// in vsetup.
|
||||
//state_container<std::vector<nl_fptype>> m_last_V;
|
||||
//state_container<std::vector<nl_fptype>> m_DD_n_m_1;
|
||||
//state_container<std::vector<nl_fptype>> m_h_n_m_1;
|
||||
|
||||
constexpr std::size_t size() const noexcept { return (SIZE > 0) ? static_cast<std::size_t>(SIZE) : m_dim; }
|
||||
|
||||
netlist_time compute_next_timestep(const nl_fptype cur_ts) override
|
||||
{
|
||||
nl_fptype new_solver_timestep = m_params.m_max_timestep;
|
||||
|
||||
if (m_params.m_dynamic_ts)
|
||||
{
|
||||
for (std::size_t k = 0; k < m_terms.size(); k++)
|
||||
{
|
||||
auto &t = m_terms[k];
|
||||
//const nl_fptype DD_n = (n->Q_Analog() - t->m_last_V);
|
||||
// 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,(t.getV() - m_last_V[k])));
|
||||
const nl_fptype hn = cur_ts;
|
||||
|
||||
//printf("%g %g %g %g\n", DD_n, hn, t.m_DD_n_m_1, t.m_h_n_m_1);
|
||||
nl_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 (std::fabs(DD2) > fp_constants<nl_fptype>::TIMESTEP_MINDIV) // avoid div-by-zero
|
||||
new_net_timestep = std::sqrt(m_params.m_dynamic_lte / std::fabs(plib::constants<nl_fptype>::cast(0.5)*DD2));
|
||||
else
|
||||
new_net_timestep = m_params.m_max_timestep;
|
||||
|
||||
if (new_net_timestep < new_solver_timestep)
|
||||
new_solver_timestep = new_net_timestep;
|
||||
|
||||
m_last_V[k] = t.getV();
|
||||
}
|
||||
if (new_solver_timestep < m_params.m_min_timestep)
|
||||
{
|
||||
new_solver_timestep = m_params.m_min_timestep;
|
||||
}
|
||||
}
|
||||
//if (new_solver_timestep > 10.0 * hn)
|
||||
// new_solver_timestep = 10.0 * hn;
|
||||
/*
|
||||
* FIXME: Factor 2 below is important. Without, we get timing issues. This must be a bug elsewhere.
|
||||
*/
|
||||
return std::max(netlist_time::from_double(new_solver_timestep), netlist_time::quantum() * 2);
|
||||
}
|
||||
|
||||
|
||||
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]);
|
||||
for (std::size_t i=0; i < n; i++)
|
||||
p[i] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
template <typename RT>
|
||||
void fill_matrix(RT &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]);
|
||||
|
||||
const std::size_t term_count = net.count();
|
||||
const std::size_t railstart = net.railstart();
|
||||
const auto &go = m_gonn[k];
|
||||
const auto > = m_gtn[k];
|
||||
const auto &Idr = m_Idrn[k];
|
||||
const auto &cnV = m_connected_net_Vn[k];
|
||||
|
||||
for (std::size_t i = 0; i < railstart; i++)
|
||||
*tcr_r[i] += go[i];
|
||||
|
||||
typename RT::value_type gtot_t = 0.0;
|
||||
typename RT::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;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
} // namespace solver
|
||||
} // namespace netlist
|
||||
|
@ -22,7 +22,7 @@ namespace solver
|
||||
{
|
||||
|
||||
template <typename FT, int SIZE>
|
||||
class matrix_solver_direct_t: public matrix_solver_t
|
||||
class matrix_solver_direct_t: public matrix_solver_ext_t<FT, SIZE>
|
||||
{
|
||||
friend class matrix_solver_t;
|
||||
public:
|
||||
@ -37,7 +37,6 @@ namespace solver
|
||||
|
||||
private:
|
||||
|
||||
const std::size_t m_dim;
|
||||
const std::size_t m_pitch;
|
||||
|
||||
protected:
|
||||
@ -47,8 +46,6 @@ namespace solver
|
||||
unsigned vsolve_non_dynamic(const bool newton_raphson) override;
|
||||
unsigned solve_non_dynamic(const bool newton_raphson);
|
||||
|
||||
constexpr std::size_t size() const noexcept { return (SIZE > 0) ? static_cast<std::size_t>(SIZE) : m_dim; }
|
||||
|
||||
void LE_solve();
|
||||
|
||||
template <typename T>
|
||||
@ -69,15 +66,15 @@ namespace solver
|
||||
template <typename FT, int SIZE>
|
||||
void matrix_solver_direct_t<FT, SIZE>::LE_solve()
|
||||
{
|
||||
const std::size_t kN = size();
|
||||
if (!m_params.m_pivot)
|
||||
const std::size_t kN = this->size();
|
||||
if (!this->m_params.m_pivot)
|
||||
{
|
||||
for (std::size_t i = 0; i < kN; i++)
|
||||
{
|
||||
/* FIXME: Singular matrix? */
|
||||
const FT f = 1.0 / m_A[i][i];
|
||||
const auto &nzrd = m_terms[i].m_nzrd;
|
||||
const auto &nzbd = m_terms[i].m_nzbd;
|
||||
const auto &nzrd = this->m_terms[i].m_nzrd;
|
||||
const auto &nzbd = this->m_terms[i].m_nzbd;
|
||||
|
||||
for (auto &j : nzbd)
|
||||
{
|
||||
@ -138,10 +135,10 @@ namespace solver
|
||||
void matrix_solver_direct_t<FT, SIZE>::LE_back_subst(
|
||||
T & x)
|
||||
{
|
||||
const std::size_t kN = size();
|
||||
const std::size_t kN = this->size();
|
||||
|
||||
/* back substitution */
|
||||
if (m_params.m_pivot)
|
||||
if (this->m_params.m_pivot)
|
||||
{
|
||||
for (std::size_t j = kN; j-- > 0; )
|
||||
{
|
||||
@ -156,7 +153,7 @@ namespace solver
|
||||
for (std::size_t j = kN; j-- > 0; )
|
||||
{
|
||||
FT tmp = 0;
|
||||
const auto &nzrd = m_terms[j].m_nzrd;
|
||||
const auto &nzrd = this->m_terms[j].m_nzrd;
|
||||
const auto e = nzrd.size(); // - 1; /* exclude RHS element */
|
||||
for ( std::size_t k = 0; k < e; k++)
|
||||
tmp += m_A[j][nzrd[k]] * x[nzrd[k]];
|
||||
@ -171,20 +168,17 @@ namespace solver
|
||||
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);
|
||||
const FT err = (newton_raphson ? this->delta(m_new_V) : 0.0);
|
||||
this->store(m_new_V);
|
||||
return (err > this->m_params.m_accuracy) ? 2 : 1;
|
||||
}
|
||||
|
||||
template <typename FT, int SIZE>
|
||||
unsigned matrix_solver_direct_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson)
|
||||
{
|
||||
const std::size_t iN = this->size();
|
||||
|
||||
/* populate matrix */
|
||||
|
||||
this->clear_square_mat(iN, m_A);
|
||||
this->fill_matrix(iN, m_RHS);
|
||||
this->clear_square_mat(m_A);
|
||||
this->fill_matrix(m_RHS);
|
||||
|
||||
this->m_stat_calculations++;
|
||||
return this->solve_non_dynamic(newton_raphson);
|
||||
@ -195,14 +189,13 @@ namespace solver
|
||||
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_pitch(m_pitch_ABS ? m_pitch_ABS : (((m_dim + 0) + 7) / 8) * 8)
|
||||
: matrix_solver_ext_t<FT, SIZE>(anetlist, name, nets, params, size)
|
||||
, m_pitch(m_pitch_ABS ? m_pitch_ABS : (((size + 0) + 7) / 8) * 8)
|
||||
, m_new_V(size)
|
||||
, m_A(size, m_pitch)
|
||||
, m_RHS(size)
|
||||
{
|
||||
this->build_mat_ptr(size, m_A);
|
||||
this->build_mat_ptr(m_A);
|
||||
}
|
||||
|
||||
} // namespace solver
|
||||
|
@ -34,8 +34,8 @@ namespace solver
|
||||
// ----------------------------------------------------------------------------------------
|
||||
unsigned vsolve_non_dynamic(const bool newton_raphson) override
|
||||
{
|
||||
this->clear_square_mat(1, this->m_A);
|
||||
this->fill_matrix(1, this->m_RHS);
|
||||
this->clear_square_mat(this->m_A);
|
||||
this->fill_matrix(this->m_RHS);
|
||||
|
||||
std::array<FT, 1> new_V = { this->m_RHS[0] / this->m_A[0][0] };
|
||||
|
||||
|
@ -34,8 +34,8 @@ namespace solver
|
||||
{}
|
||||
unsigned vsolve_non_dynamic(const bool newton_raphson) override
|
||||
{
|
||||
this->clear_square_mat(2, this->m_A);
|
||||
this->fill_matrix(2, this->m_RHS);
|
||||
this->clear_square_mat(this->m_A);
|
||||
this->fill_matrix(this->m_RHS);
|
||||
|
||||
const float_type a = this->m_A[0][0];
|
||||
const float_type b = this->m_A[0][1];
|
||||
|
@ -26,25 +26,22 @@ namespace solver
|
||||
{
|
||||
|
||||
template <typename FT, int SIZE>
|
||||
class matrix_solver_GCR_t: public matrix_solver_t
|
||||
class matrix_solver_GCR_t: public matrix_solver_ext_t<FT, SIZE>
|
||||
{
|
||||
public:
|
||||
|
||||
using mat_type = plib::pGEmatrix_cr_t<plib::pmatrix_cr_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 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)
|
||||
: matrix_solver_ext_t<FT, SIZE>(anetlist, name, nets, params, size)
|
||||
, RHS(size)
|
||||
, new_V(size)
|
||||
, mat(static_cast<typename mat_type::index_type>(size))
|
||||
, m_proc()
|
||||
{
|
||||
const std::size_t iN = this->N();
|
||||
const std::size_t iN = this->size();
|
||||
|
||||
/* build the final matrix */
|
||||
|
||||
@ -65,7 +62,7 @@ namespace solver
|
||||
|
||||
auto gr = mat.gaussian_extend_fill_mat(fill);
|
||||
|
||||
log_fill(fill, mat);
|
||||
this->log_fill(fill, mat);
|
||||
|
||||
mat.build_from_fill_mat(fill);
|
||||
|
||||
@ -79,13 +76,13 @@ namespace solver
|
||||
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];
|
||||
this->m_mat_ptr[k][j] = &mat.A[i];
|
||||
cnt++;
|
||||
break;
|
||||
}
|
||||
}
|
||||
nl_assert(cnt == this->m_terms[k].railstart());
|
||||
m_mat_ptr[k][this->m_terms[k].railstart()] = &mat.A[mat.diag[k]];
|
||||
this->m_mat_ptr[k][this->m_terms[k].railstart()] = &mat.A[mat.diag[k]];
|
||||
}
|
||||
|
||||
this->log().verbose("maximum fill: {1}", gr.first);
|
||||
@ -96,7 +93,7 @@ namespace solver
|
||||
|
||||
// FIXME: Move me
|
||||
|
||||
if (state().lib().isLoaded())
|
||||
if (this->state().lib().isLoaded())
|
||||
{
|
||||
pstring symname = static_compile_name();
|
||||
m_proc.load(this->state().lib(), symname);
|
||||
@ -111,8 +108,6 @@ namespace solver
|
||||
}
|
||||
}
|
||||
|
||||
constexpr std::size_t N() const { return m_dim; }
|
||||
|
||||
unsigned vsolve_non_dynamic(const bool newton_raphson) override;
|
||||
|
||||
std::pair<pstring, pstring> create_solver_code() override;
|
||||
@ -125,7 +120,6 @@ namespace solver
|
||||
|
||||
pstring static_compile_name();
|
||||
|
||||
const std::size_t m_dim;
|
||||
plib::parray<FT, SIZE> RHS;
|
||||
plib::parray<FT, SIZE> new_V;
|
||||
|
||||
@ -142,7 +136,7 @@ namespace solver
|
||||
template <typename FT, int SIZE>
|
||||
void matrix_solver_GCR_t<FT, SIZE>::generate_code(plib::putf8_fmt_writer &strm)
|
||||
{
|
||||
const std::size_t iN = N();
|
||||
const std::size_t iN = this->size();
|
||||
|
||||
for (std::size_t i = 0; i < mat.nz_num; i++)
|
||||
strm("double m_A{1} = m_A[{2}];\n", i, i);
|
||||
@ -232,13 +226,9 @@ namespace solver
|
||||
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, RHS);
|
||||
mat.set_scalar(0.0);
|
||||
this->fill_matrix(RHS);
|
||||
|
||||
/* now solve it */
|
||||
|
||||
@ -257,8 +247,8 @@ namespace solver
|
||||
|
||||
this->m_stat_calculations++;
|
||||
|
||||
const FT err = (newton_raphson ? delta(new_V) : 0.0);
|
||||
store(new_V);
|
||||
const FT err = (newton_raphson ? this->delta(new_V) : 0.0);
|
||||
this->store(new_V);
|
||||
return (err > this->m_params.m_accuracy) ? 2 : 1;
|
||||
}
|
||||
|
||||
|
@ -106,7 +106,7 @@ namespace solver
|
||||
m_ops.m_mat.set_scalar(0.0);
|
||||
|
||||
/* populate matrix and V for first estimate */
|
||||
this->fill_matrix(iN, RHS);
|
||||
this->fill_matrix(RHS);
|
||||
|
||||
for (std::size_t k = 0; k < iN; k++)
|
||||
{
|
||||
|
@ -45,7 +45,7 @@ namespace solver
|
||||
{
|
||||
|
||||
template <typename FT, int SIZE>
|
||||
class matrix_solver_sm_t: public matrix_solver_t
|
||||
class matrix_solver_sm_t: public matrix_solver_ext_t<FT, SIZE>
|
||||
{
|
||||
friend class matrix_solver_t;
|
||||
|
||||
@ -59,11 +59,10 @@ namespace solver
|
||||
matrix_solver_sm_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)
|
||||
: matrix_solver_ext_t<FT, SIZE>(anetlist, name, nets, params, size)
|
||||
, m_cnt(0)
|
||||
{
|
||||
this->build_mat_ptr(this->size(), m_A);
|
||||
this->build_mat_ptr(m_A);
|
||||
}
|
||||
|
||||
void reset() override { matrix_solver_t::reset(); }
|
||||
@ -72,8 +71,6 @@ namespace solver
|
||||
unsigned vsolve_non_dynamic(const bool newton_raphson) override;
|
||||
unsigned solve_non_dynamic(const bool newton_raphson);
|
||||
|
||||
constexpr std::size_t size() const { return m_dim; }
|
||||
|
||||
void LE_invert();
|
||||
|
||||
template <typename T>
|
||||
@ -109,7 +106,6 @@ namespace solver
|
||||
|
||||
//float_ext_type m_RHSx[storage_N];
|
||||
|
||||
const std::size_t m_dim;
|
||||
std::size_t m_cnt;
|
||||
|
||||
};
|
||||
@ -121,7 +117,7 @@ namespace solver
|
||||
template <typename FT, int SIZE>
|
||||
void matrix_solver_sm_t<FT, SIZE>::LE_invert()
|
||||
{
|
||||
const std::size_t kN = size();
|
||||
const std::size_t kN = this->size();
|
||||
|
||||
for (std::size_t i = 0; i < kN; i++)
|
||||
{
|
||||
@ -137,13 +133,13 @@ namespace solver
|
||||
{
|
||||
/* 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();
|
||||
const auto * const p = this->m_terms[i].m_nzrd.data();
|
||||
const std::size_t e = this->m_terms[i].m_nzrd.size();
|
||||
|
||||
/* Eliminate column i from row j */
|
||||
|
||||
const auto * const pb = m_terms[i].m_nzbd.data();
|
||||
const std::size_t eb = m_terms[i].m_nzbd.size();
|
||||
const auto * const pb = this->m_terms[i].m_nzbd.data();
|
||||
const std::size_t eb = this->m_terms[i].m_nzbd.size();
|
||||
for (std::size_t jb = 0; jb < eb; jb++)
|
||||
{
|
||||
const unsigned j = pb[jb];
|
||||
@ -186,7 +182,7 @@ namespace solver
|
||||
void matrix_solver_sm_t<FT, SIZE>::LE_compute_x(
|
||||
T & x)
|
||||
{
|
||||
const std::size_t kN = size();
|
||||
const std::size_t kN = this->size();
|
||||
|
||||
for (std::size_t i=0; i<kN; i++)
|
||||
x[i] = 0.0;
|
||||
@ -204,7 +200,7 @@ namespace solver
|
||||
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();
|
||||
const std::size_t iN = this->size();
|
||||
|
||||
// NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init)
|
||||
std::array<float_type, storage_N> new_V;
|
||||
@ -232,7 +228,7 @@ namespace solver
|
||||
{
|
||||
std::size_t colcount = 0;
|
||||
|
||||
auto &nz = m_terms[row].m_nz;
|
||||
auto &nz = this->m_terms[row].m_nz;
|
||||
for (unsigned & col : nz)
|
||||
{
|
||||
v[col] = A(row,col) - lA(row,col);
|
||||
@ -279,8 +275,8 @@ namespace solver
|
||||
|
||||
this->LE_compute_x(new_V);
|
||||
|
||||
const float_type err = (newton_raphson ? delta(new_V) : 0.0);
|
||||
store(new_V);
|
||||
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;
|
||||
}
|
||||
|
||||
@ -288,9 +284,8 @@ namespace solver
|
||||
unsigned matrix_solver_sm_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson)
|
||||
{
|
||||
|
||||
const std::size_t iN = this->size();
|
||||
this->clear_square_mat(iN, this->m_A);
|
||||
this->fill_matrix(iN, this->m_RHS);
|
||||
this->clear_square_mat(this->m_A);
|
||||
this->fill_matrix(this->m_RHS);
|
||||
|
||||
this->m_stat_calculations++;
|
||||
return this->solve_non_dynamic(newton_raphson);
|
||||
|
@ -114,11 +114,10 @@ namespace solver
|
||||
* Need something like that for gaussian elimination as well.
|
||||
*/
|
||||
|
||||
|
||||
const std::size_t iN = this->size();
|
||||
|
||||
this->clear_square_mat(iN, this->m_A);
|
||||
this->fill_matrix(iN, this->m_RHS);
|
||||
this->clear_square_mat(this->m_A);
|
||||
this->fill_matrix(this->m_RHS);
|
||||
|
||||
bool resched = false;
|
||||
|
||||
|
@ -52,7 +52,7 @@ namespace solver
|
||||
{
|
||||
|
||||
template <typename FT, int SIZE>
|
||||
class matrix_solver_w_t: public matrix_solver_t
|
||||
class matrix_solver_w_t: public matrix_solver_ext_t<FT, SIZE>
|
||||
{
|
||||
friend class matrix_solver_t;
|
||||
|
||||
@ -66,11 +66,10 @@ namespace solver
|
||||
matrix_solver_w_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)
|
||||
: matrix_solver_ext_t<FT, SIZE>(anetlist, name, nets, params, size)
|
||||
, m_cnt(0)
|
||||
, m_dim(size)
|
||||
{
|
||||
this->build_mat_ptr(this->size(), m_A);
|
||||
this->build_mat_ptr(m_A);
|
||||
}
|
||||
|
||||
void reset() override { matrix_solver_t::reset(); }
|
||||
@ -79,8 +78,6 @@ namespace solver
|
||||
unsigned vsolve_non_dynamic(const bool newton_raphson) override;
|
||||
unsigned solve_non_dynamic(const bool newton_raphson);
|
||||
|
||||
constexpr std::size_t size() const { return m_dim; }
|
||||
|
||||
void LE_invert();
|
||||
|
||||
template <typename T>
|
||||
@ -121,11 +118,6 @@ namespace solver
|
||||
std::array<unsigned, storage_N> colcount;
|
||||
|
||||
unsigned m_cnt;
|
||||
|
||||
//float_ext_type m_RHSx[storage_N];
|
||||
|
||||
const std::size_t m_dim;
|
||||
|
||||
};
|
||||
|
||||
// ----------------------------------------------------------------------------------------
|
||||
@ -135,7 +127,7 @@ namespace solver
|
||||
template <typename FT, int SIZE>
|
||||
void matrix_solver_w_t<FT, SIZE>::LE_invert()
|
||||
{
|
||||
const std::size_t kN = size();
|
||||
const std::size_t kN = this->size();
|
||||
|
||||
for (std::size_t i = 0; i < kN; i++)
|
||||
{
|
||||
@ -151,13 +143,13 @@ namespace solver
|
||||
{
|
||||
/* FIXME: Singular matrix? */
|
||||
const float_type f = 1.0 / W(i,i);
|
||||
const auto * const p = m_terms[i].m_nzrd.data();
|
||||
const size_t e = m_terms[i].m_nzrd.size();
|
||||
const auto * const p = this->m_terms[i].m_nzrd.data();
|
||||
const size_t e = this->m_terms[i].m_nzrd.size();
|
||||
|
||||
/* Eliminate column i from row j */
|
||||
|
||||
const auto * const pb = m_terms[i].m_nzbd.data();
|
||||
const size_t eb = m_terms[i].m_nzbd.size();
|
||||
const auto * const pb = this->m_terms[i].m_nzbd.data();
|
||||
const size_t eb = this->m_terms[i].m_nzbd.size();
|
||||
for (std::size_t jb = 0; jb < eb; jb++)
|
||||
{
|
||||
const auto j = pb[jb];
|
||||
@ -199,7 +191,7 @@ namespace solver
|
||||
void matrix_solver_w_t<FT, SIZE>::LE_compute_x(
|
||||
T & x)
|
||||
{
|
||||
const std::size_t kN = size();
|
||||
const std::size_t kN = this->size();
|
||||
|
||||
for (std::size_t i=0; i<kN; i++)
|
||||
x[i] = 0.0;
|
||||
@ -217,7 +209,7 @@ namespace solver
|
||||
template <typename FT, int SIZE>
|
||||
unsigned matrix_solver_w_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphson)
|
||||
{
|
||||
const auto iN = size();
|
||||
const auto iN = this->size();
|
||||
|
||||
// NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init)
|
||||
std::array<float_type, storage_N> new_V;
|
||||
@ -245,7 +237,7 @@ namespace solver
|
||||
for (unsigned row = 0; row < iN; row ++)
|
||||
{
|
||||
unsigned cc=0;
|
||||
auto &nz = m_terms[row].m_nz;
|
||||
auto &nz = this->m_terms[row].m_nz;
|
||||
for (auto & col : nz)
|
||||
{
|
||||
if (A(row,col) != lA(row,col))
|
||||
@ -349,17 +341,16 @@ namespace solver
|
||||
plib::perrlogger("{} failed on row {}: {} RHS: {}\n", this->name(), i, std::abs(tmp-RHS(i)), RHS(i));
|
||||
}
|
||||
|
||||
const float_type err = (newton_raphson ? delta(new_V) : 0.0);
|
||||
store(new_V);
|
||||
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;
|
||||
}
|
||||
|
||||
template <typename FT, int SIZE>
|
||||
unsigned matrix_solver_w_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson)
|
||||
{
|
||||
const std::size_t iN = this->size();
|
||||
this->clear_square_mat(iN, this->m_A);
|
||||
this->fill_matrix(iN, this->m_RHS);
|
||||
this->clear_square_mat(this->m_A);
|
||||
this->fill_matrix(this->m_RHS);
|
||||
|
||||
this->m_stat_calculations++;
|
||||
return this->solve_non_dynamic(newton_raphson);
|
||||
|
Loading…
Reference in New Issue
Block a user