From 63efcfcb85e2f475b26541d0973b0321c9204598 Mon Sep 17 00:00:00 2001 From: couriersud Date: Sun, 20 Sep 2020 08:04:22 +0200 Subject: netlist: code refactoring. - rename mat_cr.h to pmatrix_cr.h - Optimization to the gmres solver. - Simplifcation of vector operation code --- scripts/src/netlist.lua | 2 +- src/lib/netlist/plib/gmres.h | 189 +++++---- src/lib/netlist/plib/mat_cr.h | 656 ----------------------------- src/lib/netlist/plib/pgsl.h | 4 + src/lib/netlist/plib/pmatrix_cr.h | 635 ++++++++++++++++++++++++++++ src/lib/netlist/plib/vector_ops.h | 73 ++-- src/lib/netlist/solver/nld_matrix_solver.h | 4 +- src/lib/netlist/solver/nld_ms_gcr.h | 8 +- src/lib/netlist/solver/nld_ms_gmres.h | 4 +- 9 files changed, 786 insertions(+), 789 deletions(-) delete mode 100644 src/lib/netlist/plib/mat_cr.h create mode 100644 src/lib/netlist/plib/pmatrix_cr.h diff --git a/scripts/src/netlist.lua b/scripts/src/netlist.lua index c5dec60b6b4..6f90d32a06d 100644 --- a/scripts/src/netlist.lua +++ b/scripts/src/netlist.lua @@ -81,6 +81,7 @@ project "netlist" MAME_DIR .. "src/lib/netlist/plib/pmain.cpp", MAME_DIR .. "src/lib/netlist/plib/pmain.h", MAME_DIR .. "src/lib/netlist/plib/pmath.h", + MAME_DIR .. "src/lib/netlist/plib/pmatrix_cr.h", MAME_DIR .. "src/lib/netlist/plib/pmempool.h", MAME_DIR .. "src/lib/netlist/plib/pmulti_threading.h", MAME_DIR .. "src/lib/netlist/plib/pomp.h", @@ -123,7 +124,6 @@ project "netlist" MAME_DIR .. "src/lib/netlist/solver/nld_ms_sor.h", MAME_DIR .. "src/lib/netlist/solver/nld_ms_sor_mat.h", MAME_DIR .. "src/lib/netlist/solver/nld_ms_gmres.h", - MAME_DIR .. "src/lib/netlist/solver/mat_cr.h", MAME_DIR .. "src/lib/netlist/solver/nld_ms_sm.h", MAME_DIR .. "src/lib/netlist/solver/nld_ms_w.h", MAME_DIR .. "src/lib/netlist/solver/nld_ms_direct_lu.h", diff --git a/src/lib/netlist/plib/gmres.h b/src/lib/netlist/plib/gmres.h index df0fc787c52..465f12021e3 100644 --- a/src/lib/netlist/plib/gmres.h +++ b/src/lib/netlist/plib/gmres.h @@ -8,7 +8,7 @@ /// \file gmres.h /// -#include "mat_cr.h" +#include "pmatrix_cr.h" #include "parray.h" #include "pconfig.h" #include "vector_ops.h" @@ -33,11 +33,11 @@ namespace plib template struct mat_precondition_ILU { - using mat_type = plib::pmatrix_cr_t; - using matLU_type = plib::pLUmatrix_cr_t; + using mat_type = plib::pmatrix_cr; + using matLU_type = plib::pLUmatrix_cr; mat_precondition_ILU(std::size_t size, std::size_t ilu_scale = 4 - , std::size_t bw = plib::pmatrix_cr_t::FILL_INFINITY) + , std::size_t bw = plib::pmatrix_cr::FILL_INFINITY) : m_mat(narrow_cast(size)) , m_LU(narrow_cast(size)) , m_ILU_scale(narrow_cast(ilu_scale)) @@ -167,7 +167,7 @@ namespace plib v[i] = v[i] * m_diag[i]; } - plib::pmatrix_cr_t m_mat; + plib::pmatrix_cr m_mat; plib::parray m_diag; plib::parray, SIZE > nzcol; }; @@ -203,19 +203,23 @@ namespace plib plib::unused_var(v); } - plib::pmatrix_cr_t m_mat; + plib::pmatrix_cr m_mat; }; // FIXME: hardcoding RESTART to 20 becomes an issue on very large // systems. - template + template struct gmres_t { public: using float_type = FT; + //constexpr static int RESTART = RESTARTMAX; + constexpr static const int RESTART = (SIZE > 0) ? ((SIZE < RESTARTMAX) ? SIZE : RESTARTMAX) + : ((SIZE < 0) ? ((-SIZE < RESTARTMAX) ? -SIZE : RESTARTMAX) : RESTARTMAX); + explicit gmres_t(std::size_t size) : residual(size) , Ax(size) @@ -226,83 +230,8 @@ namespace plib { } - void givens_mult( const FT c, const FT s, FT & g0, FT & g1 ) - { - const FT g0_last(g0); - - g0 = c * g0 - s * g1; - g1 = s * g0_last + c * g1; - } - std::size_t size() const { return (SIZE<=0) ? m_size : narrow_cast(SIZE); } - template - bool do_k(OPS &ops, VT &x, std::size_t &itr_used, FT rho_delta, bool dummy) - { - plib::unused_var(dummy); - if (do_k(ops, x, itr_used, rho_delta, do_khelper::value)) - return true; - - const std::size_t kp1 = k + 1; - const std::size_t n = size(); - - ops.calc_rhs(m_v[kp1], m_v[k]); - ops.solve_inplace(m_v[kp1]); - - for (std::size_t j = 0; j <= k; j++) - { - m_ht[j][k] = vec_mult(n, m_v[kp1], m_v[j]); - vec_add_mult_scalar(n, m_v[kp1], m_v[j], -m_ht[j][k]); - } - m_ht[kp1][k] = plib::sqrt(vec_mult2(n, m_v[kp1])); - - // FIXME: comparison to zero - if (m_ht[kp1][k] != plib::constants::zero()) - vec_scale(n, m_v[kp1], reciprocal(m_ht[kp1][k])); - - for (std::size_t j = 0; j < k; j++) - givens_mult(m_c[j], m_s[j], m_ht[j][k], m_ht[j+1][k]); - - const float_type mu = reciprocal(plib::hypot(m_ht[k][k], m_ht[kp1][k])); - - m_c[k] = m_ht[k][k] * mu; - m_s[k] = -m_ht[kp1][k] * mu; - m_ht[k][k] = m_c[k] * m_ht[k][k] - m_s[k] * m_ht[kp1][k]; - m_ht[kp1][k] = plib::constants::zero(); - - givens_mult(m_c[k], m_s[k], m_g[k], m_g[kp1]); - - FT rho = plib::abs(m_g[kp1]); - - // FIXME .. - itr_used = itr_used + 1; - - if (rho <= rho_delta || k == RESTART-1) - { - // Solve the system H * y = g - // x += m_v[j] * m_y[j] - for (std::size_t i = k + 1; i-- > 0;) - { - auto tmp(m_g[i]); - for (std::size_t j = i + 1; j <= k; j++) - tmp -= m_ht[i][j] * m_y[j]; - m_y[i] = tmp / m_ht[i][i]; - } - - for (std::size_t i = 0; i <= k; i++) - vec_add_mult_scalar(n, x, m_v[i], m_y[i]); - return true; - } - return false; - } - - template - bool do_k(OPS &ops, VT &x, std::size_t &itr_used, FT rho_delta, float dummy) - { - plib::unused_var(ops, x, itr_used, rho_delta, dummy); - return false; - } - template std::size_t solve(OPS &ops, VT &x, const VRHS & rhs, const std::size_t itr_max, float_type accuracy) { @@ -348,18 +277,22 @@ namespace plib // differently: The invest doesn't pay off. // - vec_set_scalar(n, residual, accuracy); + vec_set_scalar(residual, accuracy); ops.calc_rhs(Ax, residual); ops.solve_inplace(Ax); - const float_type rho_to_accuracy = plib::sqrt(vec_mult2(n, Ax)) / accuracy; + const float_type rho_to_accuracy = plib::sqrt(vec_mult2(Ax)) / accuracy; rho_delta = accuracy * rho_to_accuracy; } else + //rho_delta = accuracy * plib::sqrt(vec_mult2(n, rhs)) + // + 1e-4 * std::sqrt(n); rho_delta = accuracy * plib::sqrt(narrow_cast(n)); + // + // LU x = b; solve for x; // // Using // @@ -368,9 +301,6 @@ namespace plib // // to get a starting point for x degrades convergence speed compared // to using the last solution for x. - // - // LU x = b; solve for x; - // while (itr_used < itr_max) { @@ -378,11 +308,11 @@ namespace plib ops.calc_rhs(Ax, x); - vec_sub(n, residual, rhs, Ax); + vec_sub(residual, rhs, Ax); ops.solve_inplace(residual); - rho = plib::sqrt(vec_mult2(n, residual)); + rho = plib::sqrt(vec_mult2(residual)); if (rho < rho_delta) return itr_used + 1; @@ -391,10 +321,10 @@ namespace plib // on some systems / compiler versions. Issue reported by // AJR, no details known yet. - vec_set_scalar(RESTART+1, m_g, +constants::zero()); + vec_set_scalar(m_g, +constants::zero()); m_g[0] = rho; - vec_mult_scalar(n, m_v[0], residual, reciprocal(rho)); + vec_mult_scalar(m_v[0], residual, plib::reciprocal(rho)); if (do_k(ops, x, itr_used, rho_delta, true)) // converged @@ -405,6 +335,83 @@ namespace plib private: + static inline void givens_mult(FT c, FT s, FT & g0, FT & g1 ) + { + const FT g0_last(g0); + + g0 = c * g0 - s * g1; + g1 = s * g0_last + c * g1; + } + + template + bool do_k(OPS &ops, VT &x, std::size_t &itr_used, FT rho_delta, bool dummy) + { + plib::unused_var(dummy); + if (do_k(ops, x, itr_used, rho_delta, do_khelper::value)) + return true; + + constexpr const std::size_t kp1 = k + 1; + //const std::size_t n = size(); + + ops.calc_rhs(m_v[kp1], m_v[k]); + ops.solve_inplace(m_v[kp1]); + + for (std::size_t j = 0; j <= k; j++) + { + m_ht[j][k] = vec_mult(m_v[kp1], m_v[j]); + vec_add_mult_scalar(m_v[kp1], m_v[j], -m_ht[j][k]); + } + m_ht[kp1][k] = plib::sqrt(vec_mult2(m_v[kp1])); + + // FIXME: comparison to zero + if (m_ht[kp1][k] != plib::constants::zero()) + vec_scale(m_v[kp1], reciprocal(m_ht[kp1][k])); + + for (std::size_t j = 0; j < k; j++) + givens_mult(m_c[j], m_s[j], m_ht[j][k], m_ht[j+1][k]); + + const float_type mu = reciprocal(plib::hypot(m_ht[k][k], m_ht[kp1][k])); + + m_c[k] = m_ht[k][k] * mu; + m_s[k] = -m_ht[kp1][k] * mu; + m_ht[k][k] = m_c[k] * m_ht[k][k] - m_s[k] * m_ht[kp1][k]; + m_ht[kp1][k] = plib::constants::zero(); + + givens_mult(m_c[k], m_s[k], m_g[k], m_g[kp1]); + + const float_type rho = plib::abs(m_g[kp1]); + + // FIXME .. + itr_used = itr_used + 1; + + if (rho <= rho_delta || k == RESTART-1) + { + // Solve the system H * y = g + // x += m_v[j] * m_y[j] + for (std::size_t i = k + 1; i-- > 0;) + { + auto tmp(m_g[i]); + const auto htii=plib::reciprocal(m_ht[i][i]); + for (std::size_t j = i + 1; j <= k; j++) + tmp -= m_ht[i][j] * m_y[j]; + m_y[i] = tmp * htii; + vec_add_mult_scalar(x, m_v[i], m_y[i]); + } + + //for (std::size_t i = 0; i <= k; i++) + // vec_add_mult_scalar(n, x, m_v[i], m_y[i]); + return true; + } + return false; + } + + template + bool do_k(OPS &ops, VT &x, std::size_t &itr_used, FT rho_delta, float dummy) + { + plib::unused_var(ops, x, itr_used, rho_delta, dummy); + return false; + } + plib::parray residual; plib::parray Ax; diff --git a/src/lib/netlist/plib/mat_cr.h b/src/lib/netlist/plib/mat_cr.h deleted file mode 100644 index 17ce140c6e5..00000000000 --- a/src/lib/netlist/plib/mat_cr.h +++ /dev/null @@ -1,656 +0,0 @@ -// license:GPL-2.0+ -// copyright-holders:Couriersud - -#ifndef MAT_CR_H_ -#define MAT_CR_H_ - -/// -/// \file mat_cr.h -/// -/// Compressed row format matrices -/// - -#include "palloc.h" -#include "parray.h" -#include "pconfig.h" -#include "pmath.h" -#include "pmatrix2d.h" -#include "pomp.h" -#include "pstate.h" -#include "ptypes.h" -#include "putil.h" - -#include -#include -#include -#include - -namespace plib -{ - - template - struct pmatrix_cr_t - { - using index_type = C; - using value_type = T; - - static constexpr const int NSQ = (N < 0 ? -N * N : N * N); - static constexpr const int Np1 = (N == 0) ? 0 : (N < 0 ? N - 1 : N + 1); - - PCOPYASSIGNMOVE(pmatrix_cr_t, default) - - enum constants_e - { - FILL_INFINITY = 9999999 - }; - - // FIXME: these should be private - // NOLINTNEXTLINE - parray diag; // diagonal index pointer n - // NOLINTNEXTLINE - parray row_idx; // row index pointer n + 1 - // NOLINTNEXTLINE - parray col_idx; // column index array nz_num, initially (n * n) - // NOLINTNEXTLINE - parray A; // Matrix elements nz_num, initially (n * n) - // NOLINTNEXTLINE - std::size_t nz_num; - - explicit pmatrix_cr_t(std::size_t n) - : diag(n) - , row_idx(n+1) - , col_idx(n*n) - , A(n*n) - , nz_num(0) - //, nzbd(n * (n+1) / 2) - , m_nzbd(n, n) - , m_size(n) - { - for (std::size_t i=0; i0) ? narrow_cast(N) : m_size; } - - void clear() noexcept - { - nz_num = 0; - for (std::size_t i=0; i < size() + 1; i++) - row_idx[i] = 0; - } - - void set_scalar(T scalar) noexcept - { - for (std::size_t i=0, e=nz_num; iri; i--) - { - A[i] = A[i-1]; - col_idx[i] = col_idx[i-1]; - } - A[ri] = val; - col_idx[ri] = c; - for (C i = r + 1; i < size() + 1; i++) - row_idx[i]++; - nz_num++; - if (c==r) - diag[r] = ri; - } - } - - template - void build_from_fill_mat(const M &f, std::size_t max_fill = FILL_INFINITY - 1, - std::size_t band_width = FILL_INFINITY) noexcept(false) - { - C nz = 0; - if (nz_num != 0) - throw pexception("build_from_mat only allowed on empty CR matrix"); - for (std::size_t k=0; k < size(); k++) - { - row_idx[k] = nz; - - for (std::size_t j=0; j < size(); j++) - if (f[k][j] <= max_fill && plib::abs(narrow_cast(k)-narrow_cast(j)) <= narrow_cast(band_width)) - { - col_idx[nz] = narrow_cast(j); - if (j == k) - diag[k] = nz; - nz++; - } - } - - row_idx[size()] = nz; - nz_num = nz; - - // build nzbd - - for (std::size_t k=0; k < size(); k++) - { -#if 0 - for (std::size_t j=k + 1; j < size(); j++) - if (f[j][k] < FILL_INFINITY) - m_nzbd[k].push_back(narrow_cast(j)); - m_nzbd[k].push_back(0); // end of sequence -#else - for (std::size_t j=k + 1; j < size(); j++) - if (f[j][k] < FILL_INFINITY) - m_nzbd.set(k, m_nzbd.colcount(k), narrow_cast(j)); - m_nzbd.set(k, m_nzbd.colcount(k), 0); // end of sequence -#endif - } - - } - - template - void mult_vec(VTR & res, const VTV & x) const noexcept - { - - // res = A * x -#if 0 - //plib::omp::set_num_threads(4); - plib::omp::for_static(0, constants::zero(), m_size, [this, &res, &x](std::size_t row) - { - T tmp(0.0); - const std::size_t e(row_idx[row+1]); - for (std::size_t k = row_idx[row]; k < e; k++) - tmp += A[k] * x[col_idx[k]]; - res[row] = tmp; - }); -#else - // this is a bit faster than the version above - std::size_t row = 0; - std::size_t k = 0; - const std::size_t oe = nz_num; - while (k < oe) - { - T tmp = plib::constants::zero(); - const std::size_t e = row_idx[row+1]; - for (; k < e; k++) - tmp += A[k] * x[col_idx[k]]; - res[row++] = tmp; - } -#endif - } - - // throws error if P(source)>P(destination) - template - void slim_copy_from(LUMAT & src) noexcept(false) - { - for (std::size_t r=0; r - void reduction_copy_from(LUMAT & src) noexcept - { - C sp(0); - for (std::size_t r=0; r - void raw_copy_from(LUMAT & src) noexcept - { - for (std::size_t k = 0; k < nz_num; k++) - A[k] = src.A[k]; - } - - index_type * nzbd(std::size_t row) { return m_nzbd[row]; } - std::size_t nzbd_count(std::size_t row) { return m_nzbd.colcount(row) - 1; } - protected: - // FIXME: this should be private - // NOLINTNEXTLINE - //parray, N > m_nzbd; // Support for gaussian elimination - pmatrix2d_vrl m_nzbd; // Support for gaussian elimination - private: - //parray nzbd; // Support for gaussian elimination - std::size_t m_size; - }; - - template - struct pGEmatrix_cr_t : public B - { - using base = B; - using index_type = typename base::index_type; - - PCOPYASSIGNMOVE(pGEmatrix_cr_t, default) - - explicit pGEmatrix_cr_t(std::size_t n) - : B(n) - { - } - - ~pGEmatrix_cr_t() = default; - - template - std::pair gaussian_extend_fill_mat(M &fill) - { - std::size_t ops = 0; - std::size_t fill_max = 0; - - for (std::size_t k = 0; k < fill.size(); k++) - { - ops++; // 1/A(k,k) - for (std::size_t row = k + 1; row < fill.size(); row++) - { - if (fill[row][k] < base::FILL_INFINITY) - { - ops++; - for (std::size_t col = k + 1; col < fill[row].size(); col++) - //if (fill[k][col] < FILL_INFINITY) - { - auto f = std::min(fill[row][col], 1 + fill[row][k] + fill[k][col]); - if (f < base::FILL_INFINITY) - { - if (f > fill_max) - fill_max = f; - ops += 2; - } - fill[row][col] = f; - } - } - } - } - build_parallel_gaussian_execution_scheme(fill); - return { fill_max, ops }; - } - - template - void gaussian_elimination(V & RHS) - { - const std::size_t iN = base::size(); - - for (std::size_t i = 0; i < iN - 1; i++) - { - std::size_t nzbdp = 0; - std::size_t pi = base::diag[i]; - auto f = reciprocal(base::A[pi++]); - const std::size_t piie = base::row_idx[i+1]; - - const auto *nz = base::m_nzbd[i]; - while (auto j = nz[nzbdp++]) // NOLINT(bugprone-infinite-loop) - { - // proceed to column i - - std::size_t pj = base::row_idx[j]; - std::size_t pje = base::row_idx[j+1]; - - while (base::col_idx[pj] < i) - pj++; - - const typename base::value_type f1 = - base::A[pj++] * f; - - // subtract row i from j - // fill-in available assumed, i.e. matrix was prepared - - for (std::size_t pii = pi; pii(i); - return -1; - } - - template - void gaussian_elimination_parallel(V & RHS) - { - //printf("omp: %ld %d %d\n", m_ge_par.size(), nz_num, (int)m_ge_par[m_ge_par.size()-2].size()); - for (auto l = 0UL; l < m_ge_par.size(); l++) - plib::omp::for_static(base::nz_num, 0UL, m_ge_par[l].size(), [this, &RHS, &l] (unsigned ll) - { - auto &i = m_ge_par[l][ll]; - { - std::size_t nzbdp = 0; - std::size_t pi = base::diag[i]; - const auto f = reciprocal(base::A[pi++]); - const std::size_t piie = base::row_idx[i+1]; - const auto &nz = base::nzbd[i]; - - while (auto j = nz[nzbdp++]) - { - // proceed to column i - - std::size_t pj = base::row_idx[j]; - - while (base::col_idx[pj] < i) - pj++; - - auto f1 = - base::A[pj++] * f; - - // subtract row i from j - // fill-in available assumed, i.e. matrix was prepared - for (std::size_t pii = pi; pii - void gaussian_back_substitution(V1 &V, const V2 &RHS) - { - const std::size_t iN = base::size(); - // row n-1 - V[iN - 1] = RHS[iN - 1] / base::A[base::diag[iN - 1]]; - - for (std::size_t j = iN - 1; j-- > 0;) - { - typename base::value_type tmp = 0; - const auto jdiag = base::diag[j]; - const std::size_t e = base::row_idx[j+1]; - for (std::size_t pk = jdiag + 1; pk < e; pk++) - tmp += base::A[pk] * V[base::col_idx[pk]]; - V[j] = (RHS[j] - tmp) / base::A[jdiag]; - } - } - - template - void gaussian_back_substitution(V1 &V) - { - const std::size_t iN = base::size(); - // row n-1 - V[iN - 1] = V[iN - 1] / base::A[base::diag[iN - 1]]; - - for (std::size_t j = iN - 1; j-- > 0;) - { - typename base::value_type tmp = 0; - const auto jdiag = base::diag[j]; - const std::size_t e = base::row_idx[j+1]; - for (std::size_t pk = jdiag + 1; pk < e; pk++) - tmp += base::A[pk] * V[base::col_idx[pk]]; - V[j] = (V[j] - tmp) / base::A[jdiag]; - } - } - - private: - template - void build_parallel_gaussian_execution_scheme(const M &fill) - { - // calculate parallel scheme for gaussian elimination - std::vector> rt(base::size()); - for (std::size_t k = 0; k < base::size(); k++) - { - for (std::size_t j = k+1; j < base::size(); j++) - { - if (fill[j][k] < base::FILL_INFINITY) - { - rt[k].push_back(j); - } - } - } - - std::vector levGE(base::size(), 0); - std::size_t cl = 0; - - for (std::size_t k = 0; k < base::size(); k++ ) - { - if (levGE[k] >= cl) - { - std::vector t = rt[k]; - for (std::size_t j = k+1; j < base::size(); j++ ) - { - bool overlap = false; - // is there overlap - if (plib::container::contains(t, j)) - overlap = true; - for (auto &x : rt[j]) - if (plib::container::contains(t, x)) - { - overlap = true; - break; - } - if (overlap) - levGE[j] = cl + 1; - else - { - t.push_back(j); - for (auto &x : rt[j]) - t.push_back(x); - } - } - cl++; - } - } - - m_ge_par.clear(); - m_ge_par.resize(cl+1); - for (std::size_t k = 0; k < base::size(); k++) - m_ge_par[levGE[k]].push_back(k); - //for (std::size_t k = 0; k < m_ge_par.size(); k++) - // printf("%d %d\n", (int) k, (int) m_ge_par[k].size()); - } - std::vector> m_ge_par; // parallel execution support for Gauss - }; - - template - struct pLUmatrix_cr_t : public B - { - using base = B; - using index_type = typename base::index_type; - - PCOPYASSIGNMOVE(pLUmatrix_cr_t, default) - - explicit pLUmatrix_cr_t(std::size_t n) - : B(n) - , ilu_rows(n+1) - , m_ILUp(0) - { - } - - ~pLUmatrix_cr_t() = default; - - template - void build(M &fill, std::size_t ilup) - { - std::size_t p(0); - // build ilu_rows - for (decltype(fill.size()) i=1; i < fill.size(); i++) - { - bool found(false); - for (decltype(fill.size()) k = 0; k < i; k++) - { - // if (fill[i][k] < base::FILL_INFINITY) - if (fill[i][k] <= ilup) - { - // assume A[k][k]!=0 - for (decltype(fill.size()) j=k+1; j < fill.size(); j++) - { - auto f = std::min(fill[i][j], 1 + fill[i][k] + fill[k][j]); - if (f <= ilup) - fill[i][j] = f; - } - found = true; - } - } - if (found) - ilu_rows[p++] = narrow_cast(i); - } - ilu_rows[p] = 0; // end of array - this->build_from_fill_mat(fill, ilup); //, m_band_width); // ILU(2) - m_ILUp = ilup; - } - - /// \brief incomplete LU Factorization. - /// - /// We are following http://de.wikipedia.org/wiki/ILU-Zerlegung here. - /// - /// The result is stored in matrix LU - /// - /// For i = 1,...,N-1 - /// For k = 0, ... , i - 1 - /// If a[i,k] != 0 - /// a[i,k] = a[i,k] / a[k,k] - /// For j = k + 1, ... , N - 1 - /// If a[i,j] != 0 - /// a[i,j] = a[i,j] - a[i,k] * a[k,j] - /// j=j+1 - /// k=k+1 - /// i=i+1 - /// - void incomplete_LU_factorization(const base &mat) - { - if (m_ILUp < 1) - this->raw_copy_from(mat); - else - this->reduction_copy_from(mat); - - std::size_t p(0); - while (auto i = ilu_rows[p++]) // NOLINT(bugprone-infinite-loop) - { - const auto p_i_end = base::row_idx[i + 1]; - // loop over all columns k left of diag in row i - //if (row_idx[i] < diag[i]) - // printf("occ %d\n", (int)i); - for (auto i_k = base::row_idx[i]; i_k < base::diag[i]; i_k++) - { - const auto k(base::col_idx[i_k]); - const auto p_k_end(base::row_idx[k + 1]); - const typename base::value_type LUp_i_k = base::A[i_k] = base::A[i_k] / base::A[base::diag[k]]; - - std::size_t k_j(base::diag[k] + 1); - std::size_t i_j(i_k + 1); - - while (i_j < p_i_end && k_j < p_k_end ) // pj = (i, j) - { - // we can assume that within a row ja increases continuously - const std::size_t c_i_j(base::col_idx[i_j]); // row i, column j - const auto c_k_j(base::col_idx[k_j]); // row k, column j - - if (c_k_j == c_i_j) - base::A[i_j] -= LUp_i_k * base::A[k_j]; - k_j += (c_k_j <= c_i_j ? 1 : 0); - i_j += (c_k_j >= c_i_j ? 1 : 0); - - } - } - } - } - - - /// \brief Solve a linear equation - /// - /// Solve a linear equation Ax = r - /// - /// where - /// A = L*U - /// - /// L unit lower triangular - /// U upper triangular - /// - /// ==> LUx = r - /// - /// ==> Ux = L⁻¹ r = w - /// - /// ==> r = Lw - /// - /// This can be solved for w using backwards elimination in L. - /// - /// Now Ux = w - /// - /// This can be solved for x using backwards elimination in U. - /// - template - void solveLU (R &r) - { - for (std::size_t i = 1; i < base::size(); ++i ) - { - typename base::value_type tmp(0); - const auto j1(base::row_idx[i]); - const auto j2(base::diag[i]); - - for (auto j = j1; j < j2; ++j ) - tmp += base::A[j] * r[base::col_idx[j]]; - r[i] -= tmp; - } - // i now is equal to n; - for (std::size_t i = base::size(); i-- > 0; ) - { - typename base::value_type tmp(0); - const auto di(base::diag[i]); - const auto j2(base::row_idx[i+1]); - for (std::size_t j = di + 1; j < j2; j++ ) - tmp += base::A[j] * r[base::col_idx[j]]; - r[i] = (r[i] - tmp) / base::A[di]; - } - } - private: - parray ilu_rows; - std::size_t m_ILUp; - }; - -} // namespace plib - -#endif // MAT_CR_H_ diff --git a/src/lib/netlist/plib/pgsl.h b/src/lib/netlist/plib/pgsl.h index ab7798071a5..590765429eb 100644 --- a/src/lib/netlist/plib/pgsl.h +++ b/src/lib/netlist/plib/pgsl.h @@ -25,6 +25,8 @@ #if defined(__has_builtin) // clang and gcc 10 #if __has_builtin(__builtin_unreachable) #define gsl_Expects(e) ((e) ? static_cast(0) : __builtin_unreachable()) + #else + #define gsl_Expects(e) ((e) ? static_cast(0) : static_cast(0)) #endif #elif defined(__GNUC__) && !(defined( __CUDACC__ ) && defined( __CUDA_ARCH__ )) #define gsl_Expects(e) ((e) ? static_cast(0) : __builtin_unreachable()) @@ -34,8 +36,10 @@ #define gsl_Expects(e) ((e) ? static_cast(0) : static_cast(0)) #endif +#if 0 #undef gsl_Expects #define gsl_Expects(e) do {} while (0) +#endif #define gsl_Ensures(e) gsl_Expects(e) namespace plib { diff --git a/src/lib/netlist/plib/pmatrix_cr.h b/src/lib/netlist/plib/pmatrix_cr.h new file mode 100644 index 00000000000..38457650e5e --- /dev/null +++ b/src/lib/netlist/plib/pmatrix_cr.h @@ -0,0 +1,635 @@ +// license:GPL-2.0+ +// copyright-holders:Couriersud + +#ifndef PMATRIX_CR_H_ +#define PMATRIX_CR_H_ + +/// +/// \file pmatrix_cr.h +/// +/// Compressed row format matrices +/// + +#include "palloc.h" +#include "parray.h" +#include "pconfig.h" +#include "pmath.h" +#include "pmatrix2d.h" +#include "pomp.h" +#include "ptypes.h" + +#include +#include +#include +#include + +namespace plib +{ + + template + struct pmatrix_cr + { + using index_type = C; + using value_type = T; + + static constexpr const int NSQ = (N < 0 ? -N * N : N * N); + static constexpr const int Np1 = (N == 0) ? 0 : (N < 0 ? N - 1 : N + 1); + + PCOPYASSIGNMOVE(pmatrix_cr, default) + + enum constants_e + { + FILL_INFINITY = 9999999 + }; + + // FIXME: these should be private + // NOLINTNEXTLINE + parray diag; // diagonal index pointer n + // NOLINTNEXTLINE + parray row_idx; // row index pointer n + 1 + // NOLINTNEXTLINE + parray col_idx; // column index array nz_num, initially (n * n) + // NOLINTNEXTLINE + parray A; // Matrix elements nz_num, initially (n * n) + // NOLINTNEXTLINE + std::size_t nz_num; + + explicit pmatrix_cr(std::size_t n) + : diag(n) + , row_idx(n+1) + , col_idx(n*n) + , A(n*n) + , nz_num(0) + //, nzbd(n * (n+1) / 2) + , m_nzbd(n, n) + , m_size(n) + { + for (std::size_t i=0; i0) ? narrow_cast(N) : m_size; } + + void clear() noexcept + { + nz_num = 0; + for (std::size_t i=0; i < size() + 1; i++) + row_idx[i] = 0; + } + + void set_scalar(T scalar) noexcept + { + for (std::size_t i=0, e=nz_num; iri; i--) + { + A[i] = A[i-1]; + col_idx[i] = col_idx[i-1]; + } + A[ri] = val; + col_idx[ri] = c; + for (C i = r + 1; i < size() + 1; i++) + row_idx[i]++; + nz_num++; + if (c==r) + diag[r] = ri; + } + } + + template + void build_from_fill_mat(const M &f, std::size_t max_fill = FILL_INFINITY - 1, + std::size_t band_width = FILL_INFINITY) noexcept(false) + { + C nz = 0; + if (nz_num != 0) + throw pexception("build_from_mat only allowed on empty CR matrix"); + for (std::size_t k=0; k < size(); k++) + { + row_idx[k] = nz; + + for (std::size_t j=0; j < size(); j++) + if (f[k][j] <= max_fill && plib::abs(narrow_cast(k)-narrow_cast(j)) <= narrow_cast(band_width)) + { + col_idx[nz] = narrow_cast(j); + if (j == k) + diag[k] = nz; + nz++; + } + } + + row_idx[size()] = nz; + nz_num = nz; + + // build nzbd + + for (std::size_t k=0; k < size(); k++) + { + for (std::size_t j=k + 1; j < size(); j++) + if (f[j][k] < FILL_INFINITY) + m_nzbd.set(k, m_nzbd.colcount(k), narrow_cast(j)); + m_nzbd.set(k, m_nzbd.colcount(k), 0); // end of sequence + } + + } + + template + void mult_vec(VTR & res, const VTV & x) const noexcept + { + + // res = A * x + // this is a bit faster than the version above + std::size_t row = 0; + std::size_t k = 0; + const std::size_t oe = nz_num; + while (k < oe) + { + T tmp = plib::constants::zero(); + const std::size_t e = row_idx[row+1]; + for (; k < e; k++) + tmp += A[k] * x[col_idx[k]]; + res[row++] = tmp; + } + } + + // throws error if P(source)>P(destination) + template + void slim_copy_from(LUMAT & src) noexcept(false) + { + for (std::size_t r=0; r + void reduction_copy_from(LUMAT & src) noexcept + { + C sp(0); + for (std::size_t r=0; r + void raw_copy_from(LUMAT & src) noexcept + { + for (std::size_t k = 0; k < nz_num; k++) + A[k] = src.A[k]; + } + + index_type * nzbd(std::size_t row) { return m_nzbd[row]; } + std::size_t nzbd_count(std::size_t row) { return m_nzbd.colcount(row) - 1; } + protected: + // FIXME: this should be private + // NOLINTNEXTLINE + //parray, N > m_nzbd; // Support for gaussian elimination + pmatrix2d_vrl m_nzbd; // Support for gaussian elimination + private: + //parray nzbd; // Support for gaussian elimination + std::size_t m_size; + }; + + template + struct pGEmatrix_cr : public B + { + using base = B; + using index_type = typename base::index_type; + + PCOPYASSIGNMOVE(pGEmatrix_cr, default) + + explicit pGEmatrix_cr(std::size_t n) + : B(n) + { + } + + ~pGEmatrix_cr() = default; + + template + std::pair gaussian_extend_fill_mat(M &fill) + { + std::size_t ops = 0; + std::size_t fill_max = 0; + + for (std::size_t k = 0; k < fill.size(); k++) + { + ops++; // 1/A(k,k) + for (std::size_t row = k + 1; row < fill.size(); row++) + { + if (fill[row][k] < base::FILL_INFINITY) + { + ops++; + for (std::size_t col = k + 1; col < fill[row].size(); col++) + //if (fill[k][col] < FILL_INFINITY) + { + auto f = std::min(fill[row][col], 1 + fill[row][k] + fill[k][col]); + if (f < base::FILL_INFINITY) + { + if (f > fill_max) + fill_max = f; + ops += 2; + } + fill[row][col] = f; + } + } + } + } + build_parallel_gaussian_execution_scheme(fill); + return { fill_max, ops }; + } + + template + void gaussian_elimination(V & RHS) + { + const std::size_t iN = base::size(); + + for (std::size_t i = 0; i < iN - 1; i++) + { + std::size_t nzbdp = 0; + std::size_t pi = base::diag[i]; + auto f = reciprocal(base::A[pi++]); + const std::size_t piie = base::row_idx[i+1]; + + const auto *nz = base::m_nzbd[i]; + while (auto j = nz[nzbdp++]) // NOLINT(bugprone-infinite-loop) + { + // proceed to column i + + std::size_t pj = base::row_idx[j]; + std::size_t pje = base::row_idx[j+1]; + + while (base::col_idx[pj] < i) + pj++; + + const typename base::value_type f1 = - base::A[pj++] * f; + + // subtract row i from j + // fill-in available assumed, i.e. matrix was prepared + + for (std::size_t pii = pi; pii(i); + return -1; + } + + template + void gaussian_elimination_parallel(V & RHS) + { + //printf("omp: %ld %d %d\n", m_ge_par.size(), nz_num, (int)m_ge_par[m_ge_par.size()-2].size()); + for (auto l = 0UL; l < m_ge_par.size(); l++) + plib::omp::for_static(base::nz_num, 0UL, m_ge_par[l].size(), [this, &RHS, &l] (unsigned ll) + { + auto &i = m_ge_par[l][ll]; + { + std::size_t nzbdp = 0; + std::size_t pi = base::diag[i]; + const auto f = reciprocal(base::A[pi++]); + const std::size_t piie = base::row_idx[i+1]; + const auto &nz = base::nzbd[i]; + + while (auto j = nz[nzbdp++]) + { + // proceed to column i + + std::size_t pj = base::row_idx[j]; + + while (base::col_idx[pj] < i) + pj++; + + auto f1 = - base::A[pj++] * f; + + // subtract row i from j + // fill-in available assumed, i.e. matrix was prepared + for (std::size_t pii = pi; pii + void gaussian_back_substitution(V1 &V, const V2 &RHS) + { + const std::size_t iN = base::size(); + // row n-1 + V[iN - 1] = RHS[iN - 1] / base::A[base::diag[iN - 1]]; + + for (std::size_t j = iN - 1; j-- > 0;) + { + typename base::value_type tmp = 0; + const auto jdiag = base::diag[j]; + const std::size_t e = base::row_idx[j+1]; + for (std::size_t pk = jdiag + 1; pk < e; pk++) + tmp += base::A[pk] * V[base::col_idx[pk]]; + V[j] = (RHS[j] - tmp) / base::A[jdiag]; + } + } + + template + void gaussian_back_substitution(V1 &V) + { + const std::size_t iN = base::size(); + // row n-1 + V[iN - 1] = V[iN - 1] / base::A[base::diag[iN - 1]]; + + for (std::size_t j = iN - 1; j-- > 0;) + { + typename base::value_type tmp = 0; + const auto jdiag = base::diag[j]; + const std::size_t e = base::row_idx[j+1]; + for (std::size_t pk = jdiag + 1; pk < e; pk++) + tmp += base::A[pk] * V[base::col_idx[pk]]; + V[j] = (V[j] - tmp) / base::A[jdiag]; + } + } + + private: + template + void build_parallel_gaussian_execution_scheme(const M &fill) + { + // calculate parallel scheme for gaussian elimination + std::vector> rt(base::size()); + for (std::size_t k = 0; k < base::size(); k++) + { + for (std::size_t j = k+1; j < base::size(); j++) + { + if (fill[j][k] < base::FILL_INFINITY) + { + rt[k].push_back(j); + } + } + } + + std::vector levGE(base::size(), 0); + std::size_t cl = 0; + + for (std::size_t k = 0; k < base::size(); k++ ) + { + if (levGE[k] >= cl) + { + std::vector t = rt[k]; + for (std::size_t j = k+1; j < base::size(); j++ ) + { + bool overlap = false; + // is there overlap + if (plib::container::contains(t, j)) + overlap = true; + for (auto &x : rt[j]) + if (plib::container::contains(t, x)) + { + overlap = true; + break; + } + if (overlap) + levGE[j] = cl + 1; + else + { + t.push_back(j); + for (auto &x : rt[j]) + t.push_back(x); + } + } + cl++; + } + } + + m_ge_par.clear(); + m_ge_par.resize(cl+1); + for (std::size_t k = 0; k < base::size(); k++) + m_ge_par[levGE[k]].push_back(k); + //for (std::size_t k = 0; k < m_ge_par.size(); k++) + // printf("%d %d\n", (int) k, (int) m_ge_par[k].size()); + } + std::vector> m_ge_par; // parallel execution support for Gauss + }; + + template + struct pLUmatrix_cr : public B + { + using base = B; + using index_type = typename base::index_type; + + PCOPYASSIGNMOVE(pLUmatrix_cr, default) + + explicit pLUmatrix_cr(std::size_t n) + : B(n) + , ilu_rows(n+1) + , m_ILUp(0) + { + } + + ~pLUmatrix_cr() = default; + + template + void build(M &fill, std::size_t ilup) + { + std::size_t p(0); + // build ilu_rows + for (decltype(fill.size()) i=1; i < fill.size(); i++) + { + bool found(false); + for (decltype(fill.size()) k = 0; k < i; k++) + { + // if (fill[i][k] < base::FILL_INFINITY) + if (fill[i][k] <= ilup) + { + // assume A[k][k]!=0 + for (decltype(fill.size()) j=k+1; j < fill.size(); j++) + { + auto f = std::min(fill[i][j], 1 + fill[i][k] + fill[k][j]); + if (f <= ilup) + fill[i][j] = f; + } + found = true; + } + } + if (found) + ilu_rows[p++] = narrow_cast(i); + } + ilu_rows[p] = 0; // end of array + this->build_from_fill_mat(fill, ilup); //, m_band_width); // ILU(2) + m_ILUp = ilup; + } + + /// \brief incomplete LU Factorization. + /// + /// We are following http://de.wikipedia.org/wiki/ILU-Zerlegung here. + /// + /// The result is stored in matrix LU + /// + /// For i = 1,...,N-1 + /// For k = 0, ... , i - 1 + /// If a[i,k] != 0 + /// a[i,k] = a[i,k] / a[k,k] + /// For j = k + 1, ... , N - 1 + /// If a[i,j] != 0 + /// a[i,j] = a[i,j] - a[i,k] * a[k,j] + /// j=j+1 + /// k=k+1 + /// i=i+1 + /// + void incomplete_LU_factorization(const base &mat) + { + if (m_ILUp < 1) + this->raw_copy_from(mat); + else + this->reduction_copy_from(mat); + + std::size_t p(0); + while (auto i = ilu_rows[p++]) // NOLINT(bugprone-infinite-loop) + { + const auto p_i_end = base::row_idx[i + 1]; + // loop over all columns k left of diag in row i + //if (row_idx[i] < diag[i]) + // printf("occ %d\n", (int)i); + for (auto i_k = base::row_idx[i]; i_k < base::diag[i]; i_k++) + { + const auto k(base::col_idx[i_k]); + const auto p_k_end(base::row_idx[k + 1]); + const typename base::value_type LUp_i_k = base::A[i_k] = base::A[i_k] / base::A[base::diag[k]]; + + std::size_t k_j(base::diag[k] + 1); + std::size_t i_j(i_k + 1); + + while (i_j < p_i_end && k_j < p_k_end ) // pj = (i, j) + { + // we can assume that within a row ja increases continuously + const std::size_t c_i_j(base::col_idx[i_j]); // row i, column j + const auto c_k_j(base::col_idx[k_j]); // row k, column j + + if (c_k_j == c_i_j) + base::A[i_j] -= LUp_i_k * base::A[k_j]; + k_j += (c_k_j <= c_i_j ? 1 : 0); + i_j += (c_k_j >= c_i_j ? 1 : 0); + + } + } + } + } + + + /// \brief Solve a linear equation + /// + /// Solve a linear equation Ax = r + /// + /// where + /// A = L*U + /// + /// L unit lower triangular + /// U upper triangular + /// + /// ==> LUx = r + /// + /// ==> Ux = L⁻¹ r = w + /// + /// ==> r = Lw + /// + /// This can be solved for w using backwards elimination in L. + /// + /// Now Ux = w + /// + /// This can be solved for x using backwards elimination in U. + /// + template + void solveLU (R &r) + { + for (std::size_t i = 1; i < base::size(); ++i ) + { + typename base::value_type tmp(0); + const auto j1(base::row_idx[i]); + const auto j2(base::diag[i]); + + for (auto j = j1; j < j2; ++j ) + tmp += base::A[j] * r[base::col_idx[j]]; + r[i] -= tmp; + } + // i now is equal to n; + for (std::size_t i = base::size(); i-- > 0; ) + { + typename base::value_type tmp(0); + const auto di(base::diag[i]); + const auto j2(base::row_idx[i+1]); + for (std::size_t j = di + 1; j < j2; j++ ) + tmp += base::A[j] * r[base::col_idx[j]]; + r[i] = (r[i] - tmp) / base::A[di]; + } + } + private: + parray ilu_rows; + std::size_t m_ILUp; + }; + +} // namespace plib + +#endif // PMATRIX_CR_H_ diff --git a/src/lib/netlist/plib/vector_ops.h b/src/lib/netlist/plib/vector_ops.h index 6d20b9213a7..c28d97a3b97 100644 --- a/src/lib/netlist/plib/vector_ops.h +++ b/src/lib/netlist/plib/vector_ops.h @@ -12,6 +12,7 @@ /// #include "pconfig.h" #include "pmath.h" +#include "pgsl.h" #include #include @@ -27,75 +28,76 @@ namespace plib { template - void vec_set_scalar(const std::size_t n, VT &v, T && scalar) noexcept + void vec_set_scalar(VT &v, T && scalar) noexcept { + const std::size_t n(v.size()); const typename std::remove_reference::type s(std::forward(scalar)); for ( std::size_t i = 0; i < n; i++ ) v[i] = s; } template - void vec_set(const std::size_t n, VT &v, const VS & source) noexcept + void vec_set(VT &v, const VS & source) noexcept { + const std::size_t n(v.size()); + gsl_Expects(n <= source.size()); for ( std::size_t i = 0; i < n; i++ ) v[i] = source[i]; } template - T vec_mult(const std::size_t n, const V1 & v1, const V2 & v2 ) noexcept + T vec_mult(const V1 & v1, const V2 & v2 ) noexcept { - using b8 = std::array; // NOLINT - PALIGNAS_VECTOROPT() b8 value = {0}; + const std::size_t n(v1.size()); + gsl_Expects(n <= v2.size()); + T ret(plib::constants::zero()); for (std::size_t i = 0; i < n ; i++ ) { - value[i & 7] += v1[i] * v2[i]; // NOLINT + ret += v1[i] * v2[i]; // NOLINT } - return value[0] + value[1] + value[2] + value[3] + value[4] + value[5] + value[6] + value[7]; // NOLINT + return ret; } template - T vec_mult2(const std::size_t n, const VT &v) noexcept + T vec_mult2(const VT &v) noexcept { - using b8 = std::array; - PALIGNAS_VECTOROPT() b8 value = {0}; + const std::size_t n(v.size()); + T ret(plib::constants::zero()); for (std::size_t i = 0; i < n ; i++ ) { - value[i & 7] += v[i] * v[i]; + ret += v[i] * v[i]; } - return value[0] + value[1] + value[2] + value[3] + value[4] + value[5] + value[6] + value[7]; + return ret; } template - T vec_sum(const std::size_t n, const VT &v) noexcept + T vec_sum(const VT &v) noexcept { - if (n<8) + const std::size_t n(v.size()); + T ret(plib::constants::zero()); + for (std::size_t i = 0; i < n ; i++ ) { - T value(0); - for (std::size_t i = 0; i < n ; i++ ) - value += v[i]; - - return value; + ret += v[i]; } - using b8 = std::array; - PALIGNAS_VECTOROPT() b8 value = {0}; - for (std::size_t i = 0; i < n ; i++ ) - value[i & 7] += v[i]; - - return ((value[0] + value[1]) + (value[2] + value[3])) + ((value[4] + value[5]) + (value[6] + value[7])); + return ret; } template - void vec_mult_scalar(const std::size_t n, VR & result, const VV & v, T && scalar) noexcept + void vec_mult_scalar(VR & result, const VV & v, T && scalar) noexcept { + const std::size_t n(result.size()); + gsl_Expects(n <= v.size()); const typename std::remove_reference::type s(std::forward(scalar)); for ( std::size_t i = 0; i < n; i++ ) result[i] = s * v[i]; } template - void vec_add_mult_scalar(const std::size_t n, VR & result, const VV & v, T && scalar) noexcept + void vec_add_mult_scalar(VR & result, const VV & v, T && scalar) noexcept { + const std::size_t n(result.size()); + gsl_Expects(n <= v.size()); const typename std::remove_reference::type s(std::forward(scalar)); for ( std::size_t i = 0; i < n; i++ ) result[i] += s * v[i]; @@ -109,31 +111,38 @@ namespace plib } template - void vec_add_ip(const std::size_t n, R & result, const V & v) noexcept + void vec_add_ip(R & result, const V & v) noexcept { + const std::size_t n(result.size()); + gsl_Expects(n <= v.size()); for ( std::size_t i = 0; i < n; i++ ) result[i] += v[i]; } template - void vec_sub(const std::size_t n, VR & result, const V1 &v1, const V2 & v2) noexcept + void vec_sub(VR & result, const V1 &v1, const V2 & v2) noexcept { + const std::size_t n(result.size()); + gsl_Expects(n <= v1.size()); + gsl_Expects(n <= v2.size()); for ( std::size_t i = 0; i < n; i++ ) result[i] = v1[i] - v2[i]; } template - void vec_scale(const std::size_t n, V & v, T &&scalar) noexcept + void vec_scale(V & v, T &&scalar) noexcept { + const std::size_t n(v.size()); const typename std::remove_reference::type s(std::forward(scalar)); for ( std::size_t i = 0; i < n; i++ ) v[i] *= s; } template - T vec_maxabs(const std::size_t n, const V & v) noexcept + T vec_maxabs(const V & v) noexcept { - T ret = 0.0; + const std::size_t n(v.size()); + T ret(plib::constants::zero()); for ( std::size_t i = 0; i < n; i++ ) ret = std::max(ret, plib::abs(v[i])); diff --git a/src/lib/netlist/solver/nld_matrix_solver.h b/src/lib/netlist/solver/nld_matrix_solver.h index 4e392c88770..94b05d42b00 100644 --- a/src/lib/netlist/solver/nld_matrix_solver.h +++ b/src/lib/netlist/solver/nld_matrix_solver.h @@ -10,7 +10,7 @@ #include "nl_base.h" #include "nl_errstr.h" -#include "plib/mat_cr.h" +#include "plib/pmatrix_cr.h" #include "plib/palloc.h" #include "plib/penum.h" #include "plib/pmatrix2d.h" @@ -74,7 +74,7 @@ namespace solver constexpr nl_fptype m_vntol() { return nlconst::magic(1e-7); } constexpr nl_fptype m_accuracy() { return nlconst::magic(1e-7); } constexpr std::size_t m_nr_loops() { return 250; } - constexpr std::size_t m_gs_loops() { return 9; } + constexpr std::size_t m_gs_loops() { return 50; } // general parameters constexpr nl_fptype m_gmin() { return nlconst::magic(1e-9); } diff --git a/src/lib/netlist/solver/nld_ms_gcr.h b/src/lib/netlist/solver/nld_ms_gcr.h index 61bfa012fa7..5d8440e5d1d 100644 --- a/src/lib/netlist/solver/nld_ms_gcr.h +++ b/src/lib/netlist/solver/nld_ms_gcr.h @@ -10,11 +10,9 @@ /// Gaussian elimination using compressed row format. /// -#include "plib/mat_cr.h" - -//#include "nld_ms_direct.h" #include "nld_matrix_solver_ext.h" #include "nld_solver.h" +#include "plib/pmatrix_cr.h" #include "plib/pdynlib.h" #include "plib/pstream.h" #include "plib/vector_ops.h" @@ -31,7 +29,7 @@ namespace solver { public: - using mat_type = plib::pGEmatrix_cr_t>; + using mat_type = plib::pGEmatrix_cr>; using base_type = matrix_solver_ext_t; using fptype = typename base_type::fptype; @@ -119,7 +117,7 @@ namespace solver private: - using mat_index_type = typename plib::pmatrix_cr_t::index_type; + using mat_index_type = typename plib::pmatrix_cr::index_type; void generate_code(plib::putf8_fmt_writer &strm); diff --git a/src/lib/netlist/solver/nld_ms_gmres.h b/src/lib/netlist/solver/nld_ms_gmres.h index c702d2e3ba6..fd84a50c33d 100644 --- a/src/lib/netlist/solver/nld_ms_gmres.h +++ b/src/lib/netlist/solver/nld_ms_gmres.h @@ -12,7 +12,7 @@ #include "nld_ms_direct.h" #include "nld_solver.h" #include "plib/gmres.h" -#include "plib/mat_cr.h" +#include "plib/pmatrix_cr.h" #include "plib/parray.h" #include "plib/vector_ops.h" @@ -83,7 +83,7 @@ namespace solver private: - using mattype = typename plib::pmatrix_cr_t::index_type; + using mattype = typename plib::pmatrix_cr::index_type; //plib::mat_precondition_none m_ops; plib::mat_precondition_ILU m_ops; -- cgit v1.2.3-70-g09d2