summaryrefslogtreecommitdiffstatshomepage
path: root/src/lib/netlist/solver/nld_ms_gcr.h
diff options
context:
space:
mode:
Diffstat (limited to 'src/lib/netlist/solver/nld_ms_gcr.h')
-rw-r--r--src/lib/netlist/solver/nld_ms_gcr.h314
1 files changed, 72 insertions, 242 deletions
diff --git a/src/lib/netlist/solver/nld_ms_gcr.h b/src/lib/netlist/solver/nld_ms_gcr.h
index 0f5c041061f..e708cdfb281 100644
--- a/src/lib/netlist/solver/nld_ms_gcr.h
+++ b/src/lib/netlist/solver/nld_ms_gcr.h
@@ -21,17 +21,23 @@
namespace netlist
{
- namespace devices
- {
-template <std::size_t m_N, std::size_t storage_N>
+namespace devices
+{
+
+template <typename FT, int SIZE>
class matrix_solver_GCR_t: public matrix_solver_t
{
public:
+ // FIXME: dirty hack to make this compile
+ static constexpr const std::size_t storage_N = 100;
+
matrix_solver_GCR_t(netlist_base_t &anetlist, const pstring &name,
const solver_parameters_t *params, const std::size_t size)
: matrix_solver_t(anetlist, name, matrix_solver_t::ASCENDING, params)
, m_dim(size)
+ , RHS(size)
+ , new_V(size)
, mat(size)
, m_proc()
{
@@ -41,7 +47,7 @@ public:
{
}
- constexpr std::size_t N() const { return (m_N == 0) ? m_dim : m_N; }
+ constexpr std::size_t N() const { return m_dim; }
virtual void vsetup(analog_net_t::list_t &nets) override;
virtual unsigned vsolve_non_dynamic(const bool newton_raphson) override;
@@ -51,7 +57,7 @@ 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 mat_index_type;
void csc_private(plib::putf8_fmt_writer &strm);
@@ -60,8 +66,12 @@ private:
pstring static_compile_name();
const std::size_t m_dim;
- std::vector<unsigned> m_term_cr[storage_N];
- mat_cr_t<storage_N> mat;
+ plib::parray<FT, SIZE> RHS;
+ plib::parray<FT, SIZE> new_V;
+
+ std::vector<FT *> m_term_cr[storage_N];
+
+ plib::mat_cr_t<FT, SIZE> mat;
//extsolver m_proc;
plib::dynproc<void, double * RESTRICT, double * RESTRICT, double * RESTRICT> m_proc;
@@ -72,166 +82,89 @@ private:
// matrix_solver - GCR
// ----------------------------------------------------------------------------------------
-template <std::size_t m_N, std::size_t storage_N>
-void matrix_solver_GCR_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets)
+template <typename FT, int SIZE>
+void matrix_solver_GCR_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets)
{
setup_base(nets);
- mattype nz = 0;
const std::size_t iN = this->N();
/* build the final matrix */
- bool touched[storage_N][storage_N] = { { false } };
+ std::vector<std::vector<unsigned>> fill(iN);
+
+ std::size_t raw_elements = 0;
+
for (std::size_t k = 0; k < iN; k++)
{
+ fill[k].resize(iN, decltype(mat)::FILL_INFINITY);
for (auto &j : this->m_terms[k]->m_nz)
- touched[k][j] = true;
- }
+ {
+ fill[k][j] = 0;
+ raw_elements++;
+ }
- unsigned fc = 0;
+ }
- unsigned ops = 0;
+ auto gr = mat.gaussian_extend_fill_mat(fill);
for (std::size_t k = 0; k < iN; k++)
{
- ops++; // 1/A(k,k)
- for (std::size_t row = k + 1; row < iN; row++)
+ unsigned fm = 0;
+ pstring ml = "";
+ for (std::size_t j = 0; j < iN; j++)
{
- if (touched[row][k])
- {
- ops++;
- fc++;
- for (std::size_t col = k + 1; col < iN; col++)
- if (touched[k][col])
- {
- touched[row][col] = true;
- ops += 2;
- }
- }
+ ml += fill[k][j] < decltype(mat)::FILL_INFINITY ? "X" : "_";
+ if (fill[k][j] < decltype(mat)::FILL_INFINITY)
+ if (fill[k][j] > fm)
+ fm = fill[k][j];
}
+ this->log().verbose("{1:4} {2} {3:4}", k, ml, fm);
}
- for (mattype k=0; k<iN; k++)
- {
- mat.ia[k] = nz;
-
- for (mattype j=0; j<iN; j++)
- {
- if (touched[k][j])
- {
- mat.ja[nz] = j;
- if (j == k)
- mat.diag[k] = nz;
- nz++;
- }
- }
+ mat.build_from_fill_mat(fill);
+ for (mat_index_type k=0; k<iN; k++)
+ {
m_term_cr[k].clear();
/* build pointers into the compressed row format matrix for each terminal */
for (std::size_t j=0; j< this->m_terms[k]->m_railstart;j++)
{
int other = this->m_terms[k]->connected_net_idx()[j];
- for (auto i = mat.ia[k]; i < nz; i++)
- if (other == static_cast<int>(mat.ja[i]))
+ for (auto i = mat.row_idx[k]; i < mat.row_idx[k+1]; i++)
+ if (other == 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;
-
- this->log().verbose("Ops: {1} Occupancy ratio: {2}\n", ops,
- static_cast<double>(nz) / static_cast<double>(iN * iN));
+ this->log().verbose("maximum fill: {1}", gr.first);
+ this->log().verbose("Post elimination occupancy ratio: {2} Ops: {1}", gr.second,
+ static_cast<double>(mat.nz_num) / static_cast<double>(iN * iN));
+ this->log().verbose(" Pre elimination occupancy ratio: {2}",
+ static_cast<double>(raw_elements) / static_cast<double>(iN * iN));
// FIXME: Move me
if (state().lib().isLoaded())
{
pstring symname = static_compile_name();
-#if 0
- m_proc = this->netlist().lib().template getsym<extsolver>(symname);
- if (m_proc != nullptr)
- this->log().verbose("External static solver {1} found ...", symname);
- else
- this->log().warning("External static solver {1} not found ...", symname);
-#else
m_proc.load(this->state().lib(), symname);
if (m_proc.resolved())
this->log().warning("External static solver {1} found ...", symname);
else
this->log().warning("External static solver {1} not found ...", symname);
-#endif
}
}
-#if 0
-template <std::size_t m_N, std::size_t storage_N>
-void matrix_solver_GCR_t<m_N, storage_N>::csc_private(plib::putf8_fmt_writer &strm)
-{
- const std::size_t iN = N();
- for (std::size_t i = 0; i < iN - 1; i++)
- {
- const auto &nzbd = this->m_terms[i]->m_nzbd;
-
- if (nzbd.size() > 0)
- {
- std::size_t pi = mat.diag[i];
- //const nl_double f = 1.0 / m_A[pi++];
- strm("const double f{1} = 1.0 / m_A[{2}];\n", i, pi);
- pi++;
- const std::size_t piie = mat.ia[i+1];
-
- //for (auto & j : nzbd)
- for (std::size_t j : nzbd)
- {
- // proceed to column i
- std::size_t pj = mat.ia[j];
-
- while (mat.ja[pj] < i)
- pj++;
-
- //const nl_double f1 = - m_A[pj++] * f;
- strm("\tconst double f{1}_{2} = -f{3} * m_A[{4}];\n", i, j, i, pj);
- pj++;
-
- // subtract row i from j */
- for (std::size_t pii = pi; pii<piie; )
- {
- while (mat.ja[pj] < mat.ja[pii])
- pj++;
- //m_A[pj++] += m_A[pii++] * f1;
- strm("\tm_A[{1}] += m_A[{2}] * f{3}_{4};\n", pj, pii, i, j);
- pj++; pii++;
- }
- //RHS[j] += f1 * RHS[i];
- strm("\tRHS[{1}] += f{2}_{3} * RHS[{4}];\n", j, i, j, i);
- }
- }
- }
-
- //new_V[iN - 1] = RHS[iN - 1] / mat.A[mat.diag[iN - 1]];
- strm("\tV[{1}] = RHS[{2}] / m_A[{3}];\n", iN - 1, iN - 1, mat.diag[iN - 1]);
- for (std::size_t j = iN - 1; j-- > 0;)
- {
- strm("\tdouble tmp{1} = 0.0;\n", j);
- const std::size_t e = mat.ia[j+1];
- for (std::size_t pk = mat.diag[j] + 1; pk < e; pk++)
- {
- strm("\ttmp{1} += m_A[{2}] * V[{3}];\n", j, pk, mat.ja[pk]);
- }
- strm("\tV[{1}] = (RHS[{1}] - tmp{1}) / m_A[{4}];\n", j, j, j, mat.diag[j]);
- }
-}
-#else
-template <std::size_t m_N, std::size_t storage_N>
-void matrix_solver_GCR_t<m_N, storage_N>::csc_private(plib::putf8_fmt_writer &strm)
+template <typename FT, int SIZE>
+void matrix_solver_GCR_t<FT, SIZE>::csc_private(plib::putf8_fmt_writer &strm)
{
const std::size_t iN = N();
@@ -246,28 +179,28 @@ void matrix_solver_GCR_t<m_N, storage_N>::csc_private(plib::putf8_fmt_writer &st
{
std::size_t pi = mat.diag[i];
- //const nl_double f = 1.0 / m_A[pi++];
+ //const FT f = 1.0 / m_A[pi++];
strm("const double f{1} = 1.0 / m_A{2};\n", i, pi);
pi++;
- const std::size_t piie = mat.ia[i+1];
+ const std::size_t piie = mat.row_idx[i+1];
//for (auto & j : nzbd)
for (std::size_t j : nzbd)
{
// proceed to column i
- std::size_t pj = mat.ia[j];
+ std::size_t pj = mat.row_idx[j];
- while (mat.ja[pj] < i)
+ while (mat.col_idx[pj] < i)
pj++;
- //const nl_double f1 = - m_A[pj++] * f;
+ //const FT f1 = - m_A[pj++] * f;
strm("\tconst double f{1}_{2} = -f{3} * m_A{4};\n", i, j, i, pj);
pj++;
// subtract row i from j */
for (std::size_t pii = pi; pii<piie; )
{
- while (mat.ja[pj] < mat.ja[pii])
+ while (mat.col_idx[pj] < mat.col_idx[pii])
pj++;
//m_A[pj++] += m_A[pii++] * f1;
strm("\tm_A{1} += m_A{2} * f{3}_{4};\n", pj, pii, i, j);
@@ -284,18 +217,17 @@ void matrix_solver_GCR_t<m_N, storage_N>::csc_private(plib::putf8_fmt_writer &st
for (std::size_t j = iN - 1; j-- > 0;)
{
strm("\tdouble tmp{1} = 0.0;\n", j);
- const std::size_t e = mat.ia[j+1];
+ const std::size_t e = mat.row_idx[j+1];
for (std::size_t pk = mat.diag[j] + 1; pk < e; pk++)
{
- strm("\ttmp{1} += m_A{2} * V[{3}];\n", j, pk, mat.ja[pk]);
+ strm("\ttmp{1} += m_A{2} * V[{3}];\n", j, pk, mat.col_idx[pk]);
}
strm("\tV[{1}] = (RHS[{1}] - tmp{1}) / m_A{4};\n", j, j, j, mat.diag[j]);
}
}
-#endif
-template <std::size_t m_N, std::size_t storage_N>
-pstring matrix_solver_GCR_t<m_N, storage_N>::static_compile_name()
+template <typename FT, int SIZE>
+pstring matrix_solver_GCR_t<FT, SIZE>::static_compile_name()
{
plib::postringstream t;
plib::putf8_fmt_writer w(&t);
@@ -305,8 +237,8 @@ pstring matrix_solver_GCR_t<m_N, storage_N>::static_compile_name()
return plib::pfmt("nl_gcr_{1:x}_{2}")(h( t.str() ))(mat.nz_num);
}
-template <std::size_t m_N, std::size_t storage_N>
-std::pair<pstring, pstring> matrix_solver_GCR_t<m_N, storage_N>::create_solver_code()
+template <typename FT, int SIZE>
+std::pair<pstring, pstring> matrix_solver_GCR_t<FT, SIZE>::create_solver_code()
{
plib::postringstream t;
plib::putf8_fmt_writer strm(&t);
@@ -320,69 +252,16 @@ std::pair<pstring, pstring> matrix_solver_GCR_t<m_N, storage_N>::create_solver_c
}
-template <std::size_t m_N, std::size_t storage_N>
-unsigned matrix_solver_GCR_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson)
+template <typename FT, int SIZE>
+unsigned matrix_solver_GCR_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson)
{
const std::size_t iN = this->N();
- nl_double RHS[storage_N];
- nl_double new_V[storage_N];
-
mat.set_scalar(0.0);
+ /* populate matrix */
for (std::size_t k = 0; k < iN; k++)
- {
- terms_for_net_t *t = this->m_terms[k].get();
- nl_double gtot_t = 0.0;
- nl_double RHS_t = 0.0;
-
- const std::size_t term_count = t->count();
- const std::size_t railstart = t->m_railstart;
- const nl_double * const RESTRICT gt = t->gt();
- const nl_double * const RESTRICT go = t->go();
- const nl_double * const RESTRICT Idr = t->Idr();
- const nl_double * const * RESTRICT other_cur_analog = t->connected_net_V();
- const unsigned * const RESTRICT tcr = m_term_cr[k].data();
-
-#if 0
- for (std::size_t i = 0; i < term_count; i++)
- {
- gtot_t += gt[i];
- RHS_t += Idr[i];
- }
-
- for (std::size_t i = railstart; i < term_count; i++)
- 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++)
- mat.A[tcr[i]] -= go[i];
- }
-#else
- for (std::size_t i = 0; i < railstart; i++)
- mat.A[tcr[i]] -= go[i];
-
- for (std::size_t i = 0; i < railstart; i++)
- {
- gtot_t += gt[i];
- RHS_t += Idr[i];
- }
-
- for (std::size_t i = railstart; i < term_count; i++)
- {
- RHS_t += (Idr[i] + go[i] * *other_cur_analog[i]);
- gtot_t += gt[i];
- }
-
- RHS[k] = RHS_t;
- mat.A[mat.diag[k]] += gtot_t;
- }
-#endif
- mat.ia[iN] = static_cast<mattype>(mat.nz_num);
+ this->m_terms[k]->fill_matrix(m_term_cr[k], RHS[k]);
/* now solve it */
@@ -394,63 +273,14 @@ unsigned matrix_solver_GCR_t<m_N, storage_N>::vsolve_non_dynamic(const bool newt
}
else
{
- for (std::size_t i = 0; i < iN - 1; i++)
- {
- const auto &nzbd = this->m_terms[i]->m_nzbd;
-
- if (nzbd.size() > 0)
- {
- std::size_t pi = mat.diag[i];
- const nl_double f = 1.0 / mat.A[pi++];
- const std::size_t piie = mat.ia[i+1];
-
- for (std::size_t j : nzbd) // for (std::size_t j = i + 1; j < iN; j++)
- {
- // proceed to column i
- //__builtin_prefetch(&m_A[mat.diag[j+1]], 1);
- std::size_t pj = mat.ia[j];
-
- while (mat.ja[pj] < i)
- pj++;
-
- const nl_double f1 = - mat.A[pj++] * f;
-
- // subtract row i from j */
- for (std::size_t pii = pi; pii<piie; )
- {
- while (mat.ja[pj] < mat.ja[pii])
- pj++;
- mat.A[pj++] += mat.A[pii++] * f1;
- }
- RHS[j] += f1 * RHS[i];
- }
- }
- }
- /* backward substitution
- *
- */
-
- /* row n-1 */
- new_V[iN - 1] = RHS[iN - 1] / mat.A[mat.diag[iN - 1]];
-
- for (std::size_t j = iN - 1; j-- > 0;)
- {
- //__builtin_prefetch(&new_V[j-1], 1);
- //if (j>0)__builtin_prefetch(&m_A[mat.diag[j-1]], 0);
- double tmp = 0;
- auto jdiag = mat.diag[j];
- const std::size_t e = mat.ia[j+1];
- for (std::size_t pk = jdiag + 1; pk < e; pk++)
- {
- tmp += mat.A[pk] * new_V[mat.ja[pk]];
- }
- new_V[j] = (RHS[j] - tmp) / mat.A[jdiag];
- }
+ mat.gaussian_elimination(RHS);
+ /* backward substitution */
+ mat.gaussian_back_substitution(new_V, RHS);
}
this->m_stat_calculations++;
- const nl_double err = (newton_raphson ? delta(new_V) : 0.0);
+ const FT err = (newton_raphson ? delta(new_V) : 0.0);
store(new_V);
return (err > this->m_params.m_accuracy) ? 2 : 1;
}