From 6fd40f98a61df329ffbe72b1f5c4c0135c9e786d Mon Sep 17 00:00:00 2001 From: couriersud Date: Wed, 17 May 2017 14:20:25 +0200 Subject: Various code alignments across solvers. (nw) --- src/lib/netlist/solver/nld_matrix_solver.h | 7 ++-- src/lib/netlist/solver/nld_ms_direct.h | 40 ++++++++++++++++--- src/lib/netlist/solver/nld_ms_gcr.h | 8 ++-- src/lib/netlist/solver/nld_ms_gmres.h | 25 ++++++------ src/lib/netlist/solver/nld_ms_sm.h | 52 ++++++++++++------------ src/lib/netlist/solver/nld_solver.cpp | 18 ++++++--- src/lib/netlist/solver/vector_base.h | 63 +++++++++++++++++++----------- src/mame/audio/nl_kidniki.cpp | 5 ++- 8 files changed, 140 insertions(+), 78 deletions(-) diff --git a/src/lib/netlist/solver/nld_matrix_solver.h b/src/lib/netlist/solver/nld_matrix_solver.h index f91f8791e9a..da44370d0d9 100644 --- a/src/lib/netlist/solver/nld_matrix_solver.h +++ b/src/lib/netlist/solver/nld_matrix_solver.h @@ -232,9 +232,10 @@ void matrix_solver_t::build_LE_A() for (std::size_t k = 0; k < iN; k++) { terms_for_net_t *terms = m_terms[k].get(); + nl_double * Ak = &child.A(k, 0); for (std::size_t i=0; i < iN; i++) - child.A(k,i) = 0.0; + Ak[i] = 0.0; const std::size_t terms_count = terms->count(); const std::size_t railstart = terms->m_railstart; @@ -245,14 +246,14 @@ void matrix_solver_t::build_LE_A() for (std::size_t i = 0; i < terms_count; i++) akk += gt[i]; - child.A(k,k) = akk; + Ak[k] = akk; } const nl_double * const RESTRICT go = terms->go(); int * RESTRICT net_other = terms->connected_net_idx(); for (std::size_t i = 0; i < railstart; i++) - child.A(k,net_other[i]) -= go[i]; + Ak[net_other[i]] -= go[i]; } } diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h index 019b0890aca..f03cb738028 100644 --- a/src/lib/netlist/solver/nld_ms_direct.h +++ b/src/lib/netlist/solver/nld_ms_direct.h @@ -14,6 +14,12 @@ #include "nld_matrix_solver.h" #include "vector_base.h" +/* Disabling dynamic allocation gives a ~10% boost in performance + * This flag has been added to support continuous storage for arrays + * going forward in case we implement cuda solvers in the future. + */ +#define NL_USE_DYNAMIC_ALLOCATION (1) + namespace netlist { namespace devices @@ -48,10 +54,17 @@ protected: template void LE_back_subst(T * RESTRICT x); +#if (NL_USE_DYNAMIC_ALLOCATION) + template + nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r * m_pitch + c]; } + template + nl_ext_double &RHS(const T1 &r) { return m_A[r * m_pitch + N()]; } +#else template nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r][c]; } template nl_ext_double &RHS(const T1 &r) { return m_A[r][N()]; } +#endif nl_double m_last_RHS[storage_N]; // right hand side - contains currents private: @@ -59,8 +72,12 @@ private: static constexpr std::size_t m_pitch = (((storage_N + 1) + 7) / 8) * 8; //static const std::size_t m_pitch = (((storage_N + 1) + 15) / 16) * 16; //static const std::size_t m_pitch = (((storage_N + 1) + 31) / 32) * 32; +#if (NL_USE_DYNAMIC_ALLOCATION) + //nl_ext_double * RESTRICT m_A; + std::vector m_A; +#else nl_ext_double m_A[storage_N][m_pitch]; - +#endif //nl_ext_double m_RHSx[storage_N]; const std::size_t m_dim; @@ -74,6 +91,9 @@ private: template matrix_solver_direct_t::~matrix_solver_direct_t() { +#if (NL_USE_DYNAMIC_ALLOCATION) + //plib::pfree_array(m_A); +#endif } template @@ -107,15 +127,17 @@ void matrix_solver_direct_t::LE_solve() { /* FIXME: Singular matrix? */ + nl_double *Ai = &A(i, 0); 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; for (std::size_t j : nzbd) { - const nl_double f1 = -f * A(j,i); + nl_double *Aj = &A(j, 0); + const nl_double f1 = -f * Aj[i]; for (std::size_t k : nzrd) - A(j,k) += A(i,k) * f1; + Aj[k] += Ai[k] * f1; //RHS(j) += RHS(i) * f1; } } @@ -195,11 +217,11 @@ void matrix_solver_direct_t::LE_back_subst( const auto *p = m_terms[j]->m_nzrd.data(); const auto e = m_terms[j]->m_nzrd.size() - 1; /* exclude RHS element */ - + T * Aj = &A(j,0); for (std::size_t k = 0; k < e; k++) { const auto pk = p[k]; - tmp += A(j,pk) * x[pk]; + tmp += Aj[pk] * x[pk]; } x[j] = (RHS(j) - tmp) / A(j,j); } @@ -239,6 +261,10 @@ matrix_solver_direct_t::matrix_solver_direct_t(netlist_t &anetli : matrix_solver_t(anetlist, name, ASCENDING, params) , m_dim(size) { +#if (NL_USE_DYNAMIC_ALLOCATION) + m_A.resize(N() * m_pitch); + //m_A = plib::palloc_array(N() * m_pitch); +#endif for (unsigned k = 0; k < N(); k++) { m_last_RHS[k] = 0.0; @@ -251,6 +277,10 @@ matrix_solver_direct_t::matrix_solver_direct_t(netlist_t &anetli : matrix_solver_t(anetlist, name, sort, params) , m_dim(size) { +#if (NL_USE_DYNAMIC_ALLOCATION) + m_A.resize(N() * m_pitch); + //m_A = plib::palloc_array(N() * m_pitch); +#endif for (unsigned k = 0; k < N(); k++) { m_last_RHS[k] = 0.0; diff --git a/src/lib/netlist/solver/nld_ms_gcr.h b/src/lib/netlist/solver/nld_ms_gcr.h index 76e3ac7ecba..a17ca1ec282 100755 --- a/src/lib/netlist/solver/nld_ms_gcr.h +++ b/src/lib/netlist/solver/nld_ms_gcr.h @@ -364,10 +364,12 @@ unsigned matrix_solver_GCR_t::vsolve_non_dynamic(const bool newt } #else for (std::size_t i = 0; i < railstart; i++) - { mat.A[tcr[i]] -= go[i]; - gtot_t = gtot_t + gt[i]; - RHS_t = RHS_t + Idr[i]; + + for (std::size_t i = 0; i < railstart; i++) + { + gtot_t += gt[i]; + RHS_t += Idr[i]; } for (std::size_t i = railstart; i < term_count; i++) diff --git a/src/lib/netlist/solver/nld_ms_gmres.h b/src/lib/netlist/solver/nld_ms_gmres.h index c49e3a97010..2e4e447d14f 100644 --- a/src/lib/netlist/solver/nld_ms_gmres.h +++ b/src/lib/netlist/solver/nld_ms_gmres.h @@ -49,7 +49,7 @@ private: //typedef typename mat_cr_t::type mattype; typedef typename mat_cr_t::index_type mattype; - unsigned solve_ilu_gmres(nl_double (& RESTRICT x)[storage_N], const nl_double * RESTRICT rhs, const unsigned restart_max, const std::size_t mr, nl_double accuracy); + unsigned solve_ilu_gmres(nl_double (& RESTRICT x)[storage_N], const nl_double (& RESTRICT rhs)[storage_N], const unsigned restart_max, std::size_t mr, nl_double accuracy); std::vector m_term_cr[storage_N]; @@ -144,8 +144,6 @@ unsigned matrix_solver_GMRES_t::vsolve_non_dynamic(const bool ne const nl_double * const RESTRICT Idr = this->m_terms[k]->Idr(); const nl_double * const * RESTRICT other_cur_analog = this->m_terms[k]->connected_net_V(); - new_V[k] = this->m_nets[k]->Q_Analog(); - for (std::size_t i = 0; i < term_count; i++) { gtot_t = gtot_t + gt[i]; @@ -165,17 +163,20 @@ unsigned matrix_solver_GMRES_t::vsolve_non_dynamic(const bool ne const std::size_t pi = m_term_cr[k][i]; mat.A[pi] -= go[i]; } + + new_V[k] = this->m_nets[k]->Q_Analog(); + } mat.ia[iN] = static_cast(mat.nz_num); const nl_double accuracy = this->m_params.m_accuracy; - std::size_t mr = iN; + unsigned mr = iN; if (iN > 3 ) mr = static_cast(std::sqrt(iN) * 2.0); unsigned iter = std::max(1u, this->m_params.m_gs_loops); - unsigned gsl = solve_ilu_gmres(new_V, RHS, iter, static_cast(mr), accuracy); - unsigned failed = static_cast(mr) * iter; + unsigned gsl = solve_ilu_gmres(new_V, RHS, iter, mr, accuracy); + unsigned failed = mr * iter; this->m_iterative_total += gsl; this->m_stat_calculations++; @@ -202,7 +203,7 @@ inline static void givens_mult( const T c, const T s, T & g0, T & g1 ) } template -unsigned matrix_solver_GMRES_t::solve_ilu_gmres (nl_double (& RESTRICT x)[storage_N], const nl_double * RESTRICT rhs, const unsigned restart_max, const std::size_t mr, nl_double accuracy) +unsigned matrix_solver_GMRES_t::solve_ilu_gmres (nl_double (& RESTRICT x)[storage_N], const nl_double (& RESTRICT rhs)[storage_N], const unsigned restart_max, std::size_t mr, nl_double accuracy) { /*------------------------------------------------------------------------- * The code below was inspired by code published by John Burkardt under @@ -231,6 +232,8 @@ unsigned matrix_solver_GMRES_t::solve_ilu_gmres (nl_double (& RE const std::size_t n = this->N(); + if (mr > n) mr = n; + if (m_use_iLU_preconditioning) mat.incomplete_LU_factorization(m_LU); @@ -254,7 +257,7 @@ unsigned matrix_solver_GMRES_t::solve_ilu_gmres (nl_double (& RE mat.solveLUx(m_LU, Ax); - const nl_double rho_to_accuracy = std::sqrt(vecmult2(n, Ax)) / accuracy; + const nl_double rho_to_accuracy = std::sqrt(vec_mult2(n, Ax)) / accuracy; rho_delta = accuracy * rho_to_accuracy; } @@ -278,7 +281,7 @@ unsigned matrix_solver_GMRES_t::solve_ilu_gmres (nl_double (& RE mat.solveLUx(m_LU, residual); } - rho = std::sqrt(vecmult2(n, residual)); + rho = std::sqrt(vec_mult2(n, residual)); if (rho < rho_delta) return itr_used + 1; @@ -302,10 +305,10 @@ unsigned matrix_solver_GMRES_t::solve_ilu_gmres (nl_double (& RE for (std::size_t j = 0; j <= k; j++) { - m_ht[j][k] = vecmult(n, m_v[k1], m_v[j]); + m_ht[j][k] = vec_mult(n, m_v[k1], m_v[j]); vec_add_mult_scalar(n, m_v[j], -m_ht[j][k], m_v[k1]); } - m_ht[k1][k] = std::sqrt(vecmult2(n, m_v[k1])); + m_ht[k1][k] = std::sqrt(vec_mult2(n, m_v[k1])); if (m_ht[k1][k] != 0.0) vec_scale(n, m_v[k1], NL_FCONST(1.0) / m_ht[k1][k]); diff --git a/src/lib/netlist/solver/nld_ms_sm.h b/src/lib/netlist/solver/nld_ms_sm.h index acafe285da3..ed95cd2dc03 100644 --- a/src/lib/netlist/solver/nld_ms_sm.h +++ b/src/lib/netlist/solver/nld_ms_sm.h @@ -75,19 +75,19 @@ protected: template - inline nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r][c]; } + nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r][c]; } template - inline nl_ext_double &W(const T1 &r, const T2 &c) { return m_W[r][c]; } + nl_ext_double &W(const T1 &r, const T2 &c) { return m_W[r][c]; } template - inline nl_ext_double &Ainv(const T1 &r, const T2 &c) { return m_Ainv[r][c]; } + nl_ext_double &Ainv(const T1 &r, const T2 &c) { return m_Ainv[r][c]; } template - inline nl_ext_double &RHS(const T1 &r) { return m_RHS[r]; } + nl_ext_double &RHS(const T1 &r) { return m_RHS[r]; } template - inline nl_ext_double &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; } + nl_ext_double &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; } template - inline nl_ext_double &lAinv(const T1 &r, const T2 &c) { return m_lAinv[r][c]; } + nl_ext_double &lAinv(const T1 &r, const T2 &c) { return m_lAinv[r][c]; } nl_double m_last_RHS[storage_N]; // right hand side - contains currents @@ -104,6 +104,7 @@ private: //nl_ext_double m_RHSx[storage_N]; const std::size_t m_dim; + std::size_t m_cnt; }; @@ -215,13 +216,12 @@ void matrix_solver_sm_t::LE_compute_x( template unsigned matrix_solver_sm_t::solve_non_dynamic(const bool newton_raphson) { - static constexpr bool incremental = true; - static unsigned cnt = 0; - const auto iN = N(); + static constexpr const bool incremental = true; + const std::size_t iN = N(); nl_double new_V[storage_N]; // = { 0.0 }; - if (0 || ((cnt % 200) == 0)) + if ((m_cnt % 50) == 0) { /* complete calculation */ this->LE_invert(); @@ -230,18 +230,18 @@ unsigned matrix_solver_sm_t::solve_non_dynamic(const bool newton { if (!incremental) { - for (unsigned row = 0; row < iN; row ++) - for (unsigned k = 0; k < iN; k++) + for (std::size_t row = 0; row < iN; row ++) + for (std::size_t k = 0; k < iN; k++) Ainv(row,k) = lAinv(row, k); } - for (unsigned row = 0; row < iN; row ++) + for (std::size_t row = 0; row < iN; row ++) { nl_double v[m_pitch] = {0}; - unsigned cols[m_pitch]; - unsigned colcount = 0; + std::size_t cols[m_pitch]; + std::size_t colcount = 0; auto &nz = m_terms[row]->m_nz; - for (auto & col : nz) + for (unsigned & col : nz) { v[col] = A(row,col) - lA(row,col); if (incremental) @@ -254,28 +254,29 @@ unsigned matrix_solver_sm_t::solve_non_dynamic(const bool newton { nl_double lamba = 0.0; nl_double w[m_pitch] = {0}; + nl_double z[m_pitch]; /* compute w and lamba */ - for (unsigned i = 0; i < iN; i++) + for (std::size_t i = 0; i < iN; i++) z[i] = Ainv(i, row); /* u is row'th column */ - for (unsigned j = 0; j < colcount; j++) + for (std::size_t j = 0; j < colcount; j++) lamba += v[cols[j]] * z[cols[j]]; - for (unsigned j=0; j::solve_non_dynamic(const bool newton } } - cnt++; + m_cnt++; this->LE_compute_x(new_V); @@ -310,6 +311,7 @@ matrix_solver_sm_t::matrix_solver_sm_t(netlist_t &anetlist, cons const solver_parameters_t *params, const std::size_t size) : matrix_solver_t(anetlist, name, NOSORT, params) , m_dim(size) +, m_cnt(0) { for (std::size_t k = 0; k < N(); k++) { diff --git a/src/lib/netlist/solver/nld_solver.cpp b/src/lib/netlist/solver/nld_solver.cpp index 2e2df3c0f05..cae76caab1a 100644 --- a/src/lib/netlist/solver/nld_solver.cpp +++ b/src/lib/netlist/solver/nld_solver.cpp @@ -10,15 +10,16 @@ */ #if 0 -#pragma GCC optimize "fast-math" -#pragma GCC optimize "strict-aliasing" -#pragma GCC optimize "tree-vectorize" +#pragma GCC optimize "-ftree-vectorize" +#pragma GCC optimize "-ffast-math" +#pragma GCC optimize "-funsafe-math-optimizations" +#pragma GCC optimize "-funroll-loops" +#pragma GCC optimize "-funswitch-loops" +#pragma GCC optimize "-fstrict-aliasing" #pragma GCC optimize "tree-vectorizer-verbose=7" #pragma GCC optimize "opt-info-vec" #pragma GCC optimize "opt-info-vec-missed" //#pragma GCC optimize "tree-parallelize-loops=4" -#pragma GCC optimize "unroll-loops" -#pragma GCC optimize "unswitch-loops" #pragma GCC optimize "variable-expansion-in-unroller" #pragma GCC optimize "unsafe-loop-optimizations" #pragma GCC optimize "vect-cost-model" @@ -301,6 +302,7 @@ void NETLIB_NAME(solver)::post_start() else ms = create_solver<2,2>(2, sname); break; +#if 0 case 3: ms = create_solver<3,3>(3, sname); break; @@ -338,11 +340,15 @@ void NETLIB_NAME(solver)::post_start() ms = create_solver<31,31>(31, sname); break; case 35: - ms = create_solver<31,31>(31, sname); + ms = create_solver<35,35>(35, sname); + break; + case 43: + ms = create_solver<43,43>(43, sname); break; case 49: ms = create_solver<49,49>(49, sname); break; +#endif #if 0 case 87: ms = create_solver<87,87>(87, sname); diff --git a/src/lib/netlist/solver/vector_base.h b/src/lib/netlist/solver/vector_base.h index b241eb77481..28f5fa0dd44 100755 --- a/src/lib/netlist/solver/vector_base.h +++ b/src/lib/netlist/solver/vector_base.h @@ -35,31 +35,40 @@ private: #pragma GCC diagnostic ignored "-Wmaybe-uninitialized" #endif -template -inline static void vec_set (const std::size_t n, const T scalar, T * RESTRICT result) +template +inline static void vec_set (const std::size_t n, const T scalar, T (& RESTRICT v)[N]) { - for ( std::size_t i = 0; i < n; i++ ) - result[i] = scalar; + if (n != N) + for ( std::size_t i = 0; i < n; i++ ) + v[i] = scalar; + else + for ( std::size_t i = 0; i < N; i++ ) + v[i] = scalar; } -template -inline static T vecmult (const std::size_t n, const T * RESTRICT a1, const T * RESTRICT a2 ) +template +inline static T vec_mult (const std::size_t n, const T (& RESTRICT v1)[N], const T (& RESTRICT v2)[N] ) { T value = 0.0; - for ( std::size_t i = 0; i < n; i++ ) - value = value + a1[i] * a2[i]; + if (n != N) + for ( std::size_t i = 0; i < n; i++ ) + value += v1[i] * v2[i]; + else + for ( std::size_t i = 0; i < N; i++ ) + value += v1[i] * v2[i]; return value; } -template -inline static T vecmult2 (const std::size_t n, const T *a1) +template +inline static T vec_mult2 (const std::size_t n, const T (& RESTRICT v)[N]) { T value = 0.0; - for ( std::size_t i = 0; i < n; i++ ) - { - const T temp = a1[i]; - value = value + temp * temp; - } + if (n != N) + for ( std::size_t i = 0; i < n; i++ ) + value += v[i] * v[i]; + else + for ( std::size_t i = 0; i < N; i++ ) + value += v[i] * v[i]; return value; } @@ -98,18 +107,26 @@ inline static void vec_add_ip(const std::size_t n, const double * RESTRICT v, do result[i] += v[i]; } -template -inline void vec_sub(const std::size_t n, const T * RESTRICT v1, const T * RESTRICT v2, T * RESTRICT result) +template +inline void vec_sub(const std::size_t n, const T (& RESTRICT v1)[N], const T (& RESTRICT v2)[N], T (& RESTRICT result)[N]) { - for ( std::size_t i = 0; i < n; i++ ) - result[i] = v1[i] - v2[i]; + if (n != N) + for ( std::size_t i = 0; i < n; i++ ) + result[i] = v1[i] - v2[i]; + else + for ( std::size_t i = 0; i < N; i++ ) + result[i] = v1[i] - v2[i]; } -template -inline void vec_scale (const std::size_t n, T * RESTRICT v, const T scalar) +template +inline void vec_scale(const std::size_t n, T (& RESTRICT v)[N], const T scalar) { - for ( std::size_t i = 0; i < n; i++ ) - v[i] = scalar * v[i]; + if (n != N) + for ( std::size_t i = 0; i < n; i++ ) + v[i] = scalar * v[i]; + else + for ( std::size_t i = 0; i < N; i++ ) + v[i] = scalar * v[i]; } inline double vec_maxabs(const std::size_t n, const double * RESTRICT v) diff --git a/src/mame/audio/nl_kidniki.cpp b/src/mame/audio/nl_kidniki.cpp index ba158d87eec..af9fc061c8c 100644 --- a/src/mame/audio/nl_kidniki.cpp +++ b/src/mame/audio/nl_kidniki.cpp @@ -310,11 +310,12 @@ NETLIST_START(kidniki) #if (1 || USE_FRONTIERS) SOLVER(Solver, 18000) - PARAM(Solver.ACCURACY, 1e-8) + PARAM(Solver.ACCURACY, 1e-7) PARAM(Solver.NR_LOOPS, 100) - PARAM(Solver.GS_LOOPS, 5) + PARAM(Solver.GS_LOOPS, 10) //PARAM(Solver.METHOD, "SOR") PARAM(Solver.METHOD, "MAT_CR") + //PARAM(Solver.METHOD, "MAT") //PARAM(Solver.METHOD, "GMRES") PARAM(Solver.SOR_FACTOR, 1.00) PARAM(Solver.DYNAMIC_TS, 0) -- cgit v1.2.3