From 4f1ca77643e461a0dfb648d77edf57faa8c4d952 Mon Sep 17 00:00:00 2001 From: couriersud Date: Mon, 11 Apr 2016 00:50:45 +0200 Subject: [PATCH] Moved solver members to proper place. Minor code changes. (nw) --- nl_examples/kidniki.c | 4 +- src/lib/netlist/nl_base.cpp | 4 -- src/lib/netlist/nl_base.h | 2 - src/lib/netlist/solver/nld_matrix_solver.h | 16 +++++-- src/lib/netlist/solver/nld_ms_direct.h | 54 +++++++++------------- src/lib/netlist/solver/nld_ms_sm.h | 4 +- src/lib/netlist/solver/nld_ms_sor_mat.h | 2 +- src/lib/netlist/solver/nld_ms_w.h | 6 +-- src/lib/netlist/solver/nld_solver.cpp | 11 +++-- 9 files changed, 46 insertions(+), 57 deletions(-) diff --git a/nl_examples/kidniki.c b/nl_examples/kidniki.c index c4ae39f68d9..49d666dbc4d 100644 --- a/nl_examples/kidniki.c +++ b/nl_examples/kidniki.c @@ -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) diff --git a/src/lib/netlist/nl_base.cpp b/src/lib/netlist/nl_base.cpp index eb83749c79b..16bba23ed3c 100644 --- a/src/lib/netlist/nl_base.cpp +++ b/src/lib/netlist/nl_base.cpp @@ -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(); } diff --git a/src/lib/netlist/nl_base.h b/src/lib/netlist/nl_base.h index 3b6eb7b75e1..b8239addcf6 100644 --- a/src/lib/netlist/nl_base.h +++ b/src/lib/netlist/nl_base.h @@ -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; diff --git a/src/lib/netlist/solver/nld_matrix_solver.h b/src/lib/netlist/solver/nld_matrix_solver.h index 58947144ff5..96013dc04de 100644 --- a/src/lib/netlist/solver/nld_matrix_solver.h +++ b/src/lib/netlist/solver/nld_matrix_solver.h @@ -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 m_nz; /* all non zero for multiplication */ - pvector_t m_nzrd; /* non zero right of the diagonal for elimination, may include RHS element */ - pvector_t m_nzbd; /* non zero below of the diagonal for elimination */ + pvector_t m_nz; /* all non zero for multiplication */ + pvector_t m_nzrd; /* non zero right of the diagonal for elimination, may include RHS element */ + pvector_t 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 m_term; pvector_t m_net_other; pvector_t m_go; pvector_t m_gt; pvector_t m_Idr; pvector_t m_other_curanalog; + pvector_t m_term; }; diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h index 2b0ae80b855..227eed3ac54 100644 --- a/src/lib/netlist/solver/nld_ms_direct.h +++ b/src/lib/netlist/solver/nld_ms_direct.h @@ -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 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::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::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::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::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::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 0) @@ -468,7 +472,6 @@ void matrix_solver_direct_t::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::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::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); diff --git a/src/lib/netlist/solver/nld_ms_sm.h b/src/lib/netlist/solver/nld_ms_sm.h index bff8197127c..826ea7c216f 100644 --- a/src/lib/netlist/solver/nld_ms_sm.h +++ b/src/lib/netlist/solver/nld_ms_sm.h @@ -263,12 +263,12 @@ void matrix_solver_sm_t::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++) { diff --git a/src/lib/netlist/solver/nld_ms_sor_mat.h b/src/lib/netlist/solver/nld_ms_sor_mat.h index 10192717a90..70367f95593 100644 --- a/src/lib/netlist/solver/nld_ms_sor_mat.h +++ b/src/lib/netlist/solver/nld_ms_sor_mat.h @@ -178,7 +178,7 @@ int matrix_solver_SOR_mat_t::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++) diff --git a/src/lib/netlist/solver/nld_ms_w.h b/src/lib/netlist/solver/nld_ms_w.h index aa6f7104c6c..43db5b55e69 100644 --- a/src/lib/netlist/solver/nld_ms_w.h +++ b/src/lib/netlist/solver/nld_ms_w.h @@ -276,16 +276,16 @@ void matrix_solver_w_t::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) { diff --git a/src/lib/netlist/solver/nld_solver.cpp b/src/lib/netlist/solver/nld_solver.cpp index 5dd0af53313..3d2f9aeff05 100644 --- a/src/lib/netlist/solver/nld_solver.cpp +++ b/src/lib/netlist/solver/nld_solver.cpp @@ -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;