diff options
| author | 2017-05-17 14:20:25 +0200 | |
|---|---|---|
| committer | 2017-05-27 00:11:16 +0200 | |
| commit | 6fd40f98a61df329ffbe72b1f5c4c0135c9e786d (patch) | |
| tree | 8c118e321c0ab3f77450f1b81d6cfa8a109d158c /src/lib | |
| parent | 6dfe04c620590bdfcddaf1c36d07832af011f2ec (diff) | |
Various code alignments across solvers. (nw)
Diffstat (limited to 'src/lib')
| -rw-r--r-- | src/lib/netlist/solver/nld_matrix_solver.h | 7 | ||||
| -rw-r--r-- | src/lib/netlist/solver/nld_ms_direct.h | 40 | ||||
| -rwxr-xr-x | src/lib/netlist/solver/nld_ms_gcr.h | 8 | ||||
| -rw-r--r-- | src/lib/netlist/solver/nld_ms_gmres.h | 25 | ||||
| -rw-r--r-- | src/lib/netlist/solver/nld_ms_sm.h | 52 | ||||
| -rw-r--r-- | src/lib/netlist/solver/nld_solver.cpp | 18 | ||||
| -rwxr-xr-x | src/lib/netlist/solver/vector_base.h | 63 |
7 files changed, 137 insertions, 76 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 <typename T> void LE_back_subst(T * RESTRICT x); +#if (NL_USE_DYNAMIC_ALLOCATION) + template <typename T1, typename T2> + nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r * m_pitch + c]; } + template <typename T1> + nl_ext_double &RHS(const T1 &r) { return m_A[r * m_pitch + N()]; } +#else template <typename T1, typename T2> nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r][c]; } template <typename T1> 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<nl_ext_double> 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 <std::size_t m_N, std::size_t storage_N> matrix_solver_direct_t<m_N, storage_N>::~matrix_solver_direct_t() { +#if (NL_USE_DYNAMIC_ALLOCATION) + //plib::pfree_array(m_A); +#endif } template <std::size_t m_N, std::size_t storage_N> @@ -107,15 +127,17 @@ void matrix_solver_direct_t<m_N, storage_N>::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<m_N, storage_N>::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<m_N, storage_N>::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<nl_ext_double>(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<m_N, storage_N>::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<nl_ext_double>(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<m_N, storage_N>::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<storage_N>::type mattype; typedef typename mat_cr_t<storage_N>::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<unsigned> m_term_cr[storage_N]; @@ -144,8 +144,6 @@ unsigned matrix_solver_GMRES_t<m_N, storage_N>::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<m_N, storage_N>::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<mattype>(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<unsigned>(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<unsigned>(mr), accuracy); - unsigned failed = static_cast<unsigned>(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 <std::size_t m_N, std::size_t storage_N> -unsigned matrix_solver_GMRES_t<m_N, storage_N>::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<m_N, storage_N>::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<m_N, storage_N>::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<m_N, storage_N>::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<m_N, storage_N>::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<m_N, storage_N>::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 <typename T1, typename T2> - 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 <typename T1, typename T2> - 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 <typename T1, typename T2> - 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 <typename T1> - inline nl_ext_double &RHS(const T1 &r) { return m_RHS[r]; } + nl_ext_double &RHS(const T1 &r) { return m_RHS[r]; } template <typename T1, typename T2> - 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 <typename T1, typename T2> - 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<m_N, storage_N>::LE_compute_x( template <std::size_t m_N, std::size_t storage_N> unsigned matrix_solver_sm_t<m_N, storage_N>::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<m_N, storage_N>::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<m_N, storage_N>::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<colcount; j++) + for (std::size_t j=0; j<colcount; j++) { - auto col = cols[j]; - auto f = v[col]; - for (unsigned k = 0; k < iN; k++) + std::size_t col = cols[j]; + nl_double f = v[col]; + for (std::size_t k = 0; k < iN; k++) w[k] += Ainv(col,k) * f; /* Transpose(Ainv) * v */ } lamba = -1.0 / (1.0 + lamba); - for (unsigned i=0; i<iN; i++) + for (std::size_t i=0; i<iN; i++) { const nl_double f = lamba * z[i]; if (f != 0.0) - for (unsigned k = 0; k < iN; k++) + for (std::size_t k = 0; k < iN; k++) Ainv(i,k) += f * w[k]; } } @@ -283,7 +284,7 @@ unsigned matrix_solver_sm_t<m_N, storage_N>::solve_non_dynamic(const bool newton } } - cnt++; + m_cnt++; this->LE_compute_x(new_V); @@ -310,6 +311,7 @@ matrix_solver_sm_t<m_N, storage_N>::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<typename T> -inline static void vec_set (const std::size_t n, const T scalar, T * RESTRICT result) +template<typename T, std::size_t N> +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<typename T> -inline static T vecmult (const std::size_t n, const T * RESTRICT a1, const T * RESTRICT a2 ) +template<typename T, std::size_t N> +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<typename T> -inline static T vecmult2 (const std::size_t n, const T *a1) +template<typename T, std::size_t N> +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<typename T> -inline void vec_sub(const std::size_t n, const T * RESTRICT v1, const T * RESTRICT v2, T * RESTRICT result) +template<typename T, std::size_t N> +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<typename T> -inline void vec_scale (const std::size_t n, T * RESTRICT v, const T scalar) +template<typename T, std::size_t N> +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) |
