summaryrefslogtreecommitdiffstats
path: root/src/lib/netlist/plib/pmatrix_cr.h
diff options
context:
space:
mode:
Diffstat (limited to 'src/lib/netlist/plib/pmatrix_cr.h')
-rw-r--r--src/lib/netlist/plib/pmatrix_cr.h161
1 files changed, 85 insertions, 76 deletions
diff --git a/src/lib/netlist/plib/pmatrix_cr.h b/src/lib/netlist/plib/pmatrix_cr.h
index 38457650e5e..05ba82f1b0e 100644
--- a/src/lib/netlist/plib/pmatrix_cr.h
+++ b/src/lib/netlist/plib/pmatrix_cr.h
@@ -35,7 +35,10 @@ 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, default)
+ pmatrix_cr(const pmatrix_cr &) = default;
+ pmatrix_cr &operator=(const pmatrix_cr &) = default;
+ pmatrix_cr(pmatrix_cr &&) noexcept(std::is_nothrow_move_constructible<parray<value_type, NSQ>>::value) = default;
+ pmatrix_cr &operator=(pmatrix_cr &&) noexcept(std::is_nothrow_move_assignable<parray<value_type, NSQ>>::value && std::is_nothrow_move_assignable<pmatrix2d_vrl<index_type>>::value) = default;
enum constants_e
{
@@ -243,10 +246,13 @@ namespace plib
template<typename B>
struct pGEmatrix_cr : public B
{
- using base = B;
- using index_type = typename base::index_type;
+ using base_type = B;
+ using index_type = typename base_type::index_type;
- PCOPYASSIGNMOVE(pGEmatrix_cr, default)
+ pGEmatrix_cr(const pGEmatrix_cr &) = default;
+ pGEmatrix_cr &operator=(const pGEmatrix_cr &) = default;
+ pGEmatrix_cr(pGEmatrix_cr &&) noexcept(std::is_nothrow_move_constructible<base_type>::value) = default;
+ pGEmatrix_cr &operator=(pGEmatrix_cr &&) noexcept(std::is_nothrow_move_assignable<base_type>::value) = default;
explicit pGEmatrix_cr(std::size_t n)
: B(n)
@@ -266,14 +272,14 @@ namespace plib
ops++; // 1/A(k,k)
for (std::size_t row = k + 1; row < fill.size(); row++)
{
- if (fill[row][k] < base::FILL_INFINITY)
+ if (fill[row][k] < base_type::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 < base_type::FILL_INFINITY)
{
if (f > fill_max)
fill_max = f;
@@ -291,37 +297,37 @@ namespace plib
template <typename V>
void gaussian_elimination(V & RHS)
{
- const std::size_t iN = base::size();
+ const std::size_t iN = base_type::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];
+ std::size_t pi = base_type::diag[i];
+ auto f = reciprocal(base_type::A[pi++]);
+ const std::size_t piie = base_type::row_idx[i+1];
- const auto *nz = base::m_nzbd[i];
+ const auto *nz = base_type::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];
+ std::size_t pj = base_type::row_idx[j];
+ std::size_t pje = base_type::row_idx[j+1];
- while (base::col_idx[pj] < i)
+ while (base_type::col_idx[pj] < i)
pj++;
- const typename base::value_type f1 = - base::A[pj++] * f;
+ const typename base_type::value_type f1 = - base_type::A[pj++] * f;
// subtract row i from j
// fill-in available assumed, i.e. matrix was prepared
for (std::size_t pii = pi; pii<piie && pj < pje; pii++)
{
- while (base::col_idx[pj] < base::col_idx[pii])
+ while (base_type::col_idx[pj] < base_type::col_idx[pii])
pj++;
- if (base::col_idx[pj] == base::col_idx[pii])
- base::A[pj++] += base::A[pii] * f1;
+ if (base_type::col_idx[pj] == base_type::col_idx[pii])
+ base_type::A[pj++] += base_type::A[pii] * f1;
}
RHS[j] += f1 * RHS[i];
@@ -342,35 +348,35 @@ namespace plib
{
//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)
+ plib::omp::for_static(base_type::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];
+ std::size_t pi = base_type::diag[i];
+ const auto f = reciprocal(base_type::A[pi++]);
+ const std::size_t piie = base_type::row_idx[i+1];
+ const auto &nz = base_type::nzbd[i];
while (auto j = nz[nzbdp++])
{
// proceed to column i
- std::size_t pj = base::row_idx[j];
+ std::size_t pj = base_type::row_idx[j];
- while (base::col_idx[pj] < i)
+ while (base_type::col_idx[pj] < i)
pj++;
- auto f1 = - base::A[pj++] * f;
+ auto f1 = - base_type::A[pj++] * f;
// subtract row i from j
// fill-in available assumed, i.e. matrix was prepared
for (std::size_t pii = pi; pii<piie; pii++)
{
- while (base::col_idx[pj] < base::col_idx[pii])
+ while (base_type::col_idx[pj] < base_type::col_idx[pii])
pj++;
- if (base::col_idx[pj] == base::col_idx[pii])
- base::A[pj++] += base::A[pii] * f1;
+ if (base_type::col_idx[pj] == base_type::col_idx[pii])
+ base_type::A[pj++] += base_type::A[pii] * f1;
}
RHS[j] += f1 * RHS[i];
}
@@ -381,36 +387,36 @@ namespace plib
template <typename V1, typename V2>
void gaussian_back_substitution(V1 &V, const V2 &RHS)
{
- const std::size_t iN = base::size();
+ const std::size_t iN = base_type::size();
// row n-1
- V[iN - 1] = RHS[iN - 1] / base::A[base::diag[iN - 1]];
+ V[iN - 1] = RHS[iN - 1] / base_type::A[base_type::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];
+ typename base_type::value_type tmp = 0;
+ const auto jdiag = base_type::diag[j];
+ const std::size_t e = base_type::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];
+ tmp += base_type::A[pk] * V[base_type::col_idx[pk]];
+ V[j] = (RHS[j] - tmp) / base_type::A[jdiag];
}
}
template <typename V1>
void gaussian_back_substitution(V1 &V)
{
- const std::size_t iN = base::size();
+ const std::size_t iN = base_type::size();
// row n-1
- V[iN - 1] = V[iN - 1] / base::A[base::diag[iN - 1]];
+ V[iN - 1] = V[iN - 1] / base_type::A[base_type::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];
+ typename base_type::value_type tmp = 0;
+ const auto jdiag = base_type::diag[j];
+ const std::size_t e = base_type::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];
+ tmp += base_type::A[pk] * V[base_type::col_idx[pk]];
+ V[j] = (V[j] - tmp) / base_type::A[jdiag];
}
}
@@ -419,27 +425,27 @@ namespace plib
void build_parallel_gaussian_execution_scheme(const M &fill)
{
// calculate parallel scheme for gaussian elimination
- std::vector<std::vector<std::size_t>> rt(base::size());
- for (std::size_t k = 0; k < base::size(); k++)
+ std::vector<std::vector<std::size_t>> rt(base_type::size());
+ for (std::size_t k = 0; k < base_type::size(); k++)
{
- for (std::size_t j = k+1; j < base::size(); j++)
+ for (std::size_t j = k+1; j < base_type::size(); j++)
{
- if (fill[j][k] < base::FILL_INFINITY)
+ if (fill[j][k] < base_type::FILL_INFINITY)
{
rt[k].push_back(j);
}
}
}
- std::vector<std::size_t> levGE(base::size(), 0);
+ std::vector<std::size_t> levGE(base_type::size(), 0);
std::size_t cl = 0;
- for (std::size_t k = 0; k < base::size(); k++ )
+ for (std::size_t k = 0; k < base_type::size(); k++ )
{
if (levGE[k] >= cl)
{
std::vector<std::size_t> t = rt[k];
- for (std::size_t j = k+1; j < base::size(); j++ )
+ for (std::size_t j = k+1; j < base_type::size(); j++ )
{
bool overlap = false;
// is there overlap
@@ -466,7 +472,7 @@ namespace plib
m_ge_par.clear();
m_ge_par.resize(cl+1);
- for (std::size_t k = 0; k < base::size(); k++)
+ for (std::size_t k = 0; k < base_type::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());
@@ -477,10 +483,13 @@ namespace plib
template<typename B>
struct pLUmatrix_cr : public B
{
- using base = B;
- using index_type = typename base::index_type;
+ using base_type = B;
+ using index_type = typename base_type::index_type;
- PCOPYASSIGNMOVE(pLUmatrix_cr, default)
+ pLUmatrix_cr(const pLUmatrix_cr &) = default;
+ pLUmatrix_cr &operator=(const pLUmatrix_cr &) = default;
+ pLUmatrix_cr(pLUmatrix_cr &&) noexcept(std::is_nothrow_move_constructible<base_type>::value) = default;
+ pLUmatrix_cr &operator=(pLUmatrix_cr &&) noexcept(std::is_nothrow_move_assignable<base_type>::value) = default;
explicit pLUmatrix_cr(std::size_t n)
: B(n)
@@ -539,7 +548,7 @@ namespace plib
/// k=k+1
/// i=i+1
///
- void incomplete_LU_factorization(const base &mat)
+ void incomplete_LU_factorization(const base_type &mat)
{
if (m_ILUp < 1)
this->raw_copy_from(mat);
@@ -549,27 +558,27 @@ namespace plib
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];
+ const auto p_i_end = base_type::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++)
+ for (auto i_k = base_type::row_idx[i]; i_k < base_type::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]];
+ const auto k(base_type::col_idx[i_k]);
+ const auto p_k_end(base_type::row_idx[k + 1]);
+ const typename base_type::value_type LUp_i_k = base_type::A[i_k] = base_type::A[i_k] / base_type::A[base_type::diag[k]];
- std::size_t k_j(base::diag[k] + 1);
+ std::size_t k_j(base_type::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
+ const std::size_t c_i_j(base_type::col_idx[i_j]); // row i, column j
+ const auto c_k_j(base_type::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];
+ base_type::A[i_j] -= LUp_i_k * base_type::A[k_j];
k_j += (c_k_j <= c_i_j ? 1 : 0);
i_j += (c_k_j >= c_i_j ? 1 : 0);
@@ -604,29 +613,29 @@ namespace plib
template <typename R>
void solveLU (R &r)
{
- for (std::size_t i = 1; i < base::size(); ++i )
+ for (std::size_t i = 1; i < base_type::size(); ++i )
{
- typename base::value_type tmp(0);
- const auto j1(base::row_idx[i]);
- const auto j2(base::diag[i]);
+ typename base_type::value_type tmp(0);
+ const auto j1(base_type::row_idx[i]);
+ const auto j2(base_type::diag[i]);
for (auto j = j1; j < j2; ++j )
- tmp += base::A[j] * r[base::col_idx[j]];
+ tmp += base_type::A[j] * r[base_type::col_idx[j]];
r[i] -= tmp;
}
// i now is equal to n;
- for (std::size_t i = base::size(); i-- > 0; )
+ for (std::size_t i = base_type::size(); i-- > 0; )
{
- typename base::value_type tmp(0);
- const auto di(base::diag[i]);
- const auto j2(base::row_idx[i+1]);
+ typename base_type::value_type tmp(0);
+ const auto di(base_type::diag[i]);
+ const auto j2(base_type::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];
+ tmp += base_type::A[j] * r[base_type::col_idx[j]];
+ r[i] = (r[i] - tmp) / base_type::A[di];
}
}
private:
- parray<index_type, base::Np1> ilu_rows;
+ parray<index_type, base_type::Np1> ilu_rows;
std::size_t m_ILUp;
};