dcalc use try/catch
This commit is contained in:
parent
f1924594de
commit
4c3e0ec1bd
263
dcalc/DmpCeff.cc
263
dcalc/DmpCeff.cc
|
|
@ -66,6 +66,17 @@ enum DmpFunc { y20, y50, ipi };
|
||||||
|
|
||||||
static const char *dmp_func_index_strings[] = {"y20", "y50", "Ipi"};
|
static const char *dmp_func_index_strings[] = {"y20", "y50", "Ipi"};
|
||||||
|
|
||||||
|
class DmpError : public Exception
|
||||||
|
{
|
||||||
|
public:
|
||||||
|
DmpError(const char *what);
|
||||||
|
virtual ~DmpError() {}
|
||||||
|
virtual const char *what() const noexcept { return what_; }
|
||||||
|
|
||||||
|
private:
|
||||||
|
const char *what_;
|
||||||
|
};
|
||||||
|
|
||||||
static double
|
static double
|
||||||
gateModelRd(const LibertyCell *cell,
|
gateModelRd(const LibertyCell *cell,
|
||||||
GateTableModel *gate_model,
|
GateTableModel *gate_model,
|
||||||
|
|
@ -75,7 +86,7 @@ gateModelRd(const LibertyCell *cell,
|
||||||
float related_out_cap,
|
float related_out_cap,
|
||||||
const Pvt *pvt,
|
const Pvt *pvt,
|
||||||
bool pocv_enabled);
|
bool pocv_enabled);
|
||||||
static bool
|
static void
|
||||||
evalDmpEqnsState(void *state);
|
evalDmpEqnsState(void *state);
|
||||||
static void
|
static void
|
||||||
evalVoEqns(void *state,
|
evalVoEqns(void *state,
|
||||||
|
|
@ -94,16 +105,15 @@ findRoot(void (*func)(void *state, double x, double &y, double &dy),
|
||||||
double x1,
|
double x1,
|
||||||
double x2,
|
double x2,
|
||||||
double x_tol,
|
double x_tol,
|
||||||
int max_iter,
|
int max_iter);
|
||||||
const char *&error);
|
static void
|
||||||
static const char *
|
|
||||||
newtonRaphson(const int max_iter,
|
newtonRaphson(const int max_iter,
|
||||||
double x[],
|
double x[],
|
||||||
const int n,
|
const int n,
|
||||||
const double x_tol,
|
const double x_tol,
|
||||||
// eval(state) is called to fill fvec and fjac.
|
// eval(state) is called to fill fvec and fjac.
|
||||||
// Returns false if fails.
|
// Returns false if fails.
|
||||||
bool (*eval)(void *state),
|
void (*eval)(void *state),
|
||||||
void *state,
|
void *state,
|
||||||
// Temporaries supplied by caller.
|
// Temporaries supplied by caller.
|
||||||
double *fvec,
|
double *fvec,
|
||||||
|
|
@ -116,7 +126,7 @@ luSolve(double **a,
|
||||||
const int size,
|
const int size,
|
||||||
const int *index,
|
const int *index,
|
||||||
double b[]);
|
double b[]);
|
||||||
static const char *
|
static void
|
||||||
luDecomp(double **a,
|
luDecomp(double **a,
|
||||||
const int size,
|
const int size,
|
||||||
int *index,
|
int *index,
|
||||||
|
|
@ -154,7 +164,7 @@ public:
|
||||||
|
|
||||||
// Given x_ as a vector of input parameters, fill fvec_ with the
|
// Given x_ as a vector of input parameters, fill fvec_ with the
|
||||||
// equations evaluated at x_ and fjac_ with the jabobian evaluated at x_.
|
// equations evaluated at x_ and fjac_ with the jabobian evaluated at x_.
|
||||||
virtual bool evalDmpEqns() = 0;
|
virtual void evalDmpEqns() = 0;
|
||||||
// Output response to vs(t) ramp driving pi model load.
|
// Output response to vs(t) ramp driving pi model load.
|
||||||
double vo(double t);
|
double vo(double t);
|
||||||
double dVoDt(double t);
|
double dVoDt(double t);
|
||||||
|
|
@ -165,7 +175,7 @@ public:
|
||||||
|
|
||||||
protected:
|
protected:
|
||||||
// Find driver parameters t0, delta_t, Ceff.
|
// Find driver parameters t0, delta_t, Ceff.
|
||||||
const char *findDriverParams(double ceff);
|
void findDriverParams(double ceff);
|
||||||
void gateCapDelaySlew(double cl,
|
void gateCapDelaySlew(double cl,
|
||||||
double &delay,
|
double &delay,
|
||||||
double &slew);
|
double &slew);
|
||||||
|
|
@ -189,13 +199,11 @@ protected:
|
||||||
void showX();
|
void showX();
|
||||||
void showFvec();
|
void showFvec();
|
||||||
void showJacobian();
|
void showJacobian();
|
||||||
const char *findDriverDelaySlew(double &delay,
|
void findDriverDelaySlew(double &delay,
|
||||||
double &slew);
|
double &slew);
|
||||||
double findVoCrossing(double vth,
|
double findVoCrossing(double vth);
|
||||||
const char *&error);
|
|
||||||
void showVo();
|
void showVo();
|
||||||
bool findVlCrossing(double vth,
|
double findVlCrossing(double vth);
|
||||||
double &t);
|
|
||||||
void showVl();
|
void showVl();
|
||||||
void fail(const char *reason);
|
void fail(const char *reason);
|
||||||
|
|
||||||
|
|
@ -334,7 +342,7 @@ DmpAlg::init(const LibertyLibrary *drvr_library,
|
||||||
|
|
||||||
// Find Ceff, delta_t and t0 for the driver.
|
// Find Ceff, delta_t and t0 for the driver.
|
||||||
// Return error msg on failure.
|
// Return error msg on failure.
|
||||||
const char *
|
void
|
||||||
DmpAlg::findDriverParams(double ceff)
|
DmpAlg::findDriverParams(double ceff)
|
||||||
{
|
{
|
||||||
x_[DmpParam::ceff] = ceff;
|
x_[DmpParam::ceff] = ceff;
|
||||||
|
|
@ -344,29 +352,23 @@ DmpAlg::findDriverParams(double ceff)
|
||||||
double t0 = t_vth + log(1.0 - vth_) * rd_ * ceff - vth_ * dt;
|
double t0 = t_vth + log(1.0 - vth_) * rd_ * ceff - vth_ * dt;
|
||||||
x_[DmpParam::dt] = dt;
|
x_[DmpParam::dt] = dt;
|
||||||
x_[DmpParam::t0] = t0;
|
x_[DmpParam::t0] = t0;
|
||||||
const char *error = newtonRaphson(100, x_, nr_order_, driver_param_tol,
|
newtonRaphson(100, x_, nr_order_, driver_param_tol, evalDmpEqnsState,
|
||||||
evalDmpEqnsState,
|
this, fvec_, fjac_, index_, p_, scale_);
|
||||||
this, fvec_, fjac_, index_, p_, scale_);
|
t0_ = x_[DmpParam::t0];
|
||||||
if (error)
|
dt_ = x_[DmpParam::dt];
|
||||||
return error;
|
debugPrint3(debug_, "dmp_ceff", 3, " t0 = %s dt = %s ceff = %s\n",
|
||||||
else {
|
units_->timeUnit()->asString(t0_),
|
||||||
t0_ = x_[DmpParam::t0];
|
units_->timeUnit()->asString(dt_),
|
||||||
dt_ = x_[DmpParam::dt];
|
units_->capacitanceUnit()->asString(x_[DmpParam::ceff]));
|
||||||
debugPrint3(debug_, "dmp_ceff", 3, " t0 = %s dt = %s ceff = %s\n",
|
if (debug_->check("dmp_ceff", 4))
|
||||||
units_->timeUnit()->asString(t0_),
|
showVo();
|
||||||
units_->timeUnit()->asString(dt_),
|
|
||||||
units_->capacitanceUnit()->asString(x_[DmpParam::ceff]));
|
|
||||||
if (debug_->check("dmp_ceff", 4))
|
|
||||||
showVo();
|
|
||||||
return nullptr;
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|
||||||
static bool
|
static void
|
||||||
evalDmpEqnsState(void *state)
|
evalDmpEqnsState(void *state)
|
||||||
{
|
{
|
||||||
DmpAlg *alg = reinterpret_cast<DmpAlg *>(state);
|
DmpAlg *alg = reinterpret_cast<DmpAlg *>(state);
|
||||||
return alg->evalDmpEqns();
|
alg->evalDmpEqns();
|
||||||
}
|
}
|
||||||
|
|
||||||
void
|
void
|
||||||
|
|
@ -427,12 +429,9 @@ DmpAlg::dy(double t,
|
||||||
double &dydcl)
|
double &dydcl)
|
||||||
{
|
{
|
||||||
double t1 = t - t0;
|
double t1 = t - t0;
|
||||||
#if 0
|
|
||||||
if (t1 <= 0.0)
|
if (t1 <= 0.0)
|
||||||
dydt0 = dyddt = dydcl = 0.0;
|
dydt0 = dyddt = dydcl = 0.0;
|
||||||
else
|
else if (t1 <= dt) {
|
||||||
#endif
|
|
||||||
if (t1 <= dt) {
|
|
||||||
dydt0 = -y0dt(t1, cl) / dt;
|
dydt0 = -y0dt(t1, cl) / dt;
|
||||||
dyddt = -y0(t1, cl) / (dt * dt);
|
dyddt = -y0(t1, cl) / (dt * dt);
|
||||||
dydcl = y0dcl(t1, cl) / dt;
|
dydcl = y0dcl(t1, cl) / dt;
|
||||||
|
|
@ -488,34 +487,24 @@ DmpAlg::showJacobian()
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
// Return error msg on failure.
|
void
|
||||||
const char *
|
|
||||||
DmpAlg::findDriverDelaySlew(double &delay,
|
DmpAlg::findDriverDelaySlew(double &delay,
|
||||||
double &slew)
|
double &slew)
|
||||||
{
|
{
|
||||||
const char *error = nullptr;
|
const char *error = nullptr;
|
||||||
delay = findVoCrossing(vth_, error);
|
delay = findVoCrossing(vth_);
|
||||||
if (error)
|
double tl = findVoCrossing(vl_);
|
||||||
return error;
|
double th = findVoCrossing(vh_);
|
||||||
double tl = findVoCrossing(vl_, error);
|
|
||||||
if (error)
|
|
||||||
return error;
|
|
||||||
double th = findVoCrossing(vh_, error);
|
|
||||||
if (error)
|
|
||||||
return error;
|
|
||||||
slew = (th - tl) / slew_derate_;
|
slew = (th - tl) / slew_derate_;
|
||||||
return nullptr;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
// Find t such that vo(t)=v.
|
// Find t such that vo(t)=v.
|
||||||
// Return true if successful.
|
|
||||||
double
|
double
|
||||||
DmpAlg::findVoCrossing(double vth,
|
DmpAlg::findVoCrossing(double vth)
|
||||||
const char *&error)
|
|
||||||
{
|
{
|
||||||
v_cross_ = vth;
|
v_cross_ = vth;
|
||||||
double ub = voCrossingUpperBound();
|
double ub = voCrossingUpperBound();
|
||||||
return findRoot(evalVoEqns, this, t0_, ub, vth_time_tol, find_root_max_iter, error);
|
return findRoot(evalVoEqns, this, t0_, ub, vth_time_tol, find_root_max_iter);
|
||||||
}
|
}
|
||||||
|
|
||||||
static void
|
static void
|
||||||
|
|
@ -570,28 +559,27 @@ DmpAlg::loadDelaySlew(const Pin *,
|
||||||
ArcDelay &delay,
|
ArcDelay &delay,
|
||||||
Slew &slew)
|
Slew &slew)
|
||||||
{
|
{
|
||||||
if (elmore == 0.0) {
|
if (!driver_valid_
|
||||||
delay = 0.0;
|
|| elmore == 0.0
|
||||||
slew = gate_slew_;
|
// Elmore delay is small compared to driver slew.
|
||||||
}
|
|| elmore < gate_slew_ * 1e-3) {
|
||||||
else if (elmore < gate_slew_ * 1e-3) {
|
|
||||||
// Elmore delay is small compared to driver slew.
|
|
||||||
delay = elmore;
|
delay = elmore;
|
||||||
|
// solve v=1-exp(-t/rc) for t, elmore_slew_factor_ = t(vh) - t(vl)
|
||||||
|
// slew = elmore * (log(vh_) - log(vl_))
|
||||||
|
// slew = elmore * elmore_slew_factor_;
|
||||||
slew = gate_slew_;
|
slew = gate_slew_;
|
||||||
}
|
}
|
||||||
else {
|
else {
|
||||||
elmore_ = elmore;
|
|
||||||
p3_ = 1.0 / elmore;
|
|
||||||
if (driver_valid_
|
|
||||||
&& debug_->check("dmp_ceff", 4))
|
|
||||||
showVl();
|
|
||||||
// Use the driver thresholds and rely on thresholdAdjust to
|
// Use the driver thresholds and rely on thresholdAdjust to
|
||||||
// convert the delay and slew to the load's thresholds.
|
// convert the delay and slew to the load's thresholds.
|
||||||
double tl, th, load_delay;
|
try {
|
||||||
if (driver_valid_
|
if (debug_->check("dmp_ceff", 4))
|
||||||
&& findVlCrossing(vth_, load_delay)
|
showVl();
|
||||||
&& findVlCrossing(vl_, tl)
|
elmore_ = elmore;
|
||||||
&& findVlCrossing(vh_, th)) {
|
p3_ = 1.0 / elmore;
|
||||||
|
double load_delay = findVlCrossing(vth_);
|
||||||
|
double tl = findVlCrossing(vl_);
|
||||||
|
double th = findVlCrossing(vh_);
|
||||||
// Measure delay from Vo, the load dependent source excitation.
|
// Measure delay from Vo, the load dependent source excitation.
|
||||||
double delay1 = load_delay - vo_delay_;
|
double delay1 = load_delay - vo_delay_;
|
||||||
double slew1 = (th - tl) / slew_derate_;
|
double slew1 = (th - tl) / slew_derate_;
|
||||||
|
|
@ -600,7 +588,7 @@ DmpAlg::loadDelaySlew(const Pin *,
|
||||||
if (-delay1 > vth_time_tol * vo_delay_)
|
if (-delay1 > vth_time_tol * vo_delay_)
|
||||||
fail("load delay less than zero");
|
fail("load delay less than zero");
|
||||||
// Use elmore delay.
|
// Use elmore delay.
|
||||||
delay1 = 1.0 / p3_;
|
delay1 = elmore;
|
||||||
}
|
}
|
||||||
if (slew1 < gate_slew_) {
|
if (slew1 < gate_slew_) {
|
||||||
// Only report a problem if the difference is significant.
|
// Only report a problem if the difference is significant.
|
||||||
|
|
@ -611,30 +599,24 @@ DmpAlg::loadDelaySlew(const Pin *,
|
||||||
delay = delay1;
|
delay = delay1;
|
||||||
slew = slew1;
|
slew = slew1;
|
||||||
}
|
}
|
||||||
else {
|
catch (DmpError &error) {
|
||||||
// Failed - use elmore delay and driver slew.
|
// Failed - use elmore delay and driver slew.
|
||||||
delay = elmore_;
|
delay = elmore_;
|
||||||
// solve v=1-exp(-t/rc) for t, elmore_slew_factor_ = t(vh) - t(vl)
|
// solve v=1-exp(-t/rc) for t, elmore_slew_factor_ = t(vh) - t(vl)
|
||||||
// slew = elmore * (log(vh_) - log(vl_))
|
// slew = elmore * (log(vh_) - log(vl_))
|
||||||
slew = gate_slew_ + elmore * elmore_slew_factor_;
|
slew = elmore * elmore_slew_factor_;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
// Find t such that vl(t)=v.
|
// Find t such that vl(t)=v.
|
||||||
// Return true if successful.
|
// Return true if successful.
|
||||||
bool
|
double
|
||||||
DmpAlg::findVlCrossing(double vth,
|
DmpAlg::findVlCrossing(double vth)
|
||||||
double &t)
|
|
||||||
{
|
{
|
||||||
v_cross_ = vth;
|
v_cross_ = vth;
|
||||||
double ub = vlCrossingUpperBound();
|
double ub = vlCrossingUpperBound();
|
||||||
const char *error;
|
return findRoot(evalVlEqns, this, t0_, ub, vth_time_tol, find_root_max_iter);
|
||||||
t = findRoot(evalVlEqns, this, t0_, ub, vth_time_tol, find_root_max_iter,
|
|
||||||
error);
|
|
||||||
if (error)
|
|
||||||
fail("findVlCrossing: Vl(t) did not cross threshold");
|
|
||||||
return (error == nullptr);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
double
|
double
|
||||||
|
|
@ -725,7 +707,7 @@ public:
|
||||||
double elmore,
|
double elmore,
|
||||||
ArcDelay &delay,
|
ArcDelay &delay,
|
||||||
Slew &slew);
|
Slew &slew);
|
||||||
virtual bool evalDmpEqns();
|
virtual void evalDmpEqns();
|
||||||
virtual double voCrossingUpperBound();
|
virtual double voCrossingUpperBound();
|
||||||
|
|
||||||
private:
|
private:
|
||||||
|
|
@ -779,10 +761,9 @@ DmpCap::loadDelaySlew(const Pin *,
|
||||||
slew = gate_slew_;
|
slew = gate_slew_;
|
||||||
}
|
}
|
||||||
|
|
||||||
bool
|
void
|
||||||
DmpCap::evalDmpEqns()
|
DmpCap::evalDmpEqns()
|
||||||
{
|
{
|
||||||
return true;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
double
|
double
|
||||||
|
|
@ -836,11 +817,11 @@ public:
|
||||||
double c1);
|
double c1);
|
||||||
virtual void gateDelaySlew(double &delay,
|
virtual void gateDelaySlew(double &delay,
|
||||||
double &slew);
|
double &slew);
|
||||||
virtual bool evalDmpEqns();
|
virtual void evalDmpEqns();
|
||||||
virtual double voCrossingUpperBound();
|
virtual double voCrossingUpperBound();
|
||||||
|
|
||||||
private:
|
private:
|
||||||
const char *findDriverParamsPi();
|
void findDriverParamsPi();
|
||||||
virtual double v0(double t);
|
virtual double v0(double t);
|
||||||
virtual double dv0dt(double t);
|
virtual double dv0dt(double t);
|
||||||
double ipiIceff(double t0,
|
double ipiIceff(double t0,
|
||||||
|
|
@ -924,14 +905,8 @@ void
|
||||||
DmpPi::gateDelaySlew(double &delay,
|
DmpPi::gateDelaySlew(double &delay,
|
||||||
double &slew)
|
double &slew)
|
||||||
{
|
{
|
||||||
const char *error = findDriverParamsPi();
|
try {
|
||||||
if (error) {
|
findDriverParamsPi();
|
||||||
fail(error);
|
|
||||||
// Driver calculation failed - use Ceff=c1+c2.
|
|
||||||
ceff_ = c1_ + c2_;
|
|
||||||
gateCapDelaySlew(ceff_, delay, slew);
|
|
||||||
}
|
|
||||||
else {
|
|
||||||
ceff_ = x_[DmpParam::ceff];
|
ceff_ = x_[DmpParam::ceff];
|
||||||
driver_valid_ = true;
|
driver_valid_ = true;
|
||||||
double table_slew;
|
double table_slew;
|
||||||
|
|
@ -941,38 +916,50 @@ DmpPi::gateDelaySlew(double &delay,
|
||||||
// Vo slew is more accurate than table
|
// Vo slew is more accurate than table
|
||||||
// (-8% max, -3% avg vs -32% max, -12% avg).
|
// (-8% max, -3% avg vs -32% max, -12% avg).
|
||||||
// Need Vo delay to measure load wire delay waveform.
|
// Need Vo delay to measure load wire delay waveform.
|
||||||
const char *error = findDriverDelaySlew(vo_delay_, slew);
|
try {
|
||||||
if (error) {
|
findDriverDelaySlew(vo_delay_, slew);
|
||||||
fail(error);
|
}
|
||||||
|
catch (DmpError &error) {
|
||||||
|
fail(error.what());
|
||||||
// Fall back to table slew.
|
// Fall back to table slew.
|
||||||
slew = table_slew;
|
slew = table_slew;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
catch (DmpError &error) {
|
||||||
|
fail(error.what());
|
||||||
|
// Driver calculation failed - use Ceff=c1+c2.
|
||||||
|
driver_valid_ = false;
|
||||||
|
ceff_ = c1_ + c2_;
|
||||||
|
gateCapDelaySlew(ceff_, delay, slew);
|
||||||
|
}
|
||||||
// Save for wire delay calc.
|
// Save for wire delay calc.
|
||||||
gate_slew_ = slew;
|
gate_slew_ = slew;
|
||||||
}
|
}
|
||||||
|
|
||||||
const char *
|
void
|
||||||
DmpPi::findDriverParamsPi()
|
DmpPi::findDriverParamsPi()
|
||||||
{
|
{
|
||||||
const char *error;
|
try {
|
||||||
error = findDriverParams(c2_ + c1_);
|
findDriverParams(c2_ + c1_);
|
||||||
if (error)
|
}
|
||||||
error = findDriverParams(c2_);
|
catch (DmpError &) {
|
||||||
return error;
|
findDriverParams(c2_);
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
// Given x_ as a vector of input parameters, fill fvec_ with the
|
// Given x_ as a vector of input parameters, fill fvec_ with the
|
||||||
// equations evaluated at x_ and fjac_ with the jacobian evaluated at x_.
|
// equations evaluated at x_ and fjac_ with the jacobian evaluated at x_.
|
||||||
bool
|
void
|
||||||
DmpPi::evalDmpEqns()
|
DmpPi::evalDmpEqns()
|
||||||
{
|
{
|
||||||
double t0 = x_[DmpParam::t0];
|
double t0 = x_[DmpParam::t0];
|
||||||
double dt = x_[DmpParam::dt];
|
double dt = x_[DmpParam::dt];
|
||||||
double ceff = x_[DmpParam::ceff];
|
double ceff = x_[DmpParam::ceff];
|
||||||
|
|
||||||
if (ceff > (c1_ + c2_) || ceff < 0.0)
|
if (ceff < 0.0)
|
||||||
return false;
|
throw DmpError("eqn eval failed: ceff < 0");
|
||||||
|
if (ceff > (c1_ + c2_))
|
||||||
|
throw DmpError("eqn eval failed: ceff > c2 + c1");
|
||||||
|
|
||||||
double t_vth, t_vl, slew;
|
double t_vth, t_vl, slew;
|
||||||
gateDelays(ceff, t_vth, t_vl, slew);
|
gateDelays(ceff, t_vth, t_vl, slew);
|
||||||
|
|
@ -981,7 +968,7 @@ DmpPi::evalDmpEqns()
|
||||||
ceff_time = 1.4 * dt;
|
ceff_time = 1.4 * dt;
|
||||||
|
|
||||||
if (dt <= 0.0)
|
if (dt <= 0.0)
|
||||||
return false;
|
throw DmpError("eqn eval failed: dt < 0");
|
||||||
|
|
||||||
double exp_p1_dt = exp(-p1_ * dt);
|
double exp_p1_dt = exp(-p1_ * dt);
|
||||||
double exp_p2_dt = exp(-p2_ * dt);
|
double exp_p2_dt = exp(-p2_ * dt);
|
||||||
|
|
@ -1021,7 +1008,6 @@ DmpPi::evalDmpEqns()
|
||||||
showJacobian();
|
showJacobian();
|
||||||
debug_->print(".................\n");
|
debug_->print(".................\n");
|
||||||
}
|
}
|
||||||
return true;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
// Eqn 13, Eqn 14.
|
// Eqn 13, Eqn 14.
|
||||||
|
|
@ -1090,7 +1076,7 @@ class DmpOnePole : public DmpAlg
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
DmpOnePole(StaState *sta);
|
DmpOnePole(StaState *sta);
|
||||||
virtual bool evalDmpEqns();
|
virtual void evalDmpEqns();
|
||||||
virtual double voCrossingUpperBound();
|
virtual double voCrossingUpperBound();
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
@ -1099,7 +1085,7 @@ DmpOnePole::DmpOnePole(StaState *sta) :
|
||||||
{
|
{
|
||||||
}
|
}
|
||||||
|
|
||||||
bool
|
void
|
||||||
DmpOnePole::evalDmpEqns()
|
DmpOnePole::evalDmpEqns()
|
||||||
{
|
{
|
||||||
double t0 = x_[DmpParam::t0];
|
double t0 = x_[DmpParam::t0];
|
||||||
|
|
@ -1133,7 +1119,6 @@ DmpOnePole::evalDmpEqns()
|
||||||
showJacobian();
|
showJacobian();
|
||||||
debug_->print(".................\n");
|
debug_->print(".................\n");
|
||||||
}
|
}
|
||||||
return true;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
double
|
double
|
||||||
|
|
@ -1291,19 +1276,15 @@ findRoot(void (*func)(void *state, double x, double &y, double &dy),
|
||||||
double x1,
|
double x1,
|
||||||
double x2,
|
double x2,
|
||||||
double x_tol,
|
double x_tol,
|
||||||
int max_iter,
|
int max_iter)
|
||||||
const char *&error)
|
|
||||||
{
|
{
|
||||||
double y1, y2, dy;
|
double y1, y2, dy;
|
||||||
func(state, x1, y1, dy);
|
func(state, x1, y1, dy);
|
||||||
func(state, x2, y2, dy);
|
func(state, x2, y2, dy);
|
||||||
|
|
||||||
if ((y1 > 0.0 && y2 > 0.0) || (y1 < 0.0 && y2 < 0.0)) {
|
if ((y1 > 0.0 && y2 > 0.0) || (y1 < 0.0 && y2 < 0.0))
|
||||||
error = "findRoot: initial bounds do not surround a root";
|
throw DmpError("findRoot: initial bounds do not surround a root");
|
||||||
return 0.0;
|
|
||||||
}
|
|
||||||
|
|
||||||
error = nullptr;
|
|
||||||
if (y1 == 0.0)
|
if (y1 == 0.0)
|
||||||
return x1;
|
return x1;
|
||||||
|
|
||||||
|
|
@ -1346,20 +1327,19 @@ findRoot(void (*func)(void *state, double x, double &y, double &dy),
|
||||||
else
|
else
|
||||||
x2 = root;
|
x2 = root;
|
||||||
}
|
}
|
||||||
error = "findRoot: max iterations exceeded";
|
throw DmpError("findRoot: max iterations exceeded");
|
||||||
return 0.0;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
// Newton-Raphson iteration to find zeros of a function.
|
// Newton-Raphson iteration to find zeros of a function.
|
||||||
// x_tol is percentage that all changes in x must be less than (1.0 = 100%).
|
// x_tol is percentage that all changes in x must be less than (1.0 = 100%).
|
||||||
// Eval(state) is called to fill fvec and fjac (returns false if fails).
|
// Eval(state) is called to fill fvec and fjac (returns false if fails).
|
||||||
// Return error msg on failure.
|
// Return error msg on failure.
|
||||||
static const char *
|
static void
|
||||||
newtonRaphson(const int max_iter,
|
newtonRaphson(const int max_iter,
|
||||||
double x[],
|
double x[],
|
||||||
const int size,
|
const int size,
|
||||||
const double x_tol,
|
const double x_tol,
|
||||||
bool (*eval)(void *state),
|
void (*eval)(void *state),
|
||||||
void *state,
|
void *state,
|
||||||
// Temporaries supplied by caller.
|
// Temporaries supplied by caller.
|
||||||
double *fvec,
|
double *fvec,
|
||||||
|
|
@ -1369,29 +1349,23 @@ newtonRaphson(const int max_iter,
|
||||||
double *scale)
|
double *scale)
|
||||||
{
|
{
|
||||||
for (int k = 0; k < max_iter; k++) {
|
for (int k = 0; k < max_iter; k++) {
|
||||||
if (!eval(state))
|
eval(state);
|
||||||
return "Newton-Raphson eval failed";
|
|
||||||
|
|
||||||
for (int i = 0; i < size; i++)
|
for (int i = 0; i < size; i++)
|
||||||
// Right-hand side of linear equations.
|
// Right-hand side of linear equations.
|
||||||
p[i] = -fvec[i];
|
p[i] = -fvec[i];
|
||||||
const char *error = luDecomp(fjac, size, index, scale);
|
luDecomp(fjac, size, index, scale);
|
||||||
if (error)
|
luSolve(fjac, size, index, p);
|
||||||
return error;
|
|
||||||
else {
|
|
||||||
luSolve(fjac, size, index, p);
|
|
||||||
|
|
||||||
bool all_under_x_tol = true;
|
bool all_under_x_tol = true;
|
||||||
for (int i = 0; i < size; i++) {
|
for (int i = 0; i < size; i++) {
|
||||||
if (abs(p[i]) > abs(x[i]) * x_tol)
|
if (abs(p[i]) > abs(x[i]) * x_tol)
|
||||||
all_under_x_tol = false;
|
all_under_x_tol = false;
|
||||||
x[i] += p[i];
|
x[i] += p[i];
|
||||||
}
|
|
||||||
if (all_under_x_tol)
|
|
||||||
return nullptr;
|
|
||||||
}
|
}
|
||||||
|
if (all_under_x_tol)
|
||||||
|
return;
|
||||||
}
|
}
|
||||||
return "Newton-Raphson max iterations exceeded";
|
throw DmpError("Newton-Raphson max iterations exceeded");
|
||||||
}
|
}
|
||||||
|
|
||||||
// luDecomp, luSolve based on MatClass from C. R. Birchenhall,
|
// luDecomp, luSolve based on MatClass from C. R. Birchenhall,
|
||||||
|
|
@ -1406,7 +1380,7 @@ newtonRaphson(const int max_iter,
|
||||||
// Replaces a[0..size-1][0..size-1] by the LU decomposition.
|
// Replaces a[0..size-1][0..size-1] by the LU decomposition.
|
||||||
// index[0..size-1] is an output vector of the row permutations.
|
// index[0..size-1] is an output vector of the row permutations.
|
||||||
// Return error msg on failure.
|
// Return error msg on failure.
|
||||||
const char *
|
void
|
||||||
luDecomp(double **a,
|
luDecomp(double **a,
|
||||||
const int size,
|
const int size,
|
||||||
int *index,
|
int *index,
|
||||||
|
|
@ -1423,7 +1397,7 @@ luDecomp(double **a,
|
||||||
big = temp;
|
big = temp;
|
||||||
}
|
}
|
||||||
if (big == 0.0)
|
if (big == 0.0)
|
||||||
return "LU decomposition: no non-zero row element";
|
throw DmpError("LU decomposition: no non-zero row element");
|
||||||
scale[i] = 1.0 / big;
|
scale[i] = 1.0 / big;
|
||||||
}
|
}
|
||||||
int size_1 = size - 1;
|
int size_1 = size - 1;
|
||||||
|
|
@ -1473,7 +1447,6 @@ luDecomp(double **a,
|
||||||
a[i][j] *= pivot;
|
a[i][j] *= pivot;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
return nullptr;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
// Solves the set of size linear equations a*x=b, assuming A is LU form
|
// Solves the set of size linear equations a*x=b, assuming A is LU form
|
||||||
|
|
@ -1789,4 +1762,10 @@ DmpCeffDelayCalc::copyState(const StaState *sta)
|
||||||
dmp_zero_c2_->copyState(sta);
|
dmp_zero_c2_->copyState(sta);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
DmpError::DmpError(const char *what) :
|
||||||
|
what_(what)
|
||||||
|
{
|
||||||
|
//printf("DmpError %s\n", what);
|
||||||
|
}
|
||||||
|
|
||||||
} // namespace
|
} // namespace
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue