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 /src/lib/netlist/solver/nld_ms_sm.h | |
parent | 97b9fc11d63d578242713e1a2e94c538766febbe (diff) |
Fix crash in stat output when no calculations were run
Diffstat (limited to 'src/lib/netlist/solver/nld_ms_sm.h')
-rw-r--r-- | src/lib/netlist/solver/nld_ms_sm.h | 91 |
1 files changed, 52 insertions, 39 deletions
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]; + } } } |