diff options
author | couriersud <couriersud@arcor.de> | 2016-03-28 03:20:59 +0200 |
---|---|---|
committer | couriersud <couriersud@arcor.de> | 2016-03-28 13:48:27 +0200 |
commit | 8c69d3ad1562d79790bd48e3fe4521db304efceb (patch) | |
tree | d93a9cf6a6978e3a0f8ab5e1b30b55642247c034 | |
parent | 97b9fc11d63d578242713e1a2e94c538766febbe (diff) |
Fix crash in stat output when no calculations were run
-rw-r--r-- | nl_examples/kidniki.c | 3 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_matrix_solver.h | 2 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_ms_direct.h | 14 | ||||
-rw-r--r-- | src/lib/netlist/solver/nld_ms_sm.h | 91 | ||||
-rw-r--r-- | 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<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; |