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, 349 insertions, 213 deletions
diff --git a/src/lib/netlist/solver/nld_ms_gcr.h b/src/lib/netlist/solver/nld_ms_gcr.h
index f3e56a342d6..a17ca1ec282 100644
--- a/src/lib/netlist/solver/nld_ms_gcr.h
+++ b/src/lib/netlist/solver/nld_ms_gcr.h
@@ -10,201 +10,390 @@
#ifndef NLD_MS_GCR_H_
#define NLD_MS_GCR_H_
-#include "plib/mat_cr.h"
+#include <algorithm>
+#include "../plib/pdynlib.h"
+#include "mat_cr.h"
#include "nld_ms_direct.h"
#include "nld_solver.h"
-#include "plib/pdynlib.h"
-#include "plib/pstream.h"
-#include "plib/vector_ops.h"
-
-#include <algorithm>
+#include "vector_base.h"
+#include "../plib/pstream.h"
namespace netlist
{
-namespace devices
+ namespace devices
+ {
+template <std::size_t m_N, std::size_t storage_N>
+class matrix_solver_GCR_t: public matrix_solver_t
{
+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()
+ {
+ }
- template <typename FT, int SIZE>
- class matrix_solver_GCR_t: public matrix_solver_t
+ virtual ~matrix_solver_GCR_t() 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()
- {
- }
+ }
+
+ 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;
- void vsetup(analog_net_t::list_t &nets) override;
- unsigned vsolve_non_dynamic(const bool newton_raphson) override;
+ virtual std::pair<pstring, pstring> create_solver_code() override;
- std::pair<pstring, pstring> create_solver_code() override;
+private:
- private:
+ //typedef typename mat_cr_t<storage_N>::type mattype;
+ typedef typename mat_cr_t<storage_N>::index_type mattype;
- using mat_index_type = typename plib::matrix_compressed_rows_t<FT, SIZE>::index_type;
+ void csc_private(plib::putf8_fmt_writer &strm);
- void csc_private(plib::putf8_fmt_writer &strm);
+ using extsolver = void (*)(double * RESTRICT m_A, double * RESTRICT RHS, double * RESTRICT V);
- using extsolver = void (*)(double * m_A, double * RHS, double * V);
+ pstring static_compile_name();
- pstring static_compile_name();
+ const std::size_t m_dim;
+ std::vector<unsigned> m_term_cr[storage_N];
+ mat_cr_t<storage_N> mat;
- const std::size_t m_dim;
- plib::parray<FT, SIZE> RHS;
- plib::parray<FT, SIZE> new_V;
+ //extsolver m_proc;
+ plib::dynproc<void, double * RESTRICT, double * RESTRICT, double * RESTRICT> m_proc;
- mat_type mat;
+};
- //extsolver m_proc;
- plib::dynproc<void, double * , double * , double * > m_proc;
+// ----------------------------------------------------------------------------------------
+// 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)
+{
+ setup_base(nets);
- };
+ mattype nz = 0;
+ const std::size_t iN = this->N();
- // ----------------------------------------------------------------------------------------
- // matrix_solver - GCR
- // ----------------------------------------------------------------------------------------
+ /* build the final matrix */
- // FIXME: namespace or static class member
- template <typename V>
- std::size_t inline get_level(const V &v, std::size_t k)
+ bool touched[storage_N][storage_N] = { { false } };
+ for (std::size_t k = 0; k < iN; k++)
{
- 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");
+ for (auto &j : this->m_terms[k]->m_nz)
+ touched[k][j] = true;
}
- template <typename FT, int SIZE>
- void matrix_solver_GCR_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets)
- {
- setup_base(nets);
+ unsigned fc = 0;
- const std::size_t iN = this->N();
+ unsigned ops = 0;
- /* build the final matrix */
+ 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;
+ }
+ }
+ }
+ }
- std::vector<std::vector<unsigned>> fill(iN);
- std::size_t raw_elements = 0;
+ for (mattype k=0; k<iN; k++)
+ {
+ mat.ia[k] = nz;
- for (std::size_t k = 0; k < iN; k++)
+ for (mattype j=0; j<iN; j++)
{
- fill[k].resize(iN, decltype(mat)::FILL_INFINITY);
- for (auto &j : this->m_terms[k]->m_nz)
+ if (touched[k][j])
{
- fill[k][j] = 0;
- raw_elements++;
+ mat.ja[nz] = j;
+ if (j == k)
+ mat.diag[k] = nz;
+ nz++;
}
+ }
+ 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);
+ }
- auto gr = mat.gaussian_extend_fill_mat(fill);
+ mat.ia[iN] = nz;
+ mat.nz_num = nz;
- /* FIXME: move this to the cr matrix class and use computed
- * parallel ordering once it makes sense.
- */
+ 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
+ }
- std::vector<unsigned> levL(iN, 0);
- std::vector<unsigned> levU(iN, 0);
+}
+#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;
- // parallel scheme for L x = y
- for (std::size_t k = 0; k < iN; k++)
+ if (nzbd.size() > 0)
{
- 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;
+ 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);
+ }
}
+ }
- // parallel scheme for U x = y
- for (std::size_t k = iN; k-- > 0; )
+ //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++)
{
- 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("\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)
+{
+ 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 k = 0; k < iN; k++)
+ for (std::size_t i = 0; i < iN - 1; i++)
+ {
+ const auto &nzbd = this->m_terms[i]->m_nzbd;
+
+ if (nzbd.size() > 0)
{
- unsigned fm = 0;
- pstring ml = "";
- for (std::size_t j = 0; j < iN; j++)
+ 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)
{
- 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];
+ // 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);
}
- 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());
+}
- mat.build_from_fill_mat(fill);
+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];
- for (mat_index_type k=0; k<iN; k++)
+ 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++)
{
- 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]];
+ gtot_t += gt[i];
+ RHS_t += Idr[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));
+ for (std::size_t i = railstart; i < term_count; i++)
+ RHS_t += go[i] * *other_cur_analog[i];
- // FIXME: Move me
+ RHS[k] = RHS_t;
- if (state().lib().isLoaded())
+ // 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++)
{
- 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);
+ 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];
+ }
- 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();
+ RHS[k] = RHS_t;
+ mat.A[mat.diag[k]] += gtot_t;
+ }
+#endif
+ mat.ia[iN] = static_cast<mattype>(mat.nz_num);
- for (std::size_t i = 0; i < mat.nz_num; i++)
- strm("double m_A{1} = m_A[{2}];\n", i, i);
+ /* 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
+ {
for (std::size_t i = 0; i < iN - 1; i++)
{
const auto &nzbd = this->m_terms[i]->m_nzbd;
@@ -212,114 +401,61 @@ namespace devices
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];
- //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)
+ for (std::size_t j : nzbd) // for (std::size_t j = i + 1; j < iN; j++)
{
// proceed to column i
- std::size_t pj = mat.row_idx[j];
+ //__builtin_prefetch(&m_A[mat.diag[j+1]], 1);
+ std::size_t pj = mat.ia[j];
- while (mat.col_idx[pj] < i)
+ while (mat.ja[pj] < i)
pj++;
- //const FT f1 = - m_A[pj++] * f;
- strm("\tconst double f{1}_{2} = -f{3} * m_A{4};\n", i, j, i, pj);
- pj++;
+ const nl_double f1 = - mat.A[pj++] * f;
// subtract row i from j */
for (std::size_t pii = pi; pii<piie; )
{
- while (mat.col_idx[pj] < mat.col_idx[pii])
+ 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++;
+ mat.A[pj++] += mat.A[pii++] * f1;
}
- //RHS[j] += f1 * RHS[i];
- strm("\tRHS[{1}] += f{2}_{3} * RHS[{4}];\n", j, i, j, i);
+ RHS[j] += f1 * RHS[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;)
{
- 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++)
+ //__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("\ttmp{1} += m_A{2} * V[{3}];\n", j, pk, mat.col_idx[pk]);
+ tmp += mat.A[pk] * new_V[mat.ja[pk]];
}
- strm("\tV[{1}] = (RHS[{1}] - tmp{1}) / m_A{4};\n", j, j, j, mat.diag[j]);
+ new_V[j] = (RHS[j] - tmp) / mat.A[jdiag];
}
}
- 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);
- }
+ this->m_stat_calculations++;
- 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;
- }
+ const nl_double 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_ */