Skip to content

Commit

Permalink
Simplify usage of OcMatrix
Browse files Browse the repository at this point in the history
  • Loading branch information
Nicolas Cornu committed Oct 29, 2024
1 parent 29322f2 commit adfdd4e
Show file tree
Hide file tree
Showing 4 changed files with 38 additions and 31 deletions.
5 changes: 5 additions & 0 deletions src/ivoc/ocmatrix.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,11 @@ class OcMatrix {
unimp();
return nullptr;
} // matrix element pointer

inline double operator()(int i, int j) const {
return getval(i, j);
};

inline double& operator()(int i, int j) {
return *mep(i, j);
};
Expand Down
11 changes: 4 additions & 7 deletions src/nrniv/matrixmap.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,13 +28,10 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) {
mmfree();

plen_ = 0;
std::vector<int> nonzero_i, nonzero_j;
m_.nonzeros(nonzero_i, nonzero_j);
pm_.resize(nonzero_i.size());
ptree_.resize(nonzero_i.size());
for (int k = 0; k < nonzero_i.size(); k++) {
const int i = nonzero_i[k];
const int j = nonzero_j[k];
std::vector<std::pair<int, int>> nzs = m_.nonzeros();
pm_.resize(nzs.size());
ptree_.resize(nzs.size());
for (const auto [i, j]: nzs) {
int it;
if (i < nnode) {
it = nodes[i]->eqn_index_ + layer[i];
Expand Down
48 changes: 24 additions & 24 deletions src/nrnpython/rxd_extracellular.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -358,35 +358,35 @@ void* ecs_do_reactions(void* dataptr) {

for (k = 0; k < react->num_species_involved; k++) {
pd = (results_array_dx[k] - results_array[k]) / dx;
*jacobian->mep(k, j) = (j == k) - dt * pd;
jacobian(k, j) = (j == k) - dt * pd;

Check failure on line 361 in src/nrnpython/rxd_extracellular.cpp

View workflow job for this annotation

GitHub Actions / ubuntu-22.04 - cmake (-DNRN_ENABLE_CORENEURON=ON -DNRN_ENABLE_INTERVIEWS=OFF -DNMODL_SANITIZERS=undefinedundefined) 

type 'std::unique_ptr<OcFullMatrix>' does not provide a call operator
}
states_cache_dx[j] -= dx;
}
// solve for x
if (react->num_species_involved == 1) {
react->species_states[0][i] += b[0] / jacobian->getval(0, 0);
react->species_states[0][i] += b[0] / jacobian(0, 0);

Check failure on line 367 in src/nrnpython/rxd_extracellular.cpp

View workflow job for this annotation

GitHub Actions / ubuntu-22.04 - cmake (-DNRN_ENABLE_CORENEURON=ON -DNRN_ENABLE_INTERVIEWS=OFF -DNMODL_SANITIZERS=undefinedundefined) 

type 'std::unique_ptr<OcFullMatrix>' does not provide a call operator
} else {
// find entry in leftmost column with largest absolute value
// Pivot
for (j = 0; j < react->num_species_involved; j++) {
for (k = j + 1; k < react->num_species_involved; k++) {
if (abs(jacobian->getval(j, j)) < abs(jacobian->getval(k, j))) {
if (abs(jacobian(j, j)) < abs(jacobian(k, j))) {

Check failure on line 373 in src/nrnpython/rxd_extracellular.cpp

View workflow job for this annotation

GitHub Actions / ubuntu-22.04 - cmake (-DNRN_ENABLE_CORENEURON=ON -DNRN_ENABLE_INTERVIEWS=OFF -DNMODL_SANITIZERS=undefinedundefined) 

type 'std::unique_ptr<OcFullMatrix>' does not provide a call operator

Check failure on line 373 in src/nrnpython/rxd_extracellular.cpp

View workflow job for this annotation

GitHub Actions / ubuntu-22.04 - cmake (-DNRN_ENABLE_CORENEURON=ON -DNRN_ENABLE_INTERVIEWS=OFF -DNMODL_SANITIZERS=undefinedundefined) 

type 'std::unique_ptr<OcFullMatrix>' does not provide a call operator
for (n = 0; n < react->num_species_involved; n++) {
temp = jacobian->getval(j, n);
*jacobian->mep(j, n) = jacobian->getval(k, n);
*jacobian->mep(k, n) = temp;
temp = jacobian(j, n);

Check failure on line 375 in src/nrnpython/rxd_extracellular.cpp

View workflow job for this annotation

GitHub Actions / ubuntu-22.04 - cmake (-DNRN_ENABLE_CORENEURON=ON -DNRN_ENABLE_INTERVIEWS=OFF -DNMODL_SANITIZERS=undefinedundefined) 

type 'std::unique_ptr<OcFullMatrix>' does not provide a call operator
jacobian(j, n) = jacobian(k, n);

Check failure on line 376 in src/nrnpython/rxd_extracellular.cpp

View workflow job for this annotation

GitHub Actions / ubuntu-22.04 - cmake (-DNRN_ENABLE_CORENEURON=ON -DNRN_ENABLE_INTERVIEWS=OFF -DNMODL_SANITIZERS=undefinedundefined) 

type 'std::unique_ptr<OcFullMatrix>' does not provide a call operator

Check failure on line 376 in src/nrnpython/rxd_extracellular.cpp

View workflow job for this annotation

GitHub Actions / ubuntu-22.04 - cmake (-DNRN_ENABLE_CORENEURON=ON -DNRN_ENABLE_INTERVIEWS=OFF -DNMODL_SANITIZERS=undefinedundefined) 

type 'std::unique_ptr<OcFullMatrix>' does not provide a call operator
jacobian(k, n) = temp;

Check failure on line 377 in src/nrnpython/rxd_extracellular.cpp

View workflow job for this annotation

GitHub Actions / ubuntu-22.04 - cmake (-DNRN_ENABLE_CORENEURON=ON -DNRN_ENABLE_INTERVIEWS=OFF -DNMODL_SANITIZERS=undefinedundefined) 

type 'std::unique_ptr<OcFullMatrix>' does not provide a call operator
}
}
}
}

for (j = 0; j < react->num_species_involved - 1; j++) {
for (k = j + 1; k < react->num_species_involved; k++) {
ge_value = jacobian->getval(k, j) / jacobian->getval(j, j);
ge_value = jacobian(k, j) / jacobian(j, j);

Check failure on line 385 in src/nrnpython/rxd_extracellular.cpp

View workflow job for this annotation

GitHub Actions / ubuntu-22.04 - cmake (-DNRN_ENABLE_CORENEURON=ON -DNRN_ENABLE_INTERVIEWS=OFF -DNMODL_SANITIZERS=undefinedundefined) 

type 'std::unique_ptr<OcFullMatrix>' does not provide a call operator

Check failure on line 385 in src/nrnpython/rxd_extracellular.cpp

View workflow job for this annotation

GitHub Actions / ubuntu-22.04 - cmake (-DNRN_ENABLE_CORENEURON=ON -DNRN_ENABLE_INTERVIEWS=OFF -DNMODL_SANITIZERS=undefinedundefined) 

type 'std::unique_ptr<OcFullMatrix>' does not provide a call operator
for (n = 0; n < react->num_species_involved; n++) {
val_to_set = jacobian->getval(k, n) -
ge_value * jacobian->getval(j, n);
*jacobian->mep(k, n) = val_to_set;
val_to_set = jacobian(k, n) -
ge_value * jacobian(j, n);
jacobian(k, n) = val_to_set;
}
b[k] = b[k] - ge_value * b[j];
}
Expand All @@ -396,10 +396,10 @@ void* ecs_do_reactions(void* dataptr) {
x[j] = b[j];
for (k = j + 1; k < react->num_species_involved; k++) {
if (k != j) {
x[j] = x[j] - jacobian->getval(j, k) * x[k];
x[j] = x[j] - jacobian(j, k) * x[k];
}
}
x[j] = x[j] / jacobian->getval(j, j);
x[j] = x[j] / jacobian(j, j);
}
for (j = 0; j < react->num_species_involved; j++) {
// I think this should be something like
Expand Down Expand Up @@ -480,35 +480,35 @@ void* ecs_do_reactions(void* dataptr) {

for (k = 0; k < react->num_species_involved; k++) {
pd = (results_array_dx[k] - results_array[k]) / dx;
*jacobian->mep(k, j) = (j == k) - dt * pd;
jacobian(k, j) = (j == k) - dt * pd;
}
states_cache_dx[j] -= dx;
}
// solve for x
if (react->num_species_involved == 1) {
react->species_states[0][i] += b[0] / jacobian->getval(0, 0);
react->species_states[0][i] += b[0] / jacobian(0, 0);
} else {
// find entry in leftmost column with largest absolute value
// Pivot
for (j = 0; j < react->num_species_involved; j++) {
for (k = j + 1; k < react->num_species_involved; k++) {
if (abs(jacobian->getval(j, j)) < abs(jacobian->getval(k, j))) {
if (abs(jacobian(j, j)) < abs(jacobian(k, j))) {
for (n = 0; n < react->num_species_involved; n++) {
temp = jacobian->getval(j, n);
*jacobian->mep(j, n) = jacobian->getval(k, n);
*jacobian->mep(k, n) = temp;
temp = jacobian(j, n);
jacobian(j, n) = jacobian(k, n);
jacobian(k, n) = temp;
}
}
}
}

for (j = 0; j < react->num_species_involved - 1; j++) {
for (k = j + 1; k < react->num_species_involved; k++) {
ge_value = jacobian->getval(k, j) / jacobian->getval(j, j);
ge_value = jacobian(k, j) / jacobian(j, j);
for (n = 0; n < react->num_species_involved; n++) {
val_to_set = jacobian->getval(k, n) -
ge_value * jacobian->getval(j, n);
*jacobian->mep(k, n) = val_to_set;
val_to_set = jacobian(k, n) -
ge_value * jacobian(j, n);
jacobian(k, n) = val_to_set;
}
b[k] = b[k] - ge_value * b[j];
}
Expand All @@ -518,10 +518,10 @@ void* ecs_do_reactions(void* dataptr) {
x[j] = b[j];
for (k = j + 1; k < react->num_species_involved; k++) {
if (k != j) {
x[j] = x[j] - jacobian->getval(j, k) * x[k];
x[j] = x[j] - jacobian(j, k) * x[k];
}
}
x[j] = x[j] / jacobian->getval(j, j);
x[j] = x[j] / jacobian(j, j);
}
for (j = 0; j < react->num_species_involved; j++) {
// I think this should be something like
Expand Down
5 changes: 5 additions & 0 deletions test/unit_tests/matrix.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,6 +94,11 @@ SCENARIO("A Matrix", "[neuron_ivoc][OcMatrix]") {
m.setrow(3, 2.0);
REQUIRE(compareMatrix(m, {{3., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}, {2., 2., 2.}}));
}
{
std::vector<std::pair<int, int>> nzs = m.nonzeros();
std::vector<std::pair<int, int>> res = {{0, 0}, {1, 1} , {2, 2}, {3, 0}, {3, 1}, {3, 2}};
REQUIRE(nzs == res);
}
{
std::vector<int> x, y;
m.nonzeros(x, y);
Expand Down

0 comments on commit adfdd4e

Please sign in to comment.