summaryrefslogtreecommitdiffstatshomepage
path: root/src/lib/netlist/plib
diff options
context:
space:
mode:
author couriersud <couriersud@gmx.org>2019-10-04 22:34:46 +0200
committer couriersud <couriersud@gmx.org>2019-10-04 22:34:46 +0200
commit43cac214a40f1d61d4a04c3727687b1f70b2c13d (patch)
treea3a97dd067534653b7e206518120169e7135da31 /src/lib/netlist/plib
parent7fc40fc9585c1488bb089f3bbe217ec996f5caed (diff)
netlist: code maintenance (nw)
- Fix SUBMODEL - move to strongly typed matrix sort constant - extend maximum matrix size to 512x512 - optionally do parallel processing based on total operations - templatize GMRES solver loops
Diffstat (limited to 'src/lib/netlist/plib')
-rw-r--r--src/lib/netlist/plib/gmres.h135
-rw-r--r--src/lib/netlist/plib/mat_cr.h80
-rw-r--r--src/lib/netlist/plib/pomp.h25
3 files changed, 152 insertions, 88 deletions
diff --git a/src/lib/netlist/plib/gmres.h b/src/lib/netlist/plib/gmres.h
index 79ea73dd44a..7d10e0d11a4 100644
--- a/src/lib/netlist/plib/gmres.h
+++ b/src/lib/netlist/plib/gmres.h
@@ -20,6 +20,18 @@
namespace plib
{
+ template <int k>
+ struct do_khelper
+ {
+ static constexpr bool value = true;
+ };
+
+ template <>
+ struct do_khelper<-1>
+ {
+ static constexpr float value = 0.0;
+ };
+
template <typename FT, int SIZE>
struct mat_precondition_ILU
{
@@ -143,6 +155,72 @@ namespace plib
std::size_t size() const { return (SIZE<=0) ? m_size : static_cast<std::size_t>(SIZE); }
+ template <int k, typename OPS, typename VT>
+ void do_k(OPS &ops, VT &x, FT rho_delta, int &last_k, bool dummy)
+ {
+ //printf("%d\n", k);
+ do_k<k-1, OPS>(ops, x, rho_delta, last_k, do_khelper<k-1>::value);
+
+ if (last_k < RESTART)
+ return;
+
+ const std::size_t kp1 = k + 1;
+ const std::size_t n = size();
+
+ 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]);
+
+ FT rho = std::abs(m_g[kp1]);
+
+ // FIXME ..
+ //itr_used = itr_used + 1;
+
+ if (rho <= rho_delta || k == RESTART-1)
+ {
+ last_k = k;
+ /* 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]);
+ }
+ }
+
+ template <int k, typename OPS, typename VT>
+ void do_k(OPS &ops, VT &x, FT rho_delta, int &last_k, float dummy)
+ {
+ //printf("here\n");
+ }
+
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)
{
@@ -215,7 +293,7 @@ namespace plib
while (itr_used < itr_max)
{
- std::size_t last_k = RESTART;
+ int last_k = RESTART;
float_type rho;
ops.calc_rhs(Ax, x);
@@ -241,64 +319,19 @@ namespace plib
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]);
+ do_k<RESTART-1>(ops, x, rho_delta, last_k, true);
- 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;
- }
- }
+ bool converged(true);
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];
+ converged = false;
}
- 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)
+ if (converged)
break;
}
diff --git a/src/lib/netlist/plib/mat_cr.h b/src/lib/netlist/plib/mat_cr.h
index e1bc28ab621..ff5e1443728 100644
--- a/src/lib/netlist/plib/mat_cr.h
+++ b/src/lib/netlist/plib/mat_cr.h
@@ -221,43 +221,42 @@ namespace plib
template <typename V>
void gaussian_elimination_parallel(V & RHS)
{
- // FIXME: move into solver creation ...
- plib::omp::set_num_threads(4);
- //printf("omp: %ld\n", m_ge_par.size());
+ //printf("omp: %ld %d %d\n", m_ge_par.size(), nz_num, (int)m_ge_par[m_ge_par.size()-2].size());
for (auto l = 0ul; l < m_ge_par.size(); l++)
- plib::omp::for_static(0ul, m_ge_par[l].size(), [this, &RHS, &l] (unsigned ll)
- {
- auto &i = m_ge_par[l][ll];
+ plib::omp::for_static(nz_num, 0ul, m_ge_par[l].size(), [this, &RHS, &l] (unsigned ll)
{
- std::size_t nzbdp = 0;
- std::size_t pi = diag[i];
- const value_type f = 1.0 / A[pi++];
- const std::size_t piie = row_idx[i+1];
-
- while (auto j = nzbd[i][nzbdp++])
+ auto &i = m_ge_par[l][ll];
{
- // proceed to column i
-
- std::size_t pj = row_idx[j];
+ std::size_t nzbdp = 0;
+ std::size_t pi = diag[i];
+ const value_type f = 1.0 / A[pi++];
+ const std::size_t piie = row_idx[i+1];
+ const auto &nz = nzbd[i];
- while (col_idx[pj] < i)
- pj++;
+ while (auto j = nz[nzbdp++])
+ {
+ // proceed to column i
- const value_type f1 = - A[pj++] * f;
+ std::size_t pj = row_idx[j];
- // subtract row i from j
- // fill-in available assumed, i.e. matrix was prepared
- for (std::size_t pii = pi; pii<piie; pii++)
- {
- while (col_idx[pj] < col_idx[pii])
+ while (col_idx[pj] < i)
pj++;
- if (col_idx[pj] == col_idx[pii])
- A[pj++] += A[pii] * f1;
+
+ const value_type f1 = - A[pj++] * f;
+
+ // subtract row i from j
+ // fill-in available assumed, i.e. matrix was prepared
+ for (std::size_t pii = pi; pii<piie; pii++)
+ {
+ while (col_idx[pj] < col_idx[pii])
+ pj++;
+ if (col_idx[pj] == col_idx[pii])
+ A[pj++] += A[pii] * f1;
+ }
+ RHS[j] += f1 * RHS[i];
}
- RHS[j] += f1 * RHS[i];
}
- }
- });
+ });
}
template <typename V1, typename V2>
@@ -303,11 +302,29 @@ namespace plib
/*
* res = A * x
*/
+#if 0
+ parray<value_type, N < 0 ? -N * N : N *N> xi(m_size*m_size);
+ //std::vector<value_type> xi(m_size*m_size);
+ plib::omp::for_static(0, constants<index_type>::zero(), row_idx[m_size], [this, &x, &xi](index_type k)
+ {
+ xi[k] = x[col_idx[k]];
+ });
+#endif
+#if 1
+ //plib::omp::set_num_threads(4);
+ plib::omp::for_static(0, constants<index_type>::zero(), m_size, [this, &res, &x](index_type row)
+ {
+ T tmp = 0.0;
+ const index_type e = row_idx[row+1];
+ for (index_type k = row_idx[row]; k < e; k++)
+ tmp += A[k] * x[col_idx[k]];
+ res[row] = tmp;
+ });
+#else
std::size_t row = 0;
std::size_t k = 0;
const std::size_t oe = nz_num;
-
while (k < oe)
{
T tmp = 0.0;
@@ -316,6 +333,7 @@ namespace plib
tmp += A[k] * x[col_idx[k]];
res[row++] = tmp;
}
+#endif
}
/* throws error if P(source)>P(destination) */
@@ -409,7 +427,7 @@ namespace plib
{
// we can assume that within a row ja increases continuously */
const std::size_t c_i_j = col_idx[i_j]; // row i, column j
- const std::size_t c_k_j = col_idx[k_j]; // row i, column j
+ const std::size_t c_k_j = col_idx[k_j]; // row k, column j
if (c_k_j < c_i_j)
k_j++;
else if (c_k_j == c_i_j)
@@ -520,6 +538,8 @@ namespace plib
m_ge_par.resize(cl+1);
for (index_type k = 0; k < size(); k++)
m_ge_par[levGE[k]].push_back(k);
+ //for (std::size_t k = 0; k < m_ge_par.size(); k++)
+ // printf("%d %d\n", (int) k, (int) m_ge_par[k].size());
}
index_type m_size;
diff --git a/src/lib/netlist/plib/pomp.h b/src/lib/netlist/plib/pomp.h
index ecd9eacc41b..400fbee9e64 100644
--- a/src/lib/netlist/plib/pomp.h
+++ b/src/lib/netlist/plib/pomp.h
@@ -22,18 +22,29 @@ namespace plib {
namespace omp {
template <typename I, class T>
-void for_static(const I start, const I end, const T &what)
+void for_static(std::size_t numops, const I start, const I end, const T &what)
{
-#if HAS_OPENMP && USE_OPENMP
- #pragma omp parallel
-#endif
+ if (numops>1000)
{
-#if HAS_OPENMP && USE_OPENMP
- #pragma omp for //schedule(static)
-#endif
+ #if HAS_OPENMP && USE_OPENMP
+ #pragma omp parallel for schedule(static)
+ #endif
for (I i = start; i < end; i++)
what(i);
}
+ else
+ for (I i = start; i < end; i++)
+ what(i);
+}
+
+template <typename I, class T>
+void for_static(const I start, const I end, const T &what)
+{
+#if HAS_OPENMP && USE_OPENMP
+ #pragma omp parallel for schedule(static)
+#endif
+ for (I i = start; i < end; i++)
+ what(i);
}
template <typename I, class T>