summaryrefslogtreecommitdiffstatshomepage
path: root/src/lib/netlist/solver/nld_ms_direct.h
diff options
context:
space:
mode:
Diffstat (limited to 'src/lib/netlist/solver/nld_ms_direct.h')
-rw-r--r--src/lib/netlist/solver/nld_ms_direct.h88
1 files changed, 37 insertions, 51 deletions
diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h
index 582a7996947..f3c31dc3371 100644
--- a/src/lib/netlist/solver/nld_ms_direct.h
+++ b/src/lib/netlist/solver/nld_ms_direct.h
@@ -35,7 +35,15 @@ namespace solver
void reset() override { matrix_solver_t::reset(); }
+ private:
+
+ const std::size_t m_dim;
+ const std::size_t m_pitch;
+
protected:
+ static constexpr const std::size_t SIZEABS = plib::parray<FT, SIZE>::SIZEABS();
+ static constexpr const std::size_t m_pitch_ABS = (((SIZEABS + 0) + 7) / 8) * 8;
+
unsigned vsolve_non_dynamic(const bool newton_raphson) override;
unsigned solve_non_dynamic(const bool newton_raphson);
@@ -46,24 +54,12 @@ namespace solver
template <typename T>
void LE_back_subst(T & x);
- const FT &A(std::size_t r, std::size_t c) const noexcept { return m_A[r][c]; }
- FT &A(std::size_t r, std::size_t c) noexcept { return m_A[r][c]; }
-
- const FT &RHS(std::size_t r) const noexcept { return m_A[r][size()]; }
- FT &RHS(std::size_t r) noexcept { return m_A[r][size()]; }
-
PALIGNAS_VECTOROPT()
plib::parray<FT, SIZE> m_new_V;
-
- private:
- static constexpr const std::size_t SIZEABS = plib::parray<FT, SIZE>::SIZEABS();
- static constexpr const std::size_t m_pitch_ABS = (((SIZEABS + 1) + 7) / 8) * 8;
-
- const std::size_t m_dim;
- const std::size_t m_pitch;
-
PALIGNAS_VECTOROPT()
plib::parray2D<FT, SIZE, m_pitch_ABS> m_A;
+ PALIGNAS_VECTOROPT()
+ plib::parray<FT, SIZE> m_RHS;
};
// ----------------------------------------------------------------------------------------
@@ -79,16 +75,16 @@ namespace solver
for (std::size_t i = 0; i < kN; i++)
{
/* FIXME: Singular matrix? */
- const FT f = 1.0 / A(i,i);
+ const FT f = 1.0 / m_A[i][i];
const auto &nzrd = m_terms[i].m_nzrd;
const auto &nzbd = m_terms[i].m_nzbd;
for (auto &j : nzbd)
{
- const FT f1 = -f * A(j, i);
+ const FT f1 = -f * m_A[j][i];
for (auto &k : nzrd)
- A(j, k) += A(i, k) * f1;
- //RHS(j) += RHS(i) * f1;
+ m_A[j][k] += m_A[i][k] * f1;
+ m_RHS[j] += m_RHS[i] * f1;
}
}
}
@@ -101,7 +97,7 @@ namespace solver
for (std::size_t j = i + 1; j < kN; j++)
{
//if (std::abs(m_A[j][i]) > std::abs(m_A[maxrow][i]))
- if (A(j,i) * A(j,i) > A(maxrow,i) * A(maxrow,i))
+ if (m_A[j][i] * m_A[j][i] > m_A[maxrow][i] * m_A[maxrow][i])
maxrow = j;
}
@@ -109,32 +105,28 @@ namespace solver
{
/* Swap the maxrow and ith row */
for (std::size_t k = 0; k < kN + 1; k++) {
- std::swap(A(i,k), A(maxrow,k));
+ std::swap(m_A[i][k], m_A[maxrow][k]);
}
//std::swap(RHS(i), RHS(maxrow));
}
/* FIXME: Singular matrix? */
- const FT f = 1.0 / A(i,i);
+ const FT f = 1.0 / m_A[i][i];
/* Eliminate column i from row j */
for (std::size_t j = i + 1; j < kN; j++)
{
- const FT f1 = - A(j,i) * f;
+ const FT f1 = - m_A[j][i] * f;
if (f1 != plib::constants<FT>::zero())
{
- const FT * pi = &A(i,i+1);
- FT * pj = &A(j,i+1);
- #if 1
- plib::vec_add_mult_scalar_p(kN-i,pj, pi,f1);
- #else
- vec_add_mult_scalar_p1(kN-i-1,pj,pi,f1);
+ const FT * pi = &(m_A[i][i+1]);
+ FT * pj = &(m_A[j][i+1]);
+ plib::vec_add_mult_scalar_p(kN-i-1,pj,pi,f1);
//for (unsigned k = i+1; k < kN; k++)
// pj[k] = pj[k] + pi[k] * f1;
//for (unsigned k = i+1; k < kN; k++)
//A(j,k) += A(i,k) * f1;
- RHS(j) += RHS(i) * f1;
- #endif
+ m_RHS[j] += m_RHS[i] * f1;
}
}
}
@@ -155,8 +147,8 @@ namespace solver
{
FT tmp = 0;
for (std::size_t k = j+1; k < kN; k++)
- tmp += A(j,k) * x[k];
- x[j] = (RHS(j) - tmp) / A(j,j);
+ tmp += m_A[j][k] * x[k];
+ x[j] = (m_RHS[j] - tmp) / m_A[j][j];
}
}
else
@@ -165,10 +157,10 @@ namespace solver
{
FT tmp = 0;
const auto &nzrd = m_terms[j].m_nzrd;
- const auto e = nzrd.size() - 1; /* exclude RHS element */
+ const auto e = nzrd.size(); // - 1; /* exclude RHS element */
for ( std::size_t k = 0; k < e; k++)
- tmp += A(j, nzrd[k]) * x[nzrd[k]];
- x[j] = (RHS(j) - tmp) / A(j,j);
+ tmp += m_A[j][nzrd[k]] * x[nzrd[k]];
+ x[j] = (m_RHS[j] - tmp) / m_A[j][j];
}
}
}
@@ -187,8 +179,12 @@ namespace solver
template <typename FT, int SIZE>
unsigned matrix_solver_direct_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson)
{
- this->build_LE_A(*this);
- this->build_LE_RHS(*this);
+ const std::size_t iN = this->size();
+
+ /* populate matrix */
+
+ this->clear_square_mat(iN, m_A);
+ this->fill_matrix(iN, m_RHS);
this->m_stat_calculations++;
return this->solve_non_dynamic(newton_raphson);
@@ -200,23 +196,13 @@ namespace solver
const solver_parameters_t *params,
const std::size_t size)
: matrix_solver_t(anetlist, name, nets, params)
- , m_new_V(size)
, m_dim(size)
- , m_pitch(m_pitch_ABS ? m_pitch_ABS : (((m_dim + 1) + 7) / 8) * 8)
+ , m_pitch(m_pitch_ABS ? m_pitch_ABS : (((m_dim + 0) + 7) / 8) * 8)
+ , m_new_V(size)
, m_A(size, m_pitch)
+ , m_RHS(size)
{
- /* add RHS element */
- for (std::size_t k = 0; k < this->size(); k++)
- {
- terms_for_net_t & t = m_terms[k];
-
- if (!plib::container::contains(t.m_nzrd, static_cast<unsigned>(this->size())))
- t.m_nzrd.push_back(static_cast<unsigned>(this->size()));
- }
-
- // FIXME: This shouldn't be necessary ...
- for (std::size_t k = 0; k < this->size(); k++)
- state().save(*this, RHS(k), this->name(), plib::pfmt("RHS.{1}")(k));
+ this->build_mat_ptr(size, m_A);
}
} // namespace solver