diff options
author | 2020-09-20 08:04:22 +0200 | |
---|---|---|
committer | 2020-09-20 17:34:54 +0200 | |
commit | 63efcfcb85e2f475b26541d0973b0321c9204598 (patch) | |
tree | 63a3862e71012a254bb17a4d8d099a7a6d16b91d | |
parent | 9e4377dc4f2f22e0aecc9daaf043eedf398e38e2 (diff) |
netlist: code refactoring.
- rename mat_cr.h to pmatrix_cr.h
- Optimization to the gmres solver.
- Simplifcation of vector operation code
-rw-r--r-- | scripts/src/netlist.lua | 2 | ||||
-rw-r--r-- | src/lib/netlist/plib/gmres.h | 189 | ||||
-rw-r--r-- | src/lib/netlist/plib/pgsl.h | 4 | ||||
-rw-r--r-- | src/lib/netlist/plib/pmatrix_cr.h (renamed from src/lib/netlist/plib/mat_cr.h) | 53 | ||||
-rw-r--r-- | src/lib/netlist/plib/vector_ops.h | 73 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_matrix_solver.h | 4 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_ms_gcr.h | 8 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_ms_gmres.h | 4 |
8 files changed, 167 insertions, 170 deletions
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 <typename FT, int SIZE> struct mat_precondition_ILU { - using mat_type = plib::pmatrix_cr_t<FT, SIZE>; - using matLU_type = plib::pLUmatrix_cr_t<mat_type>; + using mat_type = plib::pmatrix_cr<FT, SIZE>; + using matLU_type = plib::pLUmatrix_cr<mat_type>; mat_precondition_ILU(std::size_t size, std::size_t ilu_scale = 4 - , std::size_t bw = plib::pmatrix_cr_t<FT, SIZE>::FILL_INFINITY) + , std::size_t bw = plib::pmatrix_cr<FT, SIZE>::FILL_INFINITY) : m_mat(narrow_cast<typename mat_type::index_type>(size)) , m_LU(narrow_cast<typename mat_type::index_type>(size)) , m_ILU_scale(narrow_cast<std::size_t>(ilu_scale)) @@ -167,7 +167,7 @@ namespace plib v[i] = v[i] * m_diag[i]; } - plib::pmatrix_cr_t<FT, SIZE> m_mat; + plib::pmatrix_cr<FT, SIZE> m_mat; plib::parray<FT, SIZE> m_diag; plib::parray<std::vector<std::size_t>, SIZE > nzcol; }; @@ -203,19 +203,23 @@ namespace plib plib::unused_var(v); } - plib::pmatrix_cr_t<FT, SIZE> m_mat; + plib::pmatrix_cr<FT, SIZE> m_mat; }; // FIXME: hardcoding RESTART to 20 becomes an issue on very large // systems. - template <typename FT, int SIZE, int RESTART = 16> + template <typename FT, int SIZE, int RESTARTMAX = 16> 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<std::size_t>(SIZE); } - template <int k, typename OPS, typename VT> - bool do_k(OPS &ops, VT &x, std::size_t &itr_used, FT rho_delta, bool dummy) - { - plib::unused_var(dummy); - if (do_k<k-1, OPS>(ops, x, itr_used, rho_delta, do_khelper<k-1>::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<FT>(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<FT>(n, m_v[kp1])); - - // FIXME: comparison to zero - if (m_ht[kp1][k] != plib::constants<FT>::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<FT>::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 <int k, typename OPS, typename VT> - 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 <typename OPS, typename VT, typename VRHS> std::size_t solve(OPS &ops, VT &x, const VRHS & rhs, const std::size_t itr_max, float_type accuracy) { @@ -348,19 +277,23 @@ 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<FT>(n, Ax)) / accuracy; + const float_type rho_to_accuracy = plib::sqrt(vec_mult2<FT>(Ax)) / accuracy; rho_delta = accuracy * rho_to_accuracy; } else + //rho_delta = accuracy * plib::sqrt(vec_mult2<FT>(n, rhs)) + // + 1e-4 * std::sqrt(n); rho_delta = accuracy * plib::sqrt(narrow_cast<FT>(n)); // + // LU x = b; solve for x; + // // Using // // vec_set(n, x, rhs); @@ -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<FT>(n, residual)); + rho = plib::sqrt(vec_mult2<FT>(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<FT>::zero()); + vec_set_scalar(m_g, +constants<FT>::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<RESTART-1>(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 <int k, typename OPS, typename VT> + bool do_k(OPS &ops, VT &x, std::size_t &itr_used, FT rho_delta, bool dummy) + { + plib::unused_var(dummy); + if (do_k<k-1, OPS>(ops, x, itr_used, rho_delta, do_khelper<k-1>::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<FT>(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<FT>(m_v[kp1])); + + // FIXME: comparison to zero + if (m_ht[kp1][k] != plib::constants<FT>::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<FT>::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 <int k, typename OPS, typename VT> + 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<float_type, SIZE> residual; plib::parray<float_type, SIZE> Ax; 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<void>(0) : __builtin_unreachable()) + #else + #define gsl_Expects(e) ((e) ? static_cast<void>(0) : static_cast<void>(0)) #endif #elif defined(__GNUC__) && !(defined( __CUDACC__ ) && defined( __CUDA_ARCH__ )) #define gsl_Expects(e) ((e) ? static_cast<void>(0) : __builtin_unreachable()) @@ -34,8 +36,10 @@ #define gsl_Expects(e) ((e) ? static_cast<void>(0) : static_cast<void>(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/mat_cr.h b/src/lib/netlist/plib/pmatrix_cr.h index 17ce140c6e5..38457650e5e 100644 --- a/src/lib/netlist/plib/mat_cr.h +++ b/src/lib/netlist/plib/pmatrix_cr.h @@ -1,11 +1,11 @@ // license:GPL-2.0+ // copyright-holders:Couriersud -#ifndef MAT_CR_H_ -#define MAT_CR_H_ +#ifndef PMATRIX_CR_H_ +#define PMATRIX_CR_H_ /// -/// \file mat_cr.h +/// \file pmatrix_cr.h /// /// Compressed row format matrices /// @@ -16,9 +16,7 @@ #include "pmath.h" #include "pmatrix2d.h" #include "pomp.h" -#include "pstate.h" #include "ptypes.h" -#include "putil.h" #include <algorithm> #include <array> @@ -29,7 +27,7 @@ namespace plib { template<typename T, int N, typename C = uint16_t> - struct pmatrix_cr_t + struct pmatrix_cr { using index_type = C; using value_type = T; @@ -37,7 +35,7 @@ namespace plib 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) + PCOPYASSIGNMOVE(pmatrix_cr, default) enum constants_e { @@ -56,7 +54,7 @@ namespace plib // NOLINTNEXTLINE std::size_t nz_num; - explicit pmatrix_cr_t(std::size_t n) + explicit pmatrix_cr(std::size_t n) : diag(n) , row_idx(n+1) , col_idx(n*n) @@ -72,7 +70,7 @@ namespace plib } } - ~pmatrix_cr_t() = default; + ~pmatrix_cr() = default; constexpr std::size_t size() const noexcept { return (N>0) ? narrow_cast<std::size_t>(N) : m_size; } @@ -149,17 +147,10 @@ namespace plib 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<C>(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<C>(j)); m_nzbd.set(k, m_nzbd.colcount(k), 0); // end of sequence -#endif } } @@ -169,17 +160,6 @@ namespace plib { // res = A * x -#if 0 - //plib::omp::set_num_threads(4); - plib::omp::for_static(0, constants<std::size_t>::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; @@ -192,7 +172,6 @@ namespace plib tmp += A[k] * x[col_idx[k]]; res[row++] = tmp; } -#endif } // throws error if P(source)>P(destination) @@ -262,19 +241,19 @@ namespace plib }; template<typename B> - struct pGEmatrix_cr_t : public B + struct pGEmatrix_cr : public B { using base = B; using index_type = typename base::index_type; - PCOPYASSIGNMOVE(pGEmatrix_cr_t, default) + PCOPYASSIGNMOVE(pGEmatrix_cr, default) - explicit pGEmatrix_cr_t(std::size_t n) + explicit pGEmatrix_cr(std::size_t n) : B(n) { } - ~pGEmatrix_cr_t() = default; + ~pGEmatrix_cr() = default; template <typename M> std::pair<std::size_t, std::size_t> gaussian_extend_fill_mat(M &fill) @@ -496,21 +475,21 @@ namespace plib }; template<typename B> - struct pLUmatrix_cr_t : public B + struct pLUmatrix_cr : public B { using base = B; using index_type = typename base::index_type; - PCOPYASSIGNMOVE(pLUmatrix_cr_t, default) + PCOPYASSIGNMOVE(pLUmatrix_cr, default) - explicit pLUmatrix_cr_t(std::size_t n) + explicit pLUmatrix_cr(std::size_t n) : B(n) , ilu_rows(n+1) , m_ILUp(0) { } - ~pLUmatrix_cr_t() = default; + ~pLUmatrix_cr() = default; template <typename M> void build(M &fill, std::size_t ilup) @@ -653,4 +632,4 @@ namespace plib } // namespace plib -#endif // MAT_CR_H_ +#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 <algorithm> #include <array> @@ -27,75 +28,76 @@ namespace plib { template<typename VT, typename T> - 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<decltype(v[0])>::type s(std::forward<T>(scalar)); for ( std::size_t i = 0; i < n; i++ ) v[i] = s; } template<typename VT, typename VS> - 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<typename T, typename V1, typename V2> - 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<T, 8>; // NOLINT - PALIGNAS_VECTOROPT() b8 value = {0}; + const std::size_t n(v1.size()); + gsl_Expects(n <= v2.size()); + T ret(plib::constants<T>::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<typename T, typename VT> - T vec_mult2(const std::size_t n, const VT &v) noexcept + T vec_mult2(const VT &v) noexcept { - using b8 = std::array<T, 8>; - PALIGNAS_VECTOROPT() b8 value = {0}; + const std::size_t n(v.size()); + T ret(plib::constants<T>::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<typename T, typename VT> - 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<T>::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<T, 8>; - 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<typename VV, typename T, typename VR> - 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<decltype(v[0])>::type s(std::forward<T>(scalar)); for ( std::size_t i = 0; i < n; i++ ) result[i] = s * v[i]; } template<typename VR, typename VV, typename T> - 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<decltype(v[0])>::type s(std::forward<T>(scalar)); for ( std::size_t i = 0; i < n; i++ ) result[i] += s * v[i]; @@ -109,31 +111,38 @@ namespace plib } template<typename R, typename V> - 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<typename VR, typename V1, typename V2> - 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<typename V, typename T> - 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<decltype(v[0])>::type s(std::forward<T>(scalar)); for ( std::size_t i = 0; i < n; i++ ) v[i] *= s; } template<typename T, typename V> - 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<T>::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<plib::pmatrix_cr_t<FT, SIZE>>; + using mat_type = plib::pGEmatrix_cr<plib::pmatrix_cr<FT, SIZE>>; using base_type = matrix_solver_ext_t<FT, SIZE>; using fptype = typename base_type::fptype; @@ -119,7 +117,7 @@ namespace solver private: - using mat_index_type = typename plib::pmatrix_cr_t<FT, SIZE>::index_type; + using mat_index_type = typename plib::pmatrix_cr<FT, SIZE>::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<FT, SIZE>::index_type; + using mattype = typename plib::pmatrix_cr<FT, SIZE>::index_type; //plib::mat_precondition_none<FT, SIZE> m_ops; plib::mat_precondition_ILU<FT, SIZE> m_ops; |