From adfdd4ee71ea12681b954cc9d3d1d409cec27ba5 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 29 Oct 2024 16:28:56 +0100 Subject: [PATCH] Simplify usage of OcMatrix --- src/ivoc/ocmatrix.h | 5 +++ src/nrniv/matrixmap.cpp | 11 +++---- src/nrnpython/rxd_extracellular.cpp | 48 ++++++++++++++--------------- test/unit_tests/matrix.cpp | 5 +++ 4 files changed, 38 insertions(+), 31 deletions(-) diff --git a/src/ivoc/ocmatrix.h b/src/ivoc/ocmatrix.h index 3bcdb87975..17df1d3713 100644 --- a/src/ivoc/ocmatrix.h +++ b/src/ivoc/ocmatrix.h @@ -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); }; diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp index f07cfbed2b..05673ff858 100644 --- a/src/nrniv/matrixmap.cpp +++ b/src/nrniv/matrixmap.cpp @@ -28,13 +28,10 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { mmfree(); plen_ = 0; - std::vector 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> 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]; diff --git a/src/nrnpython/rxd_extracellular.cpp b/src/nrnpython/rxd_extracellular.cpp index aeb377eed4..5bb0ec8f4c 100644 --- a/src/nrnpython/rxd_extracellular.cpp +++ b/src/nrnpython/rxd_extracellular.cpp @@ -358,23 +358,23 @@ 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; } } } @@ -382,11 +382,11 @@ void* ecs_do_reactions(void* dataptr) { 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]; } @@ -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 @@ -480,23 +480,23 @@ 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; } } } @@ -504,11 +504,11 @@ void* ecs_do_reactions(void* dataptr) { 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]; } @@ -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 diff --git a/test/unit_tests/matrix.cpp b/test/unit_tests/matrix.cpp index 27386501ca..115dbfe120 100644 --- a/test/unit_tests/matrix.cpp +++ b/test/unit_tests/matrix.cpp @@ -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> nzs = m.nonzeros(); + std::vector> res = {{0, 0}, {1, 1} , {2, 2}, {3, 0}, {3, 1}, {3, 2}}; + REQUIRE(nzs == res); + } { std::vector x, y; m.nonzeros(x, y);