Skip to content

Commit

Permalink
References for 'BlueBrain/nmodl#1417'.
Browse files Browse the repository at this point in the history
  • Loading branch information
GitHub Actions Bot committed Aug 30, 2024
1 parent 95f53e2 commit 6a09015
Show file tree
Hide file tree
Showing 9 changed files with 93 additions and 1,763 deletions.
32 changes: 11 additions & 21 deletions global/neuron/thread_newton.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -248,22 +248,8 @@ namespace newton {
* @{
*/

static constexpr int MAX_ITER = 50;
static constexpr double EPS = 1e-13;

template <int N>
bool is_converged(const Eigen::Matrix<double, N, 1>& X,
const Eigen::Matrix<double, N, N>& J,
const Eigen::Matrix<double, N, 1>& F,
double eps) {
for (Eigen::Index i = 0; i < N; ++i) {
double square_error = J(i, Eigen::all).cwiseAbs2() * (eps * X).cwiseAbs2();
if (F(i) * F(i) > square_error) {
return false;
}
}
return true;
}
static constexpr int MAX_ITER = 1e3;
static constexpr double EPS = 1e-12;

/**
* \brief Newton method with user-provided Jacobian
Expand All @@ -286,14 +272,17 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
int max_iter = MAX_ITER) {
// Vector to store result of function F(X):
Eigen::Matrix<double, N, 1> F;
// Matrix to store Jacobian of F(X):
// Matrix to store jacobian of F(X):
Eigen::Matrix<double, N, N> J;
// Solver iteration count:
int iter = -1;
while (++iter < max_iter) {
// calculate F, J from X using user-supplied functor
functor(X, F, J);
if (is_converged(X, J, F, eps)) {
// get error norm: here we use sqrt(|F|^2)
double error = F.norm();
if (error < eps) {
// we have converged: return iteration count
return iter;
}
// In Eigen the default storage order is ColMajor.
Expand Down Expand Up @@ -334,7 +323,8 @@ EIGEN_DEVICE_FUNC int newton_solver_small_N(Eigen::Matrix<double, N, 1>& X,
int iter = -1;
while (++iter < max_iter) {
functor(X, F, J);
if (is_converged(X, J, F, eps)) {
double error = F.norm();
if (error < eps) {
return iter;
}
// The inverse can be called from within OpenACC regions without any issue, as opposed to
Expand Down Expand Up @@ -571,8 +561,8 @@ namespace neuron {
const double* nmodl_eigen_x = nmodl_eigen_xm.data();
double* nmodl_eigen_j = nmodl_eigen_jm.data();
double* nmodl_eigen_f = nmodl_eigen_fm.data();
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;
nmodl_eigen_f[static_cast<int>(0)] = -nmodl_eigen_x[static_cast<int>(0)] + nt->_dt * source0_ + old_X;
nmodl_eigen_j[static_cast<int>(0)] = -1.0;
}

void finalize() {
Expand Down
40 changes: 15 additions & 25 deletions kinetic/coreneuron/X2Y.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -257,22 +257,8 @@ namespace newton {
* @{
*/

static constexpr int MAX_ITER = 50;
static constexpr double EPS = 1e-13;

template <int N>
bool is_converged(const Eigen::Matrix<double, N, 1>& X,
const Eigen::Matrix<double, N, N>& J,
const Eigen::Matrix<double, N, 1>& F,
double eps) {
for (Eigen::Index i = 0; i < N; ++i) {
double square_error = J(i, Eigen::all).cwiseAbs2() * (eps * X).cwiseAbs2();
if (F(i) * F(i) > square_error) {
return false;
}
}
return true;
}
static constexpr int MAX_ITER = 1e3;
static constexpr double EPS = 1e-12;

/**
* \brief Newton method with user-provided Jacobian
Expand All @@ -295,14 +281,17 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
int max_iter = MAX_ITER) {
// Vector to store result of function F(X):
Eigen::Matrix<double, N, 1> F;
// Matrix to store Jacobian of F(X):
// Matrix to store jacobian of F(X):
Eigen::Matrix<double, N, N> J;
// Solver iteration count:
int iter = -1;
while (++iter < max_iter) {
// calculate F, J from X using user-supplied functor
functor(X, F, J);
if (is_converged(X, J, F, eps)) {
// get error norm: here we use sqrt(|F|^2)
double error = F.norm();
if (error < eps) {
// we have converged: return iteration count
return iter;
}
// In Eigen the default storage order is ColMajor.
Expand Down Expand Up @@ -343,7 +332,8 @@ EIGEN_DEVICE_FUNC int newton_solver_small_N(Eigen::Matrix<double, N, 1>& X,
int iter = -1;
while (++iter < max_iter) {
functor(X, F, J);
if (is_converged(X, J, F, eps)) {
double error = F.norm();
if (error < eps) {
return iter;
}
// The inverse can be called from within OpenACC regions without any issue, as opposed to
Expand Down Expand Up @@ -629,12 +619,12 @@ namespace coreneuron {
const double* nmodl_eigen_x = nmodl_eigen_xm.data();
double* nmodl_eigen_j = nmodl_eigen_jm.data();
double* nmodl_eigen_f = nmodl_eigen_fm.data();
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_;
nmodl_eigen_f[static_cast<int>(1)] = ( -nmodl_eigen_x[static_cast<int>(1)] + nt->_dt * (nmodl_eigen_x[static_cast<int>(0)] * kf0_ - nmodl_eigen_x[static_cast<int>(1)] * kb0_) + old_Y) / nt->_dt;
nmodl_eigen_j[static_cast<int>(1)] = kf0_;
nmodl_eigen_j[static_cast<int>(3)] = -kb0_ - 1.0 / nt->_dt;
nmodl_eigen_f[static_cast<int>(0)] = -nmodl_eigen_x[static_cast<int>(0)] * nt->_dt * kf0_ - nmodl_eigen_x[static_cast<int>(0)] + nmodl_eigen_x[static_cast<int>(1)] * nt->_dt * kb0_ + old_X;
nmodl_eigen_j[static_cast<int>(0)] = -nt->_dt * kf0_ - 1.0;
nmodl_eigen_j[static_cast<int>(2)] = nt->_dt * kb0_;
nmodl_eigen_f[static_cast<int>(1)] = nmodl_eigen_x[static_cast<int>(0)] * nt->_dt * kf0_ - nmodl_eigen_x[static_cast<int>(1)] * nt->_dt * kb0_ - nmodl_eigen_x[static_cast<int>(1)] + old_Y;
nmodl_eigen_j[static_cast<int>(1)] = nt->_dt * kf0_;
nmodl_eigen_j[static_cast<int>(3)] = -nt->_dt * kb0_ - 1.0;
}

void finalize() {
Expand Down
40 changes: 15 additions & 25 deletions kinetic/coreneuron/side_effects.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -257,22 +257,8 @@ namespace newton {
* @{
*/

static constexpr int MAX_ITER = 50;
static constexpr double EPS = 1e-13;

template <int N>
bool is_converged(const Eigen::Matrix<double, N, 1>& X,
const Eigen::Matrix<double, N, N>& J,
const Eigen::Matrix<double, N, 1>& F,
double eps) {
for (Eigen::Index i = 0; i < N; ++i) {
double square_error = J(i, Eigen::all).cwiseAbs2() * (eps * X).cwiseAbs2();
if (F(i) * F(i) > square_error) {
return false;
}
}
return true;
}
static constexpr int MAX_ITER = 1e3;
static constexpr double EPS = 1e-12;

/**
* \brief Newton method with user-provided Jacobian
Expand All @@ -295,14 +281,17 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
int max_iter = MAX_ITER) {
// Vector to store result of function F(X):
Eigen::Matrix<double, N, 1> F;
// Matrix to store Jacobian of F(X):
// Matrix to store jacobian of F(X):
Eigen::Matrix<double, N, N> J;
// Solver iteration count:
int iter = -1;
while (++iter < max_iter) {
// calculate F, J from X using user-supplied functor
functor(X, F, J);
if (is_converged(X, J, F, eps)) {
// get error norm: here we use sqrt(|F|^2)
double error = F.norm();
if (error < eps) {
// we have converged: return iteration count
return iter;
}
// In Eigen the default storage order is ColMajor.
Expand Down Expand Up @@ -343,7 +332,8 @@ EIGEN_DEVICE_FUNC int newton_solver_small_N(Eigen::Matrix<double, N, 1>& X,
int iter = -1;
while (++iter < max_iter) {
functor(X, F, J);
if (is_converged(X, J, F, eps)) {
double error = F.norm();
if (error < eps) {
return iter;
}
// The inverse can be called from within OpenACC regions without any issue, as opposed to
Expand Down Expand Up @@ -628,12 +618,12 @@ namespace coreneuron {
const double* nmodl_eigen_x = nmodl_eigen_xm.data();
double* nmodl_eigen_j = nmodl_eigen_jm.data();
double* nmodl_eigen_f = nmodl_eigen_fm.data();
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_;
nmodl_eigen_f[static_cast<int>(1)] = ( -nmodl_eigen_x[static_cast<int>(1)] + nt->_dt * (nmodl_eigen_x[static_cast<int>(0)] * kf0_ - nmodl_eigen_x[static_cast<int>(1)] * kb0_) + old_Y) / nt->_dt;
nmodl_eigen_j[static_cast<int>(1)] = kf0_;
nmodl_eigen_j[static_cast<int>(3)] = -kb0_ - 1.0 / nt->_dt;
nmodl_eigen_f[static_cast<int>(0)] = -nmodl_eigen_x[static_cast<int>(0)] * nt->_dt * kf0_ - nmodl_eigen_x[static_cast<int>(0)] + nmodl_eigen_x[static_cast<int>(1)] * nt->_dt * kb0_ + old_X;
nmodl_eigen_j[static_cast<int>(0)] = -nt->_dt * kf0_ - 1.0;
nmodl_eigen_j[static_cast<int>(2)] = nt->_dt * kb0_;
nmodl_eigen_f[static_cast<int>(1)] = nmodl_eigen_x[static_cast<int>(0)] * nt->_dt * kf0_ - nmodl_eigen_x[static_cast<int>(1)] * nt->_dt * kb0_ - nmodl_eigen_x[static_cast<int>(1)] + old_Y;
nmodl_eigen_j[static_cast<int>(1)] = nt->_dt * kf0_;
nmodl_eigen_j[static_cast<int>(3)] = -nt->_dt * kb0_ - 1.0;
}

void finalize() {
Expand Down
40 changes: 15 additions & 25 deletions kinetic/neuron/X2Y.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -248,22 +248,8 @@ namespace newton {
* @{
*/

static constexpr int MAX_ITER = 50;
static constexpr double EPS = 1e-13;

template <int N>
bool is_converged(const Eigen::Matrix<double, N, 1>& X,
const Eigen::Matrix<double, N, N>& J,
const Eigen::Matrix<double, N, 1>& F,
double eps) {
for (Eigen::Index i = 0; i < N; ++i) {
double square_error = J(i, Eigen::all).cwiseAbs2() * (eps * X).cwiseAbs2();
if (F(i) * F(i) > square_error) {
return false;
}
}
return true;
}
static constexpr int MAX_ITER = 1e3;
static constexpr double EPS = 1e-12;

/**
* \brief Newton method with user-provided Jacobian
Expand All @@ -286,14 +272,17 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
int max_iter = MAX_ITER) {
// Vector to store result of function F(X):
Eigen::Matrix<double, N, 1> F;
// Matrix to store Jacobian of F(X):
// Matrix to store jacobian of F(X):
Eigen::Matrix<double, N, N> J;
// Solver iteration count:
int iter = -1;
while (++iter < max_iter) {
// calculate F, J from X using user-supplied functor
functor(X, F, J);
if (is_converged(X, J, F, eps)) {
// get error norm: here we use sqrt(|F|^2)
double error = F.norm();
if (error < eps) {
// we have converged: return iteration count
return iter;
}
// In Eigen the default storage order is ColMajor.
Expand Down Expand Up @@ -334,7 +323,8 @@ EIGEN_DEVICE_FUNC int newton_solver_small_N(Eigen::Matrix<double, N, 1>& X,
int iter = -1;
while (++iter < max_iter) {
functor(X, F, J);
if (is_converged(X, J, F, eps)) {
double error = F.norm();
if (error < eps) {
return iter;
}
// The inverse can be called from within OpenACC regions without any issue, as opposed to
Expand Down Expand Up @@ -570,12 +560,12 @@ namespace neuron {
const double* nmodl_eigen_x = nmodl_eigen_xm.data();
double* nmodl_eigen_j = nmodl_eigen_jm.data();
double* nmodl_eigen_f = nmodl_eigen_fm.data();
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_;
nmodl_eigen_f[static_cast<int>(1)] = ( -nmodl_eigen_x[static_cast<int>(1)] + nt->_dt * (nmodl_eigen_x[static_cast<int>(0)] * kf0_ - nmodl_eigen_x[static_cast<int>(1)] * kb0_) + old_Y) / nt->_dt;
nmodl_eigen_j[static_cast<int>(1)] = kf0_;
nmodl_eigen_j[static_cast<int>(3)] = -kb0_ - 1.0 / nt->_dt;
nmodl_eigen_f[static_cast<int>(0)] = -nmodl_eigen_x[static_cast<int>(0)] * nt->_dt * kf0_ - nmodl_eigen_x[static_cast<int>(0)] + nmodl_eigen_x[static_cast<int>(1)] * nt->_dt * kb0_ + old_X;
nmodl_eigen_j[static_cast<int>(0)] = -nt->_dt * kf0_ - 1.0;
nmodl_eigen_j[static_cast<int>(2)] = nt->_dt * kb0_;
nmodl_eigen_f[static_cast<int>(1)] = nmodl_eigen_x[static_cast<int>(0)] * nt->_dt * kf0_ - nmodl_eigen_x[static_cast<int>(1)] * nt->_dt * kb0_ - nmodl_eigen_x[static_cast<int>(1)] + old_Y;
nmodl_eigen_j[static_cast<int>(1)] = nt->_dt * kf0_;
nmodl_eigen_j[static_cast<int>(3)] = -nt->_dt * kb0_ - 1.0;
}

void finalize() {
Expand Down
40 changes: 15 additions & 25 deletions kinetic/neuron/side_effects.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -248,22 +248,8 @@ namespace newton {
* @{
*/

static constexpr int MAX_ITER = 50;
static constexpr double EPS = 1e-13;

template <int N>
bool is_converged(const Eigen::Matrix<double, N, 1>& X,
const Eigen::Matrix<double, N, N>& J,
const Eigen::Matrix<double, N, 1>& F,
double eps) {
for (Eigen::Index i = 0; i < N; ++i) {
double square_error = J(i, Eigen::all).cwiseAbs2() * (eps * X).cwiseAbs2();
if (F(i) * F(i) > square_error) {
return false;
}
}
return true;
}
static constexpr int MAX_ITER = 1e3;
static constexpr double EPS = 1e-12;

/**
* \brief Newton method with user-provided Jacobian
Expand All @@ -286,14 +272,17 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
int max_iter = MAX_ITER) {
// Vector to store result of function F(X):
Eigen::Matrix<double, N, 1> F;
// Matrix to store Jacobian of F(X):
// Matrix to store jacobian of F(X):
Eigen::Matrix<double, N, N> J;
// Solver iteration count:
int iter = -1;
while (++iter < max_iter) {
// calculate F, J from X using user-supplied functor
functor(X, F, J);
if (is_converged(X, J, F, eps)) {
// get error norm: here we use sqrt(|F|^2)
double error = F.norm();
if (error < eps) {
// we have converged: return iteration count
return iter;
}
// In Eigen the default storage order is ColMajor.
Expand Down Expand Up @@ -334,7 +323,8 @@ EIGEN_DEVICE_FUNC int newton_solver_small_N(Eigen::Matrix<double, N, 1>& X,
int iter = -1;
while (++iter < max_iter) {
functor(X, F, J);
if (is_converged(X, J, F, eps)) {
double error = F.norm();
if (error < eps) {
return iter;
}
// The inverse can be called from within OpenACC regions without any issue, as opposed to
Expand Down Expand Up @@ -571,12 +561,12 @@ namespace neuron {
const double* nmodl_eigen_x = nmodl_eigen_xm.data();
double* nmodl_eigen_j = nmodl_eigen_jm.data();
double* nmodl_eigen_f = nmodl_eigen_fm.data();
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_;
nmodl_eigen_f[static_cast<int>(1)] = ( -nmodl_eigen_x[static_cast<int>(1)] + nt->_dt * (nmodl_eigen_x[static_cast<int>(0)] * kf0_ - nmodl_eigen_x[static_cast<int>(1)] * kb0_) + old_Y) / nt->_dt;
nmodl_eigen_j[static_cast<int>(1)] = kf0_;
nmodl_eigen_j[static_cast<int>(3)] = -kb0_ - 1.0 / nt->_dt;
nmodl_eigen_f[static_cast<int>(0)] = -nmodl_eigen_x[static_cast<int>(0)] * nt->_dt * kf0_ - nmodl_eigen_x[static_cast<int>(0)] + nmodl_eigen_x[static_cast<int>(1)] * nt->_dt * kb0_ + old_X;
nmodl_eigen_j[static_cast<int>(0)] = -nt->_dt * kf0_ - 1.0;
nmodl_eigen_j[static_cast<int>(2)] = nt->_dt * kb0_;
nmodl_eigen_f[static_cast<int>(1)] = nmodl_eigen_x[static_cast<int>(0)] * nt->_dt * kf0_ - nmodl_eigen_x[static_cast<int>(1)] * nt->_dt * kb0_ - nmodl_eigen_x[static_cast<int>(1)] + old_Y;
nmodl_eigen_j[static_cast<int>(1)] = nt->_dt * kf0_;
nmodl_eigen_j[static_cast<int>(3)] = -nt->_dt * kb0_ - 1.0;
}

void finalize() {
Expand Down
Loading

0 comments on commit 6a09015

Please sign in to comment.