diff options
Diffstat (limited to 'src/lib/netlist/plib/gmres.h')
-rw-r--r-- | src/lib/netlist/plib/gmres.h | 499 |
1 files changed, 287 insertions, 212 deletions
diff --git a/src/lib/netlist/plib/gmres.h b/src/lib/netlist/plib/gmres.h index 2c357e97624..6ce96ed2ae6 100644 --- a/src/lib/netlist/plib/gmres.h +++ b/src/lib/netlist/plib/gmres.h @@ -1,36 +1,50 @@ -// license:GPL-2.0+ +// license:BSD-3-Clause // copyright-holders:Couriersud -/* - * gmres.h - * - */ #ifndef PLIB_GMRES_H_ #define PLIB_GMRES_H_ -#include "mat_cr.h" +// Names +// spell-checker: words Burkardt, Saad, Yousef, Givens +// + +/// +/// \file gmres.h +/// + #include "parray.h" #include "pconfig.h" +#include "pmatrix_cr.h" #include "vector_ops.h" #include <algorithm> -#include <cmath> - namespace plib { - template <typename FT, int SIZE> + template <int k> + struct do_k_helper + { + static constexpr bool value = true; + }; + + template <> + struct do_k_helper<-1> + { + static constexpr float value = 0.0; + }; + + template <typename ARENA, typename FT, int SIZE> struct mat_precondition_ILU { - using mat_type = plib::matrix_compressed_rows_t<FT, SIZE>; - - mat_precondition_ILU(std::size_t size, int ilu_scale = 4 - , std::size_t bw = plib::matrix_compressed_rows_t<FT, SIZE>::FILL_INFINITY) - : m_mat(static_cast<typename mat_type::index_type>(size)) - , m_LU(static_cast<typename mat_type::index_type>(size)) - , m_use_iLU_preconditioning(ilu_scale >= 0) - , m_ILU_scale(static_cast<std::size_t>(ilu_scale)) + using mat_type = plib::pmatrix_cr<ARENA, FT, SIZE>; + using matLU_type = plib::pLUmatrix_cr<mat_type>; + + mat_precondition_ILU(ARENA &arena, std::size_t size, std::size_t ilu_scale = 4 + , std::size_t bw = plib::pmatrix_cr<ARENA, FT, SIZE>::FILL_INFINITY) + : m_mat(arena, narrow_cast<typename mat_type::index_type>(size)) + , m_LU(arena, narrow_cast<typename mat_type::index_type>(size)) + , m_ILU_scale(narrow_cast<std::size_t>(ilu_scale)) , m_band_width(bw) { } @@ -39,16 +53,11 @@ namespace plib void build(M &fill) { m_mat.build_from_fill_mat(fill, 0); - if (m_use_iLU_preconditioning) - { - m_LU.gaussian_extend_fill_mat(fill); - m_LU.build_from_fill_mat(fill, m_ILU_scale, m_band_width); // ILU(2) - //m_LU.build_from_fill_mat(fill, 9999, 20); // Band matrix width 20 - } + m_LU.build(fill, m_ILU_scale); } - template<typename R, typename V> + template <typename R, typename V> void calc_rhs(R &rhs, const V &v) { m_mat.mult_vec(rhs, v); @@ -56,41 +65,30 @@ namespace plib void precondition() { - if (m_use_iLU_preconditioning) - { - if (m_ILU_scale < 1) - m_LU.raw_copy_from(m_mat); - else - m_LU.reduction_copy_from(m_mat); - m_LU.incomplete_LU_factorization(); - } + m_LU.incomplete_LU_factorization(m_mat); } - template<typename V> - void solve_LU_inplace(V &v) + template <typename V> + void solve_inplace(V &v) { - if (m_use_iLU_preconditioning) - { - m_LU.solveLUx(v); - } + m_LU.solveLU(v); } PALIGNAS_VECTOROPT() mat_type m_mat; PALIGNAS_VECTOROPT() - mat_type m_LU; - bool m_use_iLU_preconditioning; + matLU_type m_LU; std::size_t m_ILU_scale; std::size_t m_band_width; }; - template <typename FT, int SIZE> - struct mat_precondition_diag + template <typename ARENA, typename FT, int SIZE> + struct mat_precondition_diagonal { - mat_precondition_diag(std::size_t size) - : m_mat(size) - , m_diag(size) - , m_use_iLU_preconditioning(true) + mat_precondition_diagonal(ARENA &arena, std::size_t size, [[maybe_unused]] int dummy = 0) + : m_mat(arena, size) + , m_diagonal(size) + , nz_col(size) { } @@ -98,9 +96,21 @@ namespace plib void build(M &fill) { m_mat.build_from_fill_mat(fill, 0); + for (std::size_t i = 0; i< m_diagonal.size(); i++) + { + for (std::size_t j = 0; j< m_diagonal.size(); j++) + { + std::size_t k=m_mat.row_idx[j]; + while (m_mat.col_idx[k] < i && k < m_mat.row_idx[j+1]) + k++; + if (m_mat.col_idx[k] == i && k < m_mat.row_idx[j+1]) + nz_col[i].push_back(k); + } + nz_col[i].push_back(narrow_cast<std::size_t>(-1)); + } } - template<typename R, typename V> + template <typename R, typename V> void calc_rhs(R &rhs, const V &v) { m_mat.mult_vec(rhs, v); @@ -108,87 +118,148 @@ namespace plib void precondition() { - if (m_use_iLU_preconditioning) + for (std::size_t i = 0; i< m_diagonal.size(); i++) { - for (std::size_t i = 0; i< m_diag.size(); i++) + // ILUT: 265% + FT v(0.0); +#if 0 + // doesn't works, Mame performance drops significantly% + // 136% + for (std::size_t j = m_mat.row_idx[i]; j< m_mat.row_idx[i+1]; j++) + v += m_mat.A[j] * m_mat.A[j]; + m_diagonal[i] = reciprocal(std::sqrt(v)); +#elif 0 + // works halfway, i.e. Mame performance 50% + // 147% - lowest average solution time with 7.094 + for (std::size_t j = m_mat.row_idx[i]; j< m_mat.row_idx[i+1]; j++) + v += m_mat.A[j] * m_mat.A[j]; + m_diagonal[i] = m_mat.A[m_mat.diagonal[i]] / v; +#elif 0 + // works halfway, i.e. Mame performance 50% + // sum over column i + // 344% - lowest average solution time with 3.06 + std::size_t nzcolp = 0; + const auto &nz = nz_col[i]; + std::size_t j; + + while ((j = nz[nzcolp++])!=narrow_cast<std::size_t>(-1)) // NOLINT(bugprone-infinite-loop) { - m_diag[i] = 1.0 / m_mat.A[m_mat.diag[i]]; + v += m_mat.A[j] * m_mat.A[j]; } + m_diagonal[i] = m_mat.A[m_mat.diagonal[i]] / v; +#elif 0 + // works halfway, i.e. Mame performance 50% + // 151% + for (std::size_t j = m_mat.row_idx[i]; j< m_mat.row_idx[i+1]; j++) + v += plib::abs(m_mat.A[j]); + m_diagonal[i] = reciprocal(v); +#else + // 124% + for (std::size_t j = m_mat.row_idx[i]; j< m_mat.row_idx[i+1]; j++) + v = std::max(v, plib::abs(m_mat.A[j])); + m_diagonal[i] = reciprocal(v); +#endif + //m_diagonal[i] = reciprocal(m_mat.A[m_mat.diagonal[i]]); } } - template<typename V> - void solve_LU_inplace(V &v) + template <typename V> + void solve_inplace(V &v) { - if (m_use_iLU_preconditioning) - { - for (std::size_t i = 0; i< m_diag.size(); i++) - v[i] = v[i] * m_diag[i]; - } + for (std::size_t i = 0; i< m_diagonal.size(); i++) + v[i] = v[i] * m_diagonal[i]; } - plib::matrix_compressed_rows_t<FT, SIZE> m_mat; - plib::parray<FT, SIZE> m_diag; - bool m_use_iLU_preconditioning; + plib::pmatrix_cr<ARENA, FT, SIZE> m_mat; + plib::parray<FT, SIZE> m_diagonal; + plib::parray<std::vector<std::size_t>, SIZE > nz_col; }; - /* FIXME: hardcoding RESTART to 20 becomes an issue on very large - * systems. - */ - template <typename FT, int SIZE, int RESTART = 20> + template <typename ARENA, typename FT, int SIZE> + struct mat_precondition_none + { + mat_precondition_none(std::size_t size, [[maybe_unused]] int dummy = 0) + : m_mat(size) + { + } + + template <typename M> + void build(M &fill) + { + m_mat.build_from_fill_mat(fill, 0); + } + + template <typename R, typename V> + void calc_rhs(R &rhs, const V &v) + { + m_mat.mult_vec(rhs, v); + } + + void precondition() + { + } + + template <typename V> + void solve_inplace([[maybe_unused]] V &v) + { + } + + plib::pmatrix_cr<ARENA, FT, SIZE> m_mat; + }; + + // FIXME: hard coding RESTART to 20 becomes an issue on very large + // systems. + + template <typename FT, int SIZE, int RESTARTMAX = 16> struct gmres_t { public: using float_type = FT; - // FIXME: dirty hack to make this compile - static constexpr const std::size_t storage_N = plib::sizeabs<FT, SIZE>::ABS(); - gmres_t(std::size_t size) + //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) + , m_ht(RESTART +1, RESTART) + , m_v(RESTART + 1, size) , m_size(size) , m_use_more_precise_stop_condition(false) { } - 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 : static_cast<std::size_t>(SIZE); } + std::size_t size() const { return (SIZE<=0) ? m_size : narrow_cast<std::size_t>(SIZE); } 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) { - /*------------------------------------------------------------------------- - * The code below was inspired by code published by John Burkardt under - * the LPGL here: - * - * http://people.sc.fsu.edu/~jburkardt/cpp_src/mgmres/mgmres.html - * - * The code below was completely written from scratch based on the pseudo code - * found here: - * - * http://de.wikipedia.org/wiki/GMRES-Verfahren - * - * The Algorithm itself is described in - * - * Yousef Saad, - * Iterative Methods for Sparse Linear Systems, - * Second Edition, - * SIAM, 20003, - * ISBN: 0898715342, - * LC: QA188.S17. - * - *------------------------------------------------------------------------*/ + // ------------------------------------------------------------------------- + // The code below was inspired by code published by John Burkardt under + // the LPGL here: + // + // http://people.sc.fsu.edu/~jburkardt/cpp_src/mgmres/mgmres.html + // + // The code below was completely written from scratch based on the pseudo code + // found here: + // + // http://de.wikipedia.org/wiki/GMRES-Verfahren + // + // The Algorithm itself is described in + // + // Yousef Saad, + // Iterative Methods for Sparse Linear Systems, + // Second Edition, + // SIAM, 20003, + // ISBN: 0898715342, + // LC: QA188.S17. + // + //------------------------------------------------------------------------ std::size_t itr_used = 0; - double rho_delta = 0.0; + float_type rho_delta(plib::constants<float_type>::zero()); const std::size_t n = size(); @@ -196,165 +267,174 @@ namespace plib if (m_use_more_precise_stop_condition) { - /* derive residual for a given delta x - * - * LU y = A dx - * - * ==> rho / accuracy = sqrt(y * y) - * - * This approach will approximate the iterative stop condition - * based |xnew - xold| pretty precisely. But it is slow, or expressed - * differently: The invest doesn't pay off. - */ - - vec_set_scalar(n, residual, accuracy); + // derive residual for a given delta x + // + // LU y = A dx + // + // ==> rho / accuracy = sqrt(y * y) + // + // This approach will approximate the iterative stop condition + // based `|xnew - xold|` pretty precisely. But it is slow, or expressed + // differently: The invest doesn't pay off. + // + + vec_set_scalar(residual, accuracy); ops.calc_rhs(Ax, residual); - ops.solve_LU_inplace(Ax); + ops.solve_inplace(Ax); - const float_type rho_to_accuracy = std::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 * std::sqrt(static_cast<FT>(n)); - - /* - * Using - * - * vec_set(n, x, rhs); - * ops.solve_LU_inplace(x); - * - * to get a starting point for x degrades convergence speed compared - * to using the last solution for x. - * - * LU x = b; solve for x; - * - */ + //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); + // ops.solve_inplace(x); + // + // to get a starting point for x degrades convergence speed compared + // to using the last solution for x. while (itr_used < itr_max) { - std::size_t last_k = RESTART; float_type rho; ops.calc_rhs(Ax, x); - vec_sub(n, residual, rhs, Ax); + vec_sub(residual, rhs, Ax); - ops.solve_LU_inplace(residual); + ops.solve_inplace(residual); - rho = std::sqrt(vec_mult2<FT>(n, residual)); + rho = plib::sqrt(vec_mult2<FT>(residual)); if (rho < rho_delta) return itr_used + 1; - /* FIXME: The "+" is necessary to avoid link issues - * on some systems / compiler versions. Issue reported by - * AJR, no details known yet. - */ - vec_set_scalar(RESTART+1, m_g, +constants<FT>::zero()); + // FIXME: The "+" is necessary to avoid link issues + // on some systems / compiler versions. Issue reported by + // AJR, no details known yet. + + vec_set_scalar(m_g, +constants<FT>::zero()); m_g[0] = rho; - //for (std::size_t i = 0; i < mr + 1; i++) - // vec_set_scalar(mr, m_ht[i], NL_FCONST(0.0)); + vec_mult_scalar(m_v[0], residual, plib::reciprocal(rho)); - vec_mult_scalar(n, m_v[0], residual, constants<FT>::one() / rho); + if (do_k<RESTART-1>(ops, x, itr_used, rho_delta, true)) + // converged + break; + } + return itr_used; + } - for (std::size_t k = 0; k < RESTART; k++) - { - const std::size_t kp1 = k + 1; + private: - ops.calc_rhs(m_v[kp1], m_v[k]); - ops.solve_LU_inplace(m_v[kp1]); + static 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; + } - 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])); + template <int k, typename OPS, typename VT> + bool do_k(OPS &ops, VT &x, std::size_t &itr_used, FT rho_delta, [[maybe_unused]] bool dummy) + { + if (do_k<k-1, OPS>(ops, x, itr_used, rho_delta, do_k_helper<k-1>::value)) + return true; + + constexpr const std::size_t kp1 = k + 1; + //const std::size_t n = size(); - if (m_ht[kp1][k] != 0.0) - vec_scale(n, m_v[kp1], constants<FT>::one() / m_ht[kp1][k]); + ops.calc_rhs(m_v[kp1], m_v[k]); + ops.solve_inplace(m_v[kp1]); - 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]); + 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])); - const float_type mu = 1.0 / std::hypot(m_ht[k][k], m_ht[kp1][k]); + // FIXME: comparison to zero + if (m_ht[kp1][k] != plib::constants<FT>::zero()) + vec_scale(m_v[kp1], reciprocal(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; + 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]); - givens_mult(m_c[k], m_s[k], m_g[k], m_g[kp1]); + const float_type mu = reciprocal(plib::hypot(m_ht[k][k], m_ht[kp1][k])); - rho = std::abs(m_g[kp1]); + 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(); - itr_used = itr_used + 1; + givens_mult(m_c[k], m_s[k], m_g[k], m_g[kp1]); - if (rho <= rho_delta) - { - last_k = k; - break; - } - } + const float_type rho = plib::abs(m_g[kp1]); - if (last_k >= RESTART) - /* didn't converge within accuracy */ - last_k = RESTART - 1; + // FIXME .. + itr_used = itr_used + 1; - /* Solve the system H * y = g */ - /* x += m_v[j] * m_y[j] */ - for (std::size_t i = last_k + 1; i-- > 0;) + 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;) { - double tmp = m_g[i]; - for (std::size_t j = i + 1; j <= last_k; j++) + auto tmp = m_g[i]; + const auto ht_i_i = 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 / m_ht[i][i]; + m_y[i] = tmp * ht_i_i; + vec_add_mult_scalar(x, m_v[i], m_y[i]); } - 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) - break; - + //for (std::size_t i = 0; i <= k; i++) + // vec_add_mult_scalar(n, x, m_v[i], m_y[i]); + return true; } - return itr_used; + return false; } - private: - - //typedef typename plib::mat_cr_t<FT, SIZE>::index_type mattype; + 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; - plib::parray<float_type, RESTART + 1> m_c; /* mr + 1 */ - plib::parray<float_type, RESTART + 1> m_g; /* mr + 1 */ - plib::parray<plib::parray<float_type, RESTART>, RESTART + 1> m_ht; /* (mr + 1), mr */ - plib::parray<float_type, RESTART + 1> m_s; /* mr + 1 */ - plib::parray<float_type, RESTART + 1> m_y; /* mr + 1 */ + plib::parray<float_type, RESTART + 1> m_c; // mr + 1 + plib::parray<float_type, RESTART + 1> m_g; // mr + 1 + plib::parray2D<float_type, RESTART + 1, RESTART> m_ht; // (mr + 1), mr + plib::parray<float_type, RESTART + 1> m_s; // mr + 1 + plib::parray<float_type, RESTART + 1> m_y; // mr + 1 - //plib::parray<float_type, SIZE> m_v[RESTART + 1]; /* mr + 1, n */ - plib::parray<plib::parray<float_type, storage_N>, RESTART + 1> m_v; /* mr + 1, n */ + plib::parray2D<float_type, RESTART + 1, SIZE> m_v; // mr + 1, n std::size_t m_size; bool m_use_more_precise_stop_condition; - - }; #if 0 - /* Example of a Chebyshev iteration solver. This one doesn't work yet, - * it needs to be extended for non-symmetric matrix operation and - * depends on spectral radius estimates - which we don't have. - * - * Left here as another example. - */ + // Example of a Chebyshev iteration solver. This one doesn't work yet, + // it needs to be extended for non-symmetric matrix operation and + // depends on spectral radius estimates - which we don't have. + // + // Left here as another example. template <typename FT, int SIZE> struct ch_t @@ -375,16 +455,11 @@ namespace plib { } - std::size_t size() const { return (SIZE<=0) ? m_size : static_cast<std::size_t>(SIZE); } + std::size_t size() const { return (SIZE<=0) ? m_size : narrow_cast<std::size_t>(SIZE); } template <typename OPS, typename VT, typename VRHS> std::size_t solve(OPS &ops, VT &x0, const VRHS & rhs, const std::size_t iter_max, float_type accuracy) { - /*------------------------------------------------------------------------- - * - * - *------------------------------------------------------------------------*/ - ops.precondition(); const FT lmax = 20.0; @@ -404,13 +479,13 @@ namespace plib ops.calc_rhs(Ax, x); vec_sub(size(), rhs, Ax, residual); - FT rho_delta = accuracy * std::sqrt(static_cast<FT>(size())); + FT rho_delta = accuracy * std::sqrt(narrow_cast<FT>(size())); rho_delta = 1e-9; for (int i = 0; i < iter_max; i++) { - ops.solve_LU_inplace(residual); + ops.solve_inplace(residual); if (i==0) { vec_set(size(), p, residual); @@ -419,7 +494,7 @@ namespace plib else { beta = alpha * ( c / 2.0)*( c / 2.0); - alpha = 1.0 / (d - beta); + alpha = reciprocal(d - beta); for (std::size_t k = 0; k < size(); k++) p[k] = residual[k] + beta * p[k]; } @@ -435,7 +510,7 @@ namespace plib } private: - //typedef typename plib::mat_cr_t<FT, SIZE>::index_type mattype; + //#typedef typename plib::mat_cr_t<FT, SIZE>::index_type mattype; plib::parray<float_type, SIZE> residual; plib::parray<float_type, SIZE> Ax; @@ -447,4 +522,4 @@ namespace plib } // namespace plib -#endif /* PLIB_GMRES_H_ */ +#endif // PLIB_GMRES_H_ |