- more code optimization
- hide matrix_solver_t implementation
- use netlist_time for time deltas
This commit is contained in:
couriersud 2016-03-27 03:15:27 +02:00
parent 2e21930b27
commit 9d2f61ee92
13 changed files with 248 additions and 256 deletions

View File

@ -8,14 +8,18 @@
#define USE_FIXED_STV 1 #define USE_FIXED_STV 1
NETLIST_START(dummy) NETLIST_START(dummy)
#if (USE_FRONTIERS) #if (1 || USE_FRONTIERS)
SOLVER(Solver, 18000) SOLVER(Solver, 18000)
PARAM(Solver.ACCURACY, 1e-8) PARAM(Solver.ACCURACY, 1e-8)
PARAM(Solver.NR_LOOPS, 300) PARAM(Solver.NR_LOOPS, 300)
PARAM(Solver.GS_LOOPS, 1) PARAM(Solver.GS_LOOPS, 1)
PARAM(Solver.GS_THRESHOLD, 6) PARAM(Solver.GS_THRESHOLD, 6)
//PARAM(Solver.ITERATIVE, "MAT") //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)
PARAM(Solver.PARALLEL, 0) PARAM(Solver.PARALLEL, 0)
PARAM(Solver.SOR_FACTOR, 1.00) PARAM(Solver.SOR_FACTOR, 1.00)
#else #else
@ -46,14 +50,15 @@ NETLIST_START(dummy)
//ANALOG_INPUT(I_V0, 0) //ANALOG_INPUT(I_V0, 0)
ALIAS(I_V0.Q, GND) ALIAS(I_V0.Q, GND)
#if 0 #if 0
ANALOG_INPUT(I_SD0, 0) TTL_INPUT(I_SD0, 0)
ANALOG_INPUT(I_BD0, 0) TTL_INPUT(I_BD0, 0)
ANALOG_INPUT(I_CH0, 0) TTL_INPUT(I_CH0, 0)
ANALOG_INPUT(I_OH0, 0) TTL_INPUT(I_OH0, 0)
ANALOG_INPUT(I_SOUNDIC0, 0) TTL_INPUT(I_SOUNDIC0, 0)
ANALOG_INPUT(I_OKI0, 0) ANALOG_INPUT(I_MSM2K0, 0)
ANALOG_INPUT(I_SOUND0, 0) ANALOG_INPUT(I_MSM3K0, 0)
ANALOG_INPUT(I_SINH0, 0) TTL_INPUT(I_SOUND0, 0)
TTL_INPUT(I_SINH0, 0)
#else #else
TTL_INPUT(I_SD0, 1) TTL_INPUT(I_SD0, 1)
//CLOCK(I_SD0, 5) //CLOCK(I_SD0, 5)

View File

@ -73,6 +73,7 @@ project "netlist"
MAME_DIR .. "src/lib/netlist/analog/nld_opamps.h", MAME_DIR .. "src/lib/netlist/analog/nld_opamps.h",
MAME_DIR .. "src/lib/netlist/solver/nld_solver.cpp", MAME_DIR .. "src/lib/netlist/solver/nld_solver.cpp",
MAME_DIR .. "src/lib/netlist/solver/nld_solver.h", MAME_DIR .. "src/lib/netlist/solver/nld_solver.h",
MAME_DIR .. "src/lib/netlist/solver/nld_matrix_solver.h",
MAME_DIR .. "src/lib/netlist/solver/nld_ms_direct.h", MAME_DIR .. "src/lib/netlist/solver/nld_ms_direct.h",
MAME_DIR .. "src/lib/netlist/solver/nld_ms_direct1.h", MAME_DIR .. "src/lib/netlist/solver/nld_ms_direct1.h",
MAME_DIR .. "src/lib/netlist/solver/nld_ms_direct2.h", MAME_DIR .. "src/lib/netlist/solver/nld_ms_direct2.h",

View File

@ -5,7 +5,7 @@
* *
*/ */
#include <solver/nld_solver.h> #include <solver/nld_matrix_solver.h>
#include "nld_mm5837.h" #include "nld_mm5837.h"
#include "nl_setup.h" #include "nl_setup.h"

View File

@ -6,6 +6,7 @@
*/ */
#include <solver/nld_solver.h> #include <solver/nld_solver.h>
#include <solver/nld_matrix_solver.h>
#include "nld_system.h" #include "nld_system.h"
NETLIB_NAMESPACE_DEVICES_START() NETLIB_NAMESPACE_DEVICES_START()

View File

@ -5,7 +5,8 @@
* *
*/ */
#include <solver/nld_solver.h> #include <solver/nld_matrix_solver.h>
#include <cstring> #include <cstring>
#include <algorithm> #include <algorithm>

View File

@ -1062,7 +1062,7 @@ namespace netlist
ATTR_HOT nl_double INPANALOG(const analog_input_t &inp) const { return inp.Q_Analog(); } ATTR_HOT nl_double INPANALOG(const analog_input_t &inp) const { return inp.Q_Analog(); }
ATTR_HOT nl_double TERMANALOG(const terminal_t &term) const { return term.net().as_analog().Q_Analog(); } ATTR_HOT nl_double TERMANALOG(const terminal_t &term) const { return term.net().Q_Analog(); }
ATTR_HOT void OUTANALOG(analog_output_t &out, const nl_double val) ATTR_HOT void OUTANALOG(analog_output_t &out, const nl_double val)
{ {

View File

@ -0,0 +1,144 @@
// license:GPL-2.0+
// copyright-holders:Couriersud
/*
* nld_matrix_solver.h
*
*/
#ifndef NLD_MATRIX_SOLVER_H_
#define NLD_MATRIX_SOLVER_H_
#include "solver/nld_solver.h"
NETLIB_NAMESPACE_DEVICES_START()
class terms_t
{
P_PREVENT_COPYING(terms_t)
public:
ATTR_COLD terms_t() : m_railstart(0)
{}
ATTR_COLD void clear()
{
m_term.clear();
m_net_other.clear();
m_gt.clear();
m_go.clear();
m_Idr.clear();
m_other_curanalog.clear();
}
ATTR_COLD void add(terminal_t *term, int net_other, bool sorted);
inline unsigned count() { return m_term.size(); }
inline terminal_t **terms() { return m_term.data(); }
inline int *net_other() { return m_net_other.data(); }
inline nl_double *gt() { return m_gt.data(); }
inline nl_double *go() { return m_go.data(); }
inline nl_double *Idr() { return m_Idr.data(); }
inline nl_double **other_curanalog() { return m_other_curanalog.data(); }
ATTR_COLD void set_pointers();
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_nzbd; /* non zero below of the diagonal for elimination */
private:
pvector_t<terminal_t *> m_term;
pvector_t<int> m_net_other;
pvector_t<nl_double> m_go;
pvector_t<nl_double> m_gt;
pvector_t<nl_double> m_Idr;
pvector_t<nl_double *> m_other_curanalog;
};
class matrix_solver_t : public device_t
{
public:
typedef pvector_t<matrix_solver_t *> list_t;
typedef core_device_t::list_t dev_list_t;
enum eSolverType
{
GAUSSIAN_ELIMINATION,
GAUSS_SEIDEL
};
matrix_solver_t(const eSolverType type, const solver_parameters_t *params);
virtual ~matrix_solver_t();
void setup(analog_net_t::list_t &nets) { vsetup(nets); }
netlist_time solve_base();
netlist_time solve();
inline bool is_dynamic() const { return m_dynamic_devices.size() > 0; }
inline bool is_timestep() const { return m_step_devices.size() > 0; }
void update_forced();
void update_after(const netlist_time after)
{
m_Q_sync.net().reschedule_in_queue(after);
}
/* netdevice functions */
virtual void update() override;
virtual void start() override;
virtual void reset() override;
ATTR_COLD int get_net_idx(net_t *net);
eSolverType type() const { return m_type; }
plog_base<NL_DEBUG> &log() { return netlist().log(); }
virtual void log_stats();
protected:
ATTR_COLD void setup_base(analog_net_t::list_t &nets);
void update_dynamic();
virtual void add_term(int net_idx, terminal_t *term) = 0;
virtual void vsetup(analog_net_t::list_t &nets) = 0;
virtual int vsolve_non_dynamic(const bool newton_raphson) = 0;
virtual netlist_time compute_next_timestep() = 0;
pvector_t<analog_net_t *> m_nets;
pvector_t<analog_output_t *> m_inps;
int m_stat_calculations;
int m_stat_newton_raphson;
int m_stat_vsolver_calls;
int m_iterative_fail;
int m_iterative_total;
const solver_parameters_t &m_params;
inline nl_double current_timestep() { return m_cur_ts; }
private:
netlist_time m_last_step;
nl_double m_cur_ts;
dev_list_t m_step_devices;
dev_list_t m_dynamic_devices;
logic_input_t m_fb_sync;
logic_output_t m_Q_sync;
void step(const netlist_time &delta);
void update_inputs();
const eSolverType m_type;
};
NETLIB_NAMESPACE_DEVICES_END()
#endif /* NLD_MS_DIRECT_H_ */

View File

@ -39,32 +39,30 @@ public:
virtual void vsetup(analog_net_t::list_t &nets) override; virtual void vsetup(analog_net_t::list_t &nets) override;
virtual void reset() override { matrix_solver_t::reset(); } virtual void reset() override { matrix_solver_t::reset(); }
inline unsigned N() const { if (m_N == 0) return m_dim; else return m_N; }
virtual int vsolve_non_dynamic(const bool newton_raphson) override;
protected: protected:
virtual void add_term(int net_idx, terminal_t *term) override; virtual void add_term(int net_idx, terminal_t *term) override;
virtual int vsolve_non_dynamic(const bool newton_raphson) override;
int solve_non_dynamic(const bool newton_raphson); int solve_non_dynamic(const bool newton_raphson);
inline const unsigned N() const { if (m_N == 0) return m_dim; else return m_N; }
void build_LE_A(); void build_LE_A();
void build_LE_RHS(); void build_LE_RHS();
void LE_solve(); void LE_solve();
void LE_back_subst(nl_double * RESTRICT x);
/* Full LU back substitution, not used currently, in for future use */ template <typename T>
void LE_back_subst(T * RESTRICT x);
void LE_back_subst_full(nl_double * RESTRICT x); template <typename T>
T delta(const T * RESTRICT V);
nl_double delta(const nl_double * RESTRICT V); template <typename T>
void store(const nl_double * RESTRICT V); void store(const T * RESTRICT V);
/* bring the whole system to the current time virtual netlist_time compute_next_timestep() override;
* Don't schedule a new calculation time. The recalculation has to be
* triggered by the caller after the netlist element was changed.
*/
virtual nl_double compute_next_timestep() override;
terms_t **m_terms;
terms_t *m_rails_temp;
#if (NL_USE_DYNAMIC_ALLOCATION) #if (NL_USE_DYNAMIC_ALLOCATION)
template <typename T1, typename T2> template <typename T1, typename T2>
@ -76,13 +74,11 @@ protected:
inline nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r][c]; } inline nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r][c]; }
template <typename T1> template <typename T1>
inline nl_ext_double &RHS(const T1 &r) { return m_A[r][N()]; } inline nl_ext_double &RHS(const T1 &r) { return m_A[r][N()]; }
#endif #endif
ATTR_ALIGN nl_double m_last_RHS[_storage_N]; // right hand side - contains currents ATTR_ALIGN nl_double m_last_RHS[_storage_N]; // right hand side - contains currents
ATTR_ALIGN nl_double m_last_V[_storage_N]; ATTR_ALIGN nl_double m_last_V[_storage_N];
terms_t **m_terms;
terms_t *m_rails_temp;
private: private:
static const unsigned m_pitch = (((_storage_N + 1) + 7) / 8) * 8; static const unsigned m_pitch = (((_storage_N + 1) + 7) / 8) * 8;
@ -94,6 +90,7 @@ private:
//ATTR_ALIGN nl_ext_double m_RHSx[_storage_N]; //ATTR_ALIGN nl_ext_double m_RHSx[_storage_N];
const unsigned m_dim; const unsigned m_dim;
}; };
// ---------------------------------------------------------------------------------------- // ----------------------------------------------------------------------------------------
@ -115,7 +112,7 @@ matrix_solver_direct_t<m_N, _storage_N>::~matrix_solver_direct_t()
} }
template <unsigned m_N, unsigned _storage_N> template <unsigned m_N, unsigned _storage_N>
nl_double matrix_solver_direct_t<m_N, _storage_N>::compute_next_timestep() netlist_time matrix_solver_direct_t<m_N, _storage_N>::compute_next_timestep()
{ {
nl_double new_solver_timestep = m_params.m_max_timestep; nl_double new_solver_timestep = m_params.m_max_timestep;
@ -152,7 +149,7 @@ nl_double matrix_solver_direct_t<m_N, _storage_N>::compute_next_timestep()
} }
//if (new_solver_timestep > 10.0 * hn) //if (new_solver_timestep > 10.0 * hn)
// new_solver_timestep = 10.0 * hn; // new_solver_timestep = 10.0 * hn;
return new_solver_timestep; return netlist_time::from_double(new_solver_timestep);
} }
template <unsigned m_N, unsigned _storage_N> template <unsigned m_N, unsigned _storage_N>
@ -284,9 +281,15 @@ ATTR_COLD void matrix_solver_direct_t<m_N, _storage_N>::vsetup(analog_net_t::lis
t->m_nz.push_back(other[i]); t->m_nz.push_back(other[i]);
} }
} }
/* Add RHS element */
if (!t->m_nzrd.contains(N()))
t->m_nzrd.push_back(N());
/* and sort */
psort_list(t->m_nzrd); psort_list(t->m_nzrd);
t->m_nz.push_back(k); // add diagonal t->m_nz.push_back(k); // add diagonal
psort_list(t->m_nz); psort_list(t->m_nz);
} }
@ -437,7 +440,7 @@ void matrix_solver_direct_t<m_N, _storage_N>::LE_solve()
if (f1 != NL_FCONST(0.0)) if (f1 != NL_FCONST(0.0))
{ {
nl_double * RESTRICT pi = &A(i,i+1); nl_double * RESTRICT pi = &A(i,i+1);
nl_double * RESTRICT pj = &A(j,i+1); const nl_double * RESTRICT pj = &A(j,i+1);
#if 1 #if 1
vec_add_mult_scalar(kN-i,pj,f1,pi); vec_add_mult_scalar(kN-i,pj,f1,pi);
#else #else
@ -466,23 +469,18 @@ void matrix_solver_direct_t<m_N, _storage_N>::LE_solve()
{ {
const unsigned j = pb[jb]; const unsigned j = pb[jb];
const nl_double f1 = - A(j,i) * f; const nl_double f1 = - A(j,i) * f;
#if 0
nl_double * RESTRICT pi = &m_A[i][i+1];
nl_double * RESTRICT pj = &m_A[j][i+1];
vec_add_mult_scalar(kN-i-1,pi,f1,pj);
#else
for (unsigned k = 0; k < e; k++) for (unsigned k = 0; k < e; k++)
A(j,p[k]) += A(i,p[k]) * f1; A(j,p[k]) += A(i,p[k]) * f1;
#endif //RHS(j) += RHS(i) * f1;
RHS(j) += RHS(i) * f1;
} }
} }
} }
} }
template <unsigned m_N, unsigned _storage_N> template <unsigned m_N, unsigned _storage_N>
template <typename T>
void matrix_solver_direct_t<m_N, _storage_N>::LE_back_subst( void matrix_solver_direct_t<m_N, _storage_N>::LE_back_subst(
nl_double * RESTRICT x) T * RESTRICT x)
{ {
const unsigned kN = N(); const unsigned kN = N();
@ -491,7 +489,7 @@ void matrix_solver_direct_t<m_N, _storage_N>::LE_back_subst(
{ {
for (int j = kN - 1; j >= 0; j--) for (int j = kN - 1; j >= 0; j--)
{ {
nl_double tmp = 0; T tmp = 0;
for (unsigned k = j+1; k < kN; k++) for (unsigned k = j+1; k < kN; k++)
tmp += A(j,k) * x[k]; tmp += A(j,k) * x[k];
x[j] = (RHS(j) - tmp) / A(j,j); x[j] = (RHS(j) - tmp) / A(j,j);
@ -501,10 +499,10 @@ void matrix_solver_direct_t<m_N, _storage_N>::LE_back_subst(
{ {
for (int j = kN - 1; j >= 0; j--) for (int j = kN - 1; j >= 0; j--)
{ {
nl_double tmp = 0; T tmp = 0;
const unsigned *p = m_terms[j]->m_nzrd.data(); const unsigned *p = m_terms[j]->m_nzrd.data();
const unsigned e = m_terms[j]->m_nzrd.size(); const unsigned e = m_terms[j]->m_nzrd.size() - 1; /* exclude RHS element */
for (unsigned k = 0; k < e; k++) for (unsigned k = 0; k < e; k++)
{ {
@ -516,43 +514,11 @@ void matrix_solver_direct_t<m_N, _storage_N>::LE_back_subst(
} }
} }
template <unsigned m_N, unsigned _storage_N>
void matrix_solver_direct_t<m_N, _storage_N>::LE_back_subst_full(
nl_double * RESTRICT x)
{
const unsigned kN = N();
/* back substitution */
// int ip;
// ii=-1
//for (int i=0; i < kN; i++)
// x[i] = m_RHS[i];
for (int i=0; i < kN; i++)
{
//ip=indx[i]; USE_PIVOT_SEARCH
//sum=b[ip];
//b[ip]=b[i];
double sum=RHS(i);//x[i];
for (int j=0; j < i; j++)
sum -= A(i,j) * x[j];
x[i]=sum;
}
for (int i=kN-1; i >= 0; i--)
{
double sum=x[i];
for (int j = i+1; j < kN; j++)
sum -= A(i,j)*x[j];
x[i] = sum / A(i,i);
}
}
template <unsigned m_N, unsigned _storage_N> template <unsigned m_N, unsigned _storage_N>
nl_double matrix_solver_direct_t<m_N, _storage_N>::delta( template <typename T>
const nl_double * RESTRICT V) T matrix_solver_direct_t<m_N, _storage_N>::delta(
const T * RESTRICT V)
{ {
/* FIXME: Ideally we should also include currents (RHS) here. This would /* FIXME: Ideally we should also include currents (RHS) here. This would
* need a revaluation of the right hand side after voltages have been updated * need a revaluation of the right hand side after voltages have been updated
@ -560,15 +526,16 @@ nl_double matrix_solver_direct_t<m_N, _storage_N>::delta(
*/ */
const unsigned iN = this->N(); const unsigned iN = this->N();
nl_double cerr = 0; T cerr = 0;
for (unsigned i = 0; i < iN; i++) 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::max(cerr, nl_math::abs(V[i] - this->m_nets[i]->m_cur_Analog));
return cerr; return cerr;
} }
template <unsigned m_N, unsigned _storage_N> template <unsigned m_N, unsigned _storage_N>
template <typename T>
void matrix_solver_direct_t<m_N, _storage_N>::store( void matrix_solver_direct_t<m_N, _storage_N>::store(
const nl_double * RESTRICT V) const T * RESTRICT V)
{ {
for (unsigned i = 0, iN=N(); i < iN; i++) for (unsigned i = 0, iN=N(); i < iN; i++)
{ {

View File

@ -30,25 +30,23 @@ public:
inline int matrix_solver_direct1_t::vsolve_non_dynamic(ATTR_UNUSED const bool newton_raphson) inline int matrix_solver_direct1_t::vsolve_non_dynamic(ATTR_UNUSED const bool newton_raphson)
{ {
analog_net_t *net = m_nets[0];
this->build_LE_A(); this->build_LE_A();
this->build_LE_RHS(); this->build_LE_RHS();
//NL_VERBOSE_OUT(("{1} {2}\n", new_val, m_RHS[0] / m_A[0][0]); //NL_VERBOSE_OUT(("{1} {2}\n", new_val, m_RHS[0] / m_A[0][0]);
nl_double new_val = RHS(0) / A(0,0); nl_double new_val[1] = { RHS(0) / A(0,0) };
nl_double e = (new_val - net->Q_Analog()); if (is_dynamic())
nl_double cerr = nl_math::abs(e);
net->m_cur_Analog = new_val;
if (is_dynamic() && (cerr > m_params.m_accuracy))
{ {
return 2; nl_double err = this->delta(new_val);
store(new_val);
if (err > m_params.m_accuracy )
return 2;
else
return 1;
} }
else store(new_val);
return 1; return 1;
} }
NETLIB_NAMESPACE_DEVICES_END() NETLIB_NAMESPACE_DEVICES_END()

View File

@ -58,7 +58,7 @@ public:
private: private:
int solve_ilu_gmres(nl_double * RESTRICT x, nl_double * RESTRICT rhs, const unsigned restart_max, const unsigned mr, nl_double accuracy); int solve_ilu_gmres(nl_double * RESTRICT x, const nl_double * RESTRICT rhs, const unsigned restart_max, const unsigned mr, nl_double accuracy);
pvector_t<int> m_term_cr[_storage_N]; pvector_t<int> m_term_cr[_storage_N];
@ -157,6 +157,7 @@ int matrix_solver_GMRES_t<m_N, _storage_N>::vsolve_non_dynamic(const bool newton
const nl_double * const * RESTRICT other_cur_analog = this->m_terms[k]->other_curanalog(); const nl_double * const * RESTRICT other_cur_analog = this->m_terms[k]->other_curanalog();
l_V[k] = new_V[k] = this->m_nets[k]->m_cur_Analog; l_V[k] = new_V[k] = this->m_nets[k]->m_cur_Analog;
for (unsigned i = 0; i < term_count; i++) for (unsigned i = 0; i < term_count; i++)
{ {
gtot_t = gtot_t + gt[i]; gtot_t = gtot_t + gt[i];
@ -180,7 +181,7 @@ int matrix_solver_GMRES_t<m_N, _storage_N>::vsolve_non_dynamic(const bool newton
mat.ia[iN] = mat.nz_num; mat.ia[iN] = mat.nz_num;
const nl_double accuracy = this->m_params.m_accuracy; const nl_double accuracy = this->m_params.m_accuracy;
#if 1
int mr = _storage_N; int mr = _storage_N;
if (_storage_N > 3 ) if (_storage_N > 3 )
mr = (int) sqrt(iN); mr = (int) sqrt(iN);
@ -188,19 +189,12 @@ int matrix_solver_GMRES_t<m_N, _storage_N>::vsolve_non_dynamic(const bool newton
int iter = 4; int iter = 4;
int gsl = solve_ilu_gmres(new_V, RHS, iter, mr, accuracy); int gsl = solve_ilu_gmres(new_V, RHS, iter, mr, accuracy);
int failed = mr * iter; int failed = mr * iter;
#else
int failed = 6;
//int gsl = tt_ilu_cr(new_V, RHS, failed, accuracy);
int gsl = tt_gs_cr(new_V, RHS, failed, accuracy);
#endif
this->m_iterative_total += gsl; this->m_iterative_total += gsl;
this->m_stat_calculations++; this->m_stat_calculations++;
if (gsl>=failed) if (gsl>=failed)
{ {
//for (int k = 0; k < iN; k++)
// this->m_nets[k]->m_cur_Analog = new_V[k];
// Fallback to direct solver ...
this->m_iterative_fail++; this->m_iterative_fail++;
return matrix_solver_direct_t<m_N, _storage_N>::vsolve_non_dynamic(newton_raphson); return matrix_solver_direct_t<m_N, _storage_N>::vsolve_non_dynamic(newton_raphson);
} }
@ -226,17 +220,18 @@ int matrix_solver_GMRES_t<m_N, _storage_N>::vsolve_non_dynamic(const bool newton
} }
} }
static inline void givens_mult( const nl_double c, const nl_double s, nl_double * RESTRICT g0, nl_double * RESTRICT g1 ) template <typename T>
inline void givens_mult( const T c, const T s, T & g0, T & g1 )
{ {
const double tg0 = c * *g0 - s * *g1; const T tg0 = c * g0 - s * g1;
const double tg1 = s * *g0 + c * *g1; const T tg1 = s * g0 + c * g1;
*g0 = tg0; g0 = tg0;
*g1 = tg1; g1 = tg1;
} }
template <unsigned m_N, unsigned _storage_N> template <unsigned m_N, unsigned _storage_N>
int matrix_solver_GMRES_t<m_N, _storage_N>::solve_ilu_gmres (nl_double * RESTRICT x, nl_double * RESTRICT rhs, const unsigned restart_max, const unsigned mr, nl_double accuracy) int matrix_solver_GMRES_t<m_N, _storage_N>::solve_ilu_gmres (nl_double * RESTRICT x, const nl_double * RESTRICT rhs, const unsigned restart_max, const unsigned mr, nl_double accuracy)
{ {
/*------------------------------------------------------------------------- /*-------------------------------------------------------------------------
* The code below was inspired by code published by John Burkardt under * The code below was inspired by code published by John Burkardt under
@ -342,7 +337,7 @@ int matrix_solver_GMRES_t<m_N, _storage_N>::solve_ilu_gmres (nl_double * RESTRIC
vec_scale(n, m_v[k1], NL_FCONST(1.0) / m_ht[k1][k]); vec_scale(n, m_v[k1], NL_FCONST(1.0) / m_ht[k1][k]);
for (unsigned j = 0; j < k; j++) for (unsigned j = 0; j < k; j++)
givens_mult(m_c[j], m_s[j], &m_ht[j][k], &m_ht[j+1][k]); givens_mult(m_c[j], m_s[j], m_ht[j][k], m_ht[j+1][k]);
mu = std::sqrt(std::pow(m_ht[k][k], 2) + std::pow(m_ht[k1][k], 2)); mu = std::sqrt(std::pow(m_ht[k][k], 2) + std::pow(m_ht[k1][k], 2));
@ -351,7 +346,7 @@ int matrix_solver_GMRES_t<m_N, _storage_N>::solve_ilu_gmres (nl_double * RESTRIC
m_ht[k][k] = m_c[k] * m_ht[k][k] - m_s[k] * m_ht[k1][k]; m_ht[k][k] = m_c[k] * m_ht[k][k] - m_s[k] * m_ht[k1][k];
m_ht[k1][k] = 0.0; m_ht[k1][k] = 0.0;
givens_mult(m_c[k], m_s[k], &m_g[k], &m_g[k1]); givens_mult(m_c[k], m_s[k], m_g[k], m_g[k1]);
rho = std::abs(m_g[k1]); rho = std::abs(m_g[k1]);

View File

@ -31,7 +31,16 @@
#include <iostream> #include <iostream>
#include <algorithm> #include <algorithm>
//#include "nld_twoterm.h"
#include "nl_lists.h"
#if HAS_OPENMP
#include "omp.h"
#endif
#include "nld_solver.h" #include "nld_solver.h"
#include "nld_matrix_solver.h"
#if 1 #if 1
#include "nld_ms_direct.h" #include "nld_ms_direct.h"
#else #else
@ -42,12 +51,6 @@
#include "nld_ms_sor.h" #include "nld_ms_sor.h"
#include "nld_ms_sor_mat.h" #include "nld_ms_sor_mat.h"
#include "nld_ms_gmres.h" #include "nld_ms_gmres.h"
//#include "nld_twoterm.h"
#include "nl_lists.h"
#if HAS_OPENMP
#include "omp.h"
#endif
NETLIB_NAMESPACE_DEVICES_START() NETLIB_NAMESPACE_DEVICES_START()
@ -222,28 +225,28 @@ ATTR_COLD void matrix_solver_t::reset()
ATTR_COLD void matrix_solver_t::update() ATTR_COLD void matrix_solver_t::update()
{ {
const nl_double new_timestep = solve(); const netlist_time new_timestep = solve();
if (m_params.m_dynamic && is_timestep() && new_timestep > 0) if (m_params.m_dynamic && is_timestep() && new_timestep > netlist_time::zero)
m_Q_sync.net().reschedule_in_queue(netlist_time::from_double(new_timestep)); m_Q_sync.net().reschedule_in_queue(new_timestep);
} }
ATTR_COLD void matrix_solver_t::update_forced() ATTR_COLD void matrix_solver_t::update_forced()
{ {
ATTR_UNUSED const nl_double new_timestep = solve(); ATTR_UNUSED const netlist_time new_timestep = solve();
if (m_params.m_dynamic && is_timestep()) if (m_params.m_dynamic && is_timestep())
m_Q_sync.net().reschedule_in_queue(netlist_time::from_double(m_params.m_min_timestep)); m_Q_sync.net().reschedule_in_queue(netlist_time::from_double(m_params.m_min_timestep));
} }
void matrix_solver_t::step(const netlist_time delta) void matrix_solver_t::step(const netlist_time &delta)
{ {
const nl_double dd = delta.as_double(); const nl_double dd = delta.as_double();
for (std::size_t k=0; k < m_step_devices.size(); k++) for (std::size_t k=0; k < m_step_devices.size(); k++)
m_step_devices[k]->step_time(dd); m_step_devices[k]->step_time(dd);
} }
nl_double matrix_solver_t::solve_base() netlist_time matrix_solver_t::solve_base()
{ {
m_stat_vsolver_calls++; m_stat_vsolver_calls++;
if (is_dynamic()) if (is_dynamic())
@ -273,7 +276,7 @@ nl_double matrix_solver_t::solve_base()
return this->compute_next_timestep(); return this->compute_next_timestep();
} }
nl_double matrix_solver_t::solve() netlist_time matrix_solver_t::solve()
{ {
const netlist_time now = netlist().time(); const netlist_time now = netlist().time();
const netlist_time delta = now - m_last_step; const netlist_time delta = now - m_last_step;
@ -281,7 +284,7 @@ nl_double matrix_solver_t::solve()
// We are already up to date. Avoid oscillations. // We are already up to date. Avoid oscillations.
// FIXME: Make this a parameter! // FIXME: Make this a parameter!
if (delta < netlist_time::from_nsec(1)) // 20000 if (delta < netlist_time::from_nsec(1)) // 20000
return -1.0; return netlist_time::from_nsec(0);
/* update all terminals for new time step */ /* update all terminals for new time step */
m_last_step = now; m_last_step = now;
@ -289,7 +292,7 @@ nl_double matrix_solver_t::solve()
step(delta); step(delta);
const nl_double next_time_step = solve_base(); const netlist_time next_time_step = solve_base();
update_inputs(); update_inputs();
return next_time_step; return next_time_step;
@ -427,7 +430,7 @@ NETLIB_UPDATE(solver)
for (auto & solver : m_mat_solvers) for (auto & solver : m_mat_solvers)
if (solver->is_timestep()) if (solver->is_timestep())
// Ignore return value // Ignore return value
ATTR_UNUSED const nl_double ts = solver->solve(); ATTR_UNUSED const netlist_time ts = solver->solve();
#endif #endif
/* step circuit */ /* step circuit */

View File

@ -48,133 +48,7 @@ struct solver_parameters_t
}; };
class terms_t class matrix_solver_t;
{
P_PREVENT_COPYING(terms_t)
public:
ATTR_COLD terms_t() : m_railstart(0)
{}
ATTR_COLD void clear()
{
m_term.clear();
m_net_other.clear();
m_gt.clear();
m_go.clear();
m_Idr.clear();
m_other_curanalog.clear();
}
ATTR_COLD void add(terminal_t *term, int net_other, bool sorted);
inline unsigned count() { return m_term.size(); }
inline terminal_t **terms() { return m_term.data(); }
inline int *net_other() { return m_net_other.data(); }
inline nl_double *gt() { return m_gt.data(); }
inline nl_double *go() { return m_go.data(); }
inline nl_double *Idr() { return m_Idr.data(); }
inline nl_double **other_curanalog() { return m_other_curanalog.data(); }
ATTR_COLD void set_pointers();
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 */
pvector_t<unsigned> m_nzbd; /* non zero below of the diagonal for elimination */
private:
pvector_t<terminal_t *> m_term;
pvector_t<int> m_net_other;
pvector_t<nl_double> m_go;
pvector_t<nl_double> m_gt;
pvector_t<nl_double> m_Idr;
pvector_t<nl_double *> m_other_curanalog;
};
class matrix_solver_t : public device_t
{
public:
typedef pvector_t<matrix_solver_t *> list_t;
typedef core_device_t::list_t dev_list_t;
enum eSolverType
{
GAUSSIAN_ELIMINATION,
GAUSS_SEIDEL
};
ATTR_COLD matrix_solver_t(const eSolverType type, const solver_parameters_t *params);
virtual ~matrix_solver_t();
void setup(analog_net_t::list_t &nets) { vsetup(nets); }
nl_double solve_base();
nl_double solve();
inline bool is_dynamic() const { return m_dynamic_devices.size() > 0; }
inline bool is_timestep() const { return m_step_devices.size() > 0; }
void update_forced();
void update_after(const netlist_time after)
{
m_Q_sync.net().reschedule_in_queue(after);
}
/* netdevice functions */
virtual void update() override;
virtual void start() override;
virtual void reset() override;
ATTR_COLD int get_net_idx(net_t *net);
eSolverType type() const { return m_type; }
plog_base<NL_DEBUG> &log() { return netlist().log(); }
virtual void log_stats();
protected:
ATTR_COLD void setup_base(analog_net_t::list_t &nets);
void update_dynamic();
virtual void add_term(int net_idx, terminal_t *term) = 0;
virtual void vsetup(analog_net_t::list_t &nets) = 0;
virtual int vsolve_non_dynamic(const bool newton_raphson) = 0;
virtual nl_double compute_next_timestep() = 0;
pvector_t<analog_net_t *> m_nets;
pvector_t<analog_output_t *> m_inps;
int m_stat_calculations;
int m_stat_newton_raphson;
int m_stat_vsolver_calls;
int m_iterative_fail;
int m_iterative_total;
const solver_parameters_t &m_params;
inline nl_double current_timestep() { return m_cur_ts; }
private:
netlist_time m_last_step;
nl_double m_cur_ts;
dev_list_t m_step_devices;
dev_list_t m_dynamic_devices;
logic_input_t m_fb_sync;
logic_output_t m_Q_sync;
void step(const netlist_time delta);
void update_inputs();
const eSolverType m_type;
};
class NETLIB_NAME(solver) : public device_t class NETLIB_NAME(solver) : public device_t
{ {
@ -216,7 +90,7 @@ protected:
param_logic_t m_log_stats; param_logic_t m_log_stats;
matrix_solver_t::list_t m_mat_solvers; pvector_t<matrix_solver_t *> m_mat_solvers;
private: private:
solver_parameters_t m_params; solver_parameters_t m_params;

View File

@ -422,8 +422,11 @@ NETLIST_START(kidniki_interface)
PARAM(Solver.ITERATIVE, "SOR") PARAM(Solver.ITERATIVE, "SOR")
//PARAM(Solver.ITERATIVE, "MAT") //PARAM(Solver.ITERATIVE, "MAT")
//PARAM(Solver.ITERATIVE, "GMRES") //PARAM(Solver.ITERATIVE, "GMRES")
PARAM(Solver.PARALLEL, 1) PARAM(Solver.PARALLEL, 0)
PARAM(Solver.SOR_FACTOR, 1.00) PARAM(Solver.SOR_FACTOR, 1.00)
PARAM(Solver.DYNAMIC_TS, 0)
PARAM(Solver.DYNAMIC_LTE, 5e-4)
PARAM(Solver.MIN_TIMESTEP, 20e-6)
#else #else
SOLVER(Solver, 12000) SOLVER(Solver, 12000)
PARAM(Solver.ACCURACY, 1e-8) PARAM(Solver.ACCURACY, 1e-8)