diff options
author | 2017-01-31 22:37:12 +0100 | |
---|---|---|
committer | 2017-01-31 22:37:25 +0100 | |
commit | e02d8cad2bc67ada9b7429c6dc8401a439c0c025 (patch) | |
tree | debc396921388c052757f5da5c020c1f122c8c7d | |
parent | 282a7e76c3a2ca3de90351f6527a5aaa8ffdbd4a (diff) |
Fixed a number of issues:
- Fixed crashes on terminals without nets (i.e. connected to a rail)
- Reviewed "FIXMEs" and corrected some minor ones.
- Made m_cur_analog protected.
- Fixed pmf delegates to work with msvc.
- More optimizations to the solver code.
- Started work on a better signal pipeline in nlwav
- Only generate documentation for entities which are documented.
[Couriersud]
26 files changed, 466 insertions, 256 deletions
diff --git a/nl_examples/congo_bongo.c b/nl_examples/congo_bongo.c index c7179847048..492865c79c4 100644 --- a/nl_examples/congo_bongo.c +++ b/nl_examples/congo_bongo.c @@ -11,6 +11,7 @@ /* set to 1 to use optimizations increasing performance significantly */ #define USE_OPTMIZATIONS 1 +#define USE_FRONTIERS 1 /* ---------------------------------------------------------------------------- * Library section header START @@ -44,12 +45,18 @@ NETLIST_START(dummy) PARAM(Solver.ACCURACY, 1e-8) PARAM(Solver.NR_LOOPS, 90) - PARAM(Solver.SOR_FACTOR, 0.001) - PARAM(Solver.GS_LOOPS, 1) + PARAM(Solver.SOR_FACTOR, 1.01) + PARAM(Solver.GS_LOOPS, 4) + //PARAM(Solver.METHOD, "GMRES") PARAM(Solver.METHOD, "MAT_CR") + //PARAM(Solver.METHOD, "SOR") #if USE_OPTMIZATIONS +#if USE_FRONTIERS SOLVER(Solver, 24000) +#else + SOLVER(Solver, 48000) +#endif PARAM(Solver.DYNAMIC_TS, 0 ) PARAM(Solver.PARALLEL, 1) #else @@ -121,7 +128,7 @@ NETLIST_START(dummy) PARAM(XU13.D.MODEL, "MB3614(TYPE=1)") #endif -#if 1 +#if USE_FRONTIERS OPTIMIZE_FRONTIER(C51.1, RES_K(20), 50) OPTIMIZE_FRONTIER(R77.2, RES_K(20), 50) diff --git a/src/lib/netlist/analog/nld_bjt.cpp b/src/lib/netlist/analog/nld_bjt.cpp index 4ce91ba72b3..e40fbfa366d 100644 --- a/src/lib/netlist/analog/nld_bjt.cpp +++ b/src/lib/netlist/analog/nld_bjt.cpp @@ -47,57 +47,55 @@ private: // nld_Q - Base classes // ----------------------------------------------------------------------------- - /* FIXME: Make table pretty */ - /*! Class representing the bjt model paramers. + * * This is the model representation of the bjt model. Typically, SPICE uses * the following parameters. A "Y" in the first column indicates that the * parameter is actually used in netlist. * - * |NL? |name |parameter |units|default| example|area | - * |:--:|:-----|:--------------------------------|:----|------:|-------:|:----:| - * | Y |IS |transport saturation current|A |1E-016|1E-015|* | - * | Y |BF |ideal maximum forward beta|- |100|100| | - * | Y |NF |forward current emission coefficient|- |1|1| | - * | |VAF |forward Early voltage|V |infinite |200| | - * | |IKF |corner for forward beta high current roll-off|A |infinite |0.01|* | - * | |ISE |B-E leakage saturation current|A |0|0.0000000000001|* | - * | |NE |B-E leakage emission coefficient|- |1.5|2| | - * | Y |BR |ideal maximum reverse beta |- |1|0.1| | - * | Y |NR |reverse current emission coefficient|- |1|1| | - * | |VAR |reverse Early voltage|V |infinite |200| | - * | |IKR |corner for reverse beta high current roll-off|A |infinite |0.01|* | - * | |ISC |leakage saturation current|A |0|8| | - * | |NC |leakage emission coefficient|- |2|1.5| | - * | |RB |zero bias base resistance| |0|100|* | - * | |IRB |current where base resistance falls halfway to its min value|A |infinte |0.1|* | - * | |RBM |minimum base resistance at high currents| |RB |10|* | - * | |RE |emitter resistance| |0|1|* | - * | |RC |collector resistance | |0|10|* | - * | |CJE |B-E zero-bias depletion capacitance|F |0|2pF |* | - * | |VJE |B-E built-in potential|V |0.75|0.6| | - * | |MJE |B-E junction exponential factor |- |0.33|0.33| | - * | |TF |ideal forward transit time |sec |0|0.1ns | | - * | |XTF|coefficient for bias dependence of TF |- |0| | | - * | |VTF |voltage describing VBC dependence of TF |V |infinite | | | - * | |ITF |high-current parameter for effect on TF |A |0| |* | - * | |PTF |excess phase at freq=1.0/(TF*2PI) Hz |deg |0| | | - * | |CJC |B-C zero-bias depletion capacitance |F |0|2pF |* | - * | |VJC |B-C built-in potential |V |0.75|0.5| | - * | |MJC |B-C junction exponential factor |- |0.33|0.5| | - * | |XCJC |fraction of B-C depletion capacitance connected to internal base node |- |1| | | - * | |TR |ideal reverse transit time |sec |0|10ns | | - * | |CJS |zero-bias collector-substrate capacitance |F |0|2pF |* | - * | |VJS |substrate junction built-in potential |V |0.75| | | - * | |MJS |substrate junction exponential factor |- |0|0.5| | - * | |XTB |forward and reverse beta temperature exponent |- |0| | | - * | |EG|energy gap for temperature effect on IS |eV |1.11| | | - * | |XTI|temperature exponent for effect on IS |- |3| | | - * | |KF |flicker-noise coefficient |- |0| | | - * | |AF |flicker-noise exponent |- |1| | | - * | |FC |coefficient for forward-bias depletion capacitance formula |- |0.5| | | - * | |TNOM |Parameter measurement temperature |C |27|50| | - */ + * | NL? | name | parameter | units | default | example | area | + * |:---:|------|-----------------------------------------------------------------------|-------|---------:|----------------:|:----:| + * | Y | IS | transport saturation current | A | 1E-016 | 1E-015 | * | + * | Y | BF | ideal maximum forward beta | - | 100 | 100 | | + * | Y | NF | forward current emission coefficient | - | 1 | 1 | | + * | | VAF | forward Early voltage | V | infinite | 200 | | + * | | IKF | corner for forward beta high current roll-off | A | infinite | 0.01 | * | + * | | ISE | B-E leakage saturation current | A | 0 | 0.0000000000001 | * | + * | | NE | B-E leakage emission coefficient | - | 1.5 | 2 | | + * | Y | BR | ideal maximum reverse beta | - | 1 | 0.1 | | + * | Y | NR | reverse current emission coefficient | - | 1 | 1 | | + * | | VAR | reverse Early voltage | V | infinite | 200 | | + * | | IKR | corner for reverse beta high current roll-off | A | infinite | 0.01 | * | + * | | ISC | leakage saturation current | A | 0 | 8 | | + * | | NC | leakage emission coefficient | - | 2 | 1.5 | | + * | | RB | zero bias base resistance | | 0 | 100 | * | + * | | IRB | current where base resistance falls halfway to its min value | A | infinte | 0.1 | * | + * | | RBM | minimum base resistance at high currents | | RB | 10 | * | + * | | RE | emitter resistance | | 0 | 1 | * | + * | | RC | collector resistance | | 0 | 10 | * | + * | | CJE | B-E zero-bias depletion capacitance | F | 0 | 2pF | * | + * | | VJE | B-E built-in potential | V | 0.75 | 0.6 | | + * | | MJE | B-E junction exponential factor | - | 0.33 | 0.33 | | + * | | TF | ideal forward transit time | sec | 0 | 0.1ns | | + * | | XTF | coefficient for bias dependence of TF | - | 0 | | | + * | | VTF | voltage describing VBC dependence of TF | V | infinite | | | + * | | ITF | high-current parameter for effect on TF | A | 0 | | * | + * | | PTF | excess phase at freq=1.0/(TF*2PI) Hz | deg | 0 | | | + * | | CJC | B-C zero-bias depletion capacitance | F | 0 | 2pF | * | + * | | VJC | B-C built-in potential | V | 0.75 | 0.5 | | + * | | MJC | B-C junction exponential factor | - | 0.33 | 0.5 | | + * | | XCJC | fraction of B-C depletion capacitance connected to internal base node | - | 1 | | | + * | | TR | ideal reverse transit time | sec | 0 | 10ns | | + * | | CJS | zero-bias collector-substrate capacitance | F | 0 | 2pF | * | + * | | VJS | substrate junction built-in potential | V | 0.75 | | | + * | | MJS | substrate junction exponential factor | - | 0 | 0.5 | | + * | | XTB | forward and reverse beta temperature exponent | - | 0 | | | + * | | EG | energy gap for temperature effect on IS | eV | 1.11 | | | + * | | XTI | temperature exponent for effect on IS | - | 3 | | | + * | | KF | flicker-noise coefficient | - | 0 | | | + * | | AF | flicker-noise exponent | - | 1 | | | + * | | FC | coefficient for forward-bias depletion capacitance formula | - | 0.5 | | | + * | | TNOM | Parameter measurement temperature | C | 27 | 50 | | */ class bjt_model_t : public param_model_t { diff --git a/src/lib/netlist/analog/nlid_twoterm.cpp b/src/lib/netlist/analog/nlid_twoterm.cpp index e16f68fa706..a9a4edb604a 100644 --- a/src/lib/netlist/analog/nlid_twoterm.cpp +++ b/src/lib/netlist/analog/nlid_twoterm.cpp @@ -27,6 +27,7 @@ generic_diode::generic_diode(device_t &dev, pstring name) , m_Vt(0.0) , m_Vmin(0.0) , m_Is(0.0) + , m_logIs(0.0) , m_n(0.0) , m_gmin(1e-15) , m_VtInv(0.0) @@ -39,6 +40,7 @@ void generic_diode::set_param(const nl_double Is, const nl_double n, nl_double g { static const double csqrt2 = std::sqrt(2.0); m_Is = Is; + m_logIs = std::log(Is); m_n = n; m_gmin = gmin; @@ -62,7 +64,7 @@ void generic_diode::update_diode(const nl_double nVd) m_Vd = nVd; //m_Vd = m_Vd + 10.0 * m_Vt * std::tanh((nVd - m_Vd) / 10.0 / m_Vt); //const double IseVDVt = m_Is * std::exp(m_Vd * m_VtInv); - const double IseVDVt = m_Is * std::exp(m_Vd * m_VtInv); + const double IseVDVt = std::exp(m_logIs + m_Vd * m_VtInv); m_Id = IseVDVt - m_Is; m_G = IseVDVt * m_VtInv + m_gmin; } @@ -70,7 +72,8 @@ void generic_diode::update_diode(const nl_double nVd) { const double a = std::max((nVd - m_Vd) * m_VtInv, NL_FCONST(-0.99)); m_Vd = m_Vd + std::log1p(a) * m_Vt; - const double IseVDVt = m_Is * std::exp(m_Vd * m_VtInv); + //const double IseVDVt = m_Is * std::exp(m_Vd * m_VtInv); + const double IseVDVt = std::exp(m_logIs + m_Vd * m_VtInv); m_Id = IseVDVt - m_Is; m_G = IseVDVt * m_VtInv + m_gmin; } diff --git a/src/lib/netlist/analog/nlid_twoterm.h b/src/lib/netlist/analog/nlid_twoterm.h index 34c62b4c345..6166c707d2d 100644 --- a/src/lib/netlist/analog/nlid_twoterm.h +++ b/src/lib/netlist/analog/nlid_twoterm.h @@ -278,10 +278,10 @@ public: void set_param(const double Is, const double n, double gmin); - inline double I() const { return m_Id; } - inline double G() const { return m_G; } - inline double Ieq() const { return (m_Id - m_Vd * m_G); } - inline double Vd() const { return m_Vd; } + double I() const { return m_Id; } + double G() const { return m_G; } + double Ieq() const { return (m_Id - m_Vd * m_G); } + double Vd() const { return m_Vd; } /* owning object must save those ... */ @@ -293,6 +293,7 @@ private: double m_Vt; double m_Vmin; double m_Is; + double m_logIs; double m_n; double m_gmin; diff --git a/src/lib/netlist/build/doxygen.conf b/src/lib/netlist/build/doxygen.conf index 762fa0477a9..9de9b691114 100644 --- a/src/lib/netlist/build/doxygen.conf +++ b/src/lib/netlist/build/doxygen.conf @@ -467,7 +467,7 @@ EXTRACT_ANON_NSPACES = NO # section is generated. This option has no effect if EXTRACT_ALL is enabled. # The default value is: NO. -HIDE_UNDOC_MEMBERS = NO +HIDE_UNDOC_MEMBERS = YES # If the HIDE_UNDOC_CLASSES tag is set to YES, doxygen will hide all # undocumented classes that are normally visible in the class hierarchy. If set @@ -475,7 +475,7 @@ HIDE_UNDOC_MEMBERS = NO # has no effect if EXTRACT_ALL is enabled. # The default value is: NO. -HIDE_UNDOC_CLASSES = NO +HIDE_UNDOC_CLASSES = YES # If the HIDE_FRIEND_COMPOUNDS tag is set to YES, doxygen will hide all friend # (class|struct|union) declarations. If set to NO, these declarations will be diff --git a/src/lib/netlist/devices/nld_74193.cpp b/src/lib/netlist/devices/nld_74193.cpp index a038f6d97ff..f5370719778 100644 --- a/src/lib/netlist/devices/nld_74193.cpp +++ b/src/lib/netlist/devices/nld_74193.cpp @@ -137,8 +137,8 @@ namespace netlist for (std::size_t i=0; i<4; i++) m_Q[i].push((m_cnt >> i) & 1, delay[i]); - m_BORROWQ.push(tBorrow, NLTIME_FROM_NS(20)); //FIXME - m_CARRYQ.push(tCarry, NLTIME_FROM_NS(20)); //FIXME + m_BORROWQ.push(tBorrow, NLTIME_FROM_NS(20)); //FIXME timing + m_CARRYQ.push(tCarry, NLTIME_FROM_NS(20)); //FIXME timing } NETLIB_DEVICE_IMPL(74193) diff --git a/src/lib/netlist/macro/nlm_ttl74xx.cpp b/src/lib/netlist/macro/nlm_ttl74xx.cpp index 762a451f835..c8fe9912911 100644 --- a/src/lib/netlist/macro/nlm_ttl74xx.cpp +++ b/src/lib/netlist/macro/nlm_ttl74xx.cpp @@ -376,7 +376,7 @@ static NETLIST_START(TTL_7427_DIP) s2.C, /* C2 |5 10| B3 */ s3.B, s2.Q, /* Y2 |6 9| A3 */ s3.A, GND.I, /* GND |7 8| Y3 */ s3.Q - /* +--------------+ */ + /* +--------------+ */ ) NETLIST_END() @@ -416,7 +416,7 @@ static NETLIST_START(TTL_7430_DIP) s1.E, /* E |5 10| NC */ NC.I, s1.F, /* F |6 9| NC */ NC.I, GND.I, /* GND |7 8| Y */ s1.Q - /* +--------------+ */ + /* +--------------+ */ ) NETLIST_END() diff --git a/src/lib/netlist/nl_base.cpp b/src/lib/netlist/nl_base.cpp index 0e53ac66744..056ed77a78b 100644 --- a/src/lib/netlist/nl_base.cpp +++ b/src/lib/netlist/nl_base.cpp @@ -342,9 +342,12 @@ void netlist_t::start() auto p = setup().m_param_values.find(d->name() + ".HINT_NO_DEACTIVATE"); if (p != setup().m_param_values.end()) { - //FIXME: Error checking - auto v = p->second.as_long(); - d->set_hint_deactivate(!v); + //FIXME: turn this into a proper function + bool error; + auto v = p->second.as_double(&error); + if (error || std::abs(v - std::floor(v)) > 1e-6 ) + log().fatal(MF_1_HND_VAL_NOT_SUPPORTED, p->second); + d->set_hint_deactivate(v == 0.0); } } else @@ -436,16 +439,16 @@ void netlist_t::reset() dev->update_param(); // Step all devices once ! + /* + * INFO: The order here affects power up of e.g. breakout. However, such + * variations are explicitly stated in the breakout manual. + */ #if 0 for (std::size_t i = 0; i < m_devices.size(); i++) { m_devices[i]->update_dev(); } #else - /* FIXME: this makes breakout attract mode working again. - * It is however not acceptable that this depends on the startup order. - * Best would be, if reset would call update_dev for devices which need it. - */ std::size_t i = m_devices.size(); while (i>0) m_devices[--i]->update_dev(); @@ -706,10 +709,9 @@ detail::net_t::net_t(netlist_t &nl, const pstring &aname, core_terminal_t *mr) , m_time(*this, "m_time", netlist_time::zero()) , m_active(*this, "m_active", 0) , m_in_queue(*this, "m_in_queue", 2) - , m_railterminal(nullptr) , m_cur_Analog(*this, "m_cur_Analog", 0.0) + , m_railterminal(mr) { - m_railterminal = mr; } detail::net_t::~net_t() @@ -940,16 +942,20 @@ terminal_t::~terminal_t() void terminal_t::schedule_solve() { - // FIXME: Remove this after we found a way to remove *ALL* twoterms connected to railnets only. - if (net().solver() != nullptr) - net().solver()->update_forced(); + // Nets may belong to railnets which do not have a solver attached + // FIXME: Enforce that all terminals get connected? + if (this->has_net()) + if (net().solver() != nullptr) + net().solver()->update_forced(); } void terminal_t::schedule_after(const netlist_time &after) { - // FIXME: Remove this after we found a way to remove *ALL* twoterms connected to railnets only. - if (net().solver() != nullptr) - net().solver()->update_after(after); + // Nets may belong to railnets which do not have a solver attached + // FIXME: Enforce that all terminals get connected? + if (this->has_net()) + if (net().solver() != nullptr) + net().solver()->update_after(after); } // ---------------------------------------------------------------------------------------- @@ -1009,7 +1015,7 @@ analog_output_t::analog_output_t(core_device_t &dev, const pstring &aname) netlist().m_nets.push_back(plib::owned_ptr<analog_net_t>(&m_my_net, false)); this->set_net(&m_my_net); - net().m_cur_Analog = NL_FCONST(0.0); + //net().m_cur_Analog = NL_FCONST(0.0); netlist().setup().register_term(*this); } @@ -1019,7 +1025,7 @@ analog_output_t::~analog_output_t() void analog_output_t::initial(const nl_double val) { - net().m_cur_Analog = val; + net().set_Q_Analog(val); } // ----------------------------------------------------------------------------- diff --git a/src/lib/netlist/nl_base.h b/src/lib/netlist/nl_base.h index ac4a4001dfa..e56a10b9e20 100644 --- a/src/lib/netlist/nl_base.h +++ b/src/lib/netlist/nl_base.h @@ -547,16 +547,12 @@ namespace netlist void set(const nl_double G) { - set_ptr(m_Idr1, 0); - set_ptr(m_go1, G); - set_ptr(m_gt1, G); + set(G,G, 0.0); } void set(const nl_double GO, const nl_double GT) { - set_ptr(m_Idr1, 0); - set_ptr(m_go1, GO); - set_ptr(m_gt1, GT); + set(GO, GT, 0.0); } void set(const nl_double GO, const nl_double GT, const nl_double I) @@ -735,13 +731,13 @@ namespace netlist state_var_s32 m_active; state_var_u8 m_in_queue; /* 0: not in queue, 1: in queue, 2: last was taken */ + state_var<nl_double> m_cur_Analog; + private: plib::linkedlist_t<core_terminal_t> m_list_active; core_terminal_t * m_railterminal; public: - // FIXME: Have to fix the public at some time - state_var<nl_double> m_cur_Analog; }; @@ -796,6 +792,7 @@ namespace netlist virtual ~analog_net_t(); nl_double Q_Analog() const { return m_cur_Analog; } + void set_Q_Analog(const nl_double &v) { m_cur_Analog = v; } nl_double *Q_Analog_state_ptr() { return m_cur_Analog.ptr(); } //FIXME: needed by current solver code @@ -1450,7 +1447,7 @@ namespace netlist { if (newQ != m_my_net.Q_Analog()) { - m_my_net.m_cur_Analog = newQ; + m_my_net.set_Q_Analog(newQ); m_my_net.toggle_new_Q(); m_my_net.push_to_queue(NLTIME_FROM_NS(1)); } diff --git a/src/lib/netlist/nl_config.h b/src/lib/netlist/nl_config.h index 94dd0dfd6b0..ded567717f9 100644 --- a/src/lib/netlist/nl_config.h +++ b/src/lib/netlist/nl_config.h @@ -30,8 +30,12 @@ * */ #define USE_MEMPOOL (0) + #define USE_TRUTHTABLE (1) +// How many times do we try to resolve links (connections) +#define NL_MAX_LINK_RESOLVE_LOOPS (100) + //============================================================ // Solver defines //============================================================ diff --git a/src/lib/netlist/nl_errstr.h b/src/lib/netlist/nl_errstr.h index 83ee4fdc4b9..ec52187558f 100644 --- a/src/lib/netlist/nl_errstr.h +++ b/src/lib/netlist/nl_errstr.h @@ -18,6 +18,7 @@ #define MF_1_UNKNOWN_PARAM_TYPE "Can not determine param_type for {1}" #define MF_2_ERROR_CONNECTING_1_TO_2 "Error connecting {1} to {2}" #define MF_0_NO_SOLVER "No solver found for this netlist although analog elements are present" +#define MF_1_HND_VAL_NOT_SUPPORTED "HINT_NO_DEACTIVATE value not supported: <{1}>" // nl_factory.cpp diff --git a/src/lib/netlist/nl_setup.cpp b/src/lib/netlist/nl_setup.cpp index 271811e95d9..0572ec1e693 100644 --- a/src/lib/netlist/nl_setup.cpp +++ b/src/lib/netlist/nl_setup.cpp @@ -143,7 +143,6 @@ pstring setup_t::termtype_as_str(detail::core_terminal_t &in) const case detail::terminal_type::OUTPUT: return pstring("OUTPUT"); } - // FIXME: noreturn log().fatal(MF_1_UNKNOWN_OBJECT_TYPE_1, static_cast<unsigned>(in.type())); return pstring("Error"); } @@ -376,12 +375,11 @@ param_t *setup_t::find_param(const pstring ¶m_in, bool required) const return (ret == m_params.end() ? nullptr : &ret->second.m_param); } -// FIXME avoid dynamic cast here devices::nld_base_proxy *setup_t::get_d_a_proxy(detail::core_terminal_t &out) { nl_assert(out.is_logic()); - logic_output_t &out_cast = dynamic_cast<logic_output_t &>(out); + logic_output_t &out_cast = static_cast<logic_output_t &>(out); devices::nld_base_proxy *proxy = out_cast.get_proxy(); if (proxy == nullptr) @@ -686,8 +684,8 @@ void setup_t::resolve_inputs() * We therefore first park connecting inputs and retry * after all other terminals were connected. */ - int tries = 100; - while (m_links.size() > 0 && tries > 0) // FIXME: convert into constant + int tries = NL_MAX_LINK_RESOLVE_LOOPS; + while (m_links.size() > 0 && tries > 0) { for (auto li = m_links.begin(); li != m_links.end(); ) diff --git a/src/lib/netlist/nl_setup.h b/src/lib/netlist/nl_setup.h index 077f5e3e673..4ee449a06a8 100644 --- a/src/lib/netlist/nl_setup.h +++ b/src/lib/netlist/nl_setup.h @@ -283,8 +283,6 @@ namespace netlist const logic_family_desc_t *family_from_model(const pstring &model); - /* FIXME: truth table trampoline */ - void tt_factory_create(tt_desc &desc, const pstring &sourcefile); /* helper - also used by nltool */ diff --git a/src/lib/netlist/plib/pconfig.h b/src/lib/netlist/plib/pconfig.h index 105463fe489..e3870206cbd 100644 --- a/src/lib/netlist/plib/pconfig.h +++ b/src/lib/netlist/plib/pconfig.h @@ -56,7 +56,7 @@ typedef __int128_t INT128; //============================================================ // This will be autodetected -// #define PPMF_TYPE 1 +// #define PPMF_TYPE 0 #define PPMF_TYPE_PMF 0 #define PPMF_TYPE_GNUC_PMF_CONV 1 diff --git a/src/lib/netlist/plib/poptions.cpp b/src/lib/netlist/plib/poptions.cpp index 7616c5993d0..775ee6d8810 100644 --- a/src/lib/netlist/plib/poptions.cpp +++ b/src/lib/netlist/plib/poptions.cpp @@ -69,6 +69,13 @@ namespace plib { return (err ? 1 : 0); } + int option_long::parse(const pstring &argument) + { + bool err = false; + m_val = argument.as_long(&err); + return (err ? 1 : 0); + } + int option_vec::parse(const pstring &argument) { bool err = false; diff --git a/src/lib/netlist/plib/poptions.h b/src/lib/netlist/plib/poptions.h index 253fcbf9b34..f859b2f5b3d 100644 --- a/src/lib/netlist/plib/poptions.h +++ b/src/lib/netlist/plib/poptions.h @@ -137,6 +137,20 @@ private: double m_val; }; +class option_long : public option +{ +public: + option_long(options &parent, pstring ashort, pstring along, long defval, pstring help) + : option(parent, ashort, along, help, true), m_val(defval) + {} + + virtual int parse(const pstring &argument) override; + + long operator ()() { return m_val; } +private: + long m_val; +}; + class option_vec : public option { public: diff --git a/src/lib/netlist/plib/ppmf.h b/src/lib/netlist/plib/ppmf.h index 7d7b4e4cb7c..03a2a65baa4 100644 --- a/src/lib/netlist/plib/ppmf.h +++ b/src/lib/netlist/plib/ppmf.h @@ -105,21 +105,27 @@ namespace plib { { public: class generic_class; +#if defined (__INTEL_COMPILER) && defined (_M_X64) // needed for "Intel(R) C++ Intel(R) 64 Compiler XE for applications running on Intel(R) 64, Version 14.0.2.176 Build 20140130" at least + using generic_function = int [((sizeof(void *) + 4 * sizeof(int)) + (sizeof(int) - 1)) / sizeof(int)]; +#elif defined(_MSC_VER)// all other cases - for MSVC maximum size is one pointer, plus 3 ints; all other implementations seem to be smaller + using generic_function = int [((sizeof(void *) + 3 * sizeof(int)) + (sizeof(int) - 1)) / sizeof(int)]; +#else using generic_function = R (generic_class::*)(Targs...); - pmfp_base() : m_func(nullptr) {} +#endif + pmfp_base() {} template<typename MemberFunctionType, typename O> void set_base(MemberFunctionType mftp, O *object) { using function_ptr = R (O::*)(Targs...); function_ptr t = mftp; - m_func = reinterpret_cast<generic_function>(t); + *reinterpret_cast<function_ptr *>(&m_func) = t; } template<typename O> - inline R call(O *obj, Targs... args) const + inline R call(O *obj, Targs... args) { using function_ptr = R (O::*)(Targs...); - function_ptr t = reinterpret_cast<function_ptr>(m_func); + function_ptr t = *reinterpret_cast<function_ptr *>(&m_func); return (obj->*t)(std::forward<Targs>(args)...); } private: @@ -166,13 +172,19 @@ namespace plib { pmfp() : pmfp_base<R, Targs...>(), m_obj(nullptr) {} template<typename MemberFunctionType, typename O> + pmfp(MemberFunctionType mftp, O *object) + { + this->set(mftp, object); + } + + template<typename MemberFunctionType, typename O> void set(MemberFunctionType mftp, O *object) { this->set_base(mftp, object); m_obj = reinterpret_cast<generic_class *>(object); } - inline R operator()(Targs... args) const + inline R operator()(Targs... args) { return this->call(m_obj, std::forward<Targs>(args)...); } diff --git a/src/lib/netlist/prg/nlwav.cpp b/src/lib/netlist/prg/nlwav.cpp index 4f19e9b5960..52f59f366f6 100644 --- a/src/lib/netlist/prg/nlwav.cpp +++ b/src/lib/netlist/prg/nlwav.cpp @@ -5,6 +5,7 @@ #include "plib/plists.h" #include "plib/pstream.h" #include "plib/poptions.h" +#include "plib/ppmf.h" #include "nl_setup.h" class nlwav_options_t : public plib::options @@ -15,6 +16,7 @@ public: opt_inp(*this, "i", "input", "-", "input file"), opt_out(*this, "o", "output", "-", "output file"), opt_amp(*this, "a", "amp", 10000.0, "amplification after mean correction"), + opt_rate(*this, "r", "rate", 48000, "sample rate of output file"), opt_verb(*this, "v", "verbose", "be verbose - this produces lots of output"), opt_quiet(*this,"q", "quiet", "be quiet - no warnings"), opt_version(*this, "", "version", "display version and exit"), @@ -23,6 +25,7 @@ public: plib::option_str opt_inp; plib::option_str opt_out; plib::option_double opt_amp; + plib::option_long opt_rate; plib::option_bool opt_verb; plib::option_bool opt_quiet; plib::option_bool opt_version; @@ -139,12 +142,114 @@ private: }; -static void convert() +class log_processor +{ +public: + typedef plib::pmfp<void, double, double> callback_type; + log_processor(plib::pistream &is, callback_type cb) : m_is(is), m_cb(cb) { } + + void process() + { + plib::putf8_reader reader(m_is); + pstring line; + + while(reader.readline(line)) + { + double t = 0.0; double v = 0.0; + sscanf(line.c_str(), "%lf %lf", &t, &v); + m_cb(t, v); + } + } + +private: + plib::pistream &m_is; + callback_type m_cb; +}; + +struct aggregator +{ + typedef plib::pmfp<void, double, double> callback_type; + + aggregator(double quantum, callback_type cb) + : m_quantum(quantum) + , m_cb(cb) + , ct(0.0) + , lt(0.0) + , outsam(0.0) + , cursam(0.0) + { } + void process(double time, double val) + { + while (time >= ct) + { + outsam += (ct - lt) * cursam; + outsam = outsam / m_quantum; + m_cb(ct, outsam); + outsam = 0.0; + lt = ct; + ct += m_quantum; + } + outsam += (time-lt)*cursam; + lt = time; + cursam = val; + } + +private: + double m_quantum; + callback_type m_cb; + double ct; + double lt; + double outsam; + double cursam; +}; + +class wavwriter +{ +public: + wavwriter(plib::postream &fo, unsigned sample_rate, double ampa) + : mean(0.0) + , means(0.0) + , maxsam(-1e9) + , minsam(1e9) + , n(0) + , m_fo(fo) + , amp(ampa) + , m_wo(m_fo, sample_rate) + { } + + void process(double time, double outsam) + { + means += outsam; + maxsam = std::max(maxsam, outsam); + minsam = std::min(minsam, outsam); + n++; + //mean = means / (double) n; + mean += 5.0 / static_cast<double>(m_wo.sample_rate()) * (outsam - mean); + + outsam = (outsam - mean) * amp; + outsam = std::max(-32000.0, outsam); + outsam = std::min(32000.0, outsam); + m_wo.write_sample(static_cast<int>(outsam)); + } + + double mean; + double means; + double maxsam; + double minsam; + std::size_t n; + +private: + plib::postream &m_fo; + double amp; + wav_t m_wo; +}; + +static void convert(long sample_rate) { plib::postream *fo = (opts.opt_out() == "-" ? &pout_strm : plib::palloc<plib::pofilestream>(opts.opt_out())); plib::pistream *fin = (opts.opt_inp() == "-" ? &pin_strm : plib::palloc<plib::pifilestream>(opts.opt_inp())); plib::putf8_reader reader(*fin); - wav_t *wo = plib::palloc<wav_t>(*fo, 48000U); + wav_t *wo = plib::palloc<wav_t>(*fo, static_cast<unsigned>(sample_rate)); double dt = 1.0 / static_cast<double>(wo->sample_rate()); double ct = dt; @@ -226,10 +331,39 @@ static void convert() } } +static void convert1(long sample_rate) +{ + plib::postream *fo = (opts.opt_out() == "-" ? &pout_strm : plib::palloc<plib::pofilestream>(opts.opt_out())); + plib::pistream *fin = (opts.opt_inp() == "-" ? &pin_strm : plib::palloc<plib::pifilestream>(opts.opt_inp())); + + double dt = 1.0 / static_cast<double>(sample_rate); + + wavwriter *wo = plib::palloc<wavwriter>(*fo, static_cast<unsigned>(sample_rate), opts.opt_amp()); + aggregator ag(dt, aggregator::callback_type(&wavwriter::process, wo)); + log_processor lp(*fin, log_processor::callback_type(&aggregator::process, &ag)); + + lp.process(); + + if (!opts.opt_quiet()) + { + perr("Mean (low freq filter): {}\n", wo->mean); + perr("Mean (static): {}\n", wo->means / static_cast<double>(wo->n)); + perr("Amp + {}\n", 32000.0 / (wo->maxsam - wo->mean)); + perr("Amp - {}\n", -32000.0 / (wo->minsam - wo->mean)); + } + + plib::pfree(wo); + if (opts.opt_inp() != "-") + plib::pfree(fin); + if (opts.opt_out() != "-") + plib::pfree(fo); + +} + static void usage(plib::putf8_fmt_writer &fw) { fw("{}\n", opts.help("Convert netlist log files into wav files.\n", - "nltool [options]").c_str()); + "nltool [options]")); } @@ -262,7 +396,10 @@ int main(int argc, char *argv[]) return 0; } - convert(); + if ((1)) + convert1(opts.opt_rate()); + else + convert(opts.opt_rate()); return 0; } diff --git a/src/lib/netlist/solver/nld_matrix_solver.cpp b/src/lib/netlist/solver/nld_matrix_solver.cpp index 26ab54a7901..f314aea2314 100644 --- a/src/lib/netlist/solver/nld_matrix_solver.cpp +++ b/src/lib/netlist/solver/nld_matrix_solver.cpp @@ -66,7 +66,7 @@ void terms_for_net_t::set_pointers() for (unsigned i = 0; i < count(); i++) { m_terms[i]->set_ptrs(&m_gt[i], &m_go[i], &m_Idr[i]); - m_connected_net_V[i] = m_terms[i]->m_otherterm->net().m_cur_Analog.ptr(); + m_connected_net_V[i] = m_terms[i]->m_otherterm->net().Q_Analog_state_ptr(); } } @@ -155,7 +155,7 @@ void matrix_solver_t::setup_base(analog_net_t::list_t &nets) net_proxy_output->m_proxied_net = static_cast<analog_net_t *>(&p->net()); } net_proxy_output->net().add_terminal(*p); - // FIXME: repeated + // FIXME: repeated calling - kind of brute force net_proxy_output->net().rebuild_list(); log().debug("Added input\n"); } @@ -208,7 +208,7 @@ void matrix_solver_t::setup_matrix() * literature but I have found no articles about Gauss Seidel. * * For Gaussian Elimination however increasing order is better suited. - * FIXME: Even better would be to sort on elements right of the matrix diagonal. + * NOTE: Even better would be to sort on elements right of the matrix diagonal. * */ @@ -489,10 +489,6 @@ netlist_time matrix_solver_t::compute_next_timestep(const double cur_ts) if (m_params.m_dynamic_ts) { - /* - * FIXME: We should extend the logic to use either all nets or - * only output nets. - */ for (std::size_t k = 0, iN=m_terms.size(); k < iN; k++) { analog_net_t *n = m_nets[k]; diff --git a/src/lib/netlist/solver/nld_matrix_solver.h b/src/lib/netlist/solver/nld_matrix_solver.h index 3ed2e9f2c31..82910e9c4d1 100644 --- a/src/lib/netlist/solver/nld_matrix_solver.h +++ b/src/lib/netlist/solver/nld_matrix_solver.h @@ -198,7 +198,7 @@ private: template <typename T> T matrix_solver_t::delta(const T * RESTRICT V) { - /* FIXME: Ideally we should also include currents (RHS) here. This would + /* NOTE: Ideally we should also include currents (RHS) here. This would * need a reevaluation of the right hand side after voltages have been updated * and thus belong into a different calculation. This applies to all solvers. */ @@ -206,7 +206,7 @@ T matrix_solver_t::delta(const T * RESTRICT V) const std::size_t iN = this->m_terms.size(); T cerr = 0; for (std::size_t i = 0; i < iN; i++) - cerr = std::max(cerr, std::abs(V[i] - static_cast<T>(this->m_nets[i]->m_cur_Analog))); + cerr = std::max(cerr, std::abs(V[i] - static_cast<T>(this->m_nets[i]->Q_Analog()))); return cerr; } @@ -215,7 +215,7 @@ void matrix_solver_t::store(const T * RESTRICT V) { const std::size_t iN = this->m_terms.size(); for (std::size_t i = 0; i < iN; i++) - this->m_nets[i]->m_cur_Analog = V[i]; + this->m_nets[i]->set_Q_Analog(V[i]); } template <typename T> diff --git a/src/lib/netlist/solver/nld_ms_direct.h b/src/lib/netlist/solver/nld_ms_direct.h index b188978d2ee..b929fb6ed22 100644 --- a/src/lib/netlist/solver/nld_ms_direct.h +++ b/src/lib/netlist/solver/nld_ms_direct.h @@ -157,9 +157,9 @@ protected: inline nl_ext_double &RHS(const T1 &r) { return m_A[r * m_pitch + N()]; } #else template <typename T1, typename T2> - inline nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r][c]; } + nl_ext_double &A(const T1 &r, const T2 &c) { return m_A[r][c]; } template <typename T1> - inline nl_ext_double &RHS(const T1 &r) { return m_A[r][N()]; } + nl_ext_double &RHS(const T1 &r) { return m_A[r][N()]; } #endif nl_double m_last_RHS[storage_N]; // right hand side - contains currents @@ -245,10 +245,54 @@ template <std::size_t m_N, std::size_t storage_N> void matrix_solver_direct_t<m_N, storage_N>::LE_solve() { const std::size_t kN = N(); + if (!(!TEST_PARALLEL && m_params.m_pivot)) + { + for (std::size_t i = 0; i < kN; i++) + { +#if TEST_PARALLEL + const unsigned eb = m_terms[i]->m_nzbd.size(); + if (eb > 0) + { + //printf("here %d\n", eb); + unsigned chunks = (eb + num_thr) / (num_thr + 1); + for (int p=0; p < num_thr + 1; p++) + { + x_i[p] = i; + x_start[p] = chunks * p; + x_stop[p] = std::min(chunks*(p+1), eb); + if (p<num_thr && x_start[p] < x_stop[p]) thr_process(p, this, nullptr); + } + if (x_start[num_thr] < x_stop[num_thr]) + do_work(num_thr, nullptr); + thr_wait(); + } + else if (eb > 0) + { + x_i[0] = i; + x_start[0] = 0; + x_stop[0] = eb; + do_work(0, nullptr); + } +#else - for (std::size_t i = 0; i < kN; i++) { - // FIXME: use a parameter to enable pivoting? m_pivot - if (!TEST_PARALLEL && m_params.m_pivot) + /* FIXME: Singular matrix? */ + const nl_double f = 1.0 / A(i,i); + const auto &nzrd = m_terms[i]->m_nzrd; + const auto &nzbd = m_terms[i]->m_nzbd; + + for (std::size_t j : nzbd) + { + const nl_double f1 = -f * A(j,i); + for (std::size_t k : nzrd) + A(j,k) += A(i,k) * f1; + //RHS(j) += RHS(i) * f1; + } +#endif + } + } + else + { + for (std::size_t i = 0; i < kN; i++) { /* Find the row with the largest first value */ std::size_t maxrow = i; @@ -292,48 +336,6 @@ void matrix_solver_direct_t<m_N, storage_N>::LE_solve() } } } - else - { -#if TEST_PARALLEL - const unsigned eb = m_terms[i]->m_nzbd.size(); - if (eb > 0) - { - //printf("here %d\n", eb); - unsigned chunks = (eb + num_thr) / (num_thr + 1); - for (int p=0; p < num_thr + 1; p++) - { - x_i[p] = i; - x_start[p] = chunks * p; - x_stop[p] = std::min(chunks*(p+1), eb); - if (p<num_thr && x_start[p] < x_stop[p]) thr_process(p, this, nullptr); - } - if (x_start[num_thr] < x_stop[num_thr]) - do_work(num_thr, nullptr); - thr_wait(); - } - else if (eb > 0) - { - x_i[0] = i; - x_start[0] = 0; - x_stop[0] = eb; - do_work(0, nullptr); - } -#else - - /* FIXME: Singular matrix? */ - const nl_double f = 1.0 / A(i,i); - const auto &nzrd = m_terms[i]->m_nzrd; - const auto &nzbd = m_terms[i]->m_nzbd; - - for (std::size_t j : nzbd) - { - const nl_double f1 = -f * A(j,i); - for (std::size_t k : nzrd) - A(j,k) += A(i,k) * f1; - //RHS(j) += RHS(i) * f1; - } -#endif - } } } diff --git a/src/lib/netlist/solver/nld_ms_gcr.h b/src/lib/netlist/solver/nld_ms_gcr.h index d524447b976..eda9b14aa74 100644 --- a/src/lib/netlist/solver/nld_ms_gcr.h +++ b/src/lib/netlist/solver/nld_ms_gcr.h @@ -5,8 +5,6 @@ * * Gaussian elimination using compressed row format. * - * Fow w==1 we will do the classic Gauss-Seidel approach - * */ #ifndef NLD_MS_GCR_H_ @@ -21,11 +19,6 @@ #include "solver/vector_base.h" #include "plib/pstream.h" -#define NL_USE_SSE 0 -#if NL_USE_SSE -#include <mmintrin.h> -#endif - namespace netlist { namespace devices @@ -62,7 +55,7 @@ private: void csc_private(plib::putf8_fmt_writer &strm); - using extsolver = void (*)(double * RESTRICT m_A, double * RESTRICT RHS); + using extsolver = void (*)(double * RESTRICT m_A, double * RESTRICT RHS, double * RESTRICT V); pstring static_compile_name(); @@ -75,7 +68,7 @@ private: }; // ---------------------------------------------------------------------------------------- -// matrix_solver - GMRES +// matrix_solver - GCR // ---------------------------------------------------------------------------------------- template <std::size_t m_N, std::size_t storage_N> @@ -168,7 +161,7 @@ void matrix_solver_GCR_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets) } } - +#if 0 template <std::size_t m_N, std::size_t storage_N> void matrix_solver_GCR_t<m_N, storage_N>::csc_private(plib::putf8_fmt_writer &strm) { @@ -213,8 +206,84 @@ void matrix_solver_GCR_t<m_N, storage_N>::csc_private(plib::putf8_fmt_writer &st } } } + + //new_V[iN - 1] = RHS[iN - 1] / mat.A[mat.diag[iN - 1]]; + strm("\tV[{1}] = RHS[{2}] / m_A[{3}];\n", iN - 1, iN - 1, mat.diag[iN - 1]); + for (std::size_t j = iN - 1; j-- > 0;) + { + strm("\tdouble tmp{1} = 0.0;\n", j); + const std::size_t e = mat.ia[j+1]; + for (std::size_t pk = mat.diag[j] + 1; pk < e; pk++) + { + strm("\ttmp{1} += m_A[{2}] * V[{3}];\n", j, pk, mat.ja[pk]); + } + strm("\tV[{1}] = (RHS[{1}] - tmp{1}) / m_A[{4}];\n", j, j, j, mat.diag[j]); + } } +#else +template <std::size_t m_N, std::size_t storage_N> +void matrix_solver_GCR_t<m_N, storage_N>::csc_private(plib::putf8_fmt_writer &strm) +{ + const std::size_t iN = N(); + for (std::size_t i = 0; i < mat.nz_num; i++) + strm("double m_A{1} = m_A[{2}];\n", i, i); + + for (std::size_t i = 0; i < iN - 1; i++) + { + const auto &nzbd = this->m_terms[i]->m_nzbd; + + if (nzbd.size() > 0) + { + std::size_t pi = mat.diag[i]; + + //const nl_double f = 1.0 / m_A[pi++]; + strm("const double f{1} = 1.0 / m_A{2};\n", i, pi); + pi++; + const std::size_t piie = mat.ia[i+1]; + + //for (auto & j : nzbd) + for (std::size_t j : nzbd) + { + // proceed to column i + std::size_t pj = mat.ia[j]; + + while (mat.ja[pj] < i) + pj++; + + //const nl_double f1 = - m_A[pj++] * f; + strm("\tconst double f{1}_{2} = -f{3} * m_A{4};\n", i, j, i, pj); + pj++; + + // subtract row i from j */ + for (std::size_t pii = pi; pii<piie; ) + { + while (mat.ja[pj] < mat.ja[pii]) + pj++; + //m_A[pj++] += m_A[pii++] * f1; + strm("\tm_A{1} += m_A{2} * f{3}_{4};\n", pj, pii, i, j); + pj++; pii++; + } + //RHS[j] += f1 * RHS[i]; + strm("\tRHS[{1}] += f{2}_{3} * RHS[{4}];\n", j, i, j, i); + } + } + } + + //new_V[iN - 1] = RHS[iN - 1] / mat.A[mat.diag[iN - 1]]; + strm("\tV[{1}] = RHS[{2}] / m_A{3};\n", iN - 1, iN - 1, mat.diag[iN - 1]); + for (std::size_t j = iN - 1; j-- > 0;) + { + strm("\tdouble tmp{1} = 0.0;\n", j); + const std::size_t e = mat.ia[j+1]; + for (std::size_t pk = mat.diag[j] + 1; pk < e; pk++) + { + strm("\ttmp{1} += m_A{2} * V[{3}];\n", j, pk, mat.ja[pk]); + } + strm("\tV[{1}] = (RHS[{1}] - tmp{1}) / m_A{4};\n", j, j, j, mat.diag[j]); + } +} +#endif template <std::size_t m_N, std::size_t storage_N> pstring matrix_solver_GCR_t<m_N, storage_N>::static_compile_name() @@ -234,7 +303,7 @@ std::pair<pstring, pstring> matrix_solver_GCR_t<m_N, storage_N>::create_solver_c plib::putf8_fmt_writer strm(t); pstring name = static_compile_name(); - strm.writeline(plib::pfmt("extern \"C\" void {1}(double * __restrict m_A, double * __restrict RHS)\n")(name)); + strm.writeline(plib::pfmt("extern \"C\" void {1}(double * __restrict m_A, double * __restrict RHS, double * __restrict V)\n")(name)); strm.writeline("{\n"); csc_private(strm); strm.writeline("}\n"); @@ -266,30 +335,13 @@ unsigned matrix_solver_GCR_t<m_N, storage_N>::vsolve_non_dynamic(const bool newt const nl_double * const * RESTRICT other_cur_analog = t->connected_net_V(); const unsigned * const RESTRICT tcr = m_term_cr[k].data(); -#if 1 -#if (0 ||NL_USE_SSE) - __m128d mg = _mm_set_pd(0.0, 0.0); - __m128d mr = _mm_set_pd(0.0, 0.0); - unsigned i = 0; - for (; i < term_count - 1; i+=2) - { - mg = _mm_add_pd(mg, _mm_loadu_pd(>[i])); - mr = _mm_add_pd(mr, _mm_loadu_pd(&Idr[i])); - } - gtot_t = _mm_cvtsd_f64(mg) + _mm_cvtsd_f64(_mm_unpackhi_pd(mg,mg)); - RHS_t = _mm_cvtsd_f64(mr) + _mm_cvtsd_f64(_mm_unpackhi_pd(mr,mr)); - for (; i < term_count; i++) - { - gtot_t += gt[i]; - RHS_t += Idr[i]; - } -#else +#if 0 for (std::size_t i = 0; i < term_count; i++) { gtot_t += gt[i]; RHS_t += Idr[i]; } -#endif + for (std::size_t i = railstart; i < term_count; i++) RHS_t += go[i] * *other_cur_analog[i]; @@ -304,7 +356,7 @@ unsigned matrix_solver_GCR_t<m_N, storage_N>::vsolve_non_dynamic(const bool newt #else for (std::size_t i = 0; i < railstart; i++) { - m_A[tcr[i]] -= go[i]; + mat.A[tcr[i]] -= go[i]; gtot_t = gtot_t + gt[i]; RHS_t = RHS_t + Idr[i]; } @@ -316,7 +368,7 @@ unsigned matrix_solver_GCR_t<m_N, storage_N>::vsolve_non_dynamic(const bool newt } RHS[k] = RHS_t; - m_A[mat.diag[k]] += gtot_t; + mat.A[mat.diag[k]] += gtot_t; } #endif mat.ia[iN] = static_cast<mattype>(mat.nz_num); @@ -326,7 +378,7 @@ unsigned matrix_solver_GCR_t<m_N, storage_N>::vsolve_non_dynamic(const bool newt if (m_proc != nullptr) { //static_solver(m_A, RHS); - m_proc(mat.A, RHS); + m_proc(mat.A, RHS, new_V); } else { @@ -362,45 +414,28 @@ unsigned matrix_solver_GCR_t<m_N, storage_N>::vsolve_non_dynamic(const bool newt } } } - } - - /* backward substitution - * - */ + /* backward substitution + * + */ - /* row n-1 */ - new_V[iN - 1] = RHS[iN - 1] / mat.A[mat.diag[iN - 1]]; + /* row n-1 */ + new_V[iN - 1] = RHS[iN - 1] / mat.A[mat.diag[iN - 1]]; - for (std::size_t j = iN - 1; j-- > 0;) - { - //__builtin_prefetch(&new_V[j-1], 1); - //if (j>0)__builtin_prefetch(&m_A[mat.diag[j-1]], 0); -#if (NL_USE_SSE) - __m128d tmp = _mm_set_pd1(0.0); - const unsigned e = mat.ia[j+1]; - unsigned pk = mat.diag[j] + 1; - for (; pk < e - 1; pk+=2) + for (std::size_t j = iN - 1; j-- > 0;) { - //tmp += m_A[pk] * new_V[mat.ja[pk]]; - tmp = _mm_add_pd(tmp, _mm_mul_pd(_mm_set_pd(mat.A[pk], mat.A[pk+1]), - _mm_set_pd(new_V[mat.ja[pk]], new_V[mat.ja[pk+1]]))); - } - double tmpx = _mm_cvtsd_f64(tmp) + _mm_cvtsd_f64(_mm_unpackhi_pd(tmp,tmp)); - for (; pk < e; pk++) - { - tmpx += mat.A[pk] * new_V[mat.ja[pk]]; - } - new_V[j] = (RHS[j] - tmpx) / mat.A[mat.diag[j]]; -#else - double tmp = 0; - const std::size_t e = mat.ia[j+1]; - for (std::size_t pk = mat.diag[j] + 1; pk < e; pk++) - { - tmp += mat.A[pk] * new_V[mat.ja[pk]]; + //__builtin_prefetch(&new_V[j-1], 1); + //if (j>0)__builtin_prefetch(&m_A[mat.diag[j-1]], 0); + double tmp = 0; + auto jdiag = mat.diag[j]; + const std::size_t e = mat.ia[j+1]; + for (std::size_t pk = jdiag + 1; pk < e; pk++) + { + tmp += mat.A[pk] * new_V[mat.ja[pk]]; + } + new_V[j] = (RHS[j] - tmp) / mat.A[jdiag]; } - new_V[j] = (RHS[j] - tmp) / mat.A[mat.diag[j]]; -#endif } + this->m_stat_calculations++; if (newton_raphson) diff --git a/src/lib/netlist/solver/nld_ms_gmres.h b/src/lib/netlist/solver/nld_ms_gmres.h index 7ae04082ff8..a97de6a2446 100644 --- a/src/lib/netlist/solver/nld_ms_gmres.h +++ b/src/lib/netlist/solver/nld_ms_gmres.h @@ -144,7 +144,7 @@ unsigned matrix_solver_GMRES_t<m_N, storage_N>::vsolve_non_dynamic(const bool ne const nl_double * const RESTRICT Idr = this->m_terms[k]->Idr(); const nl_double * const * RESTRICT other_cur_analog = this->m_terms[k]->connected_net_V(); - new_V[k] = this->m_nets[k]->m_cur_Analog; + new_V[k] = this->m_nets[k]->Q_Analog(); for (std::size_t i = 0; i < term_count; i++) { diff --git a/src/lib/netlist/solver/nld_ms_sor.h b/src/lib/netlist/solver/nld_ms_sor.h index f07e7057795..dc3a5d1821c 100644 --- a/src/lib/netlist/solver/nld_ms_sor.h +++ b/src/lib/netlist/solver/nld_ms_sor.h @@ -74,7 +74,7 @@ unsigned matrix_solver_SOR_t<m_N, storage_N>::vsolve_non_dynamic(const bool newt nl_double RHS[storage_N]; nl_double new_V[storage_N]; - for (unsigned k = 0; k < iN; k++) + for (std::size_t k = 0; k < iN; k++) { nl_double gtot_t = 0.0; nl_double gabs_t = 0.0; @@ -86,7 +86,7 @@ unsigned matrix_solver_SOR_t<m_N, storage_N>::vsolve_non_dynamic(const bool newt const nl_double * const RESTRICT Idr = this->m_terms[k]->Idr(); const nl_double * const *other_cur_analog = this->m_terms[k]->connected_net_V(); - new_V[k] = this->m_nets[k]->m_cur_Analog; + new_V[k] = this->m_nets[k]->Q_Analog(); for (std::size_t i = 0; i < term_count; i++) { @@ -125,17 +125,10 @@ unsigned matrix_solver_SOR_t<m_N, storage_N>::vsolve_non_dynamic(const bool newt const nl_double accuracy = this->m_params.m_accuracy; - /* uncommenting the line below will force dynamic updates every X iterations - * althought the system has not converged yet. This is a proof of concept, - * - */ - const bool interleaved_dynamic_updates = false; - //const bool interleaved_dynamic_updates = newton_raphson; - do { resched = false; nl_double err = 0; - for (unsigned k = 0; k < iN; k++) + for (std::size_t k = 0; k < iN; k++) { const int * RESTRICT net_other = this->m_terms[k]->connected_net_idx(); const std::size_t railstart = this->m_terms[k]->m_railstart; @@ -156,29 +149,20 @@ unsigned matrix_solver_SOR_t<m_N, storage_N>::vsolve_non_dynamic(const bool newt resched_cnt++; //} while (resched && (resched_cnt < this->m_params.m_gs_loops)); - } while (resched && ((!interleaved_dynamic_updates && resched_cnt < this->m_params.m_gs_loops) || (interleaved_dynamic_updates && resched_cnt < 5 ))); + } while (resched && ((resched_cnt < this->m_params.m_gs_loops))); this->m_iterative_total += resched_cnt; + this->m_stat_calculations++; - if (resched && !interleaved_dynamic_updates) + if (resched) { // Fallback to direct solver ... this->m_iterative_fail++; return matrix_solver_direct_t<m_N, storage_N>::vsolve_non_dynamic(newton_raphson); } - this->m_stat_calculations++; - - if (interleaved_dynamic_updates) - { - for (unsigned k = 0; k < iN; k++) - this->m_nets[k]->m_cur_Analog += 1.0 * (new_V[k] - this->m_nets[k]->m_cur_Analog); - } - else - { - for (unsigned k = 0; k < iN; k++) - this->m_nets[k]->m_cur_Analog = new_V[k]; - } + for (std::size_t k = 0; k < iN; k++) + this->m_nets[k]->set_Q_Analog(new_V[k]); return resched_cnt; } diff --git a/src/lib/netlist/solver/nld_ms_sor_mat.h b/src/lib/netlist/solver/nld_ms_sor_mat.h index e9204029abc..265ef821957 100644 --- a/src/lib/netlist/solver/nld_ms_sor_mat.h +++ b/src/lib/netlist/solver/nld_ms_sor_mat.h @@ -170,7 +170,7 @@ unsigned matrix_solver_SOR_mat_t<m_N, storage_N>::vsolve_non_dynamic(const bool #endif for (std::size_t k = 0; k < iN; k++) - new_v[k] = this->m_nets[k]->m_cur_Analog; + new_v[k] = this->m_nets[k]->Q_Analog(); do { resched = false; @@ -199,10 +199,12 @@ unsigned matrix_solver_SOR_mat_t<m_N, storage_N>::vsolve_non_dynamic(const bool } while (resched && (resched_cnt < this->m_params.m_gs_loops)); this->m_stat_calculations++; + this->m_iterative_total += resched_cnt; this->m_gs_total += resched_cnt; if (resched) { + this->m_iterative_fail++; //this->netlist().warning("Falling back to direct solver .. Consider increasing RESCHED_LOOPS"); this->m_gs_fail++; diff --git a/src/lib/netlist/solver/nld_solver.cpp b/src/lib/netlist/solver/nld_solver.cpp index a388fd19d7b..d2ecf3c5bd2 100644 --- a/src/lib/netlist/solver/nld_solver.cpp +++ b/src/lib/netlist/solver/nld_solver.cpp @@ -159,8 +159,16 @@ std::unique_ptr<matrix_solver_t> NETLIB_NAME(solver)::create_solver(std::size_t } else if (pstring("MAT_CR").equals(m_method())) { - typedef matrix_solver_GCR_t<m_N,storage_N> solver_mat; - return plib::make_unique<solver_mat>(netlist(), solvername, &m_params, size); + if (size > 0) // GCR always outperforms MAT solver + { + typedef matrix_solver_GCR_t<m_N,storage_N> solver_mat; + return plib::make_unique<solver_mat>(netlist(), solvername, &m_params, size); + } + else + { + typedef matrix_solver_direct_t<m_N,storage_N> solver_mat; + return plib::make_unique<solver_mat>(netlist(), solvername, &m_params, size); + } } else if (pstring("MAT").equals(m_method())) { |