Skip to content

Commit

Permalink
References for 'BlueBrain/nmodl#1444'.
Browse files Browse the repository at this point in the history
  • Loading branch information
GitHub Actions Bot committed Sep 20, 2024
1 parent decc09b commit 00fc10f
Show file tree
Hide file tree
Showing 11 changed files with 899 additions and 29 deletions.
11 changes: 8 additions & 3 deletions global/neuron/thread_newton.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -295,6 +295,8 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
FUNC functor,
double eps = EPS,
int max_iter = MAX_ITER) {
// If finite differences are needed, this is stores the stepwidth.
Eigen::Matrix<double, N, 1> dX;
// Vector to store result of function F(X):
Eigen::Matrix<double, N, 1> F;
// Matrix to store Jacobian of F(X):
Expand All @@ -303,7 +305,7 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
int iter = -1;
while (++iter < max_iter) {
// calculate F, J from X using user-supplied functor
functor(X, F, J);
functor(X, dX, F, J);
if (is_converged(X, J, F, eps)) {
return iter;
}
Expand Down Expand Up @@ -341,10 +343,11 @@ EIGEN_DEVICE_FUNC int newton_solver_small_N(Eigen::Matrix<double, N, 1>& X,
int max_iter) {
bool invertible;
Eigen::Matrix<double, N, 1> F;
Eigen::Matrix<double, N, 1> dX;
Eigen::Matrix<double, N, N> J, J_inv;
int iter = -1;
while (++iter < max_iter) {
functor(X, F, J);
functor(X, dX, F, J);
if (is_converged(X, J, F, eps)) {
return iter;
}
Expand Down Expand Up @@ -599,10 +602,12 @@ namespace neuron {
functor_thread_newton_0(_nrn_mechanism_cache_range& _lmc, thread_newton_Instance& inst, thread_newton_NodeData& node_data, size_t id, Datum* _ppvar, Datum* _thread, thread_newton_ThreadVariables& _thread_vars, NrnThread* nt, double v)
: _lmc(_lmc), inst(inst), node_data(node_data), id(id), _ppvar(_ppvar), _thread(_thread), _thread_vars(_thread_vars), nt(nt), v(v)
{}
void operator()(const Eigen::Matrix<double, 1, 1>& nmodl_eigen_xm, Eigen::Matrix<double, 1, 1>& nmodl_eigen_fm, Eigen::Matrix<double, 1, 1>& nmodl_eigen_jm) const {
void operator()(const Eigen::Matrix<double, 1, 1>& nmodl_eigen_xm, Eigen::Matrix<double, 1, 1>& nmodl_eigen_dxm, Eigen::Matrix<double, 1, 1>& nmodl_eigen_fm, Eigen::Matrix<double, 1, 1>& nmodl_eigen_jm) const {
const double* nmodl_eigen_x = nmodl_eigen_xm.data();
double* nmodl_eigen_dx = nmodl_eigen_dxm.data();
double* nmodl_eigen_j = nmodl_eigen_jm.data();
double* nmodl_eigen_f = nmodl_eigen_fm.data();
nmodl_eigen_dx[0] = std::max(1e-6, 0.02*std::fabs(nmodl_eigen_x[0]));
nmodl_eigen_f[static_cast<int>(0)] = ( -nmodl_eigen_x[static_cast<int>(0)] + nt->_dt * source0_ + old_X) / nt->_dt;
nmodl_eigen_j[static_cast<int>(0)] = -1.0 / nt->_dt;
}
Expand Down
12 changes: 9 additions & 3 deletions kinetic/coreneuron/X2Y.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -304,6 +304,8 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
FUNC functor,
double eps = EPS,
int max_iter = MAX_ITER) {
// If finite differences are needed, this is stores the stepwidth.
Eigen::Matrix<double, N, 1> dX;
// Vector to store result of function F(X):
Eigen::Matrix<double, N, 1> F;
// Matrix to store Jacobian of F(X):
Expand All @@ -312,7 +314,7 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
int iter = -1;
while (++iter < max_iter) {
// calculate F, J from X using user-supplied functor
functor(X, F, J);
functor(X, dX, F, J);
if (is_converged(X, J, F, eps)) {
return iter;
}
Expand Down Expand Up @@ -350,10 +352,11 @@ EIGEN_DEVICE_FUNC int newton_solver_small_N(Eigen::Matrix<double, N, 1>& X,
int max_iter) {
bool invertible;
Eigen::Matrix<double, N, 1> F;
Eigen::Matrix<double, N, 1> dX;
Eigen::Matrix<double, N, N> J, J_inv;
int iter = -1;
while (++iter < max_iter) {
functor(X, F, J);
functor(X, dX, F, J);
if (is_converged(X, J, F, eps)) {
return iter;
}
Expand Down Expand Up @@ -636,10 +639,13 @@ namespace coreneuron {
functor_X2Y_0(NrnThread* nt, X2Y_Instance* inst, int id, int pnodecount, double v, const Datum* indexes, double* data, ThreadDatum* thread)
: nt(nt), inst(inst), id(id), pnodecount(pnodecount), v(v), indexes(indexes), data(data), thread(thread)
{}
void operator()(const Eigen::Matrix<double, 2, 1>& nmodl_eigen_xm, Eigen::Matrix<double, 2, 1>& nmodl_eigen_fm, Eigen::Matrix<double, 2, 2>& nmodl_eigen_jm) const {
void operator()(const Eigen::Matrix<double, 2, 1>& nmodl_eigen_xm, Eigen::Matrix<double, 2, 1>& nmodl_eigen_dxm, Eigen::Matrix<double, 2, 1>& nmodl_eigen_fm, Eigen::Matrix<double, 2, 2>& nmodl_eigen_jm) const {
const double* nmodl_eigen_x = nmodl_eigen_xm.data();
double* nmodl_eigen_dx = nmodl_eigen_dxm.data();
double* nmodl_eigen_j = nmodl_eigen_jm.data();
double* nmodl_eigen_f = nmodl_eigen_fm.data();
nmodl_eigen_dx[0] = std::max(1e-6, 0.02*std::fabs(nmodl_eigen_x[0]));
nmodl_eigen_dx[1] = std::max(1e-6, 0.02*std::fabs(nmodl_eigen_x[1]));
nmodl_eigen_f[static_cast<int>(0)] = ( -nmodl_eigen_x[static_cast<int>(0)] + nt->_dt * ( -nmodl_eigen_x[static_cast<int>(0)] * kf0_ + nmodl_eigen_x[static_cast<int>(1)] * kb0_) + old_X) / nt->_dt;
nmodl_eigen_j[static_cast<int>(0)] = -kf0_ - 1.0 / nt->_dt;
nmodl_eigen_j[static_cast<int>(2)] = kb0_;
Expand Down
12 changes: 9 additions & 3 deletions kinetic/coreneuron/side_effects.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -304,6 +304,8 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
FUNC functor,
double eps = EPS,
int max_iter = MAX_ITER) {
// If finite differences are needed, this is stores the stepwidth.
Eigen::Matrix<double, N, 1> dX;
// Vector to store result of function F(X):
Eigen::Matrix<double, N, 1> F;
// Matrix to store Jacobian of F(X):
Expand All @@ -312,7 +314,7 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
int iter = -1;
while (++iter < max_iter) {
// calculate F, J from X using user-supplied functor
functor(X, F, J);
functor(X, dX, F, J);
if (is_converged(X, J, F, eps)) {
return iter;
}
Expand Down Expand Up @@ -350,10 +352,11 @@ EIGEN_DEVICE_FUNC int newton_solver_small_N(Eigen::Matrix<double, N, 1>& X,
int max_iter) {
bool invertible;
Eigen::Matrix<double, N, 1> F;
Eigen::Matrix<double, N, 1> dX;
Eigen::Matrix<double, N, N> J, J_inv;
int iter = -1;
while (++iter < max_iter) {
functor(X, F, J);
functor(X, dX, F, J);
if (is_converged(X, J, F, eps)) {
return iter;
}
Expand Down Expand Up @@ -635,10 +638,13 @@ namespace coreneuron {
functor_side_effects_0(NrnThread* nt, side_effects_Instance* inst, int id, int pnodecount, double v, const Datum* indexes, double* data, ThreadDatum* thread)
: nt(nt), inst(inst), id(id), pnodecount(pnodecount), v(v), indexes(indexes), data(data), thread(thread)
{}
void operator()(const Eigen::Matrix<double, 2, 1>& nmodl_eigen_xm, Eigen::Matrix<double, 2, 1>& nmodl_eigen_fm, Eigen::Matrix<double, 2, 2>& nmodl_eigen_jm) const {
void operator()(const Eigen::Matrix<double, 2, 1>& nmodl_eigen_xm, Eigen::Matrix<double, 2, 1>& nmodl_eigen_dxm, Eigen::Matrix<double, 2, 1>& nmodl_eigen_fm, Eigen::Matrix<double, 2, 2>& nmodl_eigen_jm) const {
const double* nmodl_eigen_x = nmodl_eigen_xm.data();
double* nmodl_eigen_dx = nmodl_eigen_dxm.data();
double* nmodl_eigen_j = nmodl_eigen_jm.data();
double* nmodl_eigen_f = nmodl_eigen_fm.data();
nmodl_eigen_dx[0] = std::max(1e-6, 0.02*std::fabs(nmodl_eigen_x[0]));
nmodl_eigen_dx[1] = std::max(1e-6, 0.02*std::fabs(nmodl_eigen_x[1]));
nmodl_eigen_f[static_cast<int>(0)] = ( -nmodl_eigen_x[static_cast<int>(0)] + nt->_dt * ( -nmodl_eigen_x[static_cast<int>(0)] * kf0_ + nmodl_eigen_x[static_cast<int>(1)] * kb0_) + old_X) / nt->_dt;
nmodl_eigen_j[static_cast<int>(0)] = -kf0_ - 1.0 / nt->_dt;
nmodl_eigen_j[static_cast<int>(2)] = kb0_;
Expand Down
12 changes: 9 additions & 3 deletions kinetic/neuron/X2Y.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -295,6 +295,8 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
FUNC functor,
double eps = EPS,
int max_iter = MAX_ITER) {
// If finite differences are needed, this is stores the stepwidth.
Eigen::Matrix<double, N, 1> dX;
// Vector to store result of function F(X):
Eigen::Matrix<double, N, 1> F;
// Matrix to store Jacobian of F(X):
Expand All @@ -303,7 +305,7 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
int iter = -1;
while (++iter < max_iter) {
// calculate F, J from X using user-supplied functor
functor(X, F, J);
functor(X, dX, F, J);
if (is_converged(X, J, F, eps)) {
return iter;
}
Expand Down Expand Up @@ -341,10 +343,11 @@ EIGEN_DEVICE_FUNC int newton_solver_small_N(Eigen::Matrix<double, N, 1>& X,
int max_iter) {
bool invertible;
Eigen::Matrix<double, N, 1> F;
Eigen::Matrix<double, N, 1> dX;
Eigen::Matrix<double, N, N> J, J_inv;
int iter = -1;
while (++iter < max_iter) {
functor(X, F, J);
functor(X, dX, F, J);
if (is_converged(X, J, F, eps)) {
return iter;
}
Expand Down Expand Up @@ -595,10 +598,13 @@ namespace neuron {
functor_X2Y_0(_nrn_mechanism_cache_range& _lmc, X2Y_Instance& inst, X2Y_NodeData& node_data, size_t id, Datum* _ppvar, Datum* _thread, NrnThread* nt, double v)
: _lmc(_lmc), inst(inst), node_data(node_data), id(id), _ppvar(_ppvar), _thread(_thread), nt(nt), v(v)
{}
void operator()(const Eigen::Matrix<double, 2, 1>& nmodl_eigen_xm, Eigen::Matrix<double, 2, 1>& nmodl_eigen_fm, Eigen::Matrix<double, 2, 2>& nmodl_eigen_jm) const {
void operator()(const Eigen::Matrix<double, 2, 1>& nmodl_eigen_xm, Eigen::Matrix<double, 2, 1>& nmodl_eigen_dxm, Eigen::Matrix<double, 2, 1>& nmodl_eigen_fm, Eigen::Matrix<double, 2, 2>& nmodl_eigen_jm) const {
const double* nmodl_eigen_x = nmodl_eigen_xm.data();
double* nmodl_eigen_dx = nmodl_eigen_dxm.data();
double* nmodl_eigen_j = nmodl_eigen_jm.data();
double* nmodl_eigen_f = nmodl_eigen_fm.data();
nmodl_eigen_dx[0] = std::max(1e-6, 0.02*std::fabs(nmodl_eigen_x[0]));
nmodl_eigen_dx[1] = std::max(1e-6, 0.02*std::fabs(nmodl_eigen_x[1]));
nmodl_eigen_f[static_cast<int>(0)] = ( -nmodl_eigen_x[static_cast<int>(0)] + nt->_dt * ( -nmodl_eigen_x[static_cast<int>(0)] * kf0_ + nmodl_eigen_x[static_cast<int>(1)] * kb0_) + old_X) / nt->_dt;
nmodl_eigen_j[static_cast<int>(0)] = -kf0_ - 1.0 / nt->_dt;
nmodl_eigen_j[static_cast<int>(2)] = kb0_;
Expand Down
12 changes: 9 additions & 3 deletions kinetic/neuron/side_effects.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -295,6 +295,8 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
FUNC functor,
double eps = EPS,
int max_iter = MAX_ITER) {
// If finite differences are needed, this is stores the stepwidth.
Eigen::Matrix<double, N, 1> dX;
// Vector to store result of function F(X):
Eigen::Matrix<double, N, 1> F;
// Matrix to store Jacobian of F(X):
Expand All @@ -303,7 +305,7 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
int iter = -1;
while (++iter < max_iter) {
// calculate F, J from X using user-supplied functor
functor(X, F, J);
functor(X, dX, F, J);
if (is_converged(X, J, F, eps)) {
return iter;
}
Expand Down Expand Up @@ -341,10 +343,11 @@ EIGEN_DEVICE_FUNC int newton_solver_small_N(Eigen::Matrix<double, N, 1>& X,
int max_iter) {
bool invertible;
Eigen::Matrix<double, N, 1> F;
Eigen::Matrix<double, N, 1> dX;
Eigen::Matrix<double, N, N> J, J_inv;
int iter = -1;
while (++iter < max_iter) {
functor(X, F, J);
functor(X, dX, F, J);
if (is_converged(X, J, F, eps)) {
return iter;
}
Expand Down Expand Up @@ -596,10 +599,13 @@ namespace neuron {
functor_side_effects_0(_nrn_mechanism_cache_range& _lmc, side_effects_Instance& inst, side_effects_NodeData& node_data, size_t id, Datum* _ppvar, Datum* _thread, NrnThread* nt, double v)
: _lmc(_lmc), inst(inst), node_data(node_data), id(id), _ppvar(_ppvar), _thread(_thread), nt(nt), v(v)
{}
void operator()(const Eigen::Matrix<double, 2, 1>& nmodl_eigen_xm, Eigen::Matrix<double, 2, 1>& nmodl_eigen_fm, Eigen::Matrix<double, 2, 2>& nmodl_eigen_jm) const {
void operator()(const Eigen::Matrix<double, 2, 1>& nmodl_eigen_xm, Eigen::Matrix<double, 2, 1>& nmodl_eigen_dxm, Eigen::Matrix<double, 2, 1>& nmodl_eigen_fm, Eigen::Matrix<double, 2, 2>& nmodl_eigen_jm) const {
const double* nmodl_eigen_x = nmodl_eigen_xm.data();
double* nmodl_eigen_dx = nmodl_eigen_dxm.data();
double* nmodl_eigen_j = nmodl_eigen_jm.data();
double* nmodl_eigen_f = nmodl_eigen_fm.data();
nmodl_eigen_dx[0] = std::max(1e-6, 0.02*std::fabs(nmodl_eigen_x[0]));
nmodl_eigen_dx[1] = std::max(1e-6, 0.02*std::fabs(nmodl_eigen_x[1]));
nmodl_eigen_f[static_cast<int>(0)] = ( -nmodl_eigen_x[static_cast<int>(0)] + nt->_dt * ( -nmodl_eigen_x[static_cast<int>(0)] * kf0_ + nmodl_eigen_x[static_cast<int>(1)] * kb0_) + old_X) / nt->_dt;
nmodl_eigen_j[static_cast<int>(0)] = -kf0_ - 1.0 / nt->_dt;
nmodl_eigen_j[static_cast<int>(2)] = kb0_;
Expand Down
4 changes: 4 additions & 0 deletions solve/coreneuron/finite_difference.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
[NMODL] [error] :: Code incompatibility detected
[NMODL] [error] :: Cannot translate mod file to .cpp file
[NMODL] [error] :: Fix the following errors and try again
[NMODL] [error] :: Code Incompatibility :: "a" variable found at [3.12] should be defined as a RANGE variable instead of GLOBAL to enable backend transformations
11 changes: 8 additions & 3 deletions solve/neuron/derivimplicit_array.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -295,6 +295,8 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
FUNC functor,
double eps = EPS,
int max_iter = MAX_ITER) {
// If finite differences are needed, this is stores the stepwidth.
Eigen::Matrix<double, N, 1> dX;
// Vector to store result of function F(X):
Eigen::Matrix<double, N, 1> F;
// Matrix to store Jacobian of F(X):
Expand All @@ -303,7 +305,7 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
int iter = -1;
while (++iter < max_iter) {
// calculate F, J from X using user-supplied functor
functor(X, F, J);
functor(X, dX, F, J);
if (is_converged(X, J, F, eps)) {
return iter;
}
Expand Down Expand Up @@ -341,10 +343,11 @@ EIGEN_DEVICE_FUNC int newton_solver_small_N(Eigen::Matrix<double, N, 1>& X,
int max_iter) {
bool invertible;
Eigen::Matrix<double, N, 1> F;
Eigen::Matrix<double, N, 1> dX;
Eigen::Matrix<double, N, N> J, J_inv;
int iter = -1;
while (++iter < max_iter) {
functor(X, F, J);
functor(X, dX, F, J);
if (is_converged(X, J, F, eps)) {
return iter;
}
Expand Down Expand Up @@ -581,10 +584,12 @@ namespace neuron {
functor_derivimplicit_array_0(_nrn_mechanism_cache_range& _lmc, derivimplicit_array_Instance& inst, derivimplicit_array_NodeData& node_data, size_t id, Datum* _ppvar, Datum* _thread, NrnThread* nt, double v)
: _lmc(_lmc), inst(inst), node_data(node_data), id(id), _ppvar(_ppvar), _thread(_thread), nt(nt), v(v)
{}
void operator()(const Eigen::Matrix<double, 1, 1>& nmodl_eigen_xm, Eigen::Matrix<double, 1, 1>& nmodl_eigen_fm, Eigen::Matrix<double, 1, 1>& nmodl_eigen_jm) const {
void operator()(const Eigen::Matrix<double, 1, 1>& nmodl_eigen_xm, Eigen::Matrix<double, 1, 1>& nmodl_eigen_dxm, Eigen::Matrix<double, 1, 1>& nmodl_eigen_fm, Eigen::Matrix<double, 1, 1>& nmodl_eigen_jm) const {
const double* nmodl_eigen_x = nmodl_eigen_xm.data();
double* nmodl_eigen_dx = nmodl_eigen_dxm.data();
double* nmodl_eigen_j = nmodl_eigen_jm.data();
double* nmodl_eigen_f = nmodl_eigen_fm.data();
nmodl_eigen_dx[0] = std::max(1e-6, 0.02*std::fabs(nmodl_eigen_x[0]));
nmodl_eigen_f[static_cast<int>(0)] = (nmodl_eigen_x[static_cast<int>(0)] * nt->_dt * ((inst.s+id*2)[static_cast<int>(0)] + (inst.s+id*2)[static_cast<int>(1)]) * (inst.z+id*3)[static_cast<int>(0)] * (inst.z+id*3)[static_cast<int>(1)] * (inst.z+id*3)[static_cast<int>(2)] - nmodl_eigen_x[static_cast<int>(0)] + old_x) / nt->_dt;
nmodl_eigen_j[static_cast<int>(0)] = (nt->_dt * ((inst.s+id*2)[static_cast<int>(0)] + (inst.s+id*2)[static_cast<int>(1)]) * (inst.z+id*3)[static_cast<int>(0)] * (inst.z+id*3)[static_cast<int>(1)] * (inst.z+id*3)[static_cast<int>(2)] - 1.0) / nt->_dt;
}
Expand Down
11 changes: 8 additions & 3 deletions solve/neuron/derivimplicit_scalar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -295,6 +295,8 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
FUNC functor,
double eps = EPS,
int max_iter = MAX_ITER) {
// If finite differences are needed, this is stores the stepwidth.
Eigen::Matrix<double, N, 1> dX;
// Vector to store result of function F(X):
Eigen::Matrix<double, N, 1> F;
// Matrix to store Jacobian of F(X):
Expand All @@ -303,7 +305,7 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
int iter = -1;
while (++iter < max_iter) {
// calculate F, J from X using user-supplied functor
functor(X, F, J);
functor(X, dX, F, J);
if (is_converged(X, J, F, eps)) {
return iter;
}
Expand Down Expand Up @@ -341,10 +343,11 @@ EIGEN_DEVICE_FUNC int newton_solver_small_N(Eigen::Matrix<double, N, 1>& X,
int max_iter) {
bool invertible;
Eigen::Matrix<double, N, 1> F;
Eigen::Matrix<double, N, 1> dX;
Eigen::Matrix<double, N, N> J, J_inv;
int iter = -1;
while (++iter < max_iter) {
functor(X, F, J);
functor(X, dX, F, J);
if (is_converged(X, J, F, eps)) {
return iter;
}
Expand Down Expand Up @@ -569,10 +572,12 @@ namespace neuron {
functor_derivimplicit_scalar_0(_nrn_mechanism_cache_range& _lmc, derivimplicit_scalar_Instance& inst, derivimplicit_scalar_NodeData& node_data, size_t id, Datum* _ppvar, Datum* _thread, NrnThread* nt, double v)
: _lmc(_lmc), inst(inst), node_data(node_data), id(id), _ppvar(_ppvar), _thread(_thread), nt(nt), v(v)
{}
void operator()(const Eigen::Matrix<double, 1, 1>& nmodl_eigen_xm, Eigen::Matrix<double, 1, 1>& nmodl_eigen_fm, Eigen::Matrix<double, 1, 1>& nmodl_eigen_jm) const {
void operator()(const Eigen::Matrix<double, 1, 1>& nmodl_eigen_xm, Eigen::Matrix<double, 1, 1>& nmodl_eigen_dxm, Eigen::Matrix<double, 1, 1>& nmodl_eigen_fm, Eigen::Matrix<double, 1, 1>& nmodl_eigen_jm) const {
const double* nmodl_eigen_x = nmodl_eigen_xm.data();
double* nmodl_eigen_dx = nmodl_eigen_dxm.data();
double* nmodl_eigen_j = nmodl_eigen_jm.data();
double* nmodl_eigen_f = nmodl_eigen_fm.data();
nmodl_eigen_dx[0] = std::max(1e-6, 0.02*std::fabs(nmodl_eigen_x[0]));
nmodl_eigen_f[static_cast<int>(0)] = ( -nmodl_eigen_x[static_cast<int>(0)] * nt->_dt - nmodl_eigen_x[static_cast<int>(0)] + old_x) / nt->_dt;
nmodl_eigen_j[static_cast<int>(0)] = ( -nt->_dt - 1.0) / nt->_dt;
}
Expand Down
Loading

0 comments on commit 00fc10f

Please sign in to comment.