diff options
Diffstat (limited to 'src/lib/netlist/plib/gmres.h')
-rw-r--r-- | src/lib/netlist/plib/gmres.h | 450 |
1 files changed, 0 insertions, 450 deletions
diff --git a/src/lib/netlist/plib/gmres.h b/src/lib/netlist/plib/gmres.h deleted file mode 100644 index 2c357e97624..00000000000 --- a/src/lib/netlist/plib/gmres.h +++ /dev/null @@ -1,450 +0,0 @@ -// 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 <algorithm> -#include <cmath> - - -namespace plib -{ - - template <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)) - , m_band_width(bw) - { - } - - template <typename M> - 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 - } - } - - - template<typename R, typename V> - void calc_rhs(R &rhs, const V &v) - { - m_mat.mult_vec(rhs, v); - } - - 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(); - } - } - - template<typename V> - void solve_LU_inplace(V &v) - { - if (m_use_iLU_preconditioning) - { - m_LU.solveLUx(v); - } - } - - PALIGNAS_VECTOROPT() - mat_type m_mat; - PALIGNAS_VECTOROPT() - mat_type m_LU; - bool m_use_iLU_preconditioning; - std::size_t m_ILU_scale; - std::size_t m_band_width; - }; - - template <typename FT, int SIZE> - struct mat_precondition_diag - { - mat_precondition_diag(std::size_t size) - : m_mat(size) - , m_diag(size) - , m_use_iLU_preconditioning(true) - { - } - - 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() - { - if (m_use_iLU_preconditioning) - { - for (std::size_t i = 0; i< m_diag.size(); i++) - { - m_diag[i] = 1.0 / m_mat.A[m_mat.diag[i]]; - } - } - } - - template<typename V> - void solve_LU_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]; - } - } - - plib::matrix_compressed_rows_t<FT, SIZE> m_mat; - plib::parray<FT, SIZE> m_diag; - bool m_use_iLU_preconditioning; - }; - - /* FIXME: hardcoding RESTART to 20 becomes an issue on very large - * systems. - */ - template <typename FT, int SIZE, int RESTART = 20> - 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) - : residual(size) - , Ax(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); } - - 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. - * - *------------------------------------------------------------------------*/ - - 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.calc_rhs(Ax, residual); - - ops.solve_LU_inplace(Ax); - - const float_type rho_to_accuracy = std::sqrt(vec_mult2<FT>(n, 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; - * - */ - - 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); - - ops.solve_LU_inplace(residual); - - rho = std::sqrt(vec_mult2<FT>(n, 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()); - 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(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]); - - 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; - } - } - - 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]; - } - - 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; - - } - return itr_used; - } - - private: - - //typedef typename plib::mat_cr_t<FT, SIZE>::index_type mattype; - - 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, SIZE> m_v[RESTART + 1]; /* mr + 1, n */ - plib::parray<plib::parray<float_type, storage_N>, RESTART + 1> 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. - */ - - template <typename FT, int SIZE> - struct ch_t - { - public: - - typedef FT float_type; - // FIXME: dirty hack to make this compile - static constexpr const std::size_t storage_N = plib::sizeabs<FT, SIZE>::ABS(); - - // Maximum iterations before a restart ... - static constexpr const std::size_t restart_N = (storage_N > 0 ? 20 : 0); - - ch_t(std::size_t size) - : residual(size) - , Ax(size) - , m_size(size) - { - } - - std::size_t size() const { return (SIZE<=0) ? m_size : static_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; - const FT lmin = 0.0001; - - const FT d = (lmax+lmin)/2.0; - const FT c = (lmax-lmin)/2.0; - FT alpha = 0; - FT beta = 0; - std::size_t itr_used = 0; - - plib::parray<FT, SIZE> x(size()); - plib::parray<FT, SIZE> p(size()); - - plib::vec_set(size(), x, x0); - - ops.calc_rhs(Ax, x); - vec_sub(size(), rhs, Ax, residual); - - FT rho_delta = accuracy * std::sqrt(static_cast<FT>(size())); - - rho_delta = 1e-9; - - for (int i = 0; i < iter_max; i++) - { - ops.solve_LU_inplace(residual); - if (i==0) - { - vec_set(size(), p, residual); - alpha = 2.0 / d; - } - else - { - beta = alpha * ( c / 2.0)*( c / 2.0); - alpha = 1.0 / (d - beta); - for (std::size_t k = 0; k < size(); k++) - p[k] = residual[k] + beta * p[k]; - } - plib::vec_add_mult_scalar(size(), p, alpha, x); - ops.calc_rhs(Ax, x); - plib::vec_sub(size(), rhs, Ax, residual); - FT rho = std::sqrt(plib::vec_mult2<FT>(size(), residual)); - if (rho < rho_delta) - break; - itr_used++; - } - return itr_used; - } - private: - - //typedef typename plib::mat_cr_t<FT, SIZE>::index_type mattype; - - plib::parray<float_type, SIZE> residual; - plib::parray<float_type, SIZE> Ax; - - std::size_t m_size; - - }; -#endif - -} // namespace plib - -#endif /* PLIB_GMRES_H_ */ |