// license:GPL-2.0+ // copyright-holders:Couriersud /* * gmres.h * */ #ifndef PLIB_GMRES_H_ #define PLIB_GMRES_H_ #include "mat_cr.h" #include "parray.h" #include "pconfig.h" #include "vector_ops.h" #include #include namespace plib { template struct do_khelper { static constexpr bool value = true; }; template <> struct do_khelper<-1> { static constexpr float value = 0.0; }; template struct mat_precondition_ILU { using mat_type = plib::pmatrix_cr_t; using matLU_type = plib::pLUmatrix_cr_t; mat_precondition_ILU(std::size_t size, std::size_t ilu_scale = 4 , std::size_t bw = plib::pmatrix_cr_t::FILL_INFINITY) : m_mat(static_cast(size)) , m_LU(static_cast(size)) , m_ILU_scale(static_cast(ilu_scale)) , m_band_width(bw) { } template void build(M &fill) { m_mat.build_from_fill_mat(fill, 0); m_LU.build(fill, m_ILU_scale); } template void calc_rhs(R &rhs, const V &v) { m_mat.mult_vec(rhs, v); } void precondition() { m_LU.incomplete_LU_factorization(m_mat); } template void solve_inplace(V &v) { m_LU.solveLU(v); } PALIGNAS_VECTOROPT() mat_type m_mat; PALIGNAS_VECTOROPT() matLU_type m_LU; std::size_t m_ILU_scale; std::size_t m_band_width; }; template struct mat_precondition_diag { mat_precondition_diag(std::size_t size, int dummy = 0) : m_mat(size) , m_diag(size) , nzcol(size) { plib::unused_var(dummy); } template void build(M &fill) { m_mat.build_from_fill_mat(fill, 0); for (std::size_t i = 0; i< m_diag.size(); i++) { for (std::size_t j = 0; j< m_diag.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]) nzcol[i].push_back(k); } nzcol[i].push_back(static_cast(-1)); } } template void calc_rhs(R &rhs, const V &v) { m_mat.mult_vec(rhs, v); } void precondition() { for (std::size_t i = 0; i< m_diag.size(); i++) { // ILUT: 265% FT v(0.0); #if 0 // doesn't works, Mame perforamnce 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_diag[i] = 1.0 / std::sqrt(v); #elif 0 // works halfway, i.e. Mame perforamnce 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_diag[i] = m_mat.A[m_mat.diag[i]] / v; #elif 0 // works halfway, i.e. Mame perforamnce 50% // sum over column i // 344% - lowest average solution time with 3.06 std::size_t nzcolp = 0; const auto &nz = nzcol[i]; std::size_t j; while ((j = nz[nzcolp++])!=static_cast(-1)) // NOLINT(bugprone-infinite-loop) { v += m_mat.A[j] * m_mat.A[j]; } m_diag[i] = m_mat.A[m_mat.diag[i]] / v; #elif 0 // works halfway, i.e. Mame perforamnce 50% // 151% for (std::size_t j = m_mat.row_idx[i]; j< m_mat.row_idx[i+1]; j++) v += std::abs(m_mat.A[j]); m_diag[i] = 1.0 / 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, std::abs(m_mat.A[j])); m_diag[i] = 1.0 / v; #endif //m_diag[i] = 1.0 / m_mat.A[m_mat.diag[i]]; } } template void solve_inplace(V &v) { for (std::size_t i = 0; i< m_diag.size(); i++) v[i] = v[i] * m_diag[i]; } plib::pmatrix_cr_t m_mat; plib::parray m_diag; plib::parray, SIZE > nzcol; }; template struct mat_precondition_none { mat_precondition_none(std::size_t size, int dummy = 0) : m_mat(size) { plib::unused_var(dummy); } template void build(M &fill) { m_mat.build_from_fill_mat(fill, 0); } template void calc_rhs(R &rhs, const V &v) { m_mat.mult_vec(rhs, v); } void precondition() { } template void solve_inplace(V &v) { plib::unused_var(v); } plib::pmatrix_cr_t m_mat; }; /* FIXME: hardcoding RESTART to 20 becomes an issue on very large * systems. */ template 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::ABS(); 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(SIZE); } template bool do_k(OPS &ops, VT &x, std::size_t &itr_used, FT rho_delta, bool dummy) { plib::unused_var(dummy); //printf("%d\n", k); if (do_k(ops, x, itr_used, rho_delta, do_khelper::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(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(n, m_v[kp1])); if (m_ht[kp1][k] != 0.0) vec_scale(n, m_v[kp1], constants::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) { /* 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]); return true; } else return false; } template 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 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. * *------------------------------------------------------------------------*/ std::size_t itr_used = 0; double rho_delta = 0.0; const std::size_t n = size(); ops.precondition(); 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); ops.ca