diff options
Diffstat (limited to 'src/lib/netlist/solver/nld_ms_gmres.h')
-rw-r--r-- | src/lib/netlist/solver/nld_ms_gmres.h | 425 |
1 files changed, 425 insertions, 0 deletions
diff --git a/src/lib/netlist/solver/nld_ms_gmres.h b/src/lib/netlist/solver/nld_ms_gmres.h new file mode 100644 index 00000000000..5983cb1a968 --- /dev/null +++ b/src/lib/netlist/solver/nld_ms_gmres.h @@ -0,0 +1,425 @@ +// license:GPL-2.0+ +// copyright-holders:Couriersud +/* + * nld_ms_sor.h + * + * Generic successive over relaxation solver. + * + * Fow w==1 we will do the classic Gauss-Seidel approach + * + */ + +#ifndef NLD_MS_GMRES_H_ +#define NLD_MS_GMRES_H_ + +#include <algorithm> + +#include "solver/mat_cr.h" +#include "solver/nld_ms_direct.h" +#include "solver/nld_solver.h" +#include "solver/vector_base.h" + +NETLIB_NAMESPACE_DEVICES_START() + +template <unsigned m_N, unsigned _storage_N> +class matrix_solver_GMRES_t: public matrix_solver_direct_t<m_N, _storage_N> +{ +public: + + matrix_solver_GMRES_t(const solver_parameters_t *params, int size) + : matrix_solver_direct_t<m_N, _storage_N>(matrix_solver_t::GAUSS_SEIDEL, params, size) + , m_use_iLU_preconditioning(true) + , m_use_more_precise_stop_condition(false) + , m_accuracy_mult(1.0) + { + unsigned mr=this->N(); /* FIXME: maximum iterations locked in here */ + + for (unsigned i = 0; i < mr + 1; i++) + m_ht[i] = new nl_double[mr]; + + for (unsigned i = 0; i < this->N(); i++) + m_v[i] = new nl_double[_storage_N]; + + } + + virtual ~matrix_solver_GMRES_t() + { + unsigned mr=this->N(); /* FIXME: maximum iterations locked in here */ + + for (unsigned i = 0; i < mr + 1; i++) + delete[] m_ht[i]; + + for (unsigned i = 0; i < this->N(); i++) + delete[] m_v[i]; + } + + virtual void vsetup(analog_net_t::list_t &nets); + ATTR_HOT virtual int vsolve_non_dynamic(const bool newton_raphson); +protected: + ATTR_HOT virtual nl_double vsolve(); + +private: + + int solve_ilu_gmres(nl_double * RESTRICT x, nl_double * RESTRICT rhs, const unsigned restart_max, const unsigned mr, nl_double accuracy); + + plist_t<int> m_term_cr[_storage_N]; + + bool m_use_iLU_preconditioning; + bool m_use_more_precise_stop_condition; + + mat_cr_t<_storage_N> mat; + + nl_double m_A[_storage_N * _storage_N]; + nl_double m_LU[_storage_N * _storage_N]; + + nl_double m_c[_storage_N + 1]; /* mr + 1 */ + nl_double m_g[_storage_N + 1]; /* mr + 1 */ + nl_double * RESTRICT m_ht[_storage_N + 1]; /* (mr + 1), mr */ + nl_double m_s[_storage_N]; /* mr + 1 */ + nl_double * RESTRICT m_v[_storage_N + 1]; /*(mr + 1), n */ + //double m_y[_storage_N]; /* mr + 1 */ + + nl_double m_accuracy_mult; // FXIME: Save state +}; + +// ---------------------------------------------------------------------------------------- +// matrix_solver - GMRES +// ---------------------------------------------------------------------------------------- + +template <unsigned m_N, unsigned _storage_N> +void matrix_solver_GMRES_t<m_N, _storage_N>::vsetup(analog_net_t::list_t &nets) +{ + matrix_solver_direct_t<m_N, _storage_N>::vsetup(nets); + + unsigned nz = 0; + const unsigned iN = this->N(); + + for (unsigned k=0; k<iN; k++) + { + terms_t * RESTRICT row = this->m_terms[k]; + mat.ia[k] = nz; + + for (unsigned j=0; j<row->m_nz.size(); j++) + { + mat.ja[nz] = row->m_nz[j]; + if (row->m_nz[j] == k) + mat.diag[k] = nz; + nz++; + } + + /* build pointers into the compressed row format matrix for each terminal */ + + for (unsigned j=0; j< this->m_terms[k]->m_railstart;j++) + { + for (unsigned i = mat.ia[k]; i<nz; i++) + if (this->m_terms[k]->net_other()[j] == (int) mat.ja[i]) + { + m_term_cr[k].add(i); + break; + } + nl_assert(m_term_cr[k].size() == this->m_terms[k]->m_railstart); + } + } + + mat.ia[iN] = nz; + mat.nz_num = nz; +} + +template <unsigned m_N, unsigned _storage_N> +ATTR_HOT nl_double matrix_solver_GMRES_t<m_N, _storage_N>::vsolve() +{ + this->solve_base(this); + return this->compute_next_timestep(); +} + +template <unsigned m_N, unsigned _storage_N> +ATTR_HOT inline int matrix_solver_GMRES_t<m_N, _storage_N>::vsolve_non_dynamic(const bool newton_raphson) +{ + const unsigned 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 + nl_math::sqrt(1-rho)) + */ + + //nz_num = 0; + ATTR_ALIGN nl_double RHS[_storage_N]; + ATTR_ALIGN nl_double new_V[_storage_N]; + ATTR_ALIGN nl_double l_V[_storage_N]; + + for (unsigned i=0, e=mat.nz_num; i<e; i++) + m_A[i] = 0.0; + + for (unsigned k = 0; k < iN; k++) + { + nl_double gtot_t = 0.0; + nl_double RHS_t = 0.0; + + const unsigned term_count = this->m_terms[k]->count(); + const unsigned 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]->other_curanalog(); + + l_V[k] = new_V[k] = this->m_nets[k]->m_cur_Analog; + for (unsigned i = 0; i < term_count; i++) + { + gtot_t = gtot_t + gt[i]; + RHS_t = RHS_t + Idr[i]; + } + + for (unsigned i = railstart; i < term_count; i++) + RHS_t = RHS_t + go[i] * *other_cur_analog[i]; + + RHS[k] = RHS_t; + + // add diagonal element + m_A[mat.diag[k]] = gtot_t; + + for (unsigned i = 0; i < railstart; i++) + { + const unsigned pi = m_term_cr[k][i]; + m_A[pi] -= go[i]; + } + } + mat.ia[iN] = mat.nz_num; + + const nl_double accuracy = this->m_params.m_accuracy; +#if 1 + int mr = _storage_N; + if (_storage_N > 3 ) + mr = (int) sqrt(iN); + mr = std::min(mr, this->m_params.m_gs_loops); + int iter = 4; + int gsl = solve_ilu_gmres(new_V, RHS, iter, mr, accuracy); + int failed = mr * iter; +#else + int failed = 6; + //int gsl = tt_ilu_cr(new_V, RHS, failed, accuracy); + int gsl = tt_gs_cr(new_V, RHS, failed, accuracy); +#endif + this->m_iterative_total += gsl; + this->m_stat_calculations++; + + if (gsl>=failed) + { + //for (int k = 0; k < iN; k++) + // this->m_nets[k]->m_cur_Analog = new_V[k]; + // Fallback to direct solver ... + this->m_iterative_fail++; + return matrix_solver_direct_t<m_N, _storage_N>::vsolve_non_dynamic(newton_raphson); + } + + if (newton_raphson) + { + nl_double err = 0; + for (unsigned k = 0; k < iN; k++) + err = std::max(nl_math::abs(l_V[k] - new_V[k]), err); + + for (unsigned k = 0; k < iN; k++) + this->m_nets[k]->m_cur_Analog += 1.0 * (new_V[k] - this->m_nets[k]->m_cur_Analog); + if (err > accuracy) + return 2; + else + return 1; + } + else + { + for (unsigned k = 0; k < iN; k++) + this->m_nets[k]->m_cur_Analog = new_V[k]; + return 1; + } +} + +static inline void givens_mult( const nl_double c, const nl_double s, nl_double * RESTRICT g0, nl_double * RESTRICT g1 ) +{ + const double tg0 = c * *g0 - s * *g1; + const double tg1 = s * *g0 + c * *g1; + + *g0 = tg0; + *g1 = tg1; +} + +template <unsigned m_N, unsigned _storage_N> +int matrix_solver_GMRES_t<m_N, _storage_N>::solve_ilu_gmres (nl_double * RESTRICT x, nl_double * RESTRICT rhs, const unsigned restart_max, const unsigned mr, nl_double 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. + * + *------------------------------------------------------------------------*/ + + unsigned itr_used = 0; + double rho_delta = 0.0; + + const unsigned n = this->N(); + + if (m_use_iLU_preconditioning) + mat.incomplete_LU_factorization(m_A, m_LU); + + 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. + * 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(m_A, t, Ax); + mat.solveLUx(m_LU, Ax); + + const nl_double rho_to_accuracy = std::sqrt(vecmult2(n, Ax)) / accuracy; + + rho_delta = accuracy * rho_to_accuracy; + } + else + rho_delta = accuracy * std::sqrt((double) n) * m_accuracy_mult; + + for (unsigned itr = 0; itr < restart_max; itr++) + { + unsigned last_k = mr; + nl_double mu; + nl_double rho; + + nl_double Ax[_storage_N]; + nl_double residual[_storage_N]; + + mat.mult_vec(m_A, x, Ax); + + vec_sub(n, rhs, Ax, residual); + + if (m_use_iLU_preconditioning) + { + mat.solveLUx(m_LU, residual); + } + + rho = std::sqrt(vecmult2(n, residual)); + + vec_mult_scalar(n, residual, NL_FCONST(1.0) / rho, m_v[0]); + + vec_set(mr+1, NL_FCONST(0.0), m_g); + m_g[0] = rho; + + for (unsigned i = 0; i < mr; i++) + vec_set(mr + 1, NL_FCONST(0.0), m_ht[i]); + + for (unsigned k = 0; k < mr; k++) + { + const unsigned k1 = k + 1; + + mat.mult_vec(m_A, m_v[k], m_v[k1]); + + if (m_use_iLU_preconditioning) + mat.solveLUx(m_LU, m_v[k1]); + + for (unsigned j = 0; j <= k; j++) + { + m_ht[j][k] = vecmult(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(vecmult2(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]); + + for (unsigned j = 0; j < k; j++) + givens_mult(m_c[j], m_s[j], &m_ht[j][k], &m_ht[j+1][k]); + + mu = std::sqrt(std::pow(m_ht[k][k], 2) + std::pow(m_ht[k1][k], 2)); + + m_c[k] = m_ht[k][k] / mu; + m_s[k] = -m_ht[k1][k] / mu; + m_ht[k][k] = m_c[k] * m_ht[k][k] - m_s[k] * m_ht[k1][k]; + m_ht[k1][k] = 0.0; + + givens_mult(m_c[k], m_s[k], &m_g[k], &m_g[k1]); + + rho = std::abs(m_g[k1]); + + itr_used = itr_used + 1; + + if (rho <= rho_delta) + { + last_k = k; + break; + } + } + + if (last_k >= mr) + /* didn't converge within accuracy */ + last_k = mr - 1; + + nl_double m_y[_storage_N + 1]; + + /* Solve the system H * y = g */ + /* x += m_v[j] * m_y[j] */ + for (int i = last_k; i >= 0; i--) + { + double tmp = m_g[i]; + for (unsigned j = i + 1; j <= last_k; j++) + { + tmp -= m_ht[i][j] * m_y[j]; + } + m_y[i] = tmp / m_ht[i][i]; + } + + for (unsigned 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; +} + + + +NETLIB_NAMESPACE_DEVICES_END() + +#endif /* NLD_MS_GMRES_H_ */ |