From 8c69d3ad1562d79790bd48e3fe4521db304efceb Mon Sep 17 00:00:00 2001 From: couriersud Date: Mon, 28 Mar 2016 03:20:59 +0200 Subject: Fix crash in stat output when no calculations were run --- nl_examples/kidniki.c | 3 +- src/lib/netlist/solver/nld_matrix_solver.h | 2 +- src/lib/netlist/solver/nld_ms_direct.h | 14 ++--- src/lib/netlist/solver/nld_ms_sm.h | 91 +++++++++++++++++------------- src/lib/netlist/solver/nld_solver.cpp | 8 ++- 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 m_nz; /* all non zero for multiplication */ - pvector_t m_nzrd; /* non zero right of the diagonal for elimination, includes RHS element */ + pvector_t m_nzrd; /* non zero right of the diagonal for elimination, may include RHS element */ pvector_t m_nzbd; /* non zero below of the diagonal for elimination */ private: pvector_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::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::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::LE_compute_x( for (int k=0; k::store( template int matrix_solver_sm_t::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::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; jm_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 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 solver_mat; + return palloc(solver_mat(&m_params, size)); + } else if (pstring("SOR").equals(m_iterative_solver)) { typedef matrix_solver_SOR_t solver_GS; -- cgit v1.2.3