diff options
author | 2019-10-04 22:34:46 +0200 | |
---|---|---|
committer | 2019-10-04 22:34:46 +0200 | |
commit | 43cac214a40f1d61d4a04c3727687b1f70b2c13d (patch) | |
tree | a3a97dd067534653b7e206518120169e7135da31 /src/lib/netlist/plib | |
parent | 7fc40fc9585c1488bb089f3bbe217ec996f5caed (diff) |
netlist: code maintenance (nw)
- Fix SUBMODEL
- move to strongly typed matrix sort constant
- extend maximum matrix size to 512x512
- optionally do parallel processing based on total operations
- templatize GMRES solver loops
Diffstat (limited to 'src/lib/netlist/plib')
-rw-r--r-- | src/lib/netlist/plib/gmres.h | 135 | ||||
-rw-r--r-- | src/lib/netlist/plib/mat_cr.h | 80 | ||||
-rw-r--r-- | src/lib/netlist/plib/pomp.h | 25 |
3 files changed, 152 insertions, 88 deletions
diff --git a/src/lib/netlist/plib/gmres.h b/src/lib/netlist/plib/gmres.h index 79ea73dd44a..7d10e0d11a4 100644 --- a/src/lib/netlist/plib/gmres.h +++ b/src/lib/netlist/plib/gmres.h @@ -20,6 +20,18 @@ namespace plib { + template <int k> + struct do_khelper + { + static constexpr bool value = true; + }; + + template <> + struct do_khelper<-1> + { + static constexpr float value = 0.0; + }; + template <typename FT, int SIZE> struct mat_precondition_ILU { @@ -143,6 +155,72 @@ namespace plib std::size_t size() const { return (SIZE<=0) ? m_size : static_cast<std::size_t>(SIZE); } + template <int k, typename OPS, typename VT> + void do_k(OPS &ops, VT &x, FT rho_delta, int &last_k, bool dummy) + { + //printf("%d\n", k); + do_k<k-1, OPS>(ops, x, rho_delta, last_k, do_khelper<k-1>::value); + + if (last_k < RESTART) + return; + + const std::size_t kp1 = k + 1; + const std::size_t n = size(); + + ops.calc_rhs(m_v[kp1], m_v[k]); + ops.solve_LU_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] = std::sqrt(vec_mult2<FT>(n, m_v[kp1])); + + if (m_ht[kp1][k] != 0.0) + vec_scale(n, m_v[kp1], constants<FT>::one() / 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 = 1.0 / std::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] = 0.0; + + givens_mult(m_c[k], m_s[k], m_g[k], m_g[kp1]); + + FT rho = std::abs(m_g[kp1]); + + // FIXME .. + //itr_used = itr_used + 1; + + if (rho <= rho_delta || k == RESTART-1) + { + last_k = k; + /* Solve the system H * y = g */ + /* x += m_v[j] * m_y[j] */ + for (std::size_t i = k + 1; i-- > 0;) + { + double 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]); + } + } + + template <int k, typename OPS, typename VT> + void do_k(OPS &ops, VT &x, FT rho_delta, int &last_k, float dummy) + { + //printf("here\n"); + } + 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) { @@ -215,7 +293,7 @@ namespace plib while (itr_used < itr_max) { - std::size_t last_k = RESTART; + int last_k = RESTART; float_type rho; ops.calc_rhs(Ax, x); @@ -241,64 +319,19 @@ namespace plib vec_mult_scalar(n, m_v[0], residual, constants<FT>::one() / rho); - for (std::size_t k = 0; k < RESTART; k++) - { - const std::size_t kp1 = k + 1; - - ops.calc_rhs(m_v[kp1], m_v[k]); - ops.solve_LU_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] = std::sqrt(vec_mult2<FT>(n, m_v[kp1])); - - if (m_ht[kp1][k] != 0.0) - vec_scale(n, m_v[kp1], constants<FT>::one() / 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 = 1.0 / std::hypot(m_ht[k][k], m_ht[kp1][k]); + do_k<RESTART-1>(ops, x, rho_delta, last_k, true); - 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] = 0.0; - - givens_mult(m_c[k], m_s[k], m_g[k], m_g[kp1]); - - rho = std::abs(m_g[kp1]); - - itr_used = itr_used + 1; - - if (rho <= rho_delta) - { - last_k = k; - break; - } - } + bool converged(true); if (last_k >= RESTART) + { /* didn't converge within accuracy */ last_k = RESTART - 1; - - /* Solve the system H * y = g */ - /* x += m_v[j] * m_y[j] */ - for (std::size_t i = last_k + 1; i-- > 0;) - { - double tmp = m_g[i]; - for (std::size_t j = i + 1; j <= last_k; j++) - tmp -= m_ht[i][j] * m_y[j]; - m_y[i] = tmp / m_ht[i][i]; + converged = false; } - for (std::size_t i = 0; i <= last_k; i++) - vec_add_mult_scalar(n, x, m_v[i], m_y[i]); - if (rho <= rho_delta) + if (converged) break; } diff --git a/src/lib/netlist/plib/mat_cr.h b/src/lib/netlist/plib/mat_cr.h index e1bc28ab621..ff5e1443728 100644 --- a/src/lib/netlist/plib/mat_cr.h +++ b/src/lib/netlist/plib/mat_cr.h @@ -221,43 +221,42 @@ namespace plib template <typename V> void gaussian_elimination_parallel(V & RHS) { - // FIXME: move into solver creation ... - plib::omp::set_num_threads(4); - //printf("omp: %ld\n", m_ge_par.size()); + //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(0ul, m_ge_par[l].size(), [this, &RHS, &l] (unsigned ll) - { - auto &i = m_ge_par[l][ll]; + plib::omp::for_static(nz_num, 0ul, m_ge_par[l].size(), [this, &RHS, &l] (unsigned ll) { - std::size_t nzbdp = 0; - std::size_t pi = diag[i]; - const value_type f = 1.0 / A[pi++]; - const std::size_t piie = row_idx[i+1]; - - while (auto j = nzbd[i][nzbdp++]) + auto &i = m_ge_par[l][ll]; { - // proceed to column i - - std::size_t pj = row_idx[j]; + std::size_t nzbdp = 0; + std::size_t pi = diag[i]; + const value_type f = 1.0 / A[pi++]; + const std::size_t piie = row_idx[i+1]; + const auto &nz = nzbd[i]; - while (col_idx[pj] < i) - pj++; + while (auto j = nz[nzbdp++]) + { + // proceed to column i - const value_type f1 = - A[pj++] * f; + std::size_t pj = row_idx[j]; - // subtract row i from j - // fill-in available assumed, i.e. matrix was prepared - for (std::size_t pii = pi; pii<piie; pii++) - { - while (col_idx[pj] < col_idx[pii]) + while (col_idx[pj] < i) pj++; - if (col_idx[pj] == col_idx[pii]) - A[pj++] += A[pii] * f1; + + const value_type f1 = - 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 (col_idx[pj] < col_idx[pii]) + pj++; + if (col_idx[pj] == col_idx[pii]) + A[pj++] += A[pii] * f1; + } + RHS[j] += f1 * RHS[i]; } - RHS[j] += f1 * RHS[i]; } - } - }); + }); } template <typename V1, typename V2> @@ -303,11 +302,29 @@ namespace plib /* * res = A * x */ +#if 0 + parray<value_type, N < 0 ? -N * N : N *N> xi(m_size*m_size); + //std::vector<value_type> xi(m_size*m_size); + plib::omp::for_static(0, constants<index_type>::zero(), row_idx[m_size], [this, &x, &xi](index_type k) + { + xi[k] = x[col_idx[k]]; + }); +#endif +#if 1 + //plib::omp::set_num_threads(4); + plib::omp::for_static(0, constants<index_type>::zero(), m_size, [this, &res, &x](index_type row) + { + T tmp = 0.0; + const index_type e = row_idx[row+1]; + for (index_type k = row_idx[row]; k < e; k++) + tmp += A[k] * x[col_idx[k]]; + res[row] = tmp; + }); +#else std::size_t row = 0; std::size_t k = 0; const std::size_t oe = nz_num; - while (k < oe) { T tmp = 0.0; @@ -316,6 +333,7 @@ namespace plib tmp += A[k] * x[col_idx[k]]; res[row++] = tmp; } +#endif } /* throws error if P(source)>P(destination) */ @@ -409,7 +427,7 @@ namespace plib { // we can assume that within a row ja increases continuously */ const std::size_t c_i_j = col_idx[i_j]; // row i, column j - const std::size_t c_k_j = col_idx[k_j]; // row i, column j + const std::size_t c_k_j = col_idx[k_j]; // row k, column j if (c_k_j < c_i_j) k_j++; else if (c_k_j == c_i_j) @@ -520,6 +538,8 @@ namespace plib m_ge_par.resize(cl+1); for (index_type k = 0; k < 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()); } index_type m_size; diff --git a/src/lib/netlist/plib/pomp.h b/src/lib/netlist/plib/pomp.h index ecd9eacc41b..400fbee9e64 100644 --- a/src/lib/netlist/plib/pomp.h +++ b/src/lib/netlist/plib/pomp.h @@ -22,18 +22,29 @@ namespace plib { namespace omp { template <typename I, class T> -void for_static(const I start, const I end, const T &what) +void for_static(std::size_t numops, const I start, const I end, const T &what) { -#if HAS_OPENMP && USE_OPENMP - #pragma omp parallel -#endif + if (numops>1000) { -#if HAS_OPENMP && USE_OPENMP - #pragma omp for //schedule(static) -#endif + #if HAS_OPENMP && USE_OPENMP + #pragma omp parallel for schedule(static) + #endif for (I i = start; i < end; i++) what(i); } + else + for (I i = start; i < end; i++) + what(i); +} + +template <typename I, class T> +void for_static(const I start, const I end, const T &what) +{ +#if HAS_OPENMP && USE_OPENMP + #pragma omp parallel for schedule(static) +#endif + for (I i = start; i < end; i++) + what(i); } template <typename I, class T> |