mirror of
https://github.com/holub/mame
synced 2025-04-23 00:39:36 +03:00
Solver model simplification.
This commit is contained in:
couriersud
parent
db1e7f6b07
commit
0d1e57cc40
@ -11,6 +11,7 @@
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include <cstdlib>
|
||||
#include <algorithm>
|
||||
|
||||
#include "plib/pstring.h"
|
||||
#include "plib/plists.h"
|
||||
@ -40,26 +41,31 @@ private:
|
||||
nl_math() {};
|
||||
|
||||
public:
|
||||
ATTR_HOT inline static float exp(const float x) { return std::exp(x); }
|
||||
ATTR_HOT inline static double abs(const double x) { return std::fabs(x); }
|
||||
ATTR_HOT inline static float abs(const float x) { return std::fabs(x); }
|
||||
ATTR_HOT inline static double log(const double x) { return std::log(x); }
|
||||
ATTR_HOT inline static float log(const float x) { return std::log(x); }
|
||||
ATTR_HOT inline static float exp(const float &x) { return std::exp(x); }
|
||||
|
||||
ATTR_HOT inline static double abs(const double &x) { return std::abs(x); }
|
||||
ATTR_HOT inline static float abs(const float &x) { return std::abs(x); }
|
||||
|
||||
ATTR_HOT inline static double max(const double &x, const double &y) { return std::max(x, y); }
|
||||
ATTR_HOT inline static float max(const float &x, const float &y) { return std::max(x, y); }
|
||||
|
||||
ATTR_HOT inline static double log(const double &x) { return std::log(x); }
|
||||
ATTR_HOT inline static float log(const float &x) { return std::log(x); }
|
||||
#if defined(_MSC_VER) && _MSC_VER < 1800
|
||||
ATTR_HOT inline static double e_log1p(const double x) { return nl_math::log(1.0 + x); }
|
||||
ATTR_HOT inline static float e_log1p(const float x) { return nl_math::log(1.0 + x); }
|
||||
ATTR_HOT inline static double e_log1p(const double &x) { return nl_math::log(1.0 + x); }
|
||||
ATTR_HOT inline static float e_log1p(const float &x) { return nl_math::log(1.0 + x); }
|
||||
#else
|
||||
ATTR_HOT inline static double e_log1p(const double x) { return log1p(x); }
|
||||
ATTR_HOT inline static float e_log1p(const float x) { return log1pf(x); }
|
||||
ATTR_HOT inline static double e_log1p(const double &x) { return log1p(x); }
|
||||
ATTR_HOT inline static float e_log1p(const float &x) { return log1pf(x); }
|
||||
#endif
|
||||
ATTR_HOT inline static double sqrt(const double x) { return std::sqrt(x); }
|
||||
ATTR_HOT inline static float sqrt(const float x) { return std::sqrt(x); }
|
||||
ATTR_HOT inline static double sqrt(const double &x) { return std::sqrt(x); }
|
||||
ATTR_HOT inline static float sqrt(const float &x) { return std::sqrt(x); }
|
||||
|
||||
// this one has an accuracy of better than 5%. That's enough for our purpose
|
||||
// add c3 and it'll be better than 1%
|
||||
|
||||
#if 0
|
||||
inline static double fastexp_h(const double x)
|
||||
inline static double fastexp_h(const double &x)
|
||||
{
|
||||
/* static */ const double ln2r = 1.442695040888963387;
|
||||
/* static */ const double ln2 = 0.693147180559945286;
|
||||
@ -79,7 +85,7 @@ public:
|
||||
return pow(2.0, t)*e;
|
||||
}
|
||||
|
||||
ATTR_HOT inline static double exp(const double x)
|
||||
ATTR_HOT inline static double exp(const double &x)
|
||||
{
|
||||
if (x<0)
|
||||
return 1.0 / fastexp_h(-x);
|
||||
@ -87,7 +93,7 @@ public:
|
||||
return fastexp_h(x);
|
||||
}
|
||||
#else
|
||||
ATTR_HOT inline static double exp(const double x) { return std::exp(x); }
|
||||
ATTR_HOT inline static double exp(const double &x) { return std::exp(x); }
|
||||
#endif
|
||||
|
||||
};
|
||||
|
||||
@ -16,8 +16,8 @@ class terms_t
|
||||
{
|
||||
P_PREVENT_COPYING(terms_t)
|
||||
|
||||
public:
|
||||
ATTR_COLD terms_t() : m_railstart(0)
|
||||
public:
|
||||
ATTR_COLD terms_t() : m_railstart(0), m_last_V(0)
|
||||
{}
|
||||
|
||||
ATTR_COLD void clear()
|
||||
@ -48,6 +48,10 @@ class terms_t
|
||||
pvector_t<unsigned> m_nz; /* all non zero for multiplication */
|
||||
pvector_t<unsigned> m_nzrd; /* non zero right of the diagonal for elimination, may include RHS element */
|
||||
pvector_t<unsigned> m_nzbd; /* non zero below of the diagonal for elimination */
|
||||
|
||||
/* state */
|
||||
nl_double m_last_V;
|
||||
|
||||
private:
|
||||
pvector_t<terminal_t *> m_term;
|
||||
pvector_t<int> m_net_other;
|
||||
@ -55,6 +59,7 @@ private:
|
||||
pvector_t<nl_double> m_gt;
|
||||
pvector_t<nl_double> m_Idr;
|
||||
pvector_t<nl_double *> m_other_curanalog;
|
||||
|
||||
};
|
||||
|
||||
class matrix_solver_t : public device_t
|
||||
@ -104,14 +109,23 @@ protected:
|
||||
ATTR_COLD void setup_base(analog_net_t::list_t &nets);
|
||||
void update_dynamic();
|
||||
|
||||
virtual void add_term(int net_idx, terminal_t *term) = 0;
|
||||
virtual void vsetup(analog_net_t::list_t &nets) = 0;
|
||||
virtual int vsolve_non_dynamic(const bool newton_raphson) = 0;
|
||||
virtual netlist_time compute_next_timestep() = 0;
|
||||
|
||||
/* virtual */ netlist_time compute_next_timestep();
|
||||
/* virtual */ void add_term(int net_idx, terminal_t *term);
|
||||
|
||||
template <typename T>
|
||||
void store(const T * RESTRICT V);
|
||||
template <typename T>
|
||||
T delta(const T * RESTRICT V);
|
||||
|
||||
pvector_t<terms_t *> m_terms;
|
||||
pvector_t<analog_net_t *> m_nets;
|
||||
pvector_t<analog_output_t *> m_inps;
|
||||
|
||||
pvector_t<terms_t *> m_rails_temp;
|
||||
|
||||
int m_stat_calculations;
|
||||
int m_stat_newton_raphson;
|
||||
int m_stat_vsolver_calls;
|
||||
@ -138,6 +152,28 @@ private:
|
||||
const eSolverType m_type;
|
||||
};
|
||||
|
||||
template <typename T>
|
||||
T matrix_solver_t::delta(const T * RESTRICT V)
|
||||
{
|
||||
/* FIXME: Ideally we should also include currents (RHS) here. This would
|
||||
* need a revaluation of the right hand side after voltages have been updated
|
||||
* and thus belong into a different calculation. This applies to all solvers.
|
||||
*/
|
||||
|
||||
const unsigned iN = this->m_terms.size();
|
||||
T cerr = 0;
|
||||
for (unsigned i = 0; i < iN; i++)
|
||||
cerr = nl_math::max(cerr, nl_math::abs(V[i] - (T) this->m_nets[i]->m_cur_Analog));
|
||||
return cerr;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void matrix_solver_t::store(const T * RESTRICT V)
|
||||
{
|
||||
for (unsigned i = 0, iN=m_terms.size(); i < iN; i++)
|
||||
this->m_nets[i]->m_cur_Analog = V[i];
|
||||
}
|
||||
|
||||
|
||||
NETLIB_NAMESPACE_DEVICES_END()
|
||||
|
||||
|
||||
@ -129,7 +129,6 @@ public:
|
||||
virtual void reset() override { matrix_solver_t::reset(); }
|
||||
|
||||
protected:
|
||||
virtual void add_term(int net_idx, terminal_t *term) override;
|
||||
virtual int vsolve_non_dynamic(const bool newton_raphson) override;
|
||||
int solve_non_dynamic(const bool newton_raphson);
|
||||
|
||||
@ -142,14 +141,6 @@ protected:
|
||||
template <typename T>
|
||||
void LE_back_subst(T * RESTRICT x);
|
||||
|
||||
template <typename T>
|
||||
T delta(const T * RESTRICT V);
|
||||
|
||||
template <typename T>
|
||||
void store(const T * RESTRICT V);
|
||||
|
||||
virtual netlist_time compute_next_timestep() override;
|
||||
|
||||
#if TEST_PARALLEL
|
||||
int x_i[10];
|
||||
int x_start[10];
|
||||
@ -169,10 +160,6 @@ protected:
|
||||
inline nl_ext_double &RHS(const T1 &r) { return m_A[r][N()]; }
|
||||
#endif
|
||||
ATTR_ALIGN nl_double m_last_RHS[_storage_N]; // right hand side - contains currents
|
||||
ATTR_ALIGN nl_double m_last_V[_storage_N];
|
||||
|
||||
terms_t * m_terms[_storage_N];
|
||||
terms_t *m_rails_temp;
|
||||
|
||||
private:
|
||||
static const std::size_t m_pitch = (((_storage_N + 1) + 7) / 8) * 8;
|
||||
@ -196,109 +183,36 @@ private:
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
matrix_solver_direct_t<m_N, _storage_N>::~matrix_solver_direct_t()
|
||||
{
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
pfree(m_terms[k]);
|
||||
}
|
||||
pfree_array(m_rails_temp);
|
||||
#if (NL_USE_DYNAMIC_ALLOCATION)
|
||||
pfree_array(m_A);
|
||||
#endif
|
||||
#if TEST_PARALLEL
|
||||
thr_dispose();
|
||||
#endif
|
||||
|
||||
}
|
||||
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
netlist_time matrix_solver_direct_t<m_N, _storage_N>::compute_next_timestep()
|
||||
{
|
||||
nl_double new_solver_timestep = m_params.m_max_timestep;
|
||||
|
||||
if (m_params.m_dynamic)
|
||||
{
|
||||
/*
|
||||
* FIXME: We should extend the logic to use either all nets or
|
||||
* only output nets.
|
||||
*/
|
||||
for (unsigned k = 0, iN=N(); k < iN; k++)
|
||||
{
|
||||
analog_net_t *n = m_nets[k];
|
||||
|
||||
const nl_double DD_n = (n->Q_Analog() - m_last_V[k]);
|
||||
const nl_double hn = current_timestep();
|
||||
|
||||
nl_double DD2 = (DD_n / hn - n->m_DD_n_m_1 / n->m_h_n_m_1) / (hn + n->m_h_n_m_1);
|
||||
nl_double new_net_timestep;
|
||||
|
||||
n->m_h_n_m_1 = hn;
|
||||
n->m_DD_n_m_1 = DD_n;
|
||||
if (nl_math::abs(DD2) > NL_FCONST(1e-30)) // avoid div-by-zero
|
||||
new_net_timestep = nl_math::sqrt(m_params.m_lte / nl_math::abs(NL_FCONST(0.5)*DD2));
|
||||
else
|
||||
new_net_timestep = m_params.m_max_timestep;
|
||||
|
||||
if (new_net_timestep < new_solver_timestep)
|
||||
new_solver_timestep = new_net_timestep;
|
||||
|
||||
m_last_V[k] = n->Q_Analog();
|
||||
}
|
||||
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;
|
||||
return netlist_time::from_double(new_solver_timestep);
|
||||
}
|
||||
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
ATTR_COLD void matrix_solver_direct_t<m_N, _storage_N>::add_term(int k, terminal_t *term)
|
||||
{
|
||||
if (term->m_otherterm->net().isRailNet())
|
||||
{
|
||||
m_rails_temp[k].add(term, -1, false);
|
||||
}
|
||||
else
|
||||
{
|
||||
int ot = get_net_idx(&term->m_otherterm->net());
|
||||
if (ot>=0)
|
||||
{
|
||||
m_terms[k]->add(term, ot, true);
|
||||
}
|
||||
/* Should this be allowed ? */
|
||||
else // if (ot<0)
|
||||
{
|
||||
m_rails_temp[k].add(term, ot, true);
|
||||
log().fatal("found term with missing othernet {1}\n", term->name());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
ATTR_COLD void matrix_solver_direct_t<m_N, _storage_N>::vsetup(analog_net_t::list_t &nets)
|
||||
{
|
||||
if (m_dim < nets.size())
|
||||
log().fatal("Dimension {1} less than {2}", m_dim, nets.size());
|
||||
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
m_terms[k]->clear();
|
||||
m_rails_temp[k].clear();
|
||||
}
|
||||
|
||||
matrix_solver_t::setup_base(nets);
|
||||
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
m_terms[k]->m_railstart = m_terms[k]->count();
|
||||
for (unsigned i = 0; i < m_rails_temp[k].count(); i++)
|
||||
this->m_terms[k]->add(m_rails_temp[k].terms()[i], m_rails_temp[k].net_other()[i], false);
|
||||
for (unsigned i = 0; i < m_rails_temp[k]->count(); i++)
|
||||
this->m_terms[k]->add(m_rails_temp[k]->terms()[i], m_rails_temp[k]->net_other()[i], false);
|
||||
|
||||
m_rails_temp[k].clear(); // no longer needed
|
||||
m_rails_temp[k]->clear(); // no longer needed
|
||||
m_terms[k]->set_pointers();
|
||||
}
|
||||
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
pfree(m_rails_temp[k]); // no longer needed
|
||||
|
||||
m_rails_temp.clear();
|
||||
#if 1
|
||||
|
||||
/* Sort in descending order by number of connected matrix voltages.
|
||||
@ -430,13 +344,13 @@ ATTR_COLD void matrix_solver_direct_t<m_N, _storage_N>::vsetup(analog_net_t::lis
|
||||
* save states
|
||||
*/
|
||||
save(NLNAME(m_last_RHS));
|
||||
save(NLNAME(m_last_V));
|
||||
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
pstring num = pfmt("{1}")(k);
|
||||
|
||||
save(RHS(k), "RHS" + num);
|
||||
save(m_terms[k]->m_last_V, "lastV" + num);
|
||||
|
||||
save(m_terms[k]->go(),"GO" + num, m_terms[k]->count());
|
||||
save(m_terms[k]->gt(),"GT" + num, m_terms[k]->count());
|
||||
@ -682,35 +596,6 @@ void matrix_solver_direct_t<m_N, _storage_N>::LE_back_subst(
|
||||
}
|
||||
|
||||
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
template <typename T>
|
||||
T matrix_solver_direct_t<m_N, _storage_N>::delta(
|
||||
const T * RESTRICT V)
|
||||
{
|
||||
/* FIXME: Ideally we should also include currents (RHS) here. This would
|
||||
* need a revaluation of the right hand side after voltages have been updated
|
||||
* and thus belong into a different calculation. This applies to all solvers.
|
||||
*/
|
||||
|
||||
const unsigned iN = this->N();
|
||||
T cerr = 0;
|
||||
for (unsigned i = 0; i < iN; i++)
|
||||
cerr = std::fmax(cerr, nl_math::abs(V[i] - (T) this->m_nets[i]->m_cur_Analog));
|
||||
return cerr;
|
||||
}
|
||||
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
template <typename T>
|
||||
void matrix_solver_direct_t<m_N, _storage_N>::store(
|
||||
const T * RESTRICT V)
|
||||
{
|
||||
for (unsigned i = 0, iN=N(); i < iN; i++)
|
||||
{
|
||||
this->m_nets[i]->m_cur_Analog = V[i];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
int matrix_solver_direct_t<m_N, _storage_N>::solve_non_dynamic(ATTR_UNUSED const bool newton_raphson)
|
||||
{
|
||||
@ -752,15 +637,12 @@ matrix_solver_direct_t<m_N, _storage_N>::matrix_solver_direct_t(const solver_par
|
||||
: matrix_solver_t(GAUSSIAN_ELIMINATION, params)
|
||||
, m_dim(size)
|
||||
{
|
||||
m_rails_temp = palloc_array(terms_t, N());
|
||||
#if (NL_USE_DYNAMIC_ALLOCATION)
|
||||
m_A = palloc_array(nl_ext_double, N() * m_pitch);
|
||||
#endif
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
m_terms[k] = palloc(terms_t);
|
||||
m_last_RHS[k] = 0.0;
|
||||
m_last_V[k] = 0.0;
|
||||
}
|
||||
#if TEST_PARALLEL
|
||||
thr_initialize();
|
||||
@ -772,15 +654,12 @@ matrix_solver_direct_t<m_N, _storage_N>::matrix_solver_direct_t(const eSolverTyp
|
||||
: matrix_solver_t(type, params)
|
||||
, m_dim(size)
|
||||
{
|
||||
m_rails_temp = palloc_array(terms_t, N());
|
||||
#if (NL_USE_DYNAMIC_ALLOCATION)
|
||||
m_A = palloc_array(nl_ext_double, N() * m_pitch);
|
||||
#endif
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
m_terms[k] = palloc(terms_t);
|
||||
m_last_RHS[k] = 0.0;
|
||||
m_last_V[k] = 0.0;
|
||||
}
|
||||
#if TEST_PARALLEL
|
||||
thr_initialize();
|
||||
|
||||
@ -144,7 +144,6 @@ protected:
|
||||
ATTR_ALIGN nl_double m_last_RHS[_storage_N]; // right hand side - contains currents
|
||||
ATTR_ALIGN nl_double m_last_V[_storage_N];
|
||||
|
||||
terms_t **m_terms;
|
||||
terms_t *m_rails_temp;
|
||||
|
||||
private:
|
||||
@ -161,11 +160,6 @@ private:
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
matrix_solver_direct_t<m_N, _storage_N>::~matrix_solver_direct_t()
|
||||
{
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
pfree(m_terms[k]);
|
||||
}
|
||||
pfree_array(m_terms);
|
||||
pfree_array(m_rails_temp);
|
||||
}
|
||||
|
||||
@ -608,14 +602,11 @@ matrix_solver_direct_t<m_N, _storage_N>::matrix_solver_direct_t(const solver_par
|
||||
, m_dim(size)
|
||||
, m_lp_fact(0)
|
||||
{
|
||||
m_terms = palloc_array(terms_t *, N());
|
||||
m_rails_temp = palloc_array(terms_t, N());
|
||||
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
m_terms[k] = palloc(terms_t);
|
||||
m_last_RHS[k] = 0.0;
|
||||
m_last_V[k] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
@ -625,14 +616,12 @@ matrix_solver_direct_t<m_N, _storage_N>::matrix_solver_direct_t(const eSolverTyp
|
||||
, m_dim(size)
|
||||
, m_lp_fact(0)
|
||||
{
|
||||
m_terms = palloc_array(terms_t *, N());
|
||||
m_rails_temp = palloc_array(terms_t, N());
|
||||
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
m_terms[k] = palloc(terms_t);
|
||||
m_last_RHS[k] = 0.0;
|
||||
m_last_V[k] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -59,7 +59,6 @@ public:
|
||||
virtual void reset() override { matrix_solver_t::reset(); }
|
||||
|
||||
protected:
|
||||
virtual void add_term(int net_idx, terminal_t *term) override;
|
||||
virtual int vsolve_non_dynamic(const bool newton_raphson) override;
|
||||
int solve_non_dynamic(const bool newton_raphson);
|
||||
|
||||
@ -72,13 +71,6 @@ protected:
|
||||
template <typename T>
|
||||
void LE_compute_x(T * RESTRICT x);
|
||||
|
||||
template <typename T>
|
||||
T delta(const T * RESTRICT V);
|
||||
|
||||
template <typename T>
|
||||
void store(const T * RESTRICT V);
|
||||
|
||||
virtual netlist_time compute_next_timestep() override;
|
||||
|
||||
template <typename T1, typename T2>
|
||||
inline nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r][c]; }
|
||||
@ -96,10 +88,6 @@ protected:
|
||||
inline nl_ext_double &lAinv(const T1 &r, const T2 &c) { return m_lAinv[r][c]; }
|
||||
|
||||
ATTR_ALIGN nl_double m_last_RHS[_storage_N]; // right hand side - contains currents
|
||||
ATTR_ALIGN nl_double m_last_V[_storage_N];
|
||||
|
||||
terms_t * m_terms[_storage_N];
|
||||
terms_t *m_rails_temp;
|
||||
|
||||
private:
|
||||
static const std::size_t m_pitch = ((( _storage_N) + 7) / 8) * 8;
|
||||
@ -124,102 +112,26 @@ private:
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
matrix_solver_sm_t<m_N, _storage_N>::~matrix_solver_sm_t()
|
||||
{
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
pfree(m_terms[k]);
|
||||
}
|
||||
pfree_array(m_rails_temp);
|
||||
#if (NL_USE_DYNAMIC_ALLOCATION)
|
||||
pfree_array(m_A);
|
||||
#endif
|
||||
}
|
||||
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
netlist_time matrix_solver_sm_t<m_N, _storage_N>::compute_next_timestep()
|
||||
{
|
||||
nl_double new_solver_timestep = m_params.m_max_timestep;
|
||||
|
||||
if (m_params.m_dynamic)
|
||||
{
|
||||
/*
|
||||
* FIXME: We should extend the logic to use either all nets or
|
||||
* only output nets.
|
||||
*/
|
||||
for (unsigned k = 0, iN=N(); k < iN; k++)
|
||||
{
|
||||
analog_net_t *n = m_nets[k];
|
||||
|
||||
const nl_double DD_n = (n->Q_Analog() - m_last_V[k]);
|
||||
const nl_double hn = current_timestep();
|
||||
|
||||
nl_double DD2 = (DD_n / hn - n->m_DD_n_m_1 / n->m_h_n_m_1) / (hn + n->m_h_n_m_1);
|
||||
nl_double new_net_timestep;
|
||||
|
||||
n->m_h_n_m_1 = hn;
|
||||
n->m_DD_n_m_1 = DD_n;
|
||||
if (nl_math::abs(DD2) > NL_FCONST(1e-30)) // avoid div-by-zero
|
||||
new_net_timestep = nl_math::sqrt(m_params.m_lte / nl_math::abs(NL_FCONST(0.5)*DD2));
|
||||
else
|
||||
new_net_timestep = m_params.m_max_timestep;
|
||||
|
||||
if (new_net_timestep < new_solver_timestep)
|
||||
new_solver_timestep = new_net_timestep;
|
||||
|
||||
m_last_V[k] = n->Q_Analog();
|
||||
}
|
||||
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;
|
||||
return netlist_time::from_double(new_solver_timestep);
|
||||
}
|
||||
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
ATTR_COLD void matrix_solver_sm_t<m_N, _storage_N>::add_term(int k, terminal_t *term)
|
||||
{
|
||||
if (term->m_otherterm->net().isRailNet())
|
||||
{
|
||||
m_rails_temp[k].add(term, -1, false);
|
||||
}
|
||||
else
|
||||
{
|
||||
int ot = get_net_idx(&term->m_otherterm->net());
|
||||
if (ot>=0)
|
||||
{
|
||||
m_terms[k]->add(term, ot, true);
|
||||
}
|
||||
/* Should this be allowed ? */
|
||||
else // if (ot<0)
|
||||
{
|
||||
m_rails_temp[k].add(term, ot, true);
|
||||
log().fatal("found term with missing othernet {1}\n", term->name());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
ATTR_COLD void matrix_solver_sm_t<m_N, _storage_N>::vsetup(analog_net_t::list_t &nets)
|
||||
{
|
||||
if (m_dim < nets.size())
|
||||
log().fatal("Dimension {1} less than {2}", m_dim, nets.size());
|
||||
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
m_terms[k]->clear();
|
||||
m_rails_temp[k].clear();
|
||||
}
|
||||
|
||||
matrix_solver_t::setup_base(nets);
|
||||
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
m_terms[k]->m_railstart = m_terms[k]->count();
|
||||
for (unsigned i = 0; i < m_rails_temp[k].count(); i++)
|
||||
this->m_terms[k]->add(m_rails_temp[k].terms()[i], m_rails_temp[k].net_other()[i], false);
|
||||
for (unsigned i = 0; i < m_rails_temp[k]->count(); i++)
|
||||
this->m_terms[k]->add(m_rails_temp[k]->terms()[i], m_rails_temp[k]->net_other()[i], false);
|
||||
|
||||
m_rails_temp[k].clear(); // no longer needed
|
||||
m_rails_temp[k]->clear(); // no longer needed
|
||||
m_terms[k]->set_pointers();
|
||||
}
|
||||
|
||||
@ -315,13 +227,13 @@ ATTR_COLD void matrix_solver_sm_t<m_N, _storage_N>::vsetup(analog_net_t::list_t
|
||||
* save states
|
||||
*/
|
||||
save(NLNAME(m_last_RHS));
|
||||
save(NLNAME(m_last_V));
|
||||
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
pstring num = pfmt("{1}")(k);
|
||||
|
||||
save(RHS(k), "RHS" + num);
|
||||
save(m_terms[k]->m_last_V, "lastV" + num);
|
||||
|
||||
save(m_terms[k]->go(),"GO" + num, m_terms[k]->count());
|
||||
save(m_terms[k]->gt(),"GT" + num, m_terms[k]->count());
|
||||
@ -468,35 +380,6 @@ void matrix_solver_sm_t<m_N, _storage_N>::LE_compute_x(
|
||||
}
|
||||
|
||||
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
template <typename T>
|
||||
T matrix_solver_sm_t<m_N, _storage_N>::delta(
|
||||
const T * RESTRICT V)
|
||||
{
|
||||
/* FIXME: Ideally we should also include currents (RHS) here. This would
|
||||
* need a revaluation of the right hand side after voltages have been updated
|
||||
* and thus belong into a different calculation. This applies to all solvers.
|
||||
*/
|
||||
|
||||
const unsigned iN = this->N();
|
||||
T cerr = 0;
|
||||
for (unsigned i = 0; i < iN; i++)
|
||||
cerr = std::fmax(cerr, nl_math::abs(V[i] - (T) this->m_nets[i]->m_cur_Analog));
|
||||
return cerr;
|
||||
}
|
||||
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
template <typename T>
|
||||
void matrix_solver_sm_t<m_N, _storage_N>::store(
|
||||
const T * RESTRICT V)
|
||||
{
|
||||
for (unsigned i = 0, iN=N(); i < iN; i++)
|
||||
{
|
||||
this->m_nets[i]->m_cur_Analog = V[i];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
int matrix_solver_sm_t<m_N, _storage_N>::solve_non_dynamic(ATTR_UNUSED const bool newton_raphson)
|
||||
{
|
||||
@ -605,15 +488,12 @@ matrix_solver_sm_t<m_N, _storage_N>::matrix_solver_sm_t(const solver_parameters_
|
||||
: matrix_solver_t(GAUSSIAN_ELIMINATION, params)
|
||||
, m_dim(size)
|
||||
{
|
||||
m_rails_temp = palloc_array(terms_t, N());
|
||||
#if (NL_USE_DYNAMIC_ALLOCATION)
|
||||
m_A = palloc_array(nl_ext_double, N() * m_pitch);
|
||||
#endif
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
m_terms[k] = palloc(terms_t);
|
||||
m_last_RHS[k] = 0.0;
|
||||
m_last_V[k] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
@ -622,15 +502,12 @@ matrix_solver_sm_t<m_N, _storage_N>::matrix_solver_sm_t(const eSolverType type,
|
||||
: matrix_solver_t(type, params)
|
||||
, m_dim(size)
|
||||
{
|
||||
m_rails_temp = palloc_array(terms_t, N());
|
||||
#if (NL_USE_DYNAMIC_ALLOCATION)
|
||||
m_A = palloc_array(nl_ext_double, N() * m_pitch);
|
||||
#endif
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
m_terms[k] = palloc(terms_t);
|
||||
m_last_RHS[k] = 0.0;
|
||||
m_last_V[k] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -66,7 +66,6 @@ public:
|
||||
virtual void reset() override { matrix_solver_t::reset(); }
|
||||
|
||||
protected:
|
||||
virtual void add_term(int net_idx, terminal_t *term) override;
|
||||
virtual int vsolve_non_dynamic(const bool newton_raphson) override;
|
||||
int solve_non_dynamic(const bool newton_raphson);
|
||||
|
||||
@ -79,13 +78,6 @@ protected:
|
||||
template <typename T>
|
||||
void LE_compute_x(T * RESTRICT x);
|
||||
|
||||
template <typename T>
|
||||
T delta(const T * RESTRICT V);
|
||||
|
||||
template <typename T>
|
||||
void store(const T * RESTRICT V);
|
||||
|
||||
virtual netlist_time compute_next_timestep() override;
|
||||
|
||||
template <typename T1, typename T2>
|
||||
inline nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r][c]; }
|
||||
@ -103,10 +95,6 @@ protected:
|
||||
inline nl_ext_double &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; }
|
||||
|
||||
ATTR_ALIGN nl_double m_last_RHS[_storage_N]; // right hand side - contains currents
|
||||
ATTR_ALIGN nl_double m_last_V[_storage_N];
|
||||
|
||||
terms_t * m_terms[_storage_N];
|
||||
terms_t *m_rails_temp;
|
||||
|
||||
private:
|
||||
static const std::size_t m_pitch = ((( _storage_N) + 7) / 8) * 8;
|
||||
@ -138,102 +126,26 @@ private:
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
matrix_solver_w_t<m_N, _storage_N>::~matrix_solver_w_t()
|
||||
{
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
pfree(m_terms[k]);
|
||||
}
|
||||
pfree_array(m_rails_temp);
|
||||
#if (NL_USE_DYNAMIC_ALLOCATION)
|
||||
pfree_array(m_A);
|
||||
#endif
|
||||
}
|
||||
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
netlist_time matrix_solver_w_t<m_N, _storage_N>::compute_next_timestep()
|
||||
{
|
||||
nl_double new_solver_timestep = m_params.m_max_timestep;
|
||||
|
||||
if (m_params.m_dynamic)
|
||||
{
|
||||
/*
|
||||
* FIXME: We should extend the logic to use either all nets or
|
||||
* only output nets.
|
||||
*/
|
||||
for (unsigned k = 0, iN=N(); k < iN; k++)
|
||||
{
|
||||
analog_net_t *n = m_nets[k];
|
||||
|
||||
const nl_double DD_n = (n->Q_Analog() - m_last_V[k]);
|
||||
const nl_double hn = current_timestep();
|
||||
|
||||
nl_double DD2 = (DD_n / hn - n->m_DD_n_m_1 / n->m_h_n_m_1) / (hn + n->m_h_n_m_1);
|
||||
nl_double new_net_timestep;
|
||||
|
||||
n->m_h_n_m_1 = hn;
|
||||
n->m_DD_n_m_1 = DD_n;
|
||||
if (nl_math::abs(DD2) > NL_FCONST(1e-30)) // avoid div-by-zero
|
||||
new_net_timestep = nl_math::sqrt(m_params.m_lte / nl_math::abs(NL_FCONST(0.5)*DD2));
|
||||
else
|
||||
new_net_timestep = m_params.m_max_timestep;
|
||||
|
||||
if (new_net_timestep < new_solver_timestep)
|
||||
new_solver_timestep = new_net_timestep;
|
||||
|
||||
m_last_V[k] = n->Q_Analog();
|
||||
}
|
||||
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;
|
||||
return netlist_time::from_double(new_solver_timestep);
|
||||
}
|
||||
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
ATTR_COLD void matrix_solver_w_t<m_N, _storage_N>::add_term(int k, terminal_t *term)
|
||||
{
|
||||
if (term->m_otherterm->net().isRailNet())
|
||||
{
|
||||
m_rails_temp[k].add(term, -1, false);
|
||||
}
|
||||
else
|
||||
{
|
||||
int ot = get_net_idx(&term->m_otherterm->net());
|
||||
if (ot>=0)
|
||||
{
|
||||
m_terms[k]->add(term, ot, true);
|
||||
}
|
||||
/* Should this be allowed ? */
|
||||
else // if (ot<0)
|
||||
{
|
||||
m_rails_temp[k].add(term, ot, true);
|
||||
log().fatal("found term with missing othernet {1}\n", term->name());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
ATTR_COLD void matrix_solver_w_t<m_N, _storage_N>::vsetup(analog_net_t::list_t &nets)
|
||||
{
|
||||
if (m_dim < nets.size())
|
||||
log().fatal("Dimension {1} less than {2}", m_dim, nets.size());
|
||||
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
m_terms[k]->clear();
|
||||
m_rails_temp[k].clear();
|
||||
}
|
||||
|
||||
matrix_solver_t::setup_base(nets);
|
||||
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
m_terms[k]->m_railstart = m_terms[k]->count();
|
||||
for (unsigned i = 0; i < m_rails_temp[k].count(); i++)
|
||||
this->m_terms[k]->add(m_rails_temp[k].terms()[i], m_rails_temp[k].net_other()[i], false);
|
||||
for (unsigned i = 0; i < m_rails_temp[k]->count(); i++)
|
||||
this->m_terms[k]->add(m_rails_temp[k]->terms()[i], m_rails_temp[k]->net_other()[i], false);
|
||||
|
||||
m_rails_temp[k].clear(); // no longer needed
|
||||
m_rails_temp[k]->clear(); // no longer needed
|
||||
m_terms[k]->set_pointers();
|
||||
}
|
||||
|
||||
@ -329,13 +241,13 @@ ATTR_COLD void matrix_solver_w_t<m_N, _storage_N>::vsetup(analog_net_t::list_t &
|
||||
* save states
|
||||
*/
|
||||
save(NLNAME(m_last_RHS));
|
||||
save(NLNAME(m_last_V));
|
||||
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
pstring num = pfmt("{1}")(k);
|
||||
|
||||
save(RHS(k), "RHS" + num);
|
||||
save(m_terms[k]->m_last_V, "lastV" + num);
|
||||
|
||||
save(m_terms[k]->go(),"GO" + num, m_terms[k]->count());
|
||||
save(m_terms[k]->gt(),"GT" + num, m_terms[k]->count());
|
||||
@ -383,15 +295,15 @@ void matrix_solver_w_t<m_N, _storage_N>::build_LE_RHS()
|
||||
nl_double rhsk_a = 0.0;
|
||||
nl_double rhsk_b = 0.0;
|
||||
|
||||
const int terms_count = m_terms[k]->count();
|
||||
const unsigned terms_count = m_terms[k]->count();
|
||||
const nl_double * RESTRICT go = m_terms[k]->go();
|
||||
const nl_double * RESTRICT Idr = m_terms[k]->Idr();
|
||||
const nl_double * const * RESTRICT other_cur_analog = m_terms[k]->other_curanalog();
|
||||
|
||||
for (int i = 0; i < terms_count; i++)
|
||||
for (unsigned i = 0; i < terms_count; i++)
|
||||
rhsk_a = rhsk_a + Idr[i];
|
||||
|
||||
for (int i = m_terms[k]->m_railstart; i < terms_count; i++)
|
||||
for (unsigned i = m_terms[k]->m_railstart; i < terms_count; i++)
|
||||
//rhsk = rhsk + go[i] * terms[i]->m_otherterm->net().as_analog().Q_Analog();
|
||||
rhsk_b = rhsk_b + go[i] * *other_cur_analog[i];
|
||||
|
||||
@ -481,34 +393,6 @@ void matrix_solver_w_t<m_N, _storage_N>::LE_compute_x(
|
||||
}
|
||||
|
||||
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
template <typename T>
|
||||
T matrix_solver_w_t<m_N, _storage_N>::delta(
|
||||
const T * RESTRICT V)
|
||||
{
|
||||
/* FIXME: Ideally we should also include currents (RHS) here. This would
|
||||
* need a revaluation of the right hand side after voltages have been updated
|
||||
* and thus belong into a different calculation. This applies to all solvers.
|
||||
*/
|
||||
|
||||
const unsigned iN = this->N();
|
||||
T cerr = 0;
|
||||
for (unsigned i = 0; i < iN; i++)
|
||||
cerr = std::fmax(cerr, nl_math::abs(V[i] - (T) this->m_nets[i]->m_cur_Analog));
|
||||
return cerr;
|
||||
}
|
||||
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
template <typename T>
|
||||
void matrix_solver_w_t<m_N, _storage_N>::store(
|
||||
const T * RESTRICT V)
|
||||
{
|
||||
for (unsigned i = 0, iN=N(); i < iN; i++)
|
||||
{
|
||||
this->m_nets[i]->m_cur_Analog = V[i];
|
||||
}
|
||||
}
|
||||
|
||||
template <unsigned m_N, unsigned _storage_N>
|
||||
int matrix_solver_w_t<m_N, _storage_N>::solve_non_dynamic(ATTR_UNUSED const bool newton_raphson)
|
||||
{
|
||||
@ -668,15 +552,12 @@ matrix_solver_w_t<m_N, _storage_N>::matrix_solver_w_t(const solver_parameters_t
|
||||
,m_cnt(0)
|
||||
, m_dim(size)
|
||||
{
|
||||
m_rails_temp = palloc_array(terms_t, N());
|
||||
#if (NL_USE_DYNAMIC_ALLOCATION)
|
||||
m_A = palloc_array(nl_ext_double, N() * m_pitch);
|
||||
#endif
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
m_terms[k] = palloc(terms_t);
|
||||
m_last_RHS[k] = 0.0;
|
||||
m_last_V[k] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
@ -686,12 +567,9 @@ matrix_solver_w_t<m_N, _storage_N>::matrix_solver_w_t(const eSolverType type, co
|
||||
,m_cnt(0)
|
||||
, m_dim(size)
|
||||
{
|
||||
m_rails_temp = palloc_array(terms_t, N());
|
||||
for (unsigned k = 0; k < N(); k++)
|
||||
{
|
||||
m_terms[k] = palloc(terms_t);
|
||||
m_last_RHS[k] = 0.0;
|
||||
m_last_V[k] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -108,6 +108,11 @@ ATTR_COLD matrix_solver_t::matrix_solver_t(const eSolverType type, const solver_
|
||||
ATTR_COLD matrix_solver_t::~matrix_solver_t()
|
||||
{
|
||||
m_inps.clear_and_free();
|
||||
for (unsigned k = 0; k < m_terms.size(); k++)
|
||||
{
|
||||
pfree(m_terms[k]);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
ATTR_COLD void matrix_solver_t::setup_base(analog_net_t::list_t &nets)
|
||||
@ -115,9 +120,14 @@ ATTR_COLD void matrix_solver_t::setup_base(analog_net_t::list_t &nets)
|
||||
log().debug("New solver setup\n");
|
||||
|
||||
m_nets.clear();
|
||||
m_terms.clear();
|
||||
|
||||
for (auto & net : nets)
|
||||
{
|
||||
m_nets.push_back(net);
|
||||
m_terms.push_back(palloc(terms_t));
|
||||
m_rails_temp.push_back(palloc(terms_t));
|
||||
}
|
||||
|
||||
for (std::size_t k = 0; k < nets.size(); k++)
|
||||
{
|
||||
@ -308,6 +318,70 @@ ATTR_COLD int matrix_solver_t::get_net_idx(net_t *net)
|
||||
return -1;
|
||||
}
|
||||
|
||||
ATTR_COLD void matrix_solver_t::add_term(int k, terminal_t *term)
|
||||
{
|
||||
if (term->m_otherterm->net().isRailNet())
|
||||
{
|
||||
m_rails_temp[k]->add(term, -1, false);
|
||||
}
|
||||
else
|
||||
{
|
||||
int ot = get_net_idx(&term->m_otherterm->net());
|
||||
if (ot>=0)
|
||||
{
|
||||
m_terms[k]->add(term, ot, true);
|
||||
}
|
||||
/* Should this be allowed ? */
|
||||
else // if (ot<0)
|
||||
{
|
||||
m_rails_temp[k]->add(term, ot, true);
|
||||
log().fatal("found term with missing othernet {1}\n", term->name());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
netlist_time matrix_solver_t::compute_next_timestep()
|
||||
{
|
||||
nl_double new_solver_timestep = m_params.m_max_timestep;
|
||||
|
||||
if (m_params.m_dynamic)
|
||||
{
|
||||
/*
|
||||
* FIXME: We should extend the logic to use either all nets or
|
||||
* only output nets.
|
||||
*/
|
||||
for (unsigned k = 0, iN=m_terms.size(); k < iN; k++)
|
||||
{
|
||||
analog_net_t *n = m_nets[k];
|
||||
|
||||
const nl_double DD_n = (n->Q_Analog() - m_terms[k]->m_last_V);
|
||||
const nl_double hn = current_timestep();
|
||||
|
||||
nl_double DD2 = (DD_n / hn - n->m_DD_n_m_1 / n->m_h_n_m_1) / (hn + n->m_h_n_m_1);
|
||||
nl_double new_net_timestep;
|
||||
|
||||
n->m_h_n_m_1 = hn;
|
||||
n->m_DD_n_m_1 = DD_n;
|
||||
if (nl_math::abs(DD2) > NL_FCONST(1e-30)) // avoid div-by-zero
|
||||
new_net_timestep = nl_math::sqrt(m_params.m_lte / nl_math::abs(NL_FCONST(0.5)*DD2));
|
||||
else
|
||||
new_net_timestep = m_params.m_max_timestep;
|
||||
|
||||
if (new_net_timestep < new_solver_timestep)
|
||||
new_solver_timestep = new_net_timestep;
|
||||
|
||||
m_terms[k]->m_last_V = n->Q_Analog();
|
||||
}
|
||||
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;
|
||||
return netlist_time::from_double(new_solver_timestep);
|
||||
}
|
||||
|
||||
|
||||
|
||||
void matrix_solver_t::log_stats()
|
||||
{
|
||||
if (this->m_stat_calculations != 0 && this->m_stat_vsolver_calls && this->m_params.m_log_stats)
|
||||
|
||||
Loading…
Reference in New Issue
Block a user