summaryrefslogtreecommitdiffstatshomepage
diff options
context:
space:
mode:
author couriersud <couriersud@gmx.org>2020-09-20 08:04:22 +0200
committer couriersud <couriersud@gmx.org>2020-09-20 17:34:54 +0200
commit63efcfcb85e2f475b26541d0973b0321c9204598 (patch)
tree63a3862e71012a254bb17a4d8d099a7a6d16b91d
parent9e4377dc4f2f22e0aecc9daaf043eedf398e38e2 (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.lua2
-rw-r--r--src/lib/netlist/plib/gmres.h189
-rw-r--r--src/lib/netlist/plib/pgsl.h4
-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.h73
-rw-r--r--src/lib/netlist/solver/nld_matrix_solver.h4
-rw-r--r--src/lib/netlist/solver/nld_ms_gcr.h8
-rw-r--r--src/lib/netlist/solver/nld_ms_gmres.h4
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;