diff options
Diffstat (limited to 'src/lib/netlist/solver/nld_ms_gmres.h')
-rw-r--r-- | src/lib/netlist/solver/nld_ms_gmres.h | 242 |
1 files changed, 109 insertions, 133 deletions
diff --git a/src/lib/netlist/solver/nld_ms_gmres.h b/src/lib/netlist/solver/nld_ms_gmres.h index 81a7313c154..0b08d995fd3 100644 --- a/src/lib/netlist/solver/nld_ms_gmres.h +++ b/src/lib/netlist/solver/nld_ms_gmres.h @@ -15,6 +15,7 @@ #include <algorithm> #include <cmath> +#include "../plib/parray.h" #include "mat_cr.h" #include "nld_ms_direct.h" #include "nld_solver.h" @@ -24,17 +25,34 @@ namespace netlist { namespace devices { -template <std::size_t m_N, std::size_t storage_N> -class matrix_solver_GMRES_t: public matrix_solver_direct_t<m_N, storage_N> + +template <typename FT, int SIZE> +class matrix_solver_GMRES_t: public matrix_solver_direct_t<FT, SIZE> { 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); + + /* Sort rows in ascending order. This should minimize fill-in and thus + * maximize the efficiency of the incomplete LUT. + */ matrix_solver_GMRES_t(netlist_base_t &anetlist, const pstring &name, const solver_parameters_t *params, const std::size_t size) - : matrix_solver_direct_t<m_N, storage_N>(anetlist, name, matrix_solver_t::ASCENDING, params, size) + : matrix_solver_direct_t<FT, SIZE>(anetlist, name, matrix_solver_t::ASCENDING, params, size) , m_use_iLU_preconditioning(true) - , m_use_more_precise_stop_condition(false) + , m_use_more_precise_stop_condition(true) + , m_ILU_scale(0) , m_accuracy_mult(1.0) + , m_term_cr(size) , mat(size) + , residual(size) + , Ax(size) + , m_LU(size) + //, m_v(size) { } @@ -48,26 +66,32 @@ public: private: //typedef typename mat_cr_t<storage_N>::type mattype; - typedef typename mat_cr_t<storage_N>::index_type mattype; + typedef typename plib::mat_cr_t<FT, SIZE>::index_type mattype; - unsigned solve_ilu_gmres(nl_double (& RESTRICT x)[storage_N], const nl_double (& RESTRICT rhs)[storage_N], const unsigned restart_max, std::size_t mr, nl_double accuracy); + template <typename VT, typename VRHS> + std::size_t solve_ilu_gmres(VT &x, const VRHS & rhs, const std::size_t restart_max, std::size_t mr, float_type accuracy); - std::vector<unsigned> m_term_cr[storage_N]; bool m_use_iLU_preconditioning; bool m_use_more_precise_stop_condition; - nl_double m_accuracy_mult; // FXIME: Save state + std::size_t m_ILU_scale; + float_type m_accuracy_mult; // FXIME: Save state + + plib::parray<std::vector<FT *>, SIZE> m_term_cr; + plib::mat_cr_t<float_type, SIZE> mat; + plib::parray<float_type, SIZE> residual; + plib::parray<float_type, SIZE> Ax; - mat_cr_t<storage_N> mat; + plib::mat_cr_t<float_type, SIZE> m_LU; - nl_double m_LU[storage_N * storage_N]; + float_type m_c[restart_N + 1]; /* mr + 1 */ + float_type m_g[restart_N + 1]; /* mr + 1 */ + float_type m_ht[restart_N + 1][restart_N]; /* (mr + 1), mr */ + float_type m_s[restart_N + 1]; /* mr + 1 */ + float_type m_y[restart_N + 1]; /* mr + 1 */ - nl_double m_c[storage_N + 1]; /* mr + 1 */ - nl_double m_g[storage_N + 1]; /* mr + 1 */ - nl_double m_ht[storage_N + 1][storage_N]; /* (mr + 1), mr */ - nl_double m_s[storage_N + 1]; /* mr + 1 */ - nl_double m_v[storage_N + 1][storage_N]; /*(mr + 1), n */ - nl_double m_y[storage_N + 1]; /* mr + 1 */ + //plib::parray<float_type, SIZE> m_v[restart_N + 1]; /* mr + 1, n */ + float_type m_v[restart_N + 1][storage_N]; /* mr + 1, n */ }; @@ -75,106 +99,72 @@ private: // matrix_solver - GMRES // ---------------------------------------------------------------------------------------- -template <std::size_t m_N, std::size_t storage_N> -void matrix_solver_GMRES_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets) +template <typename FT, int SIZE> +void matrix_solver_GMRES_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets) { - matrix_solver_direct_t<m_N, storage_N>::vsetup(nets); + matrix_solver_direct_t<FT, SIZE>::vsetup(nets); - mattype nz = 0; const std::size_t iN = this->N(); + std::vector<std::vector<unsigned>> fill(iN); + for (std::size_t k=0; k<iN; k++) { + fill[k].resize(iN, decltype(mat)::FILL_INFINITY); terms_for_net_t * RESTRICT row = this->m_terms[k].get(); - mat.ia[k] = nz; - for (std::size_t j=0; j<row->m_nz.size(); j++) { - mat.ja[nz] = static_cast<mattype>(row->m_nz[j]); - if (row->m_nz[j] == k) - mat.diag[k] = nz; - nz++; + fill[k][static_cast<mattype>(row->m_nz[j])] = 0; } + } + + mat.build_from_fill_mat(fill, 0); + m_LU.gaussian_extend_fill_mat(fill); + m_LU.build_from_fill_mat(fill, m_ILU_scale); // ILU(2) + //m_LU.build_from_fill_mat(fill, 9999, 20); // Band matrix width 20 - /* build pointers into the compressed row format matrix for each terminal */ + /* build pointers into the compressed row format matrix for each terminal */ - for (unsigned j=0; j< this->m_terms[k]->m_railstart;j++) + for (std::size_t k=0; k<iN; k++) + { + for (std::size_t j=0; j< this->m_terms[k]->m_railstart;j++) { - for (unsigned i = mat.ia[k]; i<nz; i++) - if (this->m_terms[k]->connected_net_idx()[j] == static_cast<int>(mat.ja[i])) + for (std::size_t i = mat.row_idx[k]; i<mat.row_idx[k+1]; i++) + if (this->m_terms[k]->connected_net_idx()[j] == static_cast<int>(mat.col_idx[i])) { - m_term_cr[k].push_back(i); + m_term_cr[k].push_back(&mat.A[i]); break; } } nl_assert(m_term_cr[k].size() == this->m_terms[k]->m_railstart); + m_term_cr[k].push_back(&mat.A[mat.diag[k]]); } - - mat.ia[iN] = nz; - mat.nz_num = nz; } -template <std::size_t m_N, std::size_t storage_N> -unsigned matrix_solver_GMRES_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson) +template <typename FT, int SIZE> +unsigned matrix_solver_GMRES_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson) { const std::size_t iN = this->N(); - /* ideally, we could get an estimate for the spectral radius of - * Inv(D - L) * U - * - * and estimate using - * - * omega = 2.0 / (1.0 + std::sqrt(1-rho)) - */ - - //nz_num = 0; - nl_double RHS[storage_N]; - nl_double new_V[storage_N]; + plib::parray<FT, SIZE> RHS(iN); + //float_type new_V[storage_N]; mat.set_scalar(0.0); + /* populate matrix and V for first estimate */ for (std::size_t k = 0; k < iN; k++) { - nl_double gtot_t = 0.0; - nl_double RHS_t = 0.0; - - const std::size_t term_count = this->m_terms[k]->count(); - const std::size_t railstart = this->m_terms[k]->m_railstart; - const nl_double * const RESTRICT gt = this->m_terms[k]->gt(); - const nl_double * const RESTRICT go = this->m_terms[k]->go(); - const nl_double * const RESTRICT Idr = this->m_terms[k]->Idr(); - const nl_double * const * RESTRICT other_cur_analog = this->m_terms[k]->connected_net_V(); - - for (std::size_t i = 0; i < term_count; i++) - { - gtot_t = gtot_t + gt[i]; - RHS_t = RHS_t + Idr[i]; - } - - for (std::size_t i = railstart; i < term_count; i++) - RHS_t = RHS_t + go[i] * *other_cur_analog[i]; - - RHS[k] = RHS_t; - - // add diagonal element - mat.A[mat.diag[k]] = gtot_t; - - for (std::size_t i = 0; i < railstart; i++) - { - const std::size_t pi = m_term_cr[k][i]; - mat.A[pi] -= go[i]; - } - - new_V[k] = this->m_nets[k]->Q_Analog(); - + this->m_terms[k]->fill_matrix(m_term_cr[k], RHS[k]); + this->m_new_V[k] = this->m_nets[k]->Q_Analog(); } - mat.ia[iN] = static_cast<mattype>(mat.nz_num); - const nl_double accuracy = this->m_params.m_accuracy; + + //mat.row_idx[iN] = static_cast<mattype>(mat.nz_num); + const float_type accuracy = this->m_params.m_accuracy; const std::size_t mr = (iN > 3 ) ? static_cast<std::size_t>(std::sqrt(iN) * 2.0) : iN; - unsigned iter = std::max(1u, this->m_params.m_gs_loops); - unsigned gsl = solve_ilu_gmres(new_V, RHS, iter, mr, accuracy); + std::size_t iter = std::max(1u, this->m_params.m_gs_loops); + std::size_t gsl = solve_ilu_gmres(this->m_new_V, RHS, iter, mr, accuracy); const std::size_t failed = mr * iter; this->m_iterative_total += gsl; @@ -183,26 +173,29 @@ unsigned matrix_solver_GMRES_t<m_N, storage_N>::vsolve_non_dynamic(const bool ne if (gsl >= failed) { this->m_iterative_fail++; - return matrix_solver_direct_t<m_N, storage_N>::vsolve_non_dynamic(newton_raphson); + return matrix_solver_direct_t<FT, SIZE>::vsolve_non_dynamic(newton_raphson); } - const nl_double err = (newton_raphson ? this->delta(new_V) : 0.0); - this->store(new_V); + //if (newton_raphson) + // printf("%e %e\n", this->delta(this->m_new_V), this->m_params.m_accuracy); + + const float_type err = (newton_raphson ? this->delta(this->m_new_V) : 0.0); + this->store(this->m_new_V); return (err > this->m_params.m_accuracy) ? 2 : 1; } template <typename T> inline void givens_mult( const T c, const T s, T & g0, T & g1 ) { - const T tg0 = c * g0 - s * g1; - const T tg1 = s * g0 + c * g1; + const T g0_last(g0); - g0 = tg0; - g1 = tg1; + g0 = c * g0 - s * g1; + g1 = s * g0_last + c * g1; } -template <std::size_t m_N, std::size_t storage_N> -unsigned matrix_solver_GMRES_t<m_N, storage_N>::solve_ilu_gmres (nl_double (& RESTRICT x)[storage_N], const nl_double (& RESTRICT rhs)[storage_N], const unsigned restart_max, std::size_t mr, nl_double accuracy) +template <typename FT, int SIZE> +template <typename VT, typename VRHS> +std::size_t matrix_solver_GMRES_t<FT, SIZE>::solve_ilu_gmres (VT &x, const VRHS &rhs, const std::size_t restart_max, std::size_t mr, float_type accuracy) { /*------------------------------------------------------------------------- * The code below was inspired by code published by John Burkardt under @@ -226,15 +219,22 @@ unsigned matrix_solver_GMRES_t<m_N, storage_N>::solve_ilu_gmres (nl_double (& RE * *------------------------------------------------------------------------*/ - unsigned itr_used = 0; + std::size_t itr_used = 0; double rho_delta = 0.0; const std::size_t n = this->N(); + if (mr > restart_N) mr = restart_N; if (mr > n) mr = n; if (m_use_iLU_preconditioning) - mat.incomplete_LU_factorization(m_LU); + { + if (m_ILU_scale < 1) + m_LU.raw_copy_from(mat); + else + m_LU.reduction_copy_from(mat); + m_LU.incomplete_LU_factorization(); + } if (m_use_more_precise_stop_condition) { @@ -249,27 +249,22 @@ unsigned matrix_solver_GMRES_t<m_N, storage_N>::solve_ilu_gmres (nl_double (& RE * differently: The invest doesn't pay off. * Therefore we use the approach in the else part. */ - nl_double t[storage_N]; - nl_double Ax[storage_N]; - vec_set(n, accuracy, t); - mat.mult_vec(t, Ax); + vec_set_scalar(n, residual, accuracy); + mat.mult_vec(residual, Ax); - mat.solveLUx(m_LU, Ax); + m_LU.solveLUx(Ax); - const nl_double rho_to_accuracy = std::sqrt(vec_mult2(n, Ax)) / accuracy; + 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(n) * m_accuracy_mult; + rho_delta = accuracy * std::sqrt(static_cast<FT>(n)) * m_accuracy_mult; - for (unsigned itr = 0; itr < restart_max; itr++) + for (std::size_t itr = 0; itr < restart_max; itr++) { std::size_t last_k = mr; - nl_double rho; - - nl_double Ax[storage_N]; - nl_double residual[storage_N]; + float_type rho; mat.mult_vec(x, Ax); @@ -277,19 +272,19 @@ unsigned matrix_solver_GMRES_t<m_N, storage_N>::solve_ilu_gmres (nl_double (& RE if (m_use_iLU_preconditioning) { - mat.solveLUx(m_LU, residual); + m_LU.solveLUx(residual); } - rho = std::sqrt(vec_mult2(n, residual)); + rho = std::sqrt(vec_mult2<FT>(n, residual)); if (rho < rho_delta) return itr_used + 1; - vec_set(mr+1, NL_FCONST(0.0), m_g); + vec_set_scalar(mr+1, m_g, NL_FCONST(0.0)); m_g[0] = rho; for (std::size_t i = 0; i < mr + 1; i++) - vec_set(mr, NL_FCONST(0.0), m_ht[i]); + vec_set_scalar(mr, m_ht[i], NL_FCONST(0.0)); vec_mult_scalar(n, residual, NL_FCONST(1.0) / rho, m_v[0]); @@ -300,14 +295,14 @@ unsigned matrix_solver_GMRES_t<m_N, storage_N>::solve_ilu_gmres (nl_double (& RE mat.mult_vec(m_v[k], m_v[k1]); if (m_use_iLU_preconditioning) - mat.solveLUx(m_LU, m_v[k1]); + m_LU.solveLUx(m_v[k1]); for (std::size_t j = 0; j <= k; j++) { - m_ht[j][k] = vec_mult(n, m_v[k1], m_v[j]); + m_ht[j][k] = vec_mult<float_type>(n, m_v[k1], m_v[j]); vec_add_mult_scalar(n, m_v[j], -m_ht[j][k], m_v[k1]); } - m_ht[k1][k] = std::sqrt(vec_mult2(n, m_v[k1])); + m_ht[k1][k] = std::sqrt(vec_mult2<FT>(n, m_v[k1])); if (m_ht[k1][k] != 0.0) vec_scale(n, m_v[k1], NL_FCONST(1.0) / m_ht[k1][k]); @@ -315,7 +310,7 @@ unsigned matrix_solver_GMRES_t<m_N, storage_N>::solve_ilu_gmres (nl_double (& RE 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 nl_double mu = 1.0 / std::hypot(m_ht[k][k], m_ht[k1][k]); + const float_type mu = 1.0 / std::hypot(m_ht[k][k], m_ht[k1][k]); m_c[k] = m_ht[k][k] * mu; m_s[k] = -m_ht[k1][k] * mu; @@ -345,36 +340,17 @@ unsigned matrix_solver_GMRES_t<m_N, storage_N>::solve_ilu_gmres (nl_double (& RE { 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, m_v[i], m_y[i], x); -#if 1 if (rho <= rho_delta) - { break; - } -#else - /* we try to approximate the x difference between to steps using m_v[last_k] */ - double xdelta = m_y[last_k] * vec_maxabs(n, m_v[last_k]); - if (xdelta < accuracy) - { - if (m_accuracy_mult < 16384.0) - m_accuracy_mult = m_accuracy_mult * 2.0; - break; - } - else - m_accuracy_mult = m_accuracy_mult / 2.0; - -#endif } - return itr_used; } |