Fix pivoting and float usage.

This commit is contained in:
couriersud 2016-03-27 17:55:07 +02:00
parent b8d22b6c90
commit 93414a8bd7
5 changed files with 18 additions and 18 deletions

View File

@ -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)

View File

@ -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

View File

@ -110,6 +110,10 @@ public:
ATTR_COLD P & e(const double x, const char *f = "") { format_element(f, "", "e", x); return static_cast<P &>(*this); }
ATTR_COLD P & g(const double x, const char *f = "") { format_element(f, "", "g", x); return static_cast<P &>(*this); }
ATTR_COLD P &operator ()(const float x, const char *f = "") { format_element(f, "", "f", x); return static_cast<P &>(*this); }
ATTR_COLD P & e(const float x, const char *f = "") { format_element(f, "", "e", x); return static_cast<P &>(*this); }
ATTR_COLD P & g(const float x, const char *f = "") { format_element(f, "", "g", x); return static_cast<P &>(*this); }
ATTR_COLD P &operator ()(const char *x, const char *f = "") { format_element(f, "", "s", x); return static_cast<P &>(*this); }
ATTR_COLD P &operator ()(char *x, const char *f = "") { format_element(f, "", "s", x); return static_cast<P &>(*this); }
ATTR_COLD P &operator ()(const void *x, const char *f = "") { format_element(f, "", "p", x); return static_cast<P &>(*this); }

View File

@ -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 <unsigned m_N, unsigned _storage_N>
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 <typename T1, typename T2>
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<m_N, _storage_N>::~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<m_N, _storage_N>::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<m_N, _storage_N>::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<m_N, _storage_N>::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<m_N, _storage_N>::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<m_N, _storage_N>::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);

View File

@ -552,7 +552,7 @@ nl_double matrix_solver_direct_t<m_N, _storage_N>::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;
}