Sprinter/mame
2
0
Fork 0
mirror of https://github.com/holub/mame synced 2025-04-23 08:49:55 +03:00
Code Issues Actions 1 Packages Projects Releases Wiki Activity

netlist: move solver stuff into separate namespace. (nw)

Browse Source 
- new namespace "solver"
- minor ptime modifications
This commit is contained in:
couriersud 2019-10-31 01:15:43 +01:00
parent 9582866bbf
commit 65fb297023
15 changed files with 499 additions and 497 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
src/lib/netlist/nl_base.h
View File

@ -186,7 +186,6 @@ namespace netlist
namespace devices
{
class matrix_solver_t;
class NETLIB_NAME(solver);
class NETLIB_NAME(mainclock);
class NETLIB_NAME(base_proxy);
@ -194,6 +193,11 @@ namespace netlist
class NETLIB_NAME(base_a_to_d_proxy);
} // namespace devices
namespace solver
{
class matrix_solver_t;
} // namespace solver
class logic_output_t;
class logic_input_t;
class analog_net_t;
@ -902,12 +906,12 @@ namespace netlist
nl_double *Q_Analog_state_ptr() NL_NOEXCEPT { return m_cur_Analog.ptr(); }
//FIXME: needed by current solver code
devices::matrix_solver_t *solver() const noexcept { return m_solver; }
void set_solver(devices::matrix_solver_t *solver) noexcept { m_solver = solver; }
solver::matrix_solver_t *solver() const noexcept { return m_solver; }
void set_solver(solver::matrix_solver_t *solver) noexcept { m_solver = solver; }
private:
state_var<nl_double> m_cur_Analog;
devices::matrix_solver_t *m_solver;
solver::matrix_solver_t *m_solver;
};
// -----------------------------------------------------------------------------

24
src/lib/netlist/plib/ptime.h
View File

@ -58,7 +58,7 @@ namespace plib
return ptime(lhs.m_time + rhs.m_time);
}
friend constexpr const ptime operator*(ptime lhs, const mult_type &factor) noexcept
friend constexpr const ptime operator*(ptime lhs, const mult_type factor) noexcept
{
return ptime(lhs.m_time * factor);
}
@ -107,18 +107,18 @@ namespace plib
// for save states ....
C14CONSTEXPR internal_type *get_internaltype_ptr() noexcept { return &m_time; }
static constexpr ptime from_nsec(const internal_type ns) noexcept { return ptime(ns, UINT64_C(1000000000)); }
static constexpr ptime from_usec(const internal_type us) noexcept { return ptime(us, UINT64_C( 1000000)); }
static constexpr ptime from_msec(const internal_type ms) noexcept { return ptime(ms, UINT64_C( 1000)); }
static constexpr ptime from_sec(const internal_type s) noexcept { return ptime(s, UINT64_C( 1)); }
static constexpr ptime from_hz(const internal_type hz) noexcept { return ptime(1 , hz); }
static constexpr ptime from_raw(const internal_type raw) noexcept { return ptime(raw); }
static constexpr ptime from_double(const double t) noexcept { return ptime(static_cast<internal_type>(std::floor(t * static_cast<double>(RES) + 0.5)), RES); }
static constexpr const ptime from_nsec(const internal_type ns) noexcept { return ptime(ns, UINT64_C(1000000000)); }
static constexpr const ptime from_usec(const internal_type us) noexcept { return ptime(us, UINT64_C( 1000000)); }
static constexpr const ptime from_msec(const internal_type ms) noexcept { return ptime(ms, UINT64_C( 1000)); }
static constexpr const ptime from_sec(const internal_type s) noexcept { return ptime(s, UINT64_C( 1)); }
static constexpr const ptime from_hz(const internal_type hz) noexcept { return ptime(1 , hz); }
static constexpr const ptime from_raw(const internal_type raw) noexcept { return ptime(raw); }
static constexpr const ptime from_double(const double t) noexcept { return ptime(static_cast<internal_type>(std::floor(t * static_cast<double>(RES) + 0.5)), RES); }
static constexpr ptime zero() noexcept { return ptime(0, RES); }
static constexpr ptime quantum() noexcept { return ptime(1, RES); }
static constexpr ptime never() noexcept { return ptime(plib::numeric_limits<internal_type>::max(), RES); }
static constexpr internal_type resolution() noexcept { return RES; }
static constexpr const ptime zero() noexcept { return ptime(0, RES); }
static constexpr const ptime quantum() noexcept { return ptime(1, RES); }
static constexpr const ptime never() noexcept { return ptime(plib::numeric_limits<internal_type>::max(), RES); }
static constexpr const internal_type resolution() noexcept { return RES; }
constexpr internal_type in_nsec() const noexcept { return m_time / (RES / UINT64_C(1000000000)); }
constexpr internal_type in_usec() const noexcept { return m_time / (RES / UINT64_C( 1000000)); }

38
src/lib/netlist/solver/nld_matrix_solver.cpp
View File

@ -12,7 +12,7 @@
namespace netlist
{
namespace devices
namespace solver
{
terms_for_net_t::terms_for_net_t(analog_net_t * net)
@ -58,11 +58,13 @@ namespace devices
{
connect_post_start(m_fb_sync, m_Q_sync);
setup_base(nets);
/* now setup the matrix */
setup_matrix();
}
void matrix_solver_t::setup_base(const analog_net_t::list_t &nets)
{
log().debug("New solver setup\n");
m_terms.clear();
@ -130,9 +132,6 @@ namespace devices
}
log().debug("added net with {1} populated connections\n", net->core_terms().size());
}
/* now setup the matrix */
setup_matrix();
}
void matrix_solver_t::sort_terms(matrix_sort_type_e sort)
@ -417,8 +416,19 @@ namespace devices
d->timestep(dd);
}
void matrix_solver_t::solve_base()
const netlist_time matrix_solver_t::solve(netlist_time now)
{
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();
/* update all terminals for new time step */
m_last_step = now;
step(delta);
++m_stat_vsolver_calls;
if (has_dynamic_devices())
{
@ -444,21 +454,7 @@ namespace devices
{
this->vsolve_non_dynamic(false);
}
}
const netlist_time matrix_solver_t::solve(netlist_time now)
{
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();
/* 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;
@ -622,6 +618,6 @@ namespace devices
}
}
} // namespace devices
} // namespace solver
} // namespace netlist

10
src/lib/netlist/solver/nld_matrix_solver.h
View File

@ -20,7 +20,7 @@
namespace netlist
{
namespace devices
namespace solver
{
P_ENUM(matrix_sort_type_e,
NOSORT,
@ -162,8 +162,6 @@ namespace devices
public:
using list_t = std::vector<matrix_solver_t *>;
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.
*/
@ -481,17 +479,17 @@ namespace devices
const float_type * const * other_cur_analog = m_connected_net_Vn[k];
for (std::size_t i = 0; i < terms_count; i++)
rhsk_a = rhsk_a + Idr[i];
rhsk_a += Idr[i];
for (std::size_t i = m_terms[k].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];
rhsk_b += - go[i] * *other_cur_analog[i];
child.RHS(k) = rhsk_a + rhsk_b;
}
}
} //namespace devices
} // namespace solver
} // namespace netlist
#endif /* NLD_MS_DIRECT_H_ */

12
src/lib/netlist/solver/nld_ms_direct.h
View File

@ -18,7 +18,7 @@
namespace netlist
{
namespace devices
namespace solver
{
template <typename FT, int SIZE>
@ -51,6 +51,8 @@ namespace devices
const FT &RHS(std::size_t r) const noexcept { return m_A[r][size()]; }
FT &RHS(std::size_t r) noexcept { return m_A[r][size()]; }
PALIGNAS_VECTOROPT()
plib::parray<FT, SIZE> m_new_V;
private:
@ -59,6 +61,8 @@ namespace devices
const std::size_t m_dim;
const std::size_t m_pitch;
PALIGNAS_VECTOROPT()
plib::parray2D<FT, SIZE, m_pitch_ABS> m_A;
};
@ -79,10 +83,10 @@ namespace devices
const auto &nzrd = m_terms[i].m_nzrd;
const auto &nzbd = m_terms[i].m_nzbd;
for (auto j : nzbd)
for (auto &j : nzbd)
{
const FT f1 = -f * A(j, i);
for (auto k : nzrd)
for (auto &k : nzrd)
A(j, k) += A(i, k) * f1;
//RHS(j) += RHS(i) * f1;
}
@ -215,7 +219,7 @@ namespace devices
state().save(*this, RHS(k), this->name(), plib::pfmt("RHS.{1}")(k));
}
} // namespace devices
} // namespace solver
} // namespace netlist
#endif /* NLD_MS_DIRECT_H_ */

4
src/lib/netlist/solver/nld_ms_direct1.h
View File

@ -13,7 +13,7 @@
namespace netlist
{
namespace devices
namespace solver
{
template <typename FT>
class matrix_solver_direct1_t: public matrix_solver_direct_t<FT, 1>
@ -49,7 +49,7 @@ namespace devices
} //namespace devices
} // namespace solver
} // namespace netlist

4
src/lib/netlist/solver/nld_ms_direct2.h
View File

@ -13,7 +13,7 @@
namespace netlist
{
namespace devices
namespace solver
{
// ----------------------------------------------------------------------------------------
@ -54,7 +54,7 @@ namespace devices
};
} //namespace devices
} // namespace solver
} // namespace netlist
#endif /* NLD_MS_DIRECT2_H_ */

4
src/lib/netlist/solver/nld_ms_gcr.h
View File

@ -22,7 +22,7 @@
namespace netlist
{
namespace devices
namespace solver
{
template <typename FT, int SIZE>
@ -262,7 +262,7 @@ namespace devices
return (err > this->m_params.m_accuracy) ? 2 : 1;
}
} // namespace devices
} // namespace solver
} // namespace netlist
#endif /* NLD_MS_GCR_H_ */

4
src/lib/netlist/solver/nld_ms_gmres.h
View File

@ -21,7 +21,7 @@
namespace netlist
{
namespace devices
namespace solver
{
template <typename FT, int SIZE>
@ -135,7 +135,7 @@ namespace devices
} // namespace devices
} // namespace solver
} // namespace netlist
#endif /* NLD_MS_GMRES_H_ */

4
src/lib/netlist/solver/nld_ms_sm.h
View File

@ -41,7 +41,7 @@
namespace netlist
{
namespace devices
namespace solver
{
template <typename FT, int SIZE>
@ -297,7 +297,7 @@ namespace devices
}
} // namespace devices
} // namespace solver
} // namespace netlist
#endif /* NLD_MS_DIRECT_H_ */

246
src/lib/netlist/solver/nld_ms_sor.h
View File

@ -19,151 +19,151 @@
namespace netlist
{
namespace devices
namespace solver
{
template <typename FT, int SIZE>
class matrix_solver_SOR_t: public matrix_solver_direct_t<FT, SIZE>
{
public:
using float_type = FT;
matrix_solver_SOR_t(netlist_state_t &anetlist, const pstring &name,
analog_net_t::list_t &nets,
const solver_parameters_t *params, const std::size_t size)
: matrix_solver_direct_t<FT, SIZE>(anetlist, name, nets, 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)
{
}
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;
};
// ----------------------------------------------------------------------------------------
// matrix_solver - Gauss - Seidel
// ----------------------------------------------------------------------------------------
template <typename FT, int SIZE>
unsigned matrix_solver_SOR_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson)
{
const std::size_t iN = this->size();
bool resched = false;
unsigned resched_cnt = 0;
/* ideally, we could get an estimate for the spectral radius of
* Inv(D - L) * U
*
* and estimate using
*
* omega = 2.0 / (1.0 + std::sqrt(1-rho))
*/
const float_type ws = this->m_params.m_gs_sor;
for (std::size_t k = 0; k < iN; k++)
template <typename FT, int SIZE>
class matrix_solver_SOR_t: public matrix_solver_direct_t<FT, SIZE>
{
float_type gtot_t = 0.0;
float_type gabs_t = 0.0;
float_type RHS_t = 0.0;
public:
const std::size_t term_count = this->m_terms[k].count();
const float_type * const gt = this->m_gtn[k];
const float_type * const go = this->m_gonn[k];
const float_type * const Idr = this->m_Idrn[k];
auto other_cur_analog = this->m_connected_net_Vn[k];
using float_type = FT;
this->m_new_V[k] = this->m_terms[k].getV();
matrix_solver_SOR_t(netlist_state_t &anetlist, const pstring &name,
analog_net_t::list_t &nets,
const solver_parameters_t *params, const std::size_t size)
: matrix_solver_direct_t<FT, SIZE>(anetlist, name, nets, 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)
{
}
for (std::size_t i = 0; i < term_count; i++)
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;
};
// ----------------------------------------------------------------------------------------
// matrix_solver - Gauss - Seidel
// ----------------------------------------------------------------------------------------
template <typename FT, int SIZE>
unsigned matrix_solver_SOR_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson)
{
const std::size_t iN = this->size();
bool resched = false;
unsigned resched_cnt = 0;
/* ideally, we could get an estimate for the spectral radius of
* Inv(D - L) * U
*
* and estimate using
*
* omega = 2.0 / (1.0 + std::sqrt(1-rho))
*/
const float_type ws = this->m_params.m_gs_sor;
for (std::size_t k = 0; k < iN; k++)
{
gtot_t = gtot_t + gt[i];
RHS_t = RHS_t + Idr[i];
}
float_type gtot_t = 0.0;
float_type gabs_t = 0.0;
float_type RHS_t = 0.0;
for (std::size_t i = this->m_terms[k].railstart(); i < term_count; i++)
RHS_t = RHS_t - go[i] * *other_cur_analog[i];
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];
RHS[k] = RHS_t;
this->m_new_V[k] = this->m_terms[k].getV();
if (this->m_params.m_use_gabs)
{
for (std::size_t i = 0; i < term_count; i++)
gabs_t = gabs_t + std::abs(go[i]);
{
gtot_t = gtot_t + gt[i];
RHS_t = RHS_t + Idr[i];
}
gabs_t *= plib::constants<nl_double>::cast(0.5); // derived by try and error
if (gabs_t <= gtot_t)
for (std::size_t i = this->m_terms[k].railstart(); i < term_count; i++)
RHS_t = RHS_t - go[i] * *other_cur_analog[i];
RHS[k] = RHS_t;
if (this->m_params.m_use_gabs)
{
for (std::size_t i = 0; i < term_count; i++)
gabs_t = gabs_t + std::abs(go[i]);
gabs_t *= plib::constants<nl_double>::cast(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;
}
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);
}
}
else
{
w[k] = ws / gtot_t;
one_m_w[k] = plib::constants<FT>::one() - ws;
}
else
}
const float_type accuracy = this->m_params.m_accuracy;
do {
resched = false;
float_type err = 0;
for (std::size_t k = 0; k < iN; k++)
{
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);
const int * net_other = this->m_terms[k].m_connected_net_idx.data();
const std::size_t railstart = this->m_terms[k].railstart();
const float_type * go = this->m_gonn[k];
float_type Idrive = 0.0;
for (std::size_t i = 0; i < railstart; i++)
Idrive = Idrive - go[i] * this->m_new_V[static_cast<std::size_t>(net_other[i])];
const float_type new_val = this->m_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;
}
}
else
if (err > accuracy)
resched = true;
resched_cnt++;
} while (resched && (resched_cnt < this->m_params.m_gs_loops));
this->m_iterative_total += resched_cnt;
this->m_stat_calculations++;
if (resched)
{
w[k] = ws / gtot_t;
one_m_w[k] = plib::constants<FT>::one() - ws;
// Fallback to direct solver ...
this->m_iterative_fail++;
return matrix_solver_direct_t<FT, SIZE>::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;
}
const float_type accuracy = this->m_params.m_accuracy;
do {
resched = false;
float_type err = 0;
for (std::size_t k = 0; k < iN; k++)
{
const int * net_other = this->m_terms[k].m_connected_net_idx.data();
const std::size_t railstart = this->m_terms[k].railstart();
const float_type * go = this->m_gonn[k];
float_type Idrive = 0.0;
for (std::size_t i = 0; i < railstart; i++)
Idrive = Idrive - go[i] * this->m_new_V[static_cast<std::size_t>(net_other[i])];
const float_type new_val = this->m_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;
}
if (err > accuracy)
resched = true;
resched_cnt++;
} while (resched && (resched_cnt < this->m_params.m_gs_loops));
this->m_iterative_total += resched_cnt;
this->m_stat_calculations++;
if (resched)
{
// Fallback to direct solver ...
this->m_iterative_fail++;
return matrix_solver_direct_t<FT, SIZE>::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;
}
} //namespace devices
} // namespace solver
} // namespace netlist
#endif /* NLD_MS_SOR_H_ */

4
src/lib/netlist/solver/nld_ms_sor_mat.h
View File

@ -20,7 +20,7 @@
namespace netlist
{
namespace devices
namespace solver
{
template <typename FT, int SIZE>
@ -218,7 +218,7 @@ namespace devices
}
} // namespace devices
} // namespace solver
} // namespace netlist
#endif /* NLD_MS_GAUSS_SEIDEL_H_ */

574
src/lib/netlist/solver/nld_ms_w.h
View File

@ -48,325 +48,325 @@
namespace netlist
{
namespace devices
{
template <typename FT, int SIZE>
class matrix_solver_w_t: public matrix_solver_t
namespace solver
{
friend class matrix_solver_t;
public:
using float_ext_type = FT;
using float_type = FT;
// FIXME: dirty hack to make this compile
static constexpr const std::size_t storage_N = 100;
matrix_solver_w_t(netlist_state_t &anetlist, const pstring &name,
const analog_net_t::list_t &nets,
const solver_parameters_t *params, const std::size_t size)
: matrix_solver_t(anetlist, name, nets, params)
, m_cnt(0)
, m_dim(size)
template <typename FT, int SIZE>
class matrix_solver_w_t: public matrix_solver_t
{
// FIXME: This shouldn't be necessary, recalculate on each entry ...
for (std::size_t k = 0; k < this->size(); k++)
state().save(*this, RHS(k), this->name(), plib::pfmt("RHS.{1}")(k));
friend class matrix_solver_t;
public:
using float_ext_type = FT;
using float_type = FT;
// FIXME: dirty hack to make this compile
static constexpr const std::size_t storage_N = 100;
matrix_solver_w_t(netlist_state_t &anetlist, const pstring &name,
const analog_net_t::list_t &nets,
const solver_parameters_t *params, const std::size_t size)
: matrix_solver_t(anetlist, name, nets, params)
, m_cnt(0)
, m_dim(size)
{
// FIXME: This shouldn't be necessary, recalculate on each entry ...
for (std::size_t k = 0; k < this->size(); k++)
state().save(*this, RHS(k), this->name(), plib::pfmt("RHS.{1}")(k));
}
void reset() override { matrix_solver_t::reset(); }
protected:
unsigned vsolve_non_dynamic(const bool newton_raphson) override;
unsigned solve_non_dynamic(const bool newton_raphson);
constexpr std::size_t size() const { return m_dim; }
void LE_invert();
template <typename T>
void LE_compute_x(T & x);
template <typename T1, typename T2>
float_ext_type &A(const T1 &r, const T2 &c) { return m_A[r][c]; }
template <typename T1, typename T2>
float_ext_type &W(const T1 &r, const T2 &c) { return m_W[r][c]; }
/* access to Ainv for fixed columns over row, there store transposed */
template <typename T1, typename T2>
float_ext_type &Ainv(const T1 &r, const T2 &c) { return m_Ainv[c][r]; }
template <typename T1>
float_ext_type &RHS(const T1 &r) { return m_RHS[r]; }
template <typename T1, typename T2>
float_ext_type &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; }
private:
template <typename T, std::size_t N, std::size_t M>
using array2D = std::array<std::array<T, M>, N>;
static constexpr std::size_t m_pitch = ((( storage_N) + 7) / 8) * 8;
array2D<float_ext_type, storage_N, m_pitch> m_A;
array2D<float_ext_type, storage_N, m_pitch> m_Ainv;
array2D<float_ext_type, storage_N, m_pitch> m_W;
std::array<float_ext_type, storage_N> m_RHS; // right hand side - contains currents
array2D<float_ext_type, storage_N, m_pitch> m_lA;
/* temporary */
array2D<float_ext_type, storage_N, m_pitch> H;
std::array<unsigned, storage_N> rows;
array2D<unsigned, storage_N, m_pitch> cols;
std::array<unsigned, storage_N> colcount;
unsigned m_cnt;
//float_ext_type m_RHSx[storage_N];
const std::size_t m_dim;
};
// ----------------------------------------------------------------------------------------
// matrix_solver_direct
// ----------------------------------------------------------------------------------------
template <typename FT, int SIZE>
void matrix_solver_w_t<FT, SIZE>::LE_invert()
{
const std::size_t kN = size();
for (std::size_t i = 0; i < kN; i++)
{
for (std::size_t j = 0; j < kN; j++)
{
W(i,j) = lA(i,j) = A(i,j);
Ainv(i,j) = 0.0;
}
Ainv(i,i) = 1.0;
}
/* down */
for (std::size_t i = 0; i < kN; i++)
{
/* FIXME: Singular matrix? */
const float_type f = 1.0 / W(i,i);
const auto * const p = m_terms[i].m_nzrd.data();
const size_t e = m_terms[i].m_nzrd.size();
/* Eliminate column i from row j */
const auto * const pb = m_terms[i].m_nzbd.data();
const size_t eb = m_terms[i].m_nzbd.size();
for (std::size_t jb = 0; jb < eb; jb++)
{
const auto 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;
}
}
}
/* up */
for (std::size_t i = kN; i-- > 0; )
{
/* FIXME: Singular matrix? */
const float_type f = 1.0 / W(i,i);
for (std::size_t j = i; j-- > 0; )
{
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;
}
}
}
void reset() override { matrix_solver_t::reset(); }
protected:
unsigned vsolve_non_dynamic(const bool newton_raphson) override;
unsigned solve_non_dynamic(const bool newton_raphson);
constexpr std::size_t size() const { return m_dim; }
void LE_invert();
template <typename FT, int SIZE>
template <typename T>
void LE_compute_x(T & x);
template <typename T1, typename T2>
float_ext_type &A(const T1 &r, const T2 &c) { return m_A[r][c]; }
template <typename T1, typename T2>
float_ext_type &W(const T1 &r, const T2 &c) { return m_W[r][c]; }
/* access to Ainv for fixed columns over row, there store transposed */
template <typename T1, typename T2>
float_ext_type &Ainv(const T1 &r, const T2 &c) { return m_Ainv[c][r]; }
template <typename T1>
float_ext_type &RHS(const T1 &r) { return m_RHS[r]; }
template <typename T1, typename T2>
float_ext_type &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; }
private:
template <typename T, std::size_t N, std::size_t M>
using array2D = std::array<std::array<T, M>, N>;
static constexpr std::size_t m_pitch = ((( storage_N) + 7) / 8) * 8;
array2D<float_ext_type, storage_N, m_pitch> m_A;
array2D<float_ext_type, storage_N, m_pitch> m_Ainv;
array2D<float_ext_type, storage_N, m_pitch> m_W;
std::array<float_ext_type, storage_N> m_RHS; // right hand side - contains currents
array2D<float_ext_type, storage_N, m_pitch> m_lA;
/* temporary */
array2D<float_ext_type, storage_N, m_pitch> H;
std::array<unsigned, storage_N> rows;
array2D<unsigned, storage_N, m_pitch> cols;
std::array<unsigned, storage_N> colcount;
unsigned m_cnt;
//float_ext_type m_RHSx[storage_N];
const std::size_t m_dim;
};
// ----------------------------------------------------------------------------------------
// matrix_solver_direct
// ----------------------------------------------------------------------------------------
template <typename FT, int SIZE>
void matrix_solver_w_t<FT, SIZE>::LE_invert()
{
const std::size_t kN = size();
for (std::size_t i = 0; i < kN; i++)
void matrix_solver_w_t<FT, SIZE>::LE_compute_x(
T & x)
{
for (std::size_t j = 0; j < kN; j++)
{
W(i,j) = lA(i,j) = A(i,j);
Ainv(i,j) = 0.0;
}
Ainv(i,i) = 1.0;
}
/* down */
for (std::size_t i = 0; i < kN; i++)
{
/* FIXME: Singular matrix? */
const float_type f = 1.0 / W(i,i);
const auto * const p = m_terms[i].m_nzrd.data();
const size_t e = m_terms[i].m_nzrd.size();
/* Eliminate column i from row j */
const auto * const pb = m_terms[i].m_nzbd.data();
const size_t eb = m_terms[i].m_nzbd.size();
for (std::size_t jb = 0; jb < eb; jb++)
{
const auto 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;
}
}
}
/* up */
for (std::size_t i = kN; i-- > 0; )
{
/* FIXME: Singular matrix? */
const float_type f = 1.0 / W(i,i);
for (std::size_t j = i; j-- > 0; )
{
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;
}
}
}
template <typename FT, int SIZE>
template <typename T>
void matrix_solver_w_t<FT, SIZE>::LE_compute_x(
T & x)
{
const std::size_t kN = size();
for (std::size_t i=0; i<kN; i++)
x[i] = 0.0;
for (std::size_t k=0; k<kN; k++)
{
const float_type f = RHS(k);
const std::size_t kN = size();
for (std::size_t i=0; i<kN; i++)
x[i] += Ainv(i,k) * f;
}
}
x[i] = 0.0;
template <typename FT, int SIZE>
unsigned matrix_solver_w_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphson)
{
const auto iN = size();
// NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init)
std::array<float_type, storage_N> new_V;
// NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init)
std::array<float_type, storage_N> t; // FIXME: convert to member
// NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init)
std::array<float_type, storage_N> w;
if ((m_cnt % 50) == 0)
{
/* complete calculation */
this->LE_invert();
this->LE_compute_x(new_V);
}
else
{
/* Solve Ay = b for y */
this->LE_compute_x(new_V);
/* determine changed rows */
unsigned rowcount=0;
#define VT(r,c) (A(r,c) - lA(r,c))
for (unsigned row = 0; row < iN; row ++)
for (std::size_t k=0; k<kN; k++)
{
unsigned cc=0;
auto &nz = m_terms[row].m_nz;
for (auto & col : nz)
{
if (A(row,col) != lA(row,col))
cols[rowcount][cc++] = col;
}
if (cc > 0)
{
colcount[rowcount] = cc;
rows[rowcount++] = row;
}
const float_type f = RHS(k);
for (std::size_t i=0; i<kN; i++)
x[i] += Ainv(i,k) * f;
}
if (rowcount > 0)
}
template <typename FT, int SIZE>
unsigned matrix_solver_w_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphson)
{
const auto iN = size();
// NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init)
std::array<float_type, storage_N> new_V;
// NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init)
std::array<float_type, storage_N> t; // FIXME: convert to member
// NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init)
std::array<float_type, storage_N> w;
if ((m_cnt % 50) == 0)
{
/* construct w = transform(V) * y
* dim: rowcount x iN
* */
for (unsigned i = 0; i < rowcount; i++)
/* complete calculation */
this->LE_invert();
this->LE_compute_x(new_V);
}
else
{
/* Solve Ay = b for y */
this->LE_compute_x(new_V);
/* determine changed rows */
unsigned rowcount=0;
#define VT(r,c) (A(r,c) - lA(r,c))
for (unsigned row = 0; row < iN; row ++)
{
const unsigned r = rows[i];
double tmp = 0.0;
for (unsigned k = 0; k < iN; k++)
tmp += VT(r,k) * new_V[k];
w[i] = tmp;
}
for (unsigned i = 0; i < rowcount; i++)
for (unsigned k=0; k< rowcount; k++)
H[i][k] = 0.0;
for (unsigned i = 0; i < rowcount; i++)
H[i][i] = 1.0;
/* Construct H = (I + VT*Z) */
for (unsigned i = 0; i < rowcount; i++)
for (unsigned k=0; k< colcount[i]; k++)
unsigned cc=0;
auto &nz = m_terms[row].m_nz;
for (auto & col : nz)
{
const unsigned col = cols[i][k];
float_type f = VT(rows[i],col);
if (f!=0.0)
for (unsigned j= 0; j < rowcount; j++)
H[i][j] += f * Ainv(col,rows[j]);
if (A(row,col) != lA(row,col))
cols[rowcount][cc++] = col;
}
if (cc > 0)
{
colcount[rowcount] = cc;
rows[rowcount++] = row;
}
}
if (rowcount > 0)
{
/* construct w = transform(V) * y
* dim: rowcount x iN
* */
for (unsigned i = 0; i < rowcount; i++)
{
const unsigned r = rows[i];
double tmp = 0.0;
for (unsigned k = 0; k < iN; k++)
tmp += VT(r,k) * new_V[k];
w[i] = tmp;
}
/* Gaussian elimination of H */
for (unsigned i = 0; i < rowcount; i++)
{
if (H[i][i] == 0.0)
plib::perrlogger("{} H singular\n", this->name());
const float_type f = 1.0 / H[i][i];
for (unsigned j = i+1; j < rowcount; j++)
{
const float_type f1 = - f * H[j][i];
for (unsigned i = 0; i < rowcount; i++)
for (unsigned k=0; k< rowcount; k++)
H[i][k] = 0.0;
if (f1!=0.0)
for (unsigned i = 0; i < rowcount; i++)
H[i][i] = 1.0;
/* Construct H = (I + VT*Z) */
for (unsigned i = 0; i < rowcount; i++)
for (unsigned k=0; k< colcount[i]; k++)
{
float_type *pj = &H[j][i+1];
const float_type *pi = &H[i][i+1];
for (unsigned k = 0; k < rowcount-i-1; k++)
pj[k] += f1 * pi[k];
//H[j][k] += f1 * H[i][k];
w[j] += f1 * w[i];
const unsigned col = cols[i][k];
float_type f = VT(rows[i],col);
if (f!=0.0)
for (unsigned j= 0; j < rowcount; j++)
H[i][j] += f * Ainv(col,rows[j]);
}
/* Gaussian elimination of H */
for (unsigned i = 0; i < rowcount; i++)
{
if (H[i][i] == 0.0)
plib::perrlogger("{} H singular\n", this->name());
const float_type f = 1.0 / H[i][i];
for (unsigned j = i+1; j < rowcount; j++)
{
const float_type f1 = - f * H[j][i];
if (f1!=0.0)
{
float_type *pj = &H[j][i+1];
const float_type *pi = &H[i][i+1];
for (unsigned k = 0; k < rowcount-i-1; k++)
pj[k] += f1 * pi[k];
//H[j][k] += f1 * H[i][k];
w[j] += f1 * w[i];
}
}
}
}
/* Back substitution */
//inv(H) w = t w = H t
for (unsigned j = rowcount; j-- > 0; )
{
float_type tmp = 0;
const float_type *pj = &H[j][j+1];
const float_type *tj = &t[j+1];
for (unsigned k = 0; k < rowcount-j-1; k++)
tmp += pj[k] * tj[k];
//tmp += H[j][k] * t[k];
t[j] = (w[j] - tmp) / H[j][j];
}
/* Back substitution */
//inv(H) w = t w = H t
for (unsigned j = rowcount; j-- > 0; )
{
float_type tmp = 0;
const float_type *pj = &H[j][j+1];
const float_type *tj = &t[j+1];
for (unsigned k = 0; k < rowcount-j-1; k++)
tmp += pj[k] * tj[k];
//tmp += H[j][k] * t[k];
t[j] = (w[j] - tmp) / H[j][j];
}
/* x = y - Zt */
/* x = y - Zt */
for (unsigned i=0; i<iN; i++)
{
float_type tmp = 0.0;
for (unsigned j=0; j<rowcount;j++)
{
const unsigned row = rows[j];
tmp += Ainv(i,row) * t[j];
}
new_V[i] -= tmp;
}
}
}
m_cnt++;
if (false)
for (unsigned i=0; i<iN; i++)
{
float_type tmp = 0.0;
for (unsigned j=0; j<rowcount;j++)
for (unsigned j=0; j<iN; j++)
{
const unsigned row = rows[j];
tmp += Ainv(i,row) * t[j];
tmp += A(i,j) * new_V[j];
}
new_V[i] -= tmp;
if (std::abs(tmp-RHS(i)) > 1e-6)
plib::perrlogger("{} failed on row {}: {} RHS: {}\n", this->name(), i, std::abs(tmp-RHS(i)), RHS(i));
}
}
const float_type err = (newton_raphson ? delta(new_V) : 0.0);
store(new_V);
return (err > this->m_params.m_accuracy) ? 2 : 1;
}
m_cnt++;
if (false)
for (unsigned i=0; i<iN; i++)
{
float_type tmp = 0.0;
for (unsigned j=0; j<iN; j++)
{
tmp += A(i,j) * new_V[j];
}
if (std::abs(tmp-RHS(i)) > 1e-6)
plib::perrlogger("{} failed on row {}: {} RHS: {}\n", this->name(), i, std::abs(tmp-RHS(i)), RHS(i));
}
template <typename FT, int SIZE>
unsigned matrix_solver_w_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson)
{
this->build_LE_A(*this);
this->build_LE_RHS(*this);
const float_type err = (newton_raphson ? delta(new_V) : 0.0);
store(new_V);
return (err > this->m_params.m_accuracy) ? 2 : 1;
}
this->m_stat_calculations++;
return this->solve_non_dynamic(newton_raphson);
}
template <typename FT, int SIZE>
unsigned matrix_solver_w_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson)
{
this->build_LE_A(*this);
this->build_LE_RHS(*this);
this->m_stat_calculations++;
return this->solve_non_dynamic(newton_raphson);
}
} //namespace devices
} // namespace solver
} // namespace netlist
#endif /* NLD_MS_DIRECT_H_ */

46
src/lib/netlist/solver/nld_solver.cpp
View File

@ -81,7 +81,7 @@ namespace devices
std::size_t nthreads = std::min(static_cast<std::size_t>(m_params.m_parallel()), plib::omp::get_max_threads());
std::vector<matrix_solver_t *> &solvers = (force_solve ? m_mat_solvers_all : m_mat_solvers_timestepping);
std::vector<solver::matrix_solver_t *> &solvers = (force_solve ? m_mat_solvers_all : m_mat_solvers_timestepping);
if (nthreads > 1 && solvers.size() > 1)
{
@ -110,45 +110,45 @@ namespace devices
}
template <class C>
plib::unique_ptr<matrix_solver_t> create_it(netlist_state_t &nl, pstring name,
plib::unique_ptr<solver::matrix_solver_t> create_it(netlist_state_t &nl, pstring name,
analog_net_t::list_t &nets,
solver_parameters_t &params, std::size_t size)
solver::solver_parameters_t &params, std::size_t size)
{
return plib::make_unique<C>(nl, name, nets, &params, size);
}
template <typename FT, int SIZE>
plib::unique_ptr<matrix_solver_t> NETLIB_NAME(solver)::create_solver(std::size_t size,
plib::unique_ptr<solver::matrix_solver_t> NETLIB_NAME(solver)::create_solver(std::size_t size,
const pstring &solvername,
analog_net_t::list_t &nets)
{
switch (m_params.m_method())
{
case matrix_type_e::MAT_CR:
case solver::matrix_type_e::MAT_CR:
if (size > 0) // GCR always outperforms MAT solver
{
return create_it<matrix_solver_GCR_t<FT, SIZE>>(state(), solvername, nets, m_params, size);
return create_it<solver::matrix_solver_GCR_t<FT, SIZE>>(state(), solvername, nets, m_params, size);
}
else
{
return create_it<matrix_solver_direct_t<FT, SIZE>>(state(), solvername, nets, m_params, size);
return create_it<solver::matrix_solver_direct_t<FT, SIZE>>(state(), solvername, nets, m_params, size);
}
case matrix_type_e::SOR_MAT:
return create_it<matrix_solver_SOR_mat_t<FT, SIZE>>(state(), solvername, nets, m_params, size);
case matrix_type_e::MAT:
return create_it<matrix_solver_direct_t<FT, SIZE>>(state(), solvername, nets, m_params, size);
case matrix_type_e::SM:
case solver::matrix_type_e::SOR_MAT:
return create_it<solver::matrix_solver_SOR_mat_t<FT, SIZE>>(state(), solvername, nets, m_params, size);
case solver::matrix_type_e::MAT:
return create_it<solver::matrix_solver_direct_t<FT, SIZE>>(state(), solvername, nets, m_params, size);
case solver::matrix_type_e::SM:
/* Sherman-Morrison Formula */
return create_it<matrix_solver_sm_t<FT, SIZE>>(state(), solvername, nets, m_params, size);
case matrix_type_e::W:
return create_it<solver::matrix_solver_sm_t<FT, SIZE>>(state(), solvername, nets, m_params, size);
case solver::matrix_type_e::W:
/* Woodbury Formula */
return create_it<matrix_solver_w_t<FT, SIZE>>(state(), solvername, nets, m_params, size);
case matrix_type_e::SOR:
return create_it<matrix_solver_SOR_t<FT, SIZE>>(state(), solvername, nets, m_params, size);
case matrix_type_e::GMRES:
return create_it<matrix_solver_GMRES_t<FT, SIZE>>(state(), solvername, nets, m_params, size);
return create_it<solver::matrix_solver_w_t<FT, SIZE>>(state(), solvername, nets, m_params, size);
case solver::matrix_type_e::SOR:
return create_it<solver::matrix_solver_SOR_t<FT, SIZE>>(state(), solvername, nets, m_params, size);
case solver::matrix_type_e::GMRES:
return create_it<solver::matrix_solver_GMRES_t<FT, SIZE>>(state(), solvername, nets, m_params, size);
}
return plib::unique_ptr<matrix_solver_t>();
return plib::unique_ptr<solver::matrix_solver_t>();
}
struct net_splitter
@ -224,17 +224,17 @@ namespace devices
log().verbose("Found {1} net groups in {2} nets\n", splitter.groups.size(), state().nets().size());
for (auto & grp : splitter.groups)
{
plib::unique_ptr<matrix_solver_t> ms;
plib::unique_ptr<solver::matrix_solver_t> ms;
std::size_t net_count = grp.size();
pstring sname = plib::pfmt("Solver_{1}")(m_mat_solvers.size());
switch (net_count)
{
case 1:
ms = plib::make_unique<matrix_solver_direct1_t<double>>(state(), sname, grp, &m_params);
ms = plib::make_unique<solver::matrix_solver_direct1_t<double>>(state(), sname, grp, &m_params);
break;
case 2:
ms = plib::make_unique<matrix_solver_direct2_t<double>>(state(), sname, grp, &m_params);
ms = plib::make_unique<solver::matrix_solver_direct2_t<double>>(state(), sname, grp, &m_params);
break;
#if 0
case 3:

10
src/lib/netlist/solver/nld_solver.h
View File

@ -51,14 +51,14 @@ namespace devices
logic_input_t m_fb_step;
logic_output_t m_Q_step;
std::vector<plib::unique_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;
std::vector<plib::unique_ptr<solver::matrix_solver_t>> m_mat_solvers;
std::vector<solver::matrix_solver_t *> m_mat_solvers_all;
std::vector<solver::matrix_solver_t *> m_mat_solvers_timestepping;
solver_parameters_t m_params;
solver::solver_parameters_t m_params;
template <typename FT, int SIZE>
plib::unique_ptr<matrix_solver_t> create_solver(std::size_t size,
plib::unique_ptr<solver::matrix_solver_t> create_solver(std::size_t size,
const pstring &solvername, analog_net_t::list_t &nets);
};