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.h562
1 files changed, 213 insertions, 349 deletions
diff --git a/src/lib/netlist/solver/nld_ms_gcr.h b/src/lib/netlist/solver/nld_ms_gcr.h
index a17ca1ec282..f3e56a342d6 100644
--- a/src/lib/netlist/solver/nld_ms_gcr.h
+++ b/src/lib/netlist/solver/nld_ms_gcr.h
@@ -10,390 +10,201 @@
#ifndef NLD_MS_GCR_H_
#define NLD_MS_GCR_H_
-#include <algorithm>
+#include "plib/mat_cr.h"
-#include "../plib/pdynlib.h"
-#include "mat_cr.h"
#include "nld_ms_direct.h"
#include "nld_solver.h"
-#include "vector_base.h"
-#include "../plib/pstream.h"
+#include "plib/pdynlib.h"
+#include "plib/pstream.h"
+#include "plib/vector_ops.h"
+
+#include <algorithm>
namespace netlist
{
- namespace devices
- {
-template <std::size_t m_N, std::size_t storage_N>
-class matrix_solver_GCR_t: public matrix_solver_t
+namespace devices
{
-public:
-
- matrix_solver_GCR_t(netlist_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)
- , mat(size)
- , m_proc()
- {
- }
- virtual ~matrix_solver_GCR_t() override
+ template <typename FT, int SIZE>
+ class matrix_solver_GCR_t: public matrix_solver_t
{
- }
-
- constexpr std::size_t N() const { return (m_N == 0) ? m_dim : m_N; }
-
- virtual void vsetup(analog_net_t::list_t &nets) override;
- virtual unsigned vsolve_non_dynamic(const bool newton_raphson) override;
+ public:
+
+ using mat_type = plib::matrix_compressed_rows_t<FT, SIZE>;
+ // FIXME: dirty hack to make this compile
+ static constexpr const std::size_t storage_N = 100;
+
+ matrix_solver_GCR_t(netlist_state_t &anetlist, const pstring &name,
+ const solver_parameters_t *params, const std::size_t size)
+ : matrix_solver_t(anetlist, name, matrix_solver_t::PREFER_IDENTITY_TOP_LEFT, params)
+ , m_dim(size)
+ , RHS(size)
+ , new_V(size)
+ , mat(static_cast<typename mat_type::index_type>(size))
+ , m_proc()
+ {
+ }
- virtual std::pair<pstring, pstring> create_solver_code() override;
+ constexpr std::size_t N() const { return m_dim; }
-private:
+ void vsetup(analog_net_t::list_t &nets) override;
+ unsigned vsolve_non_dynamic(const bool newton_raphson) override;
- //typedef typename mat_cr_t<storage_N>::type mattype;
- typedef typename mat_cr_t<storage_N>::index_type mattype;
+ std::pair<pstring, pstring> create_solver_code() override;
- void csc_private(plib::putf8_fmt_writer &strm);
+ private:
- using extsolver = void (*)(double * RESTRICT m_A, double * RESTRICT RHS, double * RESTRICT V);
+ using mat_index_type = typename plib::matrix_compressed_rows_t<FT, SIZE>::index_type;
- pstring static_compile_name();
+ void csc_private(plib::putf8_fmt_writer &strm);
- const std::size_t m_dim;
- std::vector<unsigned> m_term_cr[storage_N];
- mat_cr_t<storage_N> mat;
+ using extsolver = void (*)(double * m_A, double * RHS, double * V);
- //extsolver m_proc;
- plib::dynproc<void, double * RESTRICT, double * RESTRICT, double * RESTRICT> m_proc;
+ pstring static_compile_name();
-};
+ const std::size_t m_dim;
+ plib::parray<FT, SIZE> RHS;
+ plib::parray<FT, SIZE> new_V;
-// ----------------------------------------------------------------------------------------
-// matrix_solver - GCR
-// ----------------------------------------------------------------------------------------
+ mat_type mat;
-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)
-{
- setup_base(nets);
+ //extsolver m_proc;
+ plib::dynproc<void, double * , double * , double * > m_proc;
- mattype nz = 0;
- const std::size_t iN = this->N();
+ };
- /* build the final matrix */
+ // ----------------------------------------------------------------------------------------
+ // matrix_solver - GCR
+ // ----------------------------------------------------------------------------------------
- bool touched[storage_N][storage_N] = { { false } };
- for (std::size_t k = 0; k < iN; k++)
+ // FIXME: namespace or static class member
+ template <typename V>
+ std::size_t inline get_level(const V &v, std::size_t k)
{
- for (auto &j : this->m_terms[k]->m_nz)
- touched[k][j] = true;
+ for (std::size_t i = 0; i < v.size(); i++)
+ if (plib::container::contains(v[i], k))
+ return i;
+ throw plib::pexception("Error in get_level");
}
- unsigned fc = 0;
+ template <typename FT, int SIZE>
+ void matrix_solver_GCR_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets)
+ {
+ setup_base(nets);
- unsigned ops = 0;
+ const std::size_t iN = this->N();
- for (std::size_t k = 0; k < iN; k++)
- {
- ops++; // 1/A(k,k)
- for (std::size_t row = k + 1; row < iN; row++)
- {
- 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;
- }
- }
- }
- }
+ /* build the final matrix */
+ std::vector<std::vector<unsigned>> fill(iN);
- for (mattype k=0; k<iN; k++)
- {
- mat.ia[k] = nz;
+ std::size_t raw_elements = 0;
- for (mattype j=0; j<iN; j++)
+ for (std::size_t k = 0; k < iN; k++)
{
- if (touched[k][j])
+ fill[k].resize(iN, decltype(mat)::FILL_INFINITY);
+ for (auto &j : this->m_terms[k]->m_nz)
{
- mat.ja[nz] = j;
- if (j == k)
- mat.diag[k] = nz;
- nz++;
+ fill[k][j] = 0;
+ raw_elements++;
}
- }
- 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]))
- {
- m_term_cr[k].push_back(i);
- break;
- }
}
- nl_assert(m_term_cr[k].size() == this->m_terms[k]->m_railstart);
- }
- mat.ia[iN] = nz;
- mat.nz_num = nz;
+ auto gr = mat.gaussian_extend_fill_mat(fill);
- this->log().verbose("Ops: {1} Occupancy ratio: {2}\n", ops,
- static_cast<double>(nz) / static_cast<double>(iN * iN));
-
- // FIXME: Move me
-
- if (netlist().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->netlist().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
- }
+ /* FIXME: move this to the cr matrix class and use computed
+ * parallel ordering once it makes sense.
+ */
-}
-#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;
+ std::vector<unsigned> levL(iN, 0);
+ std::vector<unsigned> levU(iN, 0);
- if (nzbd.size() > 0)
+ // parallel scheme for L x = y
+ for (std::size_t k = 0; k < iN; k++)
{
- 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);
- }
+ unsigned lm=0;
+ for (std::size_t j = 0; j<k; j++)
+ if (fill[k][j] < decltype(mat)::FILL_INFINITY)
+ lm = std::max(lm, levL[j]);
+ levL[k] = 1+lm;
}
- }
- //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++)
+ // parallel scheme for U x = y
+ for (std::size_t k = iN; k-- > 0; )
{
- strm("\ttmp{1} += m_A[{2}] * V[{3}];\n", j, pk, mat.ja[pk]);
+ unsigned lm=0;
+ for (std::size_t j = iN; --j > k; )
+ if (fill[k][j] < decltype(mat)::FILL_INFINITY)
+ lm = std::max(lm, levU[j]);
+ levU[k] = 1+lm;
}
- 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)
-{
- const std::size_t iN = N();
- for (std::size_t i = 0; i < mat.nz_num; i++)
- strm("double m_A{1} = m_A[{2}];\n", i, i);
- for (std::size_t i = 0; i < iN - 1; i++)
- {
- const auto &nzbd = this->m_terms[i]->m_nzbd;
-
- if (nzbd.size() > 0)
+ for (std::size_t k = 0; k < iN; k++)
{
- 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)
+ unsigned fm = 0;
+ pstring ml = "";
+ for (std::size_t j = 0; j < iN; j++)
{
- // 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);
+ ml += fill[k][j] == 0 ? "X" : fill[k][j] < decltype(mat)::FILL_INFINITY ? "+" : ".";
+ 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} {4:4} {5:4} {6:4}", k, ml, levL[k], levU[k], get_level(mat.m_ge_par, k), fm);
}
- }
-
- //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]);
- }
-}
-#endif
-
-template <std::size_t m_N, std::size_t storage_N>
-pstring matrix_solver_GCR_t<m_N, storage_N>::static_compile_name()
-{
- plib::postringstream t;
- plib::putf8_fmt_writer w(t);
- csc_private(w);
- std::hash<pstring> h;
-
- 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()
-{
- plib::postringstream t;
- plib::putf8_fmt_writer strm(t);
- pstring name = static_compile_name();
-
- strm.writeline(plib::pfmt("extern \"C\" void {1}(double * __restrict m_A, double * __restrict RHS, double * __restrict V)\n")(name));
- strm.writeline("{\n");
- csc_private(strm);
- strm.writeline("}\n");
- return std::pair<pstring, pstring>(name, t.str());
-}
-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)
-{
- const std::size_t iN = this->N();
-
- nl_double RHS[storage_N];
- nl_double new_V[storage_N];
+ mat.build_from_fill_mat(fill);
- mat.set_scalar(0.0);
-
- 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++)
+ for (mat_index_type k=0; k<iN; k++)
{
- gtot_t += gt[i];
- RHS_t += Idr[i];
+ std::size_t cnt(0);
+ /* 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]->m_connected_net_idx[j];
+ for (auto i = mat.row_idx[k]; i < mat.row_idx[k+1]; i++)
+ if (other == static_cast<int>(mat.col_idx[i]))
+ {
+ m_mat_ptr[k][j] = &mat.A[i];
+ cnt++;
+ break;
+ }
+ }
+ nl_assert(cnt == this->m_terms[k]->m_railstart);
+ m_mat_ptr[k][this->m_terms[k]->m_railstart] = &mat.A[mat.diag[k]];
}
- for (std::size_t i = railstart; i < term_count; i++)
- RHS_t += go[i] * *other_cur_analog[i];
+ 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));
- RHS[k] = RHS_t;
+ // FIXME: Move me
- // 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++)
+ if (state().lib().isLoaded())
{
- gtot_t += gt[i];
- RHS_t += Idr[i];
+ pstring symname = static_compile_name();
+ 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);
}
- 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);
-
- /* now solve it */
- //if (m_proc != nullptr)
- if (m_proc.resolved())
- {
- //static_solver(m_A, RHS);
- m_proc(&mat.A[0], &RHS[0], &new_V[0]);
- }
- else
+ 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();
+
+ for (std::size_t i = 0; i < mat.nz_num; i++)
+ strm("double m_A{1} = m_A[{2}];\n", i, i);
+
for (std::size_t i = 0; i < iN - 1; i++)
{
const auto &nzbd = this->m_terms[i]->m_nzbd;
@@ -401,61 +212,114 @@ unsigned matrix_solver_GCR_t<m_N, storage_N>::vsolve_non_dynamic(const bool newt
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++)
+ //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.row_idx[i+1];
+
+ //for (auto & j : nzbd)
+ for (std::size_t j : nzbd)
{
// proceed to column i
- //__builtin_prefetch(&m_A[mat.diag[j+1]], 1);
- 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 = - mat.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++;
- mat.A[pj++] += mat.A[pii++] * f1;
+ //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];
+ //RHS[j] += f1 * RHS[i];
+ strm("\tRHS[{1}] += f{2}_{3} * RHS[{4}];\n", j, i, j, i);
}
}
}
- /* backward substitution
- *
- */
-
- /* row n-1 */
- new_V[iN - 1] = RHS[iN - 1] / mat.A[mat.diag[iN - 1]];
+ //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;)
{
- //__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++)
+ strm("\tdouble tmp{1} = 0.0;\n", j);
+ const std::size_t e = mat.row_idx[j+1];
+ for (std::size_t pk = mat.diag[j] + 1; pk < e; pk++)
{
- tmp += mat.A[pk] * new_V[mat.ja[pk]];
+ strm("\ttmp{1} += m_A{2} * V[{3}];\n", j, pk, mat.col_idx[pk]);
}
- new_V[j] = (RHS[j] - tmp) / mat.A[jdiag];
+ strm("\tV[{1}] = (RHS[{1}] - tmp{1}) / m_A{4};\n", j, j, j, mat.diag[j]);
}
}
- this->m_stat_calculations++;
+ template <typename FT, int SIZE>
+ pstring matrix_solver_GCR_t<FT, SIZE>::static_compile_name()
+ {
+ plib::postringstream t;
+ plib::putf8_fmt_writer w(&t);
+ csc_private(w);
+ std::hash<pstring> h;
+
+ return plib::pfmt("nl_gcr_{1:x}_{2}")(h( t.str() ))(mat.nz_num);
+ }
- const nl_double err = (newton_raphson ? delta(new_V) : 0.0);
- store(new_V);
- return (err > this->m_params.m_accuracy) ? 2 : 1;
-}
+ 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);
+ pstring name = static_compile_name();
+
+ strm.writeline(plib::pfmt("extern \"C\" void {1}(double * __restrict m_A, double * __restrict RHS, double * __restrict V)\n")(name));
+ strm.writeline("{\n");
+ csc_private(strm);
+ strm.writeline("}\n");
+ return std::pair<pstring, pstring>(name, t.str());
+ }
+
+ 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();
+
+ mat.set_scalar(0.0);
+
+ /* populate matrix */
+
+ this->fill_matrix(iN, m_mat_ptr, RHS);
+
+ /* now solve it */
+
+ //if (m_proc != nullptr)
+ if (m_proc.resolved())
+ {
+ //static_solver(m_A, RHS);
+ m_proc(&mat.A[0], &RHS[0], &new_V[0]);
+ }
+ else
+ {
+ // mat.gaussian_elimination_parallel(RHS);
+ mat.gaussian_elimination(RHS);
+ /* backward substitution */
+ mat.gaussian_back_substitution(new_V, RHS);
+ }
+
+ this->m_stat_calculations++;
+
+ const FT err = (newton_raphson ? delta(new_V) : 0.0);
+ store(new_V);
+ return (err > this->m_params.m_accuracy) ? 2 : 1;
+ }
- } //namespace devices
+} // namespace devices
} // namespace netlist
#endif /* NLD_MS_GCR_H_ */