From 93414a8bd75b1dc378e36db0cc86c33c8a8c4bb9 Mon Sep 17 00:00:00 2001 From: couriersud Date: Sun, 27 Mar 2016 17:55:07 +0200 Subject: [PATCH] Fix pivoting and float usage. --- nl_examples/kidniki.c | 1 + src/lib/netlist/nl_util.h | 4 ++-- src/lib/netlist/plib/pfmtlog.h | 4 ++++ src/lib/netlist/solver/nld_ms_direct.h | 25 +++++++++-------------- src/lib/netlist/solver/nld_ms_direct_lu.h | 2 +- 5 files changed, 18 insertions(+), 18 deletions(-) diff --git a/nl_examples/kidniki.c b/nl_examples/kidniki.c index 3b51e40328f..f7dbf69bdfc 100644 --- a/nl_examples/kidniki.c +++ b/nl_examples/kidniki.c @@ -22,6 +22,7 @@ NETLIST_START(dummy) PARAM(Solver.MIN_TIMESTEP, 10e-6) PARAM(Solver.PARALLEL, 0) PARAM(Solver.SOR_FACTOR, 1.00) + PARAM(Solver.PIVOT, 0) #else SOLVER(Solver, 12000) PARAM(Solver.ACCURACY, 1e-8) diff --git a/src/lib/netlist/nl_util.h b/src/lib/netlist/nl_util.h index bb5a564dd74..87ec2270231 100644 --- a/src/lib/netlist/nl_util.h +++ b/src/lib/netlist/nl_util.h @@ -41,8 +41,8 @@ private: public: ATTR_HOT inline static float exp(const float x) { return std::exp(x); } - ATTR_HOT inline static double abs(const double x) { return std::abs(x); } - ATTR_HOT inline static float abs(const float x) { return std::abs(x); } + ATTR_HOT inline static double abs(const double x) { return std::fabs(x); } + ATTR_HOT inline static float abs(const float x) { return std::fabs(x); } ATTR_HOT inline static double log(const double x) { return std::log(x); } ATTR_HOT inline static float log(const float x) { return std::log(x); } #if defined(_MSC_VER) && _MSC_VER < 1800 diff --git a/src/lib/netlist/plib/pfmtlog.h b/src/lib/netlist/plib/pfmtlog.h index 7a63d14fc3a..d72bdc6ce6b 100644 --- a/src/lib/netlist/plib/pfmtlog.h +++ b/src/lib/netlist/plib/pfmtlog.h @@ -110,6 +110,10 @@ public: ATTR_COLD P & e(const double x, const char *f = "") { format_element(f, "", "e", x); return static_cast

(*this); } ATTR_COLD P & g(const double x, const char *f = "") { format_element(f, "", "g", x); return static_cast

(*this); } + ATTR_COLD P &operator ()(const float x, const char *f = "") { format_element(f, "", "f", x); return static_cast

(*this); } + ATTR_COLD P & e(const float x, const char *f = "") { format_element(f, "", "e", x); return static_cast

(*this); } + ATTR_COLD P & g(const float x, const char *f = "") { format_element(f, "", "g", x); return static_cast

(*this); } + ATTR_COLD P &operator ()(const char *x, const char *f = "") { format_element(f, "", "s", x); return static_cast

(*this); } ATTR_COLD P &operator ()(char *x, const char *f = "") { format_element(f, "", "s", x); return static_cast

(*this); } ATTR_COLD P &operator ()(const void *x, const char *f = "") { format_element(f, "", "p", x); return static_cast

(*this); } diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h index 5d0167f8758..852ceebb9ea 100644 --- a/src/lib/netlist/solver/nld_ms_direct.h +++ b/src/lib/netlist/solver/nld_ms_direct.h @@ -24,7 +24,7 @@ NETLIB_NAMESPACE_DEVICES_START() //#define nl_ext_double __float128 // slow, very slow //#define nl_ext_double long double // slightly slower -#define nl_ext_double double +#define nl_ext_double nl_double template class matrix_solver_direct_t: public matrix_solver_t @@ -61,9 +61,6 @@ protected: virtual netlist_time compute_next_timestep() override; - terms_t **m_terms; - terms_t *m_rails_temp; - #if (NL_USE_DYNAMIC_ALLOCATION) template inline nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r * m_pitch + c]; } @@ -78,10 +75,11 @@ protected: ATTR_ALIGN nl_double m_last_RHS[_storage_N]; // right hand side - contains currents ATTR_ALIGN nl_double m_last_V[_storage_N]; - + terms_t * m_terms[_storage_N]; + terms_t *m_rails_temp; private: - static const unsigned m_pitch = (((_storage_N + 1) + 7) / 8) * 8; + static const std::size_t m_pitch = (((_storage_N + 1) + 7) / 8) * 8; #if (NL_USE_DYNAMIC_ALLOCATION) ATTR_ALIGN nl_ext_double * RESTRICT m_A; #else @@ -104,7 +102,6 @@ matrix_solver_direct_t::~matrix_solver_direct_t() { pfree(m_terms[k]); } - pfree_array(m_terms); pfree_array(m_rails_temp); #if (NL_USE_DYNAMIC_ALLOCATION) pfree_array(m_A); @@ -424,10 +421,10 @@ void matrix_solver_direct_t::LE_solve() if (maxrow != i) { /* Swap the maxrow and ith row */ - for (unsigned k = 0; k < kN; k++) { + for (unsigned k = 0; k < kN + 1; k++) { std::swap(A(i,k), A(maxrow,k)); } - std::swap(RHS(i), RHS(maxrow)); + //std::swap(RHS(i), RHS(maxrow)); } /* FIXME: Singular matrix? */ const nl_double f = 1.0 / A(i,i); @@ -439,10 +436,10 @@ void matrix_solver_direct_t::LE_solve() const nl_double f1 = - A(j,i) * f; if (f1 != NL_FCONST(0.0)) { - nl_double * RESTRICT pi = &A(i,i+1); - const nl_double * RESTRICT pj = &A(j,i+1); + const nl_double * RESTRICT pi = &A(i,i+1); + nl_double * RESTRICT pj = &A(j,i+1); #if 1 - vec_add_mult_scalar(kN-i,pj,f1,pi); + vec_add_mult_scalar(kN-i,pi,f1,pj); #else vec_add_mult_scalar(kN-i-1,pj,f1,pi); //for (unsigned k = i+1; k < kN; k++) @@ -528,7 +525,7 @@ T matrix_solver_direct_t::delta( const unsigned iN = this->N(); T cerr = 0; for (unsigned i = 0; i < iN; i++) - cerr = std::max(cerr, nl_math::abs(V[i] - this->m_nets[i]->m_cur_Analog)); + cerr = std::fmax(cerr, nl_math::abs(V[i] - (T) this->m_nets[i]->m_cur_Analog)); return cerr; } @@ -585,7 +582,6 @@ matrix_solver_direct_t::matrix_solver_direct_t(const solver_par : matrix_solver_t(GAUSSIAN_ELIMINATION, params) , m_dim(size) { - m_terms = palloc_array(terms_t *, N()); m_rails_temp = palloc_array(terms_t, N()); #if (NL_USE_DYNAMIC_ALLOCATION) m_A = palloc_array(nl_ext_double, N() * m_pitch); @@ -603,7 +599,6 @@ matrix_solver_direct_t::matrix_solver_direct_t(const eSolverTyp : matrix_solver_t(type, params) , m_dim(size) { - m_terms = palloc_array(terms_t *, N()); m_rails_temp = palloc_array(terms_t, N()); #if (NL_USE_DYNAMIC_ALLOCATION) m_A = palloc_array(nl_ext_double, N() * m_pitch); diff --git a/src/lib/netlist/solver/nld_ms_direct_lu.h b/src/lib/netlist/solver/nld_ms_direct_lu.h index 5efbe38ddab..2109b0f3e41 100644 --- a/src/lib/netlist/solver/nld_ms_direct_lu.h +++ b/src/lib/netlist/solver/nld_ms_direct_lu.h @@ -552,7 +552,7 @@ nl_double matrix_solver_direct_t::delta( const unsigned iN = this->N(); nl_double cerr = 0; for (unsigned i = 0; i < iN; i++) - cerr = std::max(cerr, nl_math::abs(V[i] - this->m_nets[i]->m_cur_Analog)); + cerr = std::fmax(cerr, nl_math::abs(V[i] - this->m_nets[i]->m_cur_Analog)); return cerr; }