summaryrefslogtreecommitdiffstatshomepage
path: root/src
diff options
context:
space:
mode:
author couriersud <couriersud@arcor.de>2015-06-10 21:49:52 +0200
committer couriersud <couriersud@arcor.de>2015-06-10 21:50:10 +0200
commitdad7a6bab5201722b258ac8913f0ee2accd1b5e5 (patch)
tree58a5f2b0d9a889bdc80f6564b78890b4a09fa196 /src
parentb2a7d75ef093a5b49e34984e34beac576a546539 (diff)
From 45% to 60% to 99%. That's the improvement achieved for a 89x89
audio matrix mostly solved by elimination. Cleaned up some code as well. [Couriersud]
Diffstat (limited to 'src')
-rw-r--r--src/emu/netlist/analog/nld_ms_direct.h169
-rw-r--r--src/emu/netlist/analog/nld_ms_direct1.h2
-rw-r--r--src/emu/netlist/analog/nld_ms_direct2.h6
-rw-r--r--src/emu/netlist/analog/nld_ms_sor.h45
-rw-r--r--src/emu/netlist/analog/nld_ms_sor_mat.h129
-rw-r--r--src/emu/netlist/analog/nld_solver.h2
-rw-r--r--src/emu/netlist/plib/plists.h18
7 files changed, 201 insertions, 170 deletions
diff --git a/src/emu/netlist/analog/nld_ms_direct.h b/src/emu/netlist/analog/nld_ms_direct.h
index 3c5d404d430..056cf774968 100644
--- a/src/emu/netlist/analog/nld_ms_direct.h
+++ b/src/emu/netlist/analog/nld_ms_direct.h
@@ -36,11 +36,11 @@ protected:
ATTR_HOT int solve_non_dynamic(const bool newton_raphson);
ATTR_HOT void build_LE_A();
- ATTR_HOT void build_LE_RHS();
+ ATTR_HOT void build_LE_RHS(nl_double * RESTRICT rhs);
ATTR_HOT void LE_solve();
ATTR_HOT void LE_back_subst(nl_double * RESTRICT x);
ATTR_HOT nl_double delta(const nl_double * RESTRICT V);
- ATTR_HOT void store(const nl_double * RESTRICT V, const bool store_RHS);
+ ATTR_HOT void store(const nl_double * RESTRICT V);
/* bring the whole system to the current time
* Don't schedule a new calculation time. The recalculation has to be
@@ -88,7 +88,7 @@ ATTR_HOT nl_double netlist_matrix_solver_direct_t<m_N, _storage_N>::compute_next
* FIXME: We should extend the logic to use either all nets or
* only output nets.
*/
- for (unsigned k = 0; k < N(); k++)
+ for (unsigned k = 0, iN=N(); k < iN; k++)
{
netlist_analog_net_t *n = m_nets[k];
@@ -209,7 +209,61 @@ ATTR_COLD void netlist_matrix_solver_direct_t<m_N, _storage_N>::vsetup(netlist_a
#endif
- //ATTR_ALIGN nl_double m_A[_storage_N][((_storage_N + 7) / 8) * 8];
+ /* create a list of non zero elements right of the diagonal
+ * These list anticipate the population of array elements by
+ * Gaussian elimination.
+ */
+ for (unsigned k = 0; k < N(); k++)
+ {
+ terms_t * t = m_terms[k];
+ /* pretty brutal */
+ int *other = t->net_other();
+
+ t->m_nz.clear();
+
+ if (k==0)
+ t->m_nzrd.clear();
+ else
+ {
+ t->m_nzrd = m_terms[k-1]->m_nzrd;
+ int j=0;
+ while(j < t->m_nzrd.size())
+ {
+ if (t->m_nzrd[j] < k + 1)
+ t->m_nzrd.remove_at(j);
+ else
+ j++;
+ }
+ }
+
+ for (unsigned j = 0; j < N(); j++)
+ {
+ for (unsigned i = 0; i < t->m_railstart; i++)
+ {
+ if (!t->m_nzrd.contains(other[i]) && other[i] >= k + 1)
+ t->m_nzrd.add(other[i]);
+ if (!t->m_nz.contains(other[i]))
+ t->m_nz.add(other[i]);
+ }
+ }
+ psort_list(t->m_nzrd);
+
+ t->m_nz.add(k); // add diagonal
+ psort_list(t->m_nz);
+ }
+
+ if(0)
+ for (unsigned k = 0; k < N(); k++)
+ {
+ printf("%3d: ", k);
+ for (unsigned j = 0; j < m_terms[k]->m_nzrd.size(); j++)
+ printf(" %3d", m_terms[k]->m_nzrd[j]);
+ printf("\n");
+ }
+
+ /*
+ * save states
+ */
save(NLNAME(m_RHS));
save(NLNAME(m_last_RHS));
save(NLNAME(m_last_V));
@@ -229,9 +283,10 @@ ATTR_COLD void netlist_matrix_solver_direct_t<m_N, _storage_N>::vsetup(netlist_a
template <unsigned m_N, unsigned _storage_N>
ATTR_HOT void netlist_matrix_solver_direct_t<m_N, _storage_N>::build_LE_A()
{
- for (unsigned k = 0; k < N(); k++)
+ const unsigned iN = N();
+ for (unsigned k = 0; k < iN; k++)
{
- for (unsigned i=0; i < N(); i++)
+ for (unsigned i=0; i < iN; i++)
m_A[k][i] = 0.0;
nl_double akk = 0.0;
@@ -247,14 +302,15 @@ ATTR_HOT void netlist_matrix_solver_direct_t<m_N, _storage_N>::build_LE_A()
m_A[k][k] += akk;
for (unsigned i = 0; i < railstart; i++)
- m_A[k][net_other[i]] += -go[i];
+ m_A[k][net_other[i]] -= go[i];
}
}
template <unsigned m_N, unsigned _storage_N>
-ATTR_HOT void netlist_matrix_solver_direct_t<m_N, _storage_N>::build_LE_RHS()
+ATTR_HOT void netlist_matrix_solver_direct_t<m_N, _storage_N>::build_LE_RHS(nl_double * RESTRICT rhs)
{
- for (unsigned k = 0; k < N(); k++)
+ const unsigned iN = N();
+ for (unsigned k = 0; k < iN; k++)
{
nl_double rhsk_a = 0.0;
nl_double rhsk_b = 0.0;
@@ -271,7 +327,7 @@ ATTR_HOT void netlist_matrix_solver_direct_t<m_N, _storage_N>::build_LE_RHS()
//rhsk = rhsk + go[i] * terms[i]->m_otherterm->net().as_analog().Q_Analog();
rhsk_b = rhsk_b + go[i] * *other_cur_analog[i];
- m_RHS[k] = rhsk_a + rhsk_b;
+ rhs[k] = rhsk_a + rhsk_b;
}
}
@@ -314,16 +370,38 @@ ATTR_HOT void netlist_matrix_solver_direct_t<m_N, _storage_N>::LE_solve()
/* FIXME: Singular matrix? */
const nl_double f = 1.0 / m_A[i][i];
+ const double * RESTRICT s = &m_A[i][0];
+ const int *p = m_terms[i]->m_nzrd.data();
+ const unsigned e = m_terms[i]->m_nzrd.size();
/* Eliminate column i from row j */
for (unsigned j = i + 1; j < kN; j++)
{
- const nl_double f1 = - m_A[j][i] * f;
+ double * RESTRICT d = &m_A[j][0];
+ const nl_double f1 = - d[i] * f;
if (f1 != NL_FCONST(0.0))
{
- for (unsigned k = i + 1; k < kN; k++)
- m_A[j][k] += m_A[i][k] * f1;
+#if 0
+ /* The code below is 30% faster than the original
+ * implementation which is given here for reference.
+ *
+ * for (unsigned k = i + 1; k < kN; k++)
+ * m_A[j][k] = m_A[j][k] + m_A[i][k] * f1;
+ */
+ double * RESTRICT d = &m_A[j][i+1];
+ const double * RESTRICT s = &m_A[i][i+1];
+ const int e = kN - i - 1;
+
+ for (int k = 0; k < e; k++)
+ d[k] = d[k] + s[k] * f1;
+#else
+ for (unsigned k = 0; k < e; k++)
+ {
+ const unsigned pk = p[k];
+ d[pk] += s[pk] * f1;
+ }
+#endif
m_RHS[j] += m_RHS[i] * f1;
}
}
@@ -341,9 +419,28 @@ ATTR_HOT void netlist_matrix_solver_direct_t<m_N, _storage_N>::LE_back_subst(
{
nl_double tmp = 0;
+#if 1
+#if 0
+ const double * RESTRICT A = &m_A[j][j+1];
+ const double * RESTRICT xp = &x[j+1];
+ const int e = kN - j - 1;
+ for (int k = 0; k < e; k++)
+ tmp += A[k] * xp[k];
+#else
+ const double * RESTRICT A = &m_A[j][0];
+ const int *p = m_terms[j]->m_nzrd.data();
+ const unsigned e = m_terms[j]->m_nzrd.size();
+
+ for (unsigned k = 0; k < e; k++)
+ {
+ const unsigned pk = p[k];
+ tmp += A[pk] * x[pk];
+ }
+#endif
+#else
for (unsigned k = j + 1; k < kN; k++)
tmp += m_A[j][k] * x[k];
-
+#endif
x[j] = (m_RHS[j] - tmp) / m_A[j][j];
}
#if 0
@@ -363,40 +460,32 @@ template <unsigned m_N, unsigned _storage_N>
ATTR_HOT nl_double netlist_matrix_solver_direct_t<m_N, _storage_N>::delta(
const nl_double * 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();
nl_double cerr = 0;
- nl_double cerr2 = 0;
- for (unsigned i = 0; i < this->N(); i++)
- {
- const nl_double e = nl_math::abs(V[i] - this->m_nets[i]->m_cur_Analog);
- const nl_double e2 = nl_math::abs(m_RHS[i] - this->m_last_RHS[i]);
- cerr = (e > cerr ? e : cerr);
- cerr2 = (e2 > cerr2 ? e2 : cerr2);
- }
- // FIXME: Review
- return cerr + cerr2*NL_FCONST(100000.0);
+ for (unsigned i = 0; i < iN; i++)
+ cerr = std::max(cerr, nl_math::abs(V[i] - this->m_nets[i]->m_cur_Analog));
+ return cerr;
}
template <unsigned m_N, unsigned _storage_N>
ATTR_HOT void netlist_matrix_solver_direct_t<m_N, _storage_N>::store(
- const nl_double * RESTRICT V, const bool store_RHS)
+ const nl_double * RESTRICT V)
{
- for (unsigned i = 0; i < this->N(); i++)
+ for (unsigned i = 0, iN=N(); i < iN; i++)
{
this->m_nets[i]->m_cur_Analog = V[i];
}
- if (store_RHS)
- {
- for (unsigned i = 0; i < this->N(); i++)
- {
- this->m_last_RHS[i] = m_RHS[i];
- }
- }
}
template <unsigned m_N, unsigned _storage_N>
ATTR_HOT nl_double netlist_matrix_solver_direct_t<m_N, _storage_N>::vsolve()
{
- solve_base<netlist_matrix_solver_direct_t>(this);
+ this->solve_base(this);
return this->compute_next_timestep();
}
@@ -408,17 +497,17 @@ ATTR_HOT int netlist_matrix_solver_direct_t<m_N, _storage_N>::solve_non_dynamic(
this->LE_back_subst(new_v);
- if (this->is_dynamic())
+ if (newton_raphson)
{
nl_double err = delta(new_v);
- store(new_v, true);
+ store(new_v);
return (err > this->m_params.m_accuracy) ? 2 : 1;
}
else
{
- store(new_v, false); // ==> No need to store RHS
+ store(new_v);
return 1;
}
}
@@ -427,7 +516,11 @@ template <unsigned m_N, unsigned _storage_N>
ATTR_HOT inline int netlist_matrix_solver_direct_t<m_N, _storage_N>::vsolve_non_dynamic(const bool newton_raphson)
{
this->build_LE_A();
- this->build_LE_RHS();
+ this->build_LE_RHS(m_last_RHS);
+
+ for (unsigned i=0, iN=N(); i < iN; i++)
+ m_RHS[i] = m_last_RHS[i];
+
this->LE_solve();
return this->solve_non_dynamic(newton_raphson);
diff --git a/src/emu/netlist/analog/nld_ms_direct1.h b/src/emu/netlist/analog/nld_ms_direct1.h
index 991da037591..cbd2d284c11 100644
--- a/src/emu/netlist/analog/nld_ms_direct1.h
+++ b/src/emu/netlist/analog/nld_ms_direct1.h
@@ -38,7 +38,7 @@ ATTR_HOT inline int netlist_matrix_solver_direct1_t::vsolve_non_dynamic(ATTR_UNU
{
netlist_analog_net_t *net = m_nets[0];
this->build_LE_A();
- this->build_LE_RHS();
+ this->build_LE_RHS(m_RHS);
//NL_VERBOSE_OUT(("%f %f\n", new_val, m_RHS[0] / m_A[0][0]);
nl_double new_val = m_RHS[0] / m_A[0][0];
diff --git a/src/emu/netlist/analog/nld_ms_direct2.h b/src/emu/netlist/analog/nld_ms_direct2.h
index 69df8085c57..2996b9fbf01 100644
--- a/src/emu/netlist/analog/nld_ms_direct2.h
+++ b/src/emu/netlist/analog/nld_ms_direct2.h
@@ -39,7 +39,7 @@ ATTR_HOT nl_double netlist_matrix_solver_direct2_t::vsolve()
ATTR_HOT inline int netlist_matrix_solver_direct2_t::vsolve_non_dynamic(ATTR_UNUSED const bool newton_raphson)
{
build_LE_A();
- build_LE_RHS();
+ build_LE_RHS(m_RHS);
const nl_double a = m_A[0][0];
const nl_double b = m_A[0][1];
@@ -53,13 +53,13 @@ ATTR_HOT inline int netlist_matrix_solver_direct2_t::vsolve_non_dynamic(ATTR_UNU
if (is_dynamic())
{
nl_double err = this->delta(new_val);
- store(new_val, true);
+ store(new_val);
if (err > m_params.m_accuracy )
return 2;
else
return 1;
}
- store(new_val, false);
+ store(new_val);
return 1;
}
diff --git a/src/emu/netlist/analog/nld_ms_sor.h b/src/emu/netlist/analog/nld_ms_sor.h
index 4219363e229..d20b1d3ade1 100644
--- a/src/emu/netlist/analog/nld_ms_sor.h
+++ b/src/emu/netlist/analog/nld_ms_sor.h
@@ -100,8 +100,7 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_t<m_N, _storage_N>::vsolve_non_dyn
* omega = 2.0 / (1.0 + nl_math::sqrt(1-rho))
*/
- const nl_double ws = this->m_params.m_sor; //1.045; //2.0 / (1.0 + /*sin*/(3.14159 * 5.5 / (double) (m_nets.count()+1)));
- //const nl_double ws = 2.0 / (1.0 + sin(3.14159 * 4 / (double) (this->N())));
+ const nl_double ws = this->m_params.m_sor;
ATTR_ALIGN nl_double w[_storage_N];
ATTR_ALIGN nl_double one_m_w[_storage_N];
@@ -114,34 +113,31 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_t<m_N, _storage_N>::vsolve_non_dyn
nl_double gabs_t = 0.0;
nl_double RHS_t = 0.0;
+ const int term_count = this->m_terms[k]->count();
+ const nl_double * const RESTRICT gt = this->m_terms[k]->gt();
+ const nl_double * const RESTRICT go = this->m_terms[k]->go();
+ const nl_double * const RESTRICT Idr = this->m_terms[k]->Idr();
+ const nl_double * const *other_cur_analog = this->m_terms[k]->other_curanalog();
+
new_V[k] = this->m_nets[k]->m_cur_Analog;
+ for (unsigned i = 0; i < term_count; i++)
{
- const int term_count = this->m_terms[k]->count();
- const nl_double * const RESTRICT gt = this->m_terms[k]->gt();
- const nl_double * const RESTRICT go = this->m_terms[k]->go();
- const nl_double * const RESTRICT Idr = this->m_terms[k]->Idr();
- const nl_double * const *other_cur_analog = this->m_terms[k]->other_curanalog();
-
- for (int i = 0; i < term_count; i++)
- {
- gtot_t = gtot_t + gt[i];
- RHS_t = RHS_t + Idr[i];
- }
-
- if (USE_GABS)
- for (int i = 0; i < term_count; i++)
- gabs_t = gabs_t + nl_math::abs(go[i]);
-
- for (int i = this->m_terms[k]->m_railstart; i < term_count; i++)
- RHS_t = RHS_t + go[i] * *other_cur_analog[i];
+ gtot_t = gtot_t + gt[i];
+ RHS_t = RHS_t + Idr[i];
}
+ for (int i = this->m_terms[k]->m_railstart; i < term_count; i++)
+ RHS_t = RHS_t + go[i] * *other_cur_analog[i];
+
RHS[k] = RHS_t;
if (USE_GABS)
{
- gabs_t *= NL_FCONST(0.95); // avoid rounding issues
+ for (int i = 0; i < term_count; i++)
+ gabs_t = gabs_t + nl_math::abs(go[i]);
+
+ gabs_t *= NL_FCONST(0.5); // derived by try and error
if (gabs_t <= gtot_t)
{
w[k] = ws / gtot_t;
@@ -164,7 +160,7 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_t<m_N, _storage_N>::vsolve_non_dyn
do {
resched = false;
-
+ double err = 0;
for (int k = 0; k < iN; k++)
{
const int * RESTRICT net_other = this->m_terms[k]->net_other();
@@ -177,10 +173,13 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_t<m_N, _storage_N>::vsolve_non_dyn
const nl_double new_val = new_V[k] * one_m_w[k] + (Idrive + RHS[k]) * w[k];
- resched = resched || (nl_math::abs(new_val - new_V[k]) > accuracy);
+ err = std::max(nl_math::abs(new_val - new_V[k]), err);
new_V[k] = new_val;
}
+ if (err > accuracy)
+ resched = true;
+
resched_cnt++;
} while (resched && (resched_cnt < this->m_params.m_gs_loops));
diff --git a/src/emu/netlist/analog/nld_ms_sor_mat.h b/src/emu/netlist/analog/nld_ms_sor_mat.h
index 4a788f8f403..01a0d16a229 100644
--- a/src/emu/netlist/analog/nld_ms_sor_mat.h
+++ b/src/emu/netlist/analog/nld_ms_sor_mat.h
@@ -30,11 +30,6 @@ public:
, m_gs_fail(0)
, m_gs_total(0)
{
- pstring p = nl_util::environment("NETLIST_STATS");
- if (p != "")
- m_log_stats = (bool) p.as_long();
- else
- m_log_stats = false;
}
virtual ~netlist_matrix_solver_SOR_mat_t() {}
@@ -53,8 +48,6 @@ private:
nl_double m_lp_fact;
int m_gs_fail;
int m_gs_total;
- bool m_log_stats;
-
};
// ----------------------------------------------------------------------------------------
@@ -64,15 +57,15 @@ private:
template <unsigned m_N, unsigned _storage_N>
void netlist_matrix_solver_SOR_mat_t<m_N, _storage_N>::log_stats()
{
- if (this->m_stat_calculations != 0 && m_log_stats)
+ if (this->m_stat_calculations != 0 && this->m_params.m_log_stats)
{
- this->netlist().log("==============================================\n");
- this->netlist().log("Solver %s\n", this->name().cstr());
- this->netlist().log(" ==> %d nets\n", this->N()); //, (*(*groups[i].first())->m_core_terms.first())->name().cstr());
- this->netlist().log(" has %s elements\n", this->is_dynamic() ? "dynamic" : "no dynamic");
- this->netlist().log(" has %s elements\n", this->is_timestep() ? "timestep" : "no timestep");
- this->netlist().log(" %6.3f average newton raphson loops\n", (double) this->m_stat_newton_raphson / (double) this->m_stat_vsolver_calls);
- this->netlist().log(" %10d invocations (%6d Hz) %10d gs fails (%6.2f%%) %6.3f average\n",
+ this->netlist().log("==============================================");
+ this->netlist().log("Solver %s", this->name().cstr());
+ this->netlist().log(" ==> %d nets", this->N()); //, (*(*groups[i].first())->m_core_terms.first())->name().cstr());
+ this->netlist().log(" has %s elements", this->is_dynamic() ? "dynamic" : "no dynamic");
+ this->netlist().log(" has %s elements", this->is_timestep() ? "timestep" : "no timestep");
+ this->netlist().log(" %6.3f average newton raphson loops", (double) this->m_stat_newton_raphson / (double) this->m_stat_vsolver_calls);
+ this->netlist().log(" %10d invocations (%6d Hz) %10d gs fails (%6.2f%%) %6.3f average",
this->m_stat_calculations,
this->m_stat_calculations * 10 / (int) (this->netlist().time().as_double() * 10.0),
this->m_gs_fail,
@@ -149,20 +142,19 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_mat_t<m_N, _storage_N>::vsolve_non
ATTR_ALIGN nl_double new_v[_storage_N] = { 0.0 };
- const int iN = this->N();
+ const unsigned iN = this->N();
bool resched = false;
int resched_cnt = 0;
this->build_LE_A();
- this->build_LE_RHS();
+ this->build_LE_RHS(this->m_RHS);
-#if 1
#if 0
static int ws_cnt = 0;
ws_cnt++;
- if (0 && ws_cnt % 100 == 0)
+ if (1 && ws_cnt % 200 == 0)
{
// update omega
nl_double lambdaN = 0;
@@ -181,8 +173,8 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_mat_t<m_N, _storage_N>::vsolve_non
for (int i=0; i<iN; i++)
s = s + nl_math::abs(this->m_A[k][i]);
akk = s / akk - 1.0;
- //if ( akk > lambdaN)
- // lambdaN = akk;
+ if ( akk > lambdaN)
+ lambdaN = akk;
if (akk < lambda1)
lambda1 = akk;
#endif
@@ -191,99 +183,31 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_mat_t<m_N, _storage_N>::vsolve_non
//ws = 2.0 / (2.0 - lambdaN - lambda1);
m_omega = 2.0 / (2.0 - lambda1);
- //printf("%f\n", ws);
+ //printf("%f %f %f\n", m_omega, lambda1, lambdaN);
}
#endif
for (int k = 0; k < iN; k++)
new_v[k] = this->m_nets[k]->m_cur_Analog;
-#else
- {
- nl_double frob;
- frob = 0;
- nl_double rmin = 1e99, rmax = -1e99;
- for (int k = 0; k < iN; k++)
- {
- new_v[k] = this->m_nets[k]->m_cur_Analog;
- nl_double s=0.0;
- for (int i = 0; i < iN; i++)
- {
- frob += this->m_A[k][i] * this->m_A[k][i];
- s = s + nl_math::abs(this->m_A[k][i]);
- }
-
- if (s<rmin)
- rmin = s;
- if (s>rmax)
- rmax = s;
- }
-#if 0
- nl_double frobA = nl_math::sqrt(frob /(iN));
- if (1 &&frobA < 1.0)
- //ws = 2.0 / (1.0 + nl_math::sqrt(1.0-frobA));
- ws = 2.0 / (2.0 - frobA);
- else
- ws = 1.0;
- ws = 0.9;
-#else
- // calculate an estimate for rho.
- // This is based on the Perron???Frobenius theorem for positive matrices.
- // No mathematical proof here. The following estimates the
- // optimal relaxation parameter pretty well. Unfortunately, the
- // overhead is bigger than the gain. Consequently the fast GS below
- // uses a fixed GS. One can however use this here to determine a
- // suitable parameter.
- nl_double rm = (rmax + rmin) * 0.5;
- if (rm < 1.0)
- ws = 2.0 / (1.0 + nl_math::sqrt(1.0-rm));
- else
- ws = 1.0;
- if (ws > 1.02 && rmax > 1.001)
- printf("rmin %f rmax %f ws %f\n", rmin, rmax, ws);
-#endif
- }
-#endif
- // Frobenius norm for (D-L)^(-1)U
- //nl_double frobU;
- //nl_double frobL;
- //nl_double norm;
do {
resched = false;
nl_double cerr = 0.0;
- //frobU = 0;
- //frobL = 0;
- //norm = 0;
for (int k = 0; k < iN; k++)
{
nl_double Idrive = 0;
- //nl_double norm_t = 0;
- // Reduction loops need -ffast-math
- for (int i = 0; i < iN; i++)
- Idrive += this->m_A[k][i] * new_v[i];
-#if 0
- for (int i = 0; i < iN; i++)
- {
- if (i < k) frobL += this->m_A[k][i] * this->m_A[k][i] / this->m_A[k][k] /this-> m_A[k][k];
- if (i > k) frobU += this->m_A[k][i] * this->m_A[k][i] / this->m_A[k][k] / this->m_A[k][k];
- norm_t += nl_math::abs(this->m_A[k][i]);
- }
-#endif
- //if (norm_t > norm) norm = norm_t;
-#if 0
- const nl_double new_val = (1.0-ws) * new_v[k] + ws * (this->m_RHS[k] - Idrive + this->m_A[k][k] * new_v[k]) / this->m_A[k][k];
-
- const nl_double e = nl_math::abs(new_val - new_v[k]);
- cerr = (e > cerr ? e : cerr);
- new_v[k] = new_val;
-#else
- const nl_double delta = m_omega * (this->m_RHS[k] - Idrive) / this->m_A[k][k];
- const nl_double adelta = nl_math::abs(delta);
- cerr = (adelta > cerr ? adelta : cerr);
+ const double * RESTRICT A = &this->m_A[k][0];
+ const int *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++)
+ Idrive = Idrive + A[p[i]] * new_v[p[i]];
+
+ const nl_double delta = m_omega * (this->m_RHS[k] - Idrive) / A[k];
+ cerr = std::max(cerr, nl_math::abs(delta));
new_v[k] += delta;
-#endif
}
if (cerr > this->m_params.m_accuracy)
@@ -291,14 +215,9 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_mat_t<m_N, _storage_N>::vsolve_non
resched = true;
}
resched_cnt++;
- //ATTR_UNUSED nl_double frobUL = nl_math::sqrt((frobU + frobL) / (double) (iN) / (double) (iN));
} while (resched && (resched_cnt < this->m_params.m_gs_loops));
- //printf("Frobenius %f %f %f %f %f\n", nl_math::sqrt(frobU), nl_math::sqrt(frobL), frobUL, frobA, norm);
- //printf("Omega Estimate1 %f %f\n", 2.0 / (1.0 + nl_math::sqrt(1-frobUL)), 2.0 / (1.0 + nl_math::sqrt(1-frobA)) ); // printf("Frobenius %f\n", sqrt(frob / (double) (iN * iN) ));
- //printf("Omega Estimate2 %f %f\n", 2.0 / (2.0 - frobUL), 2.0 / (2.0 - frobA) ); // printf("Frobenius %f\n", sqrt(frob / (double) (iN * iN) ));
-
- this->store(new_v, false);
+ this->store(new_v);
this->m_gs_total += resched_cnt;
if (resched)
{
diff --git a/src/emu/netlist/analog/nld_solver.h b/src/emu/netlist/analog/nld_solver.h
index a774912ac00..c46249adc47 100644
--- a/src/emu/netlist/analog/nld_solver.h
+++ b/src/emu/netlist/analog/nld_solver.h
@@ -81,6 +81,8 @@ class terms_t
unsigned m_railstart;
+ plist_t<int> m_nzrd; /* non zero right of the diagonal for elimination */
+ plist_t<int> m_nz; /* all non zero for multiplication */
private:
plist_t<netlist_terminal_t *> m_term;
plist_t<int> m_net_other;
diff --git a/src/emu/netlist/plib/plists.h b/src/emu/netlist/plib/plists.h
index 4c223ab9150..fbc839b3757 100644
--- a/src/emu/netlist/plib/plists.h
+++ b/src/emu/netlist/plib/plists.h
@@ -11,6 +11,7 @@
#define PLISTS_H_
#include <cstring>
+#include <algorithm>
#include "palloc.h"
#include "pstring.h"
@@ -582,4 +583,21 @@ public:
}
};
+// ----------------------------------------------------------------------------------------
+// sort a list ... slow, I am lazy
+// elements must support ">" operator.
+// ----------------------------------------------------------------------------------------
+
+template<typename Class>
+static inline void psort_list(Class &sl)
+{
+ for(int i = 0; i < (int) sl.size() - 1; i++)
+ {
+ for(int j = i + 1; j < sl.size(); j++)
+ if(sl[i] > sl[j])
+ std::swap(sl[i], sl[j]);
+
+ }
+}
+
#endif /* PLISTS_H_ */