summaryrefslogtreecommitdiffstatshomepage
path: root/src/lib/netlist/solver/nld_ms_gmres.h
diff options
context:
space:
mode:
Diffstat (limited to 'src/lib/netlist/solver/nld_ms_gmres.h')
-rw-r--r--src/lib/netlist/solver/nld_ms_gmres.h242
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;
}