summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--nl_examples/kidniki.c3
-rw-r--r--src/lib/netlist/solver/nld_matrix_solver.h2
-rw-r--r--src/lib/netlist/solver/nld_ms_direct.h14
-rw-r--r--src/lib/netlist/solver/nld_ms_sm.h91
-rw-r--r--src/lib/netlist/solver/nld_solver.cpp8
5 files changed, 68 insertions, 50 deletions
diff --git a/nl_examples/kidniki.c b/nl_examples/kidniki.c
index f7dbf69bdfc..263e0d3eb7d 100644
--- a/nl_examples/kidniki.c
+++ b/nl_examples/kidniki.c
@@ -14,9 +14,10 @@ NETLIST_START(dummy)
PARAM(Solver.NR_LOOPS, 300)
PARAM(Solver.GS_LOOPS, 1)
PARAM(Solver.GS_THRESHOLD, 6)
+ PARAM(Solver.ITERATIVE, "SM")
//PARAM(Solver.ITERATIVE, "MAT")
//PARAM(Solver.ITERATIVE, "GMRES")
- PARAM(Solver.ITERATIVE, "SOR")
+ //PARAM(Solver.ITERATIVE, "SOR")
PARAM(Solver.DYNAMIC_TS, 0)
PARAM(Solver.DYNAMIC_LTE, 5e-3)
PARAM(Solver.MIN_TIMESTEP, 10e-6)
diff --git a/src/lib/netlist/solver/nld_matrix_solver.h b/src/lib/netlist/solver/nld_matrix_solver.h
index 40896d431e5..715eab6756a 100644
--- a/src/lib/netlist/solver/nld_matrix_solver.h
+++ b/src/lib/netlist/solver/nld_matrix_solver.h
@@ -46,7 +46,7 @@ class terms_t
unsigned m_railstart;
pvector_t<unsigned> m_nz; /* all non zero for multiplication */
- pvector_t<unsigned> m_nzrd; /* non zero right of the diagonal for elimination, includes RHS element */
+ pvector_t<unsigned> m_nzrd; /* non zero right of the diagonal for elimination, may include RHS element */
pvector_t<unsigned> m_nzbd; /* non zero below of the diagonal for elimination */
private:
pvector_t<terminal_t *> m_term;
diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h
index 76b448d8e8f..6911e6aa360 100644
--- a/src/lib/netlist/solver/nld_ms_direct.h
+++ b/src/lib/netlist/solver/nld_ms_direct.h
@@ -269,16 +269,14 @@ ATTR_COLD void matrix_solver_direct_t<m_N, _storage_N>::vsetup(analog_net_t::lis
}
}
- for (unsigned j = 0; j < N(); j++)
+ for (unsigned i = 0; i < t->m_railstart; i++)
{
- for (unsigned i = 0; i < t->m_railstart; i++)
- {
- if (!t->m_nzrd.contains(other[i]) && other[i] >= (int) (k + 1))
- t->m_nzrd.push_back(other[i]);
- if (!t->m_nz.contains(other[i]))
- t->m_nz.push_back(other[i]);
- }
+ if (!t->m_nzrd.contains(other[i]) && other[i] >= (int) (k + 1))
+ t->m_nzrd.push_back(other[i]);
+ if (!t->m_nz.contains(other[i]))
+ t->m_nz.push_back(other[i]);
}
+
/* Add RHS element */
if (!t->m_nzrd.contains(N()))
t->m_nzrd.push_back(N());
diff --git a/src/lib/netlist/solver/nld_ms_sm.h b/src/lib/netlist/solver/nld_ms_sm.h
index c9ee083c54a..81395a11492 100644
--- a/src/lib/netlist/solver/nld_ms_sm.h
+++ b/src/lib/netlist/solver/nld_ms_sm.h
@@ -13,11 +13,6 @@
#include "solver/nld_solver.h"
#include "solver/vector_base.h"
-/* Disabling dynamic allocation gives a ~10% boost in performance
- * This flag has been added to support continuous storage for arrays
- * going forward in case we implement cuda solvers in the future.
- */
-
NETLIB_NAMESPACE_DEVICES_START()
//#define nl_ext_double __float128 // slow, very slow
@@ -273,16 +268,14 @@ ATTR_COLD void matrix_solver_sm_t<m_N, _storage_N>::vsetup(analog_net_t::list_t
}
}
- for (unsigned j = 0; j < N(); j++)
+ for (unsigned i = 0; i < t->m_railstart; i++)
{
- for (unsigned i = 0; i < t->m_railstart; i++)
- {
- if (!t->m_nzrd.contains(other[i]) && other[i] >= (int) (k + 1))
- t->m_nzrd.push_back(other[i]);
- if (!t->m_nz.contains(other[i]))
- t->m_nz.push_back(other[i]);
- }
+ if (!t->m_nzrd.contains(other[i]) && other[i] >= (int) (k + 1))
+ t->m_nzrd.push_back(other[i]);
+ if (!t->m_nz.contains(other[i]))
+ t->m_nz.push_back(other[i]);
}
+
/* and sort */
psort_list(t->m_nzrd);
@@ -476,8 +469,10 @@ void matrix_solver_sm_t<m_N, _storage_N>::LE_compute_x(
for (int k=0; k<kN; k++)
{
+ const nl_double f = RHS(k);
+
for (int i=0; i<kN; i++)
- x[i] += Ainv(i,k) * RHS(k);
+ x[i] += Ainv(i,k) * f;
}
}
@@ -514,11 +509,12 @@ void matrix_solver_sm_t<m_N, _storage_N>::store(
template <unsigned m_N, unsigned _storage_N>
int matrix_solver_sm_t<m_N, _storage_N>::solve_non_dynamic(ATTR_UNUSED const bool newton_raphson)
{
+ static const bool incremental = true;
static uint cnt = 0;
nl_double new_V[_storage_N]; // = { 0.0 };
- if (0 || (cnt % 100 == 0))
+ if (0 || ((cnt % 200) == 0))
{
/* complete calculation */
this->LE_invert();
@@ -526,40 +522,57 @@ int matrix_solver_sm_t<m_N, _storage_N>::solve_non_dynamic(ATTR_UNUSED const boo
else
{
const auto iN = N();
- for (int row = 0; row < iN; row ++)
- for (int k = 0; k < iN; k++)
- Ainv(row,k) = lAinv(row, k);
- for (int row = 0; row < N(); row ++)
+ if (not incremental)
{
- nl_double v[m_pitch];
- bool changed = false;
- for (int k = 0; k < N(); k++)
- v[k] = A(row,k) - lA(row,k);
- for (int k = 0; k < N(); k++)
- if (v[k] != 0.0)
- {
- changed = true;
- break;
- }
+ for (int row = 0; row < iN; row ++)
+ for (int k = 0; k < iN; k++)
+ Ainv(row,k) = lAinv(row, k);
+ }
+ for (int row = 0; row < iN; row ++)
+ {
+ nl_double v[m_pitch] = {0};
+ unsigned cols[m_pitch];
+ unsigned colcount = 0;
+
+ auto &nz = m_terms[row]->m_nz;
+ for (auto & col : nz)
+ {
+ v[col] = A(row,col) - lA(row,col);
+ if (incremental)
+ lA(row,col) = A(row,col);
+ if (v[col] != 0.0)
+ cols[colcount++] = col;
+ }
- if (changed)
+ if (colcount > 0)
{
nl_double lamba = 0.0;
nl_double w[m_pitch] = {0};
- nl_double z[m_pitch] = {0};
+ nl_double z[m_pitch];
/* compute w and lamba */
- for (int k = 0; k < N(); k++)
+ for (unsigned i = 0; i < iN; i++)
+ z[i] = Ainv(i, row); /* u is row'th column */
+
+ for (unsigned j = 0; j < colcount; j++)
+ lamba += v[cols[j]] * z[cols[j]];
+
+ for (unsigned j=0; j<colcount; j++)
{
- for (int j=0; j<N(); j++)
- w[k] += Ainv(j,k) * v[j]; /* Transpose(Ainv) * v */
- z[k] = Ainv(k, row); /* u is row'th column */
- lamba += v[k] * z[k];
+ auto col = cols[j];
+ auto f = v[col];
+ for (unsigned k = 0; k < iN; k++)
+ w[k] += Ainv(col,k) * f; /* Transpose(Ainv) * v */
}
+
lamba = -1.0 / (1.0 + lamba);
- for (int i=0; i<N(); i++)
- for (int k = 0; k < N(); k++)
- Ainv(i,k) += lamba * z[i] * w[k];
+ for (int i=0; i<iN; i++)
+ {
+ const nl_double f = lamba * z[i];
+ if (f != 0.0)
+ for (int k = 0; k < iN; k++)
+ Ainv(i,k) += f * w[k];
+ }
}
}
diff --git a/src/lib/netlist/solver/nld_solver.cpp b/src/lib/netlist/solver/nld_solver.cpp
index 56033623adf..388f21bda77 100644
--- a/src/lib/netlist/solver/nld_solver.cpp
+++ b/src/lib/netlist/solver/nld_solver.cpp
@@ -309,7 +309,7 @@ ATTR_COLD int matrix_solver_t::get_net_idx(net_t *net)
void matrix_solver_t::log_stats()
{
- //if (this->m_stat_calculations != 0 && this->m_params.m_log_stats)
+ if (this->m_stat_calculations != 0 && this->m_stat_vsolver_calls && this->m_params.m_log_stats)
{
log().verbose("==============================================");
log().verbose("Solver {1}", this->name());
@@ -462,6 +462,12 @@ matrix_solver_t * NETLIB_NAME(solver)::create_solver(int size, const bool use_sp
typedef matrix_solver_sm_t<m_N,_storage_N> solver_mat;
return palloc(solver_mat(&m_params, size));
}
+ else if (pstring("SM").equals(m_iterative_solver))
+ {
+ /* Sherman-Morrison Formula */
+ typedef matrix_solver_sm_t<m_N,_storage_N> solver_mat;
+ return palloc(solver_mat(&m_params, size));
+ }
else if (pstring("SOR").equals(m_iterative_solver))
{
typedef matrix_solver_SOR_t<m_N,_storage_N> solver_GS;