diff options
Diffstat (limited to 'src/lib/netlist/plib/pmatrix_cr.h')
-rw-r--r-- | src/lib/netlist/plib/pmatrix_cr.h | 161 |
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; }; |