Added two more models to netlist:

OPAMP: Generic opamp model. This does all the annoying calculations.
Just pass the the datasheet values.

LVCCS: A limited current voltage control current source. This will allow
slew rate limiting going forward. 

In addition:

- add a (small) parallel conductance to all capacitors to improve
convergence.
- some initial work to use "long double".
This commit is contained in:
couriersud 2015-07-11 16:09:08 +02:00
parent 82735c0e36
commit 3a8d682827
15 changed files with 406 additions and 94 deletions

View File

@ -49,9 +49,9 @@ NETLIST_START(dummy)
// .END
SOLVER(Solver, 24000)
PARAM(Solver.ACCURACY, 1e-8)
PARAM(Solver.ACCURACY, 1e-7)
PARAM(Solver.NR_LOOPS, 90)
PARAM(Solver.SOR_FACTOR, 0.00001)
PARAM(Solver.SOR_FACTOR, 0.001)
PARAM(Solver.GS_LOOPS, 1)
//PARAM(Solver.GS_THRESHOLD, 99)
PARAM(Solver.ITERATIVE, "SOR")
@ -334,7 +334,7 @@ NETLIST_START(CongoBongo_schematics)
NET_C(C61.1, R94.2)
NETLIST_END()
NETLIST_START(opamp)
NETLIST_START(opamp_mod)
/* Opamp model from
*
@ -374,13 +374,18 @@ NETLIST_START(opamp)
NET_C(EBUF.ON, fVREF)
/* The opamp model */
VCCS(G1)
LVCCS(G1)
PARAM(G1.RI, RES_K(1000))
#if 0
PARAM(G1.G, 0.0022)
RES(RP1, 1e6)
CAP(CP1, 0.0318e-6)
#else
PARAM(G1.G, 0.0002)
PARAM(G1.CURLIM, 0.002)
RES(RP1, 9.5e6)
CAP(CP1, 0.0033e-6)
#endif
VCVS(EBUF)
PARAM(EBUF.RO, 50)
PARAM(EBUF.G, 1)
@ -395,21 +400,43 @@ NETLIST_START(opamp)
DIODE(DP,".model tt D(IS=1e-15 N=1)")
DIODE(DN,".model tt D(IS=1e-15 N=1)")
#if 1
NET_C(DP.K, fUH.Q)
NET_C(fUL.Q, DN.A)
NET_C(DP.A, DN.K, RP1.1)
#else
/*
* This doesn't add any performance by decreasing iteration loops.
* To the contrary, it significantly decreases iterations
*/
RES(RH1, 0.1)
RES(RL1, 0.1)
NET_C(DP.K, RH1.1)
NET_C(RH1.2, fUH.Q)
NET_C(fUL.Q, RL1.1)
NET_C(RL1.2, DN.A)
NET_C(DP.A, DN.K, RP1.1)
#endif
NET_C(EBUF.IP, RP1.1)
NETLIST_END()
NETLIST_START(MB3614_DIP)
SUBMODEL(opamp, op1)
SUBMODEL(opamp, op2)
SUBMODEL(opamp, op3)
SUBMODEL(opamp, op4)
#if 1
SUBMODEL(opamp_mod, op1)
SUBMODEL(opamp_mod, op2)
SUBMODEL(opamp_mod, op3)
SUBMODEL(opamp_mod, op4)
#else
/* The opamp actually has an FPF of about 500k. This doesn't work here and causes oscillations.
* FPF here therefore about half the Solver clock.
*/
OPAMP(op1, ".model MB3614 OPAMP(TYPE=3 VLH=2.0 VLL=0.2 FPF=5 UGF=11k SLEW=0.6M RI=1000k RO=50 DAB=0.002)")
OPAMP(op2, ".model MB3614 OPAMP(TYPE=3 VLH=2.0 VLL=0.2 FPF=5 UGF=11k SLEW=0.6M RI=1000k RO=50 DAB=0.002)")
OPAMP(op3, ".model MB3614 OPAMP(TYPE=3 VLH=2.0 VLL=0.2 FPF=5 UGF=11k SLEW=0.6M RI=1000k RO=50 DAB=0.002)")
OPAMP(op4, ".model MB3614 OPAMP(TYPE=3 VLH=2.0 VLL=0.2 FPF=5 UGF=11k SLEW=0.6M RI=1000k RO=50 DAB=0.002)")
#endif
ALIAS( 1, op1.OUT)
ALIAS( 2, op1.MINUS)
ALIAS( 3, op1.PLUS)

View File

@ -11,10 +11,12 @@ NETLIST_START(main)
/* Standard stuff */
CLOCK(clk, 1000) // 1000 Hz
SOLVER(Solver, 48000)
//PARAM(Solver.ACCURACY, 1e-3)
SOLVER(Solver, 480000)
PARAM(Solver.ACCURACY, 1e-10)
PARAM(Solver.NR_LOOPS, 30000 )
//PARAM(Solver.CONVERG, 1.0)
//PARAM(Solver.GS_LOOPS, 30)
PARAM(Solver.GS_LOOPS, 30)
PARAM(Solver.GS_THRESHOLD, 99)
// Tie up +5 to opamps thought it's not currently needed
// Stay compatible
@ -47,8 +49,8 @@ NETLIST_START(main)
NET_C(RL.2, GND)
NET_C(RL.1, op1.OUT)
//LOG(logX, op1.OUT)
//LOG(logY, clk)
LOG(logX, op1.OUT)
LOG(logY, clk)
NETLIST_END()
NETLIST_START(opamp)
@ -73,10 +75,11 @@ NETLIST_START(opamp)
/* The opamp model */
VCCS(G1)
PARAM(G1.G, 100) // typical OP-AMP amplification 100 * 1000 = 100000
RES(RP1, 1000)
CAP(CP1, 1.59e-6) // <== change to 1.59e-3 for 10Khz bandwidth
LVCCS(G1)
PARAM(G1.G, 0.0021)
PARAM(G1.CURLIM, 0.002)
RES(RP1, 1e7)
CAP(CP1, 0.00333e-6)
VCVS(EBUF)
PARAM(EBUF.RO, 50)
PARAM(EBUF.G, 1)

View File

@ -65,6 +65,7 @@ NETLIB_RESET(VCCS)
NETLIB_UPDATE_PARAM(VCCS)
{
NETLIB_NAME(VCCS)::reset();
}
NETLIB_UPDATE(VCCS)
@ -80,6 +81,56 @@ NETLIB_UPDATE(VCCS)
m_ON.schedule_solve();
}
// ----------------------------------------------------------------------------------------
// nld_LVCCS
// ----------------------------------------------------------------------------------------
NETLIB_START(LVCCS)
{
NETLIB_NAME(VCCS)::start();
register_param("CURLIM", m_cur_limit, 1000.0);
}
NETLIB_RESET(LVCCS)
{
NETLIB_NAME(VCCS)::reset();
}
NETLIB_UPDATE_PARAM(LVCCS)
{
NETLIB_NAME(VCCS)::update_param();
}
NETLIB_UPDATE(LVCCS)
{
NETLIB_NAME(VCCS)::update();
}
NETLIB_UPDATE_TERMINALS(LVCCS)
{
const nl_double m_mult = m_G.Value() * m_gfac; // 1.0 ==> 1V ==> 1A
const nl_double vi = m_IP.net().m_cur_Analog - m_IN.net().m_cur_Analog;
if (std::abs(m_mult / m_cur_limit * vi) > 0.5)
m_vi = m_vi + 0.2*std::tanh((vi - m_vi)/0.2);
else
m_vi = vi;
const nl_double x = m_mult / m_cur_limit * m_vi;
const nl_double X = std::tanh(x);
const nl_double beta = m_mult * (1.0 - X*X);
const nl_double I = m_cur_limit * X - beta * m_vi;
m_OP.set(beta, NL_FCONST(0.0), I);
m_OP1.set(-beta, NL_FCONST(0.0));
m_ON.set(-beta, NL_FCONST(0.0), -I);
m_ON1.set(beta, NL_FCONST(0.0));
//printf("vi %f beta %f I %f\n", vi, beta, I);
}
// ----------------------------------------------------------------------------------------
// nld_CCCS
// ----------------------------------------------------------------------------------------

View File

@ -16,15 +16,17 @@
// Macros
// ----------------------------------------------------------------------------------------
#define VCCS(_name) \
#define VCCS(_name) \
NET_REGISTER_DEV(VCCS, _name)
#define CCCS(_name) \
#define CCCS(_name) \
NET_REGISTER_DEV(CCCS, _name)
#define VCVS(_name) \
#define VCVS(_name) \
NET_REGISTER_DEV(VCVS, _name)
#define LVCCS(_name) \
NET_REGISTER_DEV(LVCCS, _name)
NETLIB_NAMESPACE_DEVICES_START()
@ -57,11 +59,14 @@ public:
ATTR_COLD NETLIB_NAME(VCCS)(const family_t afamily)
: device_t(afamily), m_gfac(1.0) { }
param_double_t m_G;
param_double_t m_RI;
protected:
virtual void start();
virtual void reset();
virtual void update_param();
ATTR_HOT void update();
ATTR_HOT virtual void update();
ATTR_COLD void start_internal(const nl_double def_RI);
@ -74,12 +79,31 @@ protected:
terminal_t m_OP1;
terminal_t m_ON1;
param_double_t m_G;
param_double_t m_RI;
nl_double m_gfac;
};
/* Limited Current source*/
class NETLIB_NAME(LVCCS) : public NETLIB_NAME(VCCS)
{
public:
ATTR_COLD NETLIB_NAME(LVCCS)()
: NETLIB_NAME(VCCS)(LVCCS), m_vi(0.0) { }
ATTR_COLD NETLIB_NAME(LVCCS)(const family_t afamily)
: NETLIB_NAME(VCCS)(afamily), m_vi(0.0) { }
param_double_t m_cur_limit; /* current limit */
protected:
virtual void start();
virtual void reset();
virtual void update_param();
ATTR_HOT virtual void update();
NETLIB_UPDATE_TERMINALSI();
nl_double m_vi;
};
// ----------------------------------------------------------------------------------------
// nld_CCCS
// ----------------------------------------------------------------------------------------
@ -153,6 +177,8 @@ public:
ATTR_COLD NETLIB_NAME(VCVS)()
: NETLIB_NAME(VCCS)(VCVS) { }
param_double_t m_RO;
protected:
virtual void start();
virtual void reset();
@ -162,7 +188,6 @@ protected:
terminal_t m_OP2;
terminal_t m_ON2;
param_double_t m_RO;
};
NETLIB_NAMESPACE_DEVICES_END()

View File

@ -37,3 +37,121 @@ NETLIST_START(opamp_lm3900)
PARAM(G1.RO, RES_K(8))
NETLIST_END()
NETLIB_NAMESPACE_DEVICES_START()
/*
* Type = 0: Impedance changer
* 1; Ideal opamp
* 2; opamp with first pole
* 3: opamp with first pole + output limit
* 4: opamp with input stage, first pole + output limit
*/
NETLIB_START(OPAMP)
{
register_sub("RP1", m_RP);
register_sub("CP1", m_CP);
register_sub("G1", m_G1);
register_sub("EBUF", m_EBUF);
register_sub("DN", m_DN);
register_sub("DP", m_DP);
register_input("VCC", m_VCC);
register_input("GND", m_GND);
register_output("VL", m_VL);
register_output("VH", m_VH);
register_output("VREF", m_VREF);
register_param("model", m_model, "");
m_type = m_model.model_value("TYPE");
if (m_type == 3)
{
register_subalias("PLUS", "G1.IP");
register_subalias("MINUS", "G1.IN");
register_subalias("OUT", "EBUF.OP");
connect_late("EBUF.ON", "VREF");
connect_late("G1.ON", "VREF");
connect_late("RP1.2", "VREF");
connect_late("CP1.2", "VREF");
connect_late("EBUF.IN", "VREF");
connect_late("RP1.1", "G1.OP");
connect_late("CP1.1", "RP1.1");
connect_late("DP.K", "VH");
connect_late("VL", "DN.A");
connect_late("DP.A", "DN.K");
connect_late("DN.K", "RP1.1");
connect_late("EBUF.IP", "RP1.1");
}
else
netlist().error("Unknown opamp type: %d", m_type);
}
/* .model abc OPAMP(VLH=2.0 VLL=0.2 FPF=5 UGF=10k SLEW=0.6u RI=1000k RO=50 DAB=0.002)
*
* Differential Amp Bias ~ op amp's total quiescent current.
* */
NETLIB_UPDATE(OPAMP)
{
const double cVt = 0.0258 * 1.0; // * m_n;
const double cId = m_model.model_value("DAB"); // 3 mA
const double cVd = cVt * nl_math::log(cId / 1e-15 + 1.0);
m_VH.set_Q(INPANALOG(m_VCC) - m_model.model_value("VLH") - cVd);
m_VL.set_Q(INPANALOG(m_GND) + m_model.model_value("VLL") + cVd);
m_VREF.set_Q((INPANALOG(m_VCC) + INPANALOG(m_GND)) / 2.0);
}
NETLIB_RESET(OPAMP)
{
m_EBUF.do_reset();
m_G1.do_reset();
m_DP.do_reset();
m_DN.do_reset();
m_CP.do_reset();
m_RP.do_reset();
m_EBUF.m_G.setTo(1.0);
m_G1.m_RI.setTo(m_model.model_value("RI"));
m_EBUF.m_RO.setTo(m_model.model_value("RO"));
m_DP.m_model.setTo(".model tt D(IS=1e-15 N=1)");
m_DN.m_model.setTo(".model tt D(IS=1e-15 N=1)");
double CP = m_model.model_value("DAB") / m_model.model_value("SLEW");
double RP = 0.5 / 3.1459 / CP / m_model.model_value("FPF");
double G = m_model.model_value("UGF") / m_model.model_value("FPF") / RP;
//printf("CP=%e RP=%f G=%f\n", CP, RP, G);
m_CP.m_C.setTo(CP);
m_RP.set_R(RP);
m_G1.m_G.setTo(G);
}
NETLIB_UPDATE_PARAM(OPAMP)
{
}
/*
NETLIB_DEVICE_WITH_PARAMS(OPAMPx,
NETLIB_NAME(R) m_RP;
NETLIB_NAME(C) m_CP;
NETLIB_NAME(VCCS) m_G1;
NETLIB_NAME(VCVS) m_EBUF;
param_model_t m_model;
analog_input_t m_VH;
analog_input_t m_VL;
analog_input_t m_VREF;
);
*/
NETLIB_NAMESPACE_DEVICES_END()

View File

@ -5,7 +5,7 @@
*
*/
#pragma once
//#pragma once
#ifndef NLD_OPAMPS_H_
#define NLD_OPAMPS_H_
@ -13,11 +13,16 @@
#include "../nl_base.h"
#include "../nl_setup.h"
#include "nld_twoterm.h"
#include "nld_fourterm.h"
// ----------------------------------------------------------------------------------------
// Macros
// ----------------------------------------------------------------------------------------
#define OPAMP(_name, _model) \
NET_REGISTER_DEV(OPAMP, _name) \
NETDEV_PARAMI(_name, model, _model)
#define LM3900(_name) \
SUBMODEL(opamp_lm3900, name)
@ -27,6 +32,28 @@
NETLIST_EXTERNAL(opamp_lm3900)
NETLIB_NAMESPACE_DEVICES_START()
NETLIB_DEVICE_WITH_PARAMS(OPAMP,
NETLIB_NAME(R) m_RP;
NETLIB_NAME(C) m_CP;
NETLIB_NAME(VCCS) m_G1;
NETLIB_NAME(VCVS) m_EBUF;
NETLIB_NAME(D) m_DP;
NETLIB_NAME(D) m_DN;
analog_input_t m_VCC;
analog_input_t m_GND;
param_model_t m_model;
analog_output_t m_VH;
analog_output_t m_VL;
analog_output_t m_VREF;
/* state */
unsigned m_type;
);
NETLIB_NAMESPACE_DEVICES_END()
#endif /* NLD_OPAMPS_H_ */

View File

@ -240,6 +240,7 @@ NETLIB_RESET(C)
NETLIB_UPDATE_PARAM(C)
{
//step_time(1.0/48000.0);
m_GParallel = netlist().gmin() * m_C.Value();
}
NETLIB_UPDATE(C)
@ -247,6 +248,14 @@ NETLIB_UPDATE(C)
NETLIB_NAME(twoterm)::update();
}
ATTR_HOT void NETLIB_NAME(C)::step_time(const nl_double st)
{
/* Gpar should support convergence */
const nl_double G = m_C.Value() / st + m_GParallel;
const nl_double I = -G * deltaV();
set(G, 0.0, I);
}
// ----------------------------------------------------------------------------------------
// nld_D
// ----------------------------------------------------------------------------------------

View File

@ -194,14 +194,11 @@ NETLIB_DEVICE_WITH_PARAMS(POT2,
class NETLIB_NAME(C) : public NETLIB_NAME(twoterm)
{
public:
ATTR_COLD NETLIB_NAME(C)() : NETLIB_NAME(twoterm)(CAPACITOR) { }
ATTR_COLD NETLIB_NAME(C)() : NETLIB_NAME(twoterm)(CAPACITOR), m_GParallel(0.0) { }
ATTR_HOT void step_time(const nl_double st)
{
const nl_double G = m_C.Value() / st;
const nl_double I = -G * deltaV();
set(G, 0.0, I);
}
ATTR_HOT void step_time(const nl_double st);
param_double_t m_C;
protected:
virtual void start();
@ -209,7 +206,8 @@ protected:
virtual void update_param();
ATTR_HOT void update();
param_double_t m_C;
private:
nl_double m_GParallel;
};
@ -225,6 +223,7 @@ public:
ATTR_HOT inline void update_diode(const nl_double nVd)
{
#if 1
if (nVd < NL_FCONST(-5.0) * m_Vt)
{
m_Vd = nVd;
@ -233,23 +232,32 @@ public:
}
else if (nVd < m_Vcrit)
{
const nl_double eVDVt = nl_math::exp(nVd * m_VtInv);
m_Vd = nVd;
//m_Vd = m_Vd + 10.0 * m_Vt * std::tanh((nVd - m_Vd) / 10.0 / m_Vt);
const nl_double eVDVt = nl_math::exp(m_Vd * m_VtInv);
m_Id = m_Is * (eVDVt - NL_FCONST(1.0));
m_G = m_Is * m_VtInv * eVDVt + m_gmin;
}
else
{
const nl_double a = std::max((nVd - m_Vd) * m_VtInv, NL_FCONST(1e-12) - NL_FCONST(1.0));
#if 1
const nl_double a = std::max((nVd - m_Vd) * m_VtInv, NL_FCONST(1e-1) - NL_FCONST(1.0));
m_Vd = m_Vd + nl_math::e_log1p(a) * m_Vt;
#else
m_Vd = m_Vd + 10.0 * m_Vt * std::tanh((nVd - m_Vd) / 10.0 / m_Vt);
#endif
const nl_double eVDVt = nl_math::exp(m_Vd * m_VtInv);
m_Id = m_Is * (eVDVt - NL_FCONST(1.0));
m_G = m_Is * m_VtInv * eVDVt + m_gmin;
}
//printf("nVd %f m_Vd %f Vcrit %f\n", nVd, m_Vd, m_Vcrit);
#else
m_Vd = m_Vd + 20.0 * m_Vt * std::tanh((nVd - m_Vd) / 20.0 / m_Vt);
const nl_double eVDVt = nl_math::exp(m_Vd * m_VtInv);
m_Id = m_Is * (eVDVt - NL_FCONST(1.0));
m_G = m_Is * m_VtInv * eVDVt + m_gmin;
#endif
//printf("%p nVd %f m_Vd %f Vcrit %f\n", this, nVd, m_Vd, m_Vcrit);
}
ATTR_COLD void set_param(const nl_double Is, const nl_double n, nl_double gmin);
@ -288,13 +296,13 @@ public:
NETLIB_UPDATE_TERMINALSI();
param_model_t m_model;
protected:
virtual void start();
virtual void update_param();
ATTR_HOT void update();
param_model_t m_model;
generic_diode m_D;
};

View File

@ -58,8 +58,10 @@ void initialize_factory(factory_list_t &factory)
ENTRY(VCVS, VCVS, "-")
ENTRY(VCCS, VCCS, "-")
ENTRY(CCCS, CCCS, "-")
ENTRY(LVCCS, LVCCS, "-")
ENTRY(VS, VS, "V")
ENTRY(CS, CS, "I")
ENTRY(OPAMP, OPAMP, "model")
ENTRY(dummy_input, DUMMY_INPUT, "-")
ENTRY(frontier, FRONTIER_DEV, "+I,G,Q") // not intended to be used directly
ENTRY(function, AFUNC, "N,FUNC") // only for macro devices - NO FEEDBACK loops

View File

@ -555,6 +555,19 @@ ATTR_COLD void device_t::register_subalias(const pstring &name, core_terminal_t
m_terminals.add(alias);
}
ATTR_COLD void device_t::register_subalias(const pstring &name, const pstring &aliased)
{
pstring alias = this->name() + "." + name;
pstring aliased_fqn = this->name() + "." + aliased;
// everything already fully qualified
setup().register_alias_nofqn(alias, aliased_fqn);
// FIXME: make this working again
//if (term.isType(terminal_t::INPUT) || term.isType(terminal_t::TERMINAL))
// m_terminals.add(name);
}
ATTR_COLD void device_t::register_terminal(const pstring &name, terminal_t &port)
{
setup().register_object(*this, name, port);
@ -588,13 +601,12 @@ ATTR_COLD void device_t::register_input(const pstring &name, analog_input_t &inp
ATTR_COLD void device_t::connect_late(core_terminal_t &t1, core_terminal_t &t2)
{
#if 1
//printf("device %s: connect %s to %s\n", name().cstr(), t1.name().cstr(), t2.name().cstr());
setup().register_link_fqn(t1.name(), t2.name());
#else
if (!setup().connect(t1, t2))
netlist().error("Error connecting %s to %s\n", t1.name().cstr(), t2.name().cstr());
#endif
}
ATTR_COLD void device_t::connect_late(const pstring &t1, const pstring &t2)
{
setup().register_link_fqn(name() + "." + t1, name() + "." + t2);
}
ATTR_COLD void device_t::connect_direct(core_terminal_t &t1, core_terminal_t &t2)

View File

@ -388,6 +388,7 @@ namespace netlist
BJT_SWITCH, // BJT(Switch)
VCVS, // Voltage controlled voltage source
VCCS, // Voltage controlled current source
LVCCS, // Voltage controlled current source (Current limited)
CCCS, // Current controlled current source
VS, // Voltage Source
CS, // Current Source
@ -1087,12 +1088,14 @@ namespace netlist
ATTR_COLD void register_sub(const pstring &name, device_t &dev);
ATTR_COLD void register_subalias(const pstring &name, core_terminal_t &term);
ATTR_COLD void register_subalias(const pstring &name, const pstring &aliased);
ATTR_COLD void register_terminal(const pstring &name, terminal_t &port);
ATTR_COLD void register_output(const pstring &name, analog_output_t &out);
ATTR_COLD void register_output(const pstring &name, logic_output_t &out);
ATTR_COLD void register_input(const pstring &name, analog_input_t &in);
ATTR_COLD void register_input(const pstring &name, logic_input_t &in);
ATTR_COLD void connect_late(const pstring &t1, const pstring &t2);
ATTR_COLD void connect_late(core_terminal_t &t1, core_terminal_t &t2);
ATTR_COLD void connect_direct(core_terminal_t &t1, core_terminal_t &t2);

View File

@ -919,6 +919,7 @@ nl_double setup_t::model_value(const pstring &model_str, const pstring &entity,
char numfac = *(tmp.right(1).cstr());
switch (numfac)
{
case 'M': factor = 1e6; break;
case 'k': factor = 1e3; break;
case 'm': factor = 1e-3; break;
case 'u': factor = 1e-6; break;

View File

@ -14,6 +14,10 @@
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
template <unsigned m_N, unsigned _storage_N>
class matrix_solver_direct_t: public matrix_solver_t
{
@ -50,7 +54,8 @@ protected:
*/
ATTR_HOT nl_double compute_next_timestep();
ATTR_ALIGN nl_double m_A[_storage_N][((_storage_N + 7) / 8) * 8];
ATTR_ALIGN nl_ext_double m_A[_storage_N][((_storage_N + 7) / 8) * 8];
//ATTR_ALIGN nl_double m_A[_storage_N][((_storage_N + 7) / 8) * 8];
ATTR_ALIGN nl_double m_RHS[_storage_N];
ATTR_ALIGN nl_double m_last_RHS[_storage_N]; // right hand side - contains currents
ATTR_ALIGN nl_double m_last_V[_storage_N];
@ -356,7 +361,8 @@ ATTR_HOT void matrix_solver_direct_t<m_N, _storage_N>::LE_solve()
unsigned maxrow = i;
for (unsigned j = i + 1; j < kN; j++)
{
if (nl_math::abs(m_A[j][i]) > nl_math::abs(m_A[maxrow][i]))
//if (std::abs(m_A[j][i]) > std::abs(m_A[maxrow][i]))
if (m_A[j][i] * m_A[j][i] > m_A[maxrow][i] * m_A[maxrow][i])
maxrow = j;
}
@ -368,43 +374,63 @@ ATTR_HOT void matrix_solver_direct_t<m_N, _storage_N>::LE_solve()
}
std::swap(m_RHS[i], m_RHS[maxrow]);
}
}
/* FIXME: Singular matrix? */
const nl_double f = 1.0 / m_A[i][i];
const nl_ext_double * RESTRICT s = &m_A[i][i+1];
/* FIXME: Singular matrix? */
const nl_double f = 1.0 / m_A[i][i];
const double * RESTRICT s = &m_A[i][0];
const unsigned *p = m_terms[i]->m_nzrd.data();
const unsigned e = m_terms[i]->m_nzrd.size();
/* Eliminate column i from row j */
/* Eliminate column i from row j */
for (unsigned j = i + 1; j < kN; j++)
{
double * RESTRICT d = &m_A[j][0];
const nl_double f1 = - d[i] * f;
if (f1 != NL_FCONST(0.0))
for (unsigned j = i + 1; j < kN; j++)
{
#if 0
/* The code below is 30% faster than the original
* implementation which is given here for reference.
*
* for (unsigned k = i + 1; k < kN; k++)
* m_A[j][k] = m_A[j][k] + m_A[i][k] * f1;
*/
double * RESTRICT d = &m_A[j][i+1];
const double * RESTRICT s = &m_A[i][i+1];
const int e = kN - i - 1;
for (int k = 0; k < e; k++)
d[k] = d[k] + s[k] * f1;
#else
for (unsigned k = 0; k < e; k++)
nl_ext_double * RESTRICT d = &m_A[j][i+1];
const nl_double f1 = - m_A[j][i] * f;
if (f1 != NL_FCONST(0.0))
{
const unsigned pk = p[k];
d[pk] += s[pk] * f1;
const unsigned e = kN - i - 1;
for (unsigned k = 0; k < e; k++)
d[k] = d[k] + s[k] * f1;
m_RHS[j] += m_RHS[i] * f1;
}
}
}
else
{
/* FIXME: Singular matrix? */
const nl_double f = 1.0 / m_A[i][i];
const nl_ext_double * RESTRICT s = &m_A[i][0];
const unsigned *p = m_terms[i]->m_nzrd.data();
const unsigned e = m_terms[i]->m_nzrd.size();
/* Eliminate column i from row j */
for (unsigned j = i + 1; j < kN; j++)
{
nl_ext_double * RESTRICT d = &m_A[j][0];
const nl_double f1 = - d[i] * f;
if (f1 != NL_FCONST(0.0))
{
#if 0
/* The code below is 30% faster than the original
* implementation which is given here for reference.
*
* for (unsigned k = i + 1; k < kN; k++)
* m_A[j][k] = m_A[j][k] + m_A[i][k] * f1;
*/
double * RESTRICT d = &m_A[j][i+1];
const double * RESTRICT s = &m_A[i][i+1];
const int e = kN - i - 1;
for (int k = 0; k < e; k++)
d[k] = d[k] + s[k] * f1;
#else
for (unsigned k = 0; k < e; k++)
{
const unsigned pk = p[k];
d[pk] += s[pk] * f1;
}
#endif
m_RHS[j] += m_RHS[i] * f1;
}
#endif
m_RHS[j] += m_RHS[i] * f1;
}
}
}
@ -422,14 +448,14 @@ ATTR_HOT void matrix_solver_direct_t<m_N, _storage_N>::LE_back_subst(
nl_double tmp = 0;
#if 1
#if 0
const double * RESTRICT A = &m_A[j][j+1];
const double * RESTRICT xp = &x[j+1];
const int e = kN - j - 1;
for (int k = 0; k < e; k++)
#if (USE_PIVOT_SEARCH)
const nl_ext_double * RESTRICT A = &m_A[j][j+1];
const nl_double * RESTRICT xp = &x[j+1];
const unsigned e = kN - j - 1;
for (unsigned k = 0; k < e; k++)
tmp += A[k] * xp[k];
#else
const double * RESTRICT A = &m_A[j][0];
const nl_ext_double * RESTRICT A = &m_A[j][0];
const unsigned *p = m_terms[j]->m_nzrd.data();
const unsigned e = m_terms[j]->m_nzrd.size();

View File

@ -199,7 +199,7 @@ ATTR_HOT inline int matrix_solver_SOR_mat_t<m_N, _storage_N>::vsolve_non_dynamic
{
nl_double Idrive = 0;
const double * RESTRICT A = &this->m_A[k][0];
const nl_ext_double * RESTRICT A = &this->m_A[k][0];
const unsigned *p = this->m_terms[k]->m_nz.data();
const unsigned e = this->m_terms[k]->m_nz.size();

View File

@ -134,9 +134,9 @@ ATTR_COLD void matrix_solver_t::setup(analog_net_t::list_t &nets)
break;
case device_t::BJT_EB:
case device_t::DIODE:
//case device_t::VCVS:
case device_t::LVCCS:
case device_t::BJT_SWITCH:
NL_VERBOSE_OUT(("found BJT/Diode\n"));
NL_VERBOSE_OUT(("found BJT/Diode/LVCCS\n"));
if (!m_dynamic_devices.contains(&p->device()))
m_dynamic_devices.add(&p->device());
break;