Moved solver members to proper place. Minor code changes. (nw)

This commit is contained in:
couriersud 2016-04-11 00:50:45 +02:00
parent 65e139be2f
commit 4f1ca77643
9 changed files with 46 additions and 57 deletions

View File

@ -4,7 +4,7 @@
#include "netlist/devices/nld_system.h"
#include "netlist/analog/nld_bjt.h"
#define USE_FRONTIERS 1
#define USE_FRONTIERS 0
#define USE_FIXED_STV 1
NETLIST_START(dummy)
@ -23,7 +23,7 @@ NETLIST_START(dummy)
PARAM(Solver.MIN_TIMESTEP, 10e-6)
PARAM(Solver.PARALLEL, 0)
PARAM(Solver.SOR_FACTOR, 1.00)
PARAM(Solver.PIVOT, 0)
PARAM(Solver.PIVOT, 1)
#else
SOLVER(Solver, 12000)
PARAM(Solver.ACCURACY, 1e-8)

View File

@ -784,8 +784,6 @@ ATTR_COLD void logic_net_t::save_register()
ATTR_COLD analog_net_t::analog_net_t()
: net_t(ANALOG)
, m_DD_n_m_1(0.0)
, m_h_n_m_1(1e-6)
, m_solver(NULL)
{
}
@ -797,8 +795,6 @@ ATTR_COLD void analog_net_t::reset()
ATTR_COLD void analog_net_t::save_register()
{
save(NLNAME(m_DD_n_m_1));
save(NLNAME(m_h_n_m_1));
net_t::save_register();
}

View File

@ -866,8 +866,6 @@ namespace netlist
private:
public:
nl_double m_DD_n_m_1;
nl_double m_h_n_m_1;
//FIXME: needed by current solver code
devices::matrix_solver_t *m_solver;

View File

@ -19,7 +19,11 @@ class terms_t
P_PREVENT_COPYING(terms_t)
public:
ATTR_COLD terms_t() : m_railstart(0), m_last_V(0)
ATTR_COLD terms_t()
: m_railstart(0)
, m_last_V(0.0)
, m_DD_n_m_1(0.0)
, m_h_n_m_1(1e-6)
{}
ATTR_COLD void clear()
@ -47,20 +51,22 @@ public:
unsigned m_railstart;
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 */
pvector_t<int> m_nz; /* all non zero for multiplication */
pvector_t<int> m_nzrd; /* non zero right of the diagonal for elimination, may include RHS element */
pvector_t<int> m_nzbd; /* non zero below of the diagonal for elimination */
/* state */
nl_double m_last_V;
nl_double m_DD_n_m_1;
nl_double m_h_n_m_1;
private:
pvector_t<terminal_t *> m_term;
pvector_t<int> m_net_other;
pvector_t<nl_double> m_go;
pvector_t<nl_double> m_gt;
pvector_t<nl_double> m_Idr;
pvector_t<nl_double *> m_other_curanalog;
pvector_t<terminal_t *> m_term;
};

View File

@ -34,7 +34,7 @@ NETLIB_NAMESPACE_DEVICES_START()
#if TEST_PARALLEL
#define MAXTHR 10
static const int num_thr = 2;
static const int num_thr = 1;
struct thr_intf
{
@ -46,7 +46,7 @@ struct ti_t
volatile std::atomic<int> lo;
thr_intf *intf;
void *params;
int _block[29]; /* make it 256 bytes */
// int _block[29]; /* make it 256 bytes */
};
static ti_t ti[MAXTHR];
@ -315,20 +315,27 @@ ATTR_COLD void matrix_solver_direct_t<m_N, _storage_N>::vsetup(analog_net_t::lis
touched[k][m_terms[k]->m_nz[j]] = true;
}
unsigned ops = 0;
for (unsigned k = 0; k < N(); k++)
{
ops++; // 1/A(k,k)
for (unsigned row = k + 1; row < N(); row++)
{
if (touched[row][k])
{
ops++;
if (!m_terms[k]->m_nzbd.contains(row))
m_terms[k]->m_nzbd.push_back(row);
for (unsigned col = k; col < N(); col++)
if (touched[k][col])
{
touched[row][col] = true;
ops += 2;
}
}
}
}
log().verbose("Number of mults/adds for {1}: {2}", name(), ops);
if (0)
for (unsigned k = 0; k < N(); k++)
@ -348,8 +355,10 @@ ATTR_COLD void matrix_solver_direct_t<m_N, _storage_N>::vsetup(analog_net_t::lis
{
pstring num = pfmt("{1}")(k);
save(RHS(k), "RHS" + num);
save(m_terms[k]->m_last_V, "lastV" + num);
save(RHS(k), "RHS." + num);
save(m_terms[k]->m_last_V, "lastV." + num);
save(m_terms[k]->m_DD_n_m_1, "m_DD_n_m_1." + num);
save(m_terms[k]->m_h_n_m_1, "m_h_n_m_1." + num);
save(m_terms[k]->go(),"GO" + num, m_terms[k]->count());
save(m_terms[k]->gt(),"GT" + num, m_terms[k]->count());
@ -369,14 +378,10 @@ void matrix_solver_direct_t<m_N, _storage_N>::do_work(const int id, void *param)
const nl_double f = 1.0 / A(i,i);
const unsigned * RESTRICT const p = m_terms[i]->m_nzrd.data();
const unsigned e = m_terms[i]->m_nzrd.size();
//nl_double A_cache[128];
//for (unsigned k = 0; k < e; k++)
// A_cache[k] = A(i,p[k]);
/* Eliminate column i from row j */
const unsigned * RESTRICT const pb = m_terms[i]->m_nzbd.data();
//const unsigned eb = m_terms[i]->m_nzbd.size();
const unsigned sj = x_start[id];
const unsigned se = x_stop[id];
for (unsigned jb = sj; jb < se; jb++)
@ -384,9 +389,7 @@ void matrix_solver_direct_t<m_N, _storage_N>::do_work(const int id, void *param)
const unsigned j = pb[jb];
const nl_double f1 = - A(j,i) * f;
for (unsigned k = 0; k < e; k++)
//A(j,p[k]) += A_cache[k] * f1;
A(j,p[k]) += A(i,p[k]) * f1;
//RHS(j) += RHS(i) * f1;
}
}
#endif
@ -449,15 +452,16 @@ void matrix_solver_direct_t<m_N, _storage_N>::LE_solve()
if (eb > 0)
{
//printf("here %d\n", eb);
unsigned chunks = (eb) / (num_thr + 1);
unsigned chunks = (eb + num_thr) / (num_thr + 1);
for (int p=0; p < num_thr + 1; p++)
{
x_i[p] = i;
x_start[p] = chunks * p;
x_stop[p] = std::min(chunks*(p+1), eb);
if (p<num_thr) thr_process(p, this, NULL);
if (p<num_thr && x_start[p] < x_stop[p]) thr_process(p, this, NULL);
}
do_work(num_thr, NULL);
if (x_start[num_thr] < x_stop[num_thr])
do_work(num_thr, NULL);
thr_wait();
}
else if (eb > 0)
@ -468,7 +472,6 @@ void matrix_solver_direct_t<m_N, _storage_N>::LE_solve()
do_work(0, NULL);
}
#else
#if 0
/* FIXME: Singular matrix? */
const nl_double f = 1.0 / A(i,i);
const auto &nzrd = m_terms[i]->m_nzrd;
@ -478,27 +481,12 @@ void matrix_solver_direct_t<m_N, _storage_N>::LE_solve()
for (auto & j : nzbd)
{
//__builtin_prefetch((const void*)(&A(j+2,0)),0,0);
const nl_double f1 = - A(j,i) * f;
for (auto & k : nzrd)
A(j,k) += A(i,k) * f1;
//RHS(j) += RHS(i) * f1;
}
#else
/* FIXME: Singular matrix? */
const nl_double f = 1.0 / A(i,i);
const auto &nzrd = m_terms[i]->m_nzrd;
const auto &nzbd = m_terms[i]->m_nzbd;
/* Eliminate column i from row j */
for (auto & j : nzbd)
{
const nl_double f1 = - A(j,i) * f;
for (auto & k : nzrd)
A(j,k) += A(i,k) * f1;
//RHS(j) += RHS(i) * f1;
}
#endif
#endif
}
}
@ -528,12 +516,12 @@ void matrix_solver_direct_t<m_N, _storage_N>::LE_back_subst(
{
T tmp = 0;
const unsigned *p = m_terms[j]->m_nzrd.data();
const unsigned e = m_terms[j]->m_nzrd.size() - 1; /* exclude RHS element */
const auto *p = m_terms[j]->m_nzrd.data();
const auto e = m_terms[j]->m_nzrd.size() - 1; /* exclude RHS element */
for (unsigned k = 0; k < e; k++)
{
const unsigned pk = p[k];
const auto pk = p[k];
tmp += A(j,pk) * x[pk];
}
x[j] = (RHS(j) - tmp) / A(j,j);

View File

@ -263,12 +263,12 @@ void matrix_solver_sm_t<m_N, _storage_N>::LE_invert()
{
/* FIXME: Singular matrix? */
const nl_double f = 1.0 / W(i,i);
const unsigned * RESTRICT const p = m_terms[i]->m_nzrd.data();
const auto * RESTRICT const p = m_terms[i]->m_nzrd.data();
const unsigned e = m_terms[i]->m_nzrd.size();
/* Eliminate column i from row j */
const unsigned * RESTRICT const pb = m_terms[i]->m_nzbd.data();
const auto * RESTRICT const pb = m_terms[i]->m_nzbd.data();
const unsigned eb = m_terms[i]->m_nzbd.size();
for (unsigned jb = 0; jb < eb; jb++)
{

View File

@ -178,7 +178,7 @@ int matrix_solver_SOR_mat_t<m_N, _storage_N>::vsolve_non_dynamic(const bool newt
{
nl_double Idrive = 0;
const unsigned *p = this->m_terms[k]->m_nz.data();
const auto *p = this->m_terms[k]->m_nz.data();
const unsigned e = this->m_terms[k]->m_nz.size();
for (unsigned i = 0; i < e; i++)

View File

@ -276,16 +276,16 @@ void matrix_solver_w_t<m_N, _storage_N>::LE_invert()
{
/* FIXME: Singular matrix? */
const nl_double f = 1.0 / W(i,i);
const unsigned * RESTRICT const p = m_terms[i]->m_nzrd.data();
const auto * RESTRICT const p = m_terms[i]->m_nzrd.data();
const unsigned e = m_terms[i]->m_nzrd.size();
/* Eliminate column i from row j */
const unsigned * RESTRICT const pb = m_terms[i]->m_nzbd.data();
const auto * RESTRICT const pb = m_terms[i]->m_nzbd.data();
const unsigned eb = m_terms[i]->m_nzbd.size();
for (unsigned jb = 0; jb < eb; jb++)
{
const unsigned j = pb[jb];
const auto j = pb[jb];
const nl_double f1 = - W(j,i) * f;
if (f1 != 0.0)
{

View File

@ -356,15 +356,16 @@ netlist_time matrix_solver_t::compute_next_timestep()
for (unsigned k = 0, iN=m_terms.size(); k < iN; k++)
{
analog_net_t *n = m_nets[k];
terms_t *t = m_terms[k];
const nl_double DD_n = (n->Q_Analog() - m_terms[k]->m_last_V);
const nl_double DD_n = (n->Q_Analog() - t->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 DD2 = (DD_n / hn - t->m_DD_n_m_1 / t->m_h_n_m_1) / (hn + t->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;
t->m_h_n_m_1 = hn;
t->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
@ -373,7 +374,7 @@ netlist_time matrix_solver_t::compute_next_timestep()
if (new_net_timestep < new_solver_timestep)
new_solver_timestep = new_net_timestep;
m_terms[k]->m_last_V = n->Q_Analog();
t->m_last_V = n->Q_Analog();
}
if (new_solver_timestep < m_params.m_min_timestep)
new_solver_timestep = m_params.m_min_timestep;