From 50866931374ebdf8e24c71d33aef0a5ba1f79bcc Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 22 Oct 2024 18:38:28 +0200 Subject: [PATCH 01/23] treeset / sparse13: try minimal --- src/nrncvode/occvode.cpp | 2 +- src/nrniv/matrixmap.cpp | 7 ++- src/nrnoc/fadvance.cpp | 12 ++--- src/nrnoc/multicore.cpp | 6 +-- src/nrnoc/multicore.h | 3 +- src/nrnoc/solve.cpp | 23 +++++----- src/nrnoc/treeset.cpp | 97 +++++++++++++++++++++++++--------------- 7 files changed, 90 insertions(+), 60 deletions(-) diff --git a/src/nrncvode/occvode.cpp b/src/nrncvode/occvode.cpp index a2f6bd24eb..e9f7f2f18f 100644 --- a/src/nrncvode/occvode.cpp +++ b/src/nrncvode/occvode.cpp @@ -349,7 +349,7 @@ void Cvode::daspk_init_eqn() { if (use_sparse13 == 0 || diam_changed != 0) { recalc_diam(); } - zneq = spGetSize(_nt->_sp13mat, 0); + zneq = _nt->_sp13mat->ncol(); z.neq_v_ = z.nonvint_offset_ = zneq; // now add the membrane mechanism ode's to the count for (cml = z.cv_memb_list_; cml; cml = cml->next) { diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp index f07cfbed2b..12b608e529 100644 --- a/src/nrniv/matrixmap.cpp +++ b/src/nrniv/matrixmap.cpp @@ -24,6 +24,7 @@ void MatrixMap::add(double fac) { } void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { + static double place_holder = 0; NrnThread* _nt = nrn_threads; mmfree(); @@ -54,7 +55,11 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { } else { jt = start + j - nnode; } - ptree_[plen_] = spGetElement(_nt->_sp13mat, it, jt); + if (it == 0 || jt == 0) { + ptree_[plen_] = &place_holder; + } else { + ptree_[plen_] = _nt->_sp13mat->mep(it - 1, jt - 1); + } ++plen_; } } diff --git a/src/nrnoc/fadvance.cpp b/src/nrnoc/fadvance.cpp index a8a06f314e..f51d0b9958 100644 --- a/src/nrnoc/fadvance.cpp +++ b/src/nrnoc/fadvance.cpp @@ -672,13 +672,13 @@ void nrn_print_matrix(NrnThread* _nt) { Node* nd; if (use_sparse13) { if (ifarg(1) && chkarg(1, 0., 1.) == 0.) { - spPrint(_nt->_sp13mat, 1, 0, 1); + // spPrint(_nt->_sp13mat, 1, 0, 1); } else { - int i, n = spGetSize(_nt->_sp13mat, 0); - spPrint(_nt->_sp13mat, 1, 1, 1); - for (i = 1; i <= n; ++i) { - Printf("%d %g\n", i, _nt->actual_rhs(i)); - } + // int i, n = spGetSize(_nt->_sp13mat, 0); + // spPrint(_nt->_sp13mat, 1, 1, 1); + // for (i = 1; i <= n; ++i) { + // Printf("%d %g\n", i, _nt->actual_rhs(i)); + // } } } else if (_nt) { for (inode = 0; inode < _nt->end; ++inode) { diff --git a/src/nrnoc/multicore.cpp b/src/nrnoc/multicore.cpp index b443e6c4b7..a4a97c75e4 100644 --- a/src/nrnoc/multicore.cpp +++ b/src/nrnoc/multicore.cpp @@ -345,7 +345,7 @@ void nrn_threads_create(int n, bool parallel) { nt->_ecell_memb_list = 0; nt->_ecell_child_cnt = 0; nt->_ecell_children = NULL; - nt->_sp13mat = 0; + nt->_sp13mat = nullptr; nt->_ctime = 0.0; nt->_vcv = 0; nt->_node_data_offset = 0; @@ -440,8 +440,8 @@ void nrn_threads_free() { nt->_ecell_children = NULL; } if (nt->_sp13mat) { - spDestroy(nt->_sp13mat); - nt->_sp13mat = 0; + delete(nt->_sp13mat); + nt->_sp13mat = nullptr; } nt->end = 0; nt->ncell = 0; diff --git a/src/nrnoc/multicore.h b/src/nrnoc/multicore.h index fb115b195c..ae9fe61e5e 100644 --- a/src/nrnoc/multicore.h +++ b/src/nrnoc/multicore.h @@ -29,6 +29,7 @@ actual_v, etc. #include "membfunc.h" #include +#include typedef struct NrnThreadMembList { /* patterned after CvMembList in cvodeobj.h */ struct NrnThreadMembList* next; @@ -91,7 +92,7 @@ struct NrnThread { Node** _v_parent; double* _sp13_rhs; /* rhs matrix for sparse13 solver. updates to and from this vector need to be transfered to actual_rhs */ - char* _sp13mat; /* handle to general sparse matrix */ + OcSparseMatrix* _sp13mat; /* handle to general sparse matrix */ Memb_list* _ecell_memb_list; /* normally nullptr */ Node** _ecell_children; /* nodes with no extcell but parent has it */ void* _vcv; /* replaces old cvode_instance and nrn_cvode_ */ diff --git a/src/nrnoc/solve.cpp b/src/nrnoc/solve.cpp index 38aa3ebc53..3853d41df0 100644 --- a/src/nrnoc/solve.cpp +++ b/src/nrnoc/solve.cpp @@ -59,6 +59,7 @@ node.v + extnode.v[0] #include "section.h" #include "spmatrix.h" #include "treeset.h" +#include "ivocvect.h" #include #include @@ -367,19 +368,17 @@ void nrn_solve(NrnThread* _nt) { int e; nrn_thread_error("solve use_sparse13"); update_sp13_mat_based_on_actual_d(_nt); - e = spFactor(_nt->_sp13mat); - if (e != spOKAY) { - switch (e) { - case spZERO_DIAG: - hoc_execerror("spFactor error:", "Zero Diagonal"); - case spNO_MEMORY: - hoc_execerror("spFactor error:", "No Memory"); - case spSINGULAR: - hoc_execerror("spFactor error:", "Singular"); - } - } + update_sp13_rhs_based_on_actual_rhs(_nt); - spSolve(_nt->_sp13mat, _nt->_sp13_rhs, _nt->_sp13_rhs); + int size = _nt->_sp13mat->ncol(); + std::vector vector_in(_nt->_sp13_rhs + 1, _nt->_sp13_rhs + size + 1); + IvocVect in; in.vec() = vector_in; + IvocVect out(size); + _nt->_sp13mat->solv(&in, &out, false); + for (int i = 0; i < size; ++i) { + _nt->_sp13_rhs[i + 1] = out[i]; + + } update_actual_d_based_on_sp13_mat(_nt); update_actual_rhs_based_on_sp13_rhs(_nt); } else { diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index 2a1c0abf97..51f12081df 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -394,7 +394,7 @@ void nrn_rhs(neuron::model_sorted_token const& cache_token, NrnThread& nt) { if (use_sparse13) { int i, neqn; nrn_thread_error("nrn_rhs use_sparse13"); - neqn = spGetSize(_nt->_sp13mat, 0); + neqn = _nt->_sp13mat->nrow(); for (i = 1; i <= neqn; ++i) { _nt->_sp13_rhs[i] = 0.; } @@ -500,7 +500,7 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) { if (use_sparse13) { // Zero the sparse13 matrix - spClear(_nt->_sp13mat); + _nt->_sp13mat->zero(); } // Make sure the SoA node diagonals are also zeroed (is this needed?) @@ -1789,8 +1789,8 @@ void nrn_matrix_node_free() { free(std::exchange(nt->_sp13_rhs, nullptr)); } if (nt->_sp13mat) { - spDestroy(nt->_sp13mat); - nt->_sp13mat = (char*) 0; + delete nt->_sp13mat; + nt->_sp13mat = nullptr; } } diam_changed = 1; @@ -1862,7 +1862,7 @@ and therefore is passed to spSolve as actual_rhs intead of actual_rhs-1. */ static void nrn_matrix_node_alloc(void) { - int i, b; + int b; Node* nd; NrnThread* nt; @@ -1887,55 +1887,80 @@ static void nrn_matrix_node_alloc(void) { } ++nrn_matrix_cnt_; if (use_sparse13) { - int in, err, extn, neqn, j; nt = nrn_threads; - neqn = nt->end + nrndae_extra_eqn_count(); - extn = 0; + int neqn = nt->end + nrndae_extra_eqn_count(); + int extn = 0; if (nt->_ecell_memb_list) { extn = nt->_ecell_memb_list->nodecount * nlayer; } /*printf(" %d extracellular nodes\n", extn);*/ neqn += extn; nt->_sp13_rhs = (double*) ecalloc(neqn + 1, sizeof(double)); - nt->_sp13mat = spCreate(neqn, 0, &err); - if (err != spOKAY) { - hoc_execerror("Couldn't create sparse matrix", (char*) 0); - } - for (in = 0, i = 1; in < nt->end; ++in, ++i) { + nt->_sp13mat = new OcSparseMatrix(neqn, neqn); + for (int in = 0, i = 1; in < nt->end; ++in, ++i) { nt->_v_node[in]->eqn_index_ = i; if (nt->_v_node[in]->extnode) { i += nlayer; } } - for (in = 0; in < nt->end; ++in) { - int ie, k; - Node *nd, *pnd; - Extnode* nde; - nd = nt->_v_node[in]; - nde = nd->extnode; - pnd = nt->_v_parent[in]; - i = nd->eqn_index_; + // Creating + for (int in = 0; in < nt->end; ++in) { + Node *nd = nt->_v_node[in]; + Extnode* nde = nd->extnode; + Node *pnd = nt->_v_parent[in]; + int i = nd->eqn_index_; + nt->_sp13mat->mep(i - 1, i - 1); + if (nde) { + for (int ie = 0; ie < nlayer; ++ie) { + int k = i + ie; + nt->_sp13mat->mep(k, k); + nt->_sp13mat->mep(k, k - 1); + nt->_sp13mat->mep(k - 1, k); + } + } + if (pnd) { + int j = pnd->eqn_index_; + nt->_sp13mat->mep(j - 1, i - 1); + nt->_sp13mat->mep(i - 1, j - 1); + if (nde && pnd->extnode) + for (int ie = 0; ie < nlayer; ++ie) { + int kp = j + ie; + int k = i + ie; + nt->_sp13mat->mep(kp, k); + nt->_sp13mat->mep(k, kp); + } + } else { /* not needed if index starts at 1 */ + nd->_a_matelm = nullptr; + nd->_b_matelm = nullptr; + } + } + // Collecting + for (int in = 0; in < nt->end; ++in) { + Node *nd = nt->_v_node[in]; + Extnode* nde = nd->extnode; + Node *pnd = nt->_v_parent[in]; + int i = nd->eqn_index_; nt->_sp13_rhs[i] = nt->actual_rhs(in); - nd->_d_matelm = spGetElement(nt->_sp13mat, i, i); + nd->_d_matelm = nt->_sp13mat->mep(i - 1, i - 1); if (nde) { - for (ie = 0; ie < nlayer; ++ie) { - k = i + ie + 1; - nde->_d[ie] = spGetElement(nt->_sp13mat, k, k); - nde->_rhs[ie] = nt->_sp13_rhs + k; - nde->_x21[ie] = spGetElement(nt->_sp13mat, k, k - 1); - nde->_x12[ie] = spGetElement(nt->_sp13mat, k - 1, k); + for (int ie = 0; ie < nlayer; ++ie) { + int k = i + ie; + nde->_d[ie] = nt->_sp13mat->mep(k, k); + nde->_rhs[ie] = nt->_sp13_rhs + k + 1; + nde->_x21[ie] = nt->_sp13mat->mep(k, k - 1); + nde->_x12[ie] = nt->_sp13mat->mep(k - 1, k); } } if (pnd) { - j = pnd->eqn_index_; - nd->_a_matelm = spGetElement(nt->_sp13mat, j, i); - nd->_b_matelm = spGetElement(nt->_sp13mat, i, j); + int j = pnd->eqn_index_; + nd->_a_matelm = nt->_sp13mat->mep(j - 1, i - 1); + nd->_b_matelm = nt->_sp13mat->mep(i - 1, j - 1); if (nde && pnd->extnode) - for (ie = 0; ie < nlayer; ++ie) { - int kp = j + ie + 1; - k = i + ie + 1; - nde->_a_matelm[ie] = spGetElement(nt->_sp13mat, kp, k); - nde->_b_matelm[ie] = spGetElement(nt->_sp13mat, k, kp); + for (int ie = 0; ie < nlayer; ++ie) { + int kp = j + ie; + int k = i + ie; + nde->_a_matelm[ie] = nt->_sp13mat->mep(kp, k); + nde->_b_matelm[ie] = nt->_sp13mat->mep(k, kp); } } else { /* not needed if index starts at 1 */ nd->_a_matelm = nullptr; From 1168bfdd92b6fd694c290da86a47a3d94ddecb4c Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 23 Oct 2024 19:44:34 +0200 Subject: [PATCH 02/23] Compress the matrix to keep the value stable --- src/ivoc/ocmatrix.h | 4 ++++ src/nrnoc/treeset.cpp | 1 + 2 files changed, 5 insertions(+) diff --git a/src/ivoc/ocmatrix.h b/src/ivoc/ocmatrix.h index 15960fe51f..a3ab1b40c1 100644 --- a/src/ivoc/ocmatrix.h +++ b/src/ivoc/ocmatrix.h @@ -219,6 +219,10 @@ class OcSparseMatrix final: public OcMatrix { // type 2 void zero() override; + void compress() { + m_.makeCompressed(); + } + private: Eigen::SparseMatrix m_{}; std::unique_ptr> lu_{}; diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index 51f12081df..e9a7ea6768 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -1934,6 +1934,7 @@ static void nrn_matrix_node_alloc(void) { nd->_b_matelm = nullptr; } } + nt->_sp13mat->compress(); // Collecting for (int in = 0; in < nt->end; ++in) { Node *nd = nt->_v_node[in]; From f49d7abf0a4e6046dccad6668ee6662a33e01a91 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 24 Oct 2024 10:19:12 +0200 Subject: [PATCH 03/23] forward definition because we cannot include this file to install --- src/nrnoc/multicore.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/nrnoc/multicore.h b/src/nrnoc/multicore.h index ae9fe61e5e..ac8d27b718 100644 --- a/src/nrnoc/multicore.h +++ b/src/nrnoc/multicore.h @@ -29,7 +29,7 @@ actual_v, etc. #include "membfunc.h" #include -#include +class OcSparseMatrix; typedef struct NrnThreadMembList { /* patterned after CvMembList in cvodeobj.h */ struct NrnThreadMembList* next; From e90ede272673a18ef6c2686d5e3bf9dfd40d2b30 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 24 Oct 2024 10:27:14 +0200 Subject: [PATCH 04/23] Change includes --- src/nrncvode/occvode.cpp | 2 +- src/nrnoc/solve.cpp | 2 +- src/nrnoc/treeset.cpp | 4 +--- 3 files changed, 3 insertions(+), 5 deletions(-) diff --git a/src/nrncvode/occvode.cpp b/src/nrncvode/occvode.cpp index e9f7f2f18f..e81513a61c 100644 --- a/src/nrncvode/occvode.cpp +++ b/src/nrncvode/occvode.cpp @@ -18,7 +18,7 @@ #include -#include "spmatrix.h" +#include "ocmatrix.h" extern double* sp13mat; #if 1 || NRNMPI diff --git a/src/nrnoc/solve.cpp b/src/nrnoc/solve.cpp index 3853d41df0..8c4b0286af 100644 --- a/src/nrnoc/solve.cpp +++ b/src/nrnoc/solve.cpp @@ -57,7 +57,7 @@ node.v + extnode.v[0] #include "nrnmpiuse.h" #include "ocnotify.h" #include "section.h" -#include "spmatrix.h" +#include "ocmatrix.h" #include "treeset.h" #include "ivocvect.h" diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index e9a7ea6768..54fc83dbf9 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -18,7 +18,7 @@ #include "ocnotify.h" #include "partrans.h" #include "section.h" -#include "spmatrix.h" +#include "ocmatrix.h" #include "utils/profile/profiler_interface.h" #include "multicore.h" @@ -32,8 +32,6 @@ #include -extern spREAL* spGetElement(char*, int, int); - int nrn_shape_changed_; /* for notifying Shape class in nrniv */ double* nrn_mech_wtime_; From 5d4fff5ee1ed87cf8ba9fd2c1dd0bb4f2c6edcba Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 28 Oct 2024 10:35:26 +0100 Subject: [PATCH 05/23] Replace _d_matelm by direct access --- src/nrnoc/section.h | 1 - src/nrnoc/treeset.cpp | 23 +++++++++++++++++------ 2 files changed, 17 insertions(+), 7 deletions(-) diff --git a/src/nrnoc/section.h b/src/nrnoc/section.h index 1e26dc8824..f66cca844d 100644 --- a/src/nrnoc/section.h +++ b/src/nrnoc/section.h @@ -183,7 +183,6 @@ struct Node { double _rinv{}; /* conductance uS from node to parent */ double* _a_matelm; double* _b_matelm; - double* _d_matelm; int eqn_index_; /* sparse13 matrix row/col index */ /* if no extnodes then = v_node_index +1*/ /* each extnode adds nlayer more equations after this */ diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index 54fc83dbf9..10a8df88e1 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -21,6 +21,7 @@ #include "ocmatrix.h" #include "utils/profile/profiler_interface.h" #include "multicore.h" +#include "ocmatrix.h" #include #include @@ -351,7 +352,9 @@ void update_sp13_rhs_based_on_actual_rhs(NrnThread* nt) { */ void update_actual_d_based_on_sp13_mat(NrnThread* nt) { for (int i = 0; i < nt->end; ++i) { - nt->actual_d(i) = *nt->_v_node[i]->_d_matelm; + const int index = nt->_v_node[i]->eqn_index_; + OcSparseMatrix& m = *nt->_sp13mat; + nt->actual_d(i) = m.getval(index - 1, index - 1); } } @@ -360,7 +363,9 @@ void update_actual_d_based_on_sp13_mat(NrnThread* nt) { */ void update_sp13_mat_based_on_actual_d(NrnThread* nt) { for (int i = 0; i < nt->end; ++i) { - *nt->_v_node[i]->_d_matelm = nt->actual_d(i); + const int index = nt->_v_node[i]->eqn_index_; + OcSparseMatrix& m = *nt->_sp13mat; + *m.mep(index - 1, index - 1) = nt->actual_d(i); } } @@ -569,16 +574,23 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) { update_sp13_mat_based_on_actual_d(_nt); // just because of activclamp_lhs for (i = i2; i < i3; ++i) { // note i2 Node* nd = _nt->_v_node[i]; + OcSparseMatrix& m = *_nt->_sp13mat; auto const parent_i = _nt->_v_parent_index[i]; - auto* const parent_nd = _nt->_v_node[parent_i]; auto const nd_a = NODEA(nd); auto const nd_b = NODEB(nd); // Update entries in sp13_mat *nd->_a_matelm += nd_a; *nd->_b_matelm += nd_b; /* b may have value from lincir */ - *nd->_d_matelm -= nd_b; + { + const int index = nd->eqn_index_; + *m.mep(index - 1, index - 1) -= nd_b; + } // used to update NODED (sparse13 matrix) using NODEA and NODEB ("SoA") - *parent_nd->_d_matelm -= nd_a; + { + Node* const parent_nd = _nt->_v_node[parent_i]; + const int parent_index = parent_nd->eqn_index_; + *m.mep(parent_index - 1, parent_index - 1) -= nd_a; + } // Also update the Node's d value in the SoA storage (is this needed?) vec_d[i] -= nd_b; vec_d[parent_i] -= nd_a; @@ -1940,7 +1952,6 @@ static void nrn_matrix_node_alloc(void) { Node *pnd = nt->_v_parent[in]; int i = nd->eqn_index_; nt->_sp13_rhs[i] = nt->actual_rhs(in); - nd->_d_matelm = nt->_sp13mat->mep(i - 1, i - 1); if (nde) { for (int ie = 0; ie < nlayer; ++ie) { int k = i + ie; From f4b836dcb8a9dbac653a24ec68a060db93c3a43e Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 28 Oct 2024 11:43:05 +0100 Subject: [PATCH 06/23] Replace _[ab]_matelm by direct access --- src/nrnoc/extcelln.cpp | 17 ++++++++--------- src/nrnoc/section.h | 2 -- src/nrnoc/section_fwd.hpp | 2 -- src/nrnoc/treeset.cpp | 20 +++++--------------- 4 files changed, 13 insertions(+), 28 deletions(-) diff --git a/src/nrnoc/extcelln.cpp b/src/nrnoc/extcelln.cpp index b6facf2b76..c24dbaf07f 100644 --- a/src/nrnoc/extcelln.cpp +++ b/src/nrnoc/extcelln.cpp @@ -213,14 +213,12 @@ static void extcell_init(neuron::model_sorted_token const&, void extnode_free_elements(Extnode* nde) { if (nde->v) { free(nde->v); /* along with _a and _b */ - free(nde->_d); /* along with _rhs, _a_matelm, _b_matelm, _x12, and _x21 */ + free(nde->_d); /* along with _rhs, _x12, and _x21 */ nde->v = NULL; nde->_a = NULL; nde->_b = NULL; nde->_d = NULL; nde->_rhs = NULL; - nde->_a_matelm = NULL; - nde->_b_matelm = NULL; nde->_x12 = NULL; nde->_x21 = NULL; } @@ -294,11 +292,9 @@ static void extnode_alloc_elements(Extnode* nde) { nde->_a = nde->v + nlayer; nde->_b = nde->_a + nlayer; - nde->_d = (double**) ecalloc(nlayer * 6, sizeof(double*)); + nde->_d = (double**) ecalloc(nlayer * 4, sizeof(double*)); nde->_rhs = nde->_d + nlayer; - nde->_a_matelm = nde->_rhs + nlayer; - nde->_b_matelm = nde->_a_matelm + nlayer; - nde->_x12 = nde->_b_matelm + nlayer; + nde->_x12 = nde->_rhs + nlayer; nde->_x21 = nde->_x12 + nlayer; } } @@ -477,11 +473,14 @@ void nrn_setup_ext(NrnThread* _nt) { pnde = pnd->extnode; /* axial connections */ if (pnde) { /* parent sec may not be extracellular */ + OcSparseMatrix& m = *_nt->_sp13mat; + int index = nd->eqn_index_; + int parent_index = pnd->eqn_index_; for (j = 0; j < nrn_nlayer_extracellular; ++j) { *nde->_d[j] -= nde->_b[j]; *pnde->_d[j] -= nde->_a[j]; - *nde->_a_matelm[j] += nde->_a[j]; - *nde->_b_matelm[j] += nde->_b[j]; + *m.mep(parent_index + j, index + j) += nde->_a[j]; + *m.mep(index + j, parent_index + j) += nde->_b[j]; } } } diff --git a/src/nrnoc/section.h b/src/nrnoc/section.h index f66cca844d..760b2b8764 100644 --- a/src/nrnoc/section.h +++ b/src/nrnoc/section.h @@ -181,8 +181,6 @@ struct Node { return _node_handle.non_owning_handle(); } double _rinv{}; /* conductance uS from node to parent */ - double* _a_matelm; - double* _b_matelm; int eqn_index_; /* sparse13 matrix row/col index */ /* if no extnodes then = v_node_index +1*/ /* each extnode adds nlayer more equations after this */ diff --git a/src/nrnoc/section_fwd.hpp b/src/nrnoc/section_fwd.hpp index fe1a955f6c..b39576575b 100644 --- a/src/nrnoc/section_fwd.hpp +++ b/src/nrnoc/section_fwd.hpp @@ -42,8 +42,6 @@ struct Extnode { double* _b; double** _d; double** _rhs; /* d, rhs, a, and b are analogous to those in node */ - double** _a_matelm; - double** _b_matelm; double** _x12; /* effect of v[layer] on eqn layer-1 (or internal)*/ double** _x21; /* effect of v[layer-1 or internal] on eqn layer*/ }; diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index 10a8df88e1..907066b21d 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -579,17 +579,17 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) { auto const nd_a = NODEA(nd); auto const nd_b = NODEB(nd); // Update entries in sp13_mat - *nd->_a_matelm += nd_a; - *nd->_b_matelm += nd_b; /* b may have value from lincir */ { const int index = nd->eqn_index_; *m.mep(index - 1, index - 1) -= nd_b; - } - // used to update NODED (sparse13 matrix) using NODEA and NODEB ("SoA") - { + + // used to update NODED (sparse13 matrix) using NODEA and NODEB ("SoA") Node* const parent_nd = _nt->_v_node[parent_i]; const int parent_index = parent_nd->eqn_index_; *m.mep(parent_index - 1, parent_index - 1) -= nd_a; + + *m.mep(parent_index - 1, index - 1) += nd_a; + *m.mep(index - 1, parent_index - 1) += nd_b; /* b may have value from lincir */ } // Also update the Node's d value in the SoA storage (is this needed?) vec_d[i] -= nd_b; @@ -1939,9 +1939,6 @@ static void nrn_matrix_node_alloc(void) { nt->_sp13mat->mep(kp, k); nt->_sp13mat->mep(k, kp); } - } else { /* not needed if index starts at 1 */ - nd->_a_matelm = nullptr; - nd->_b_matelm = nullptr; } } nt->_sp13mat->compress(); @@ -1963,18 +1960,11 @@ static void nrn_matrix_node_alloc(void) { } if (pnd) { int j = pnd->eqn_index_; - nd->_a_matelm = nt->_sp13mat->mep(j - 1, i - 1); - nd->_b_matelm = nt->_sp13mat->mep(i - 1, j - 1); if (nde && pnd->extnode) for (int ie = 0; ie < nlayer; ++ie) { int kp = j + ie; int k = i + ie; - nde->_a_matelm[ie] = nt->_sp13mat->mep(kp, k); - nde->_b_matelm[ie] = nt->_sp13mat->mep(k, kp); } - } else { /* not needed if index starts at 1 */ - nd->_a_matelm = nullptr; - nd->_b_matelm = nullptr; } } nrndae_alloc(); From 9b2c6eb2a5e448ad9ee5f6ba81453c49b1847c60 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 28 Oct 2024 12:12:25 +0100 Subject: [PATCH 07/23] Replace _12 and _x21 by direct access --- src/nrnoc/extcelln.cpp | 50 ++++++++++++++++++--------------------- src/nrnoc/section_fwd.hpp | 2 -- src/nrnoc/treeset.cpp | 10 -------- 3 files changed, 23 insertions(+), 39 deletions(-) diff --git a/src/nrnoc/extcelln.cpp b/src/nrnoc/extcelln.cpp index c24dbaf07f..0292efd591 100644 --- a/src/nrnoc/extcelln.cpp +++ b/src/nrnoc/extcelln.cpp @@ -213,14 +213,12 @@ static void extcell_init(neuron::model_sorted_token const&, void extnode_free_elements(Extnode* nde) { if (nde->v) { free(nde->v); /* along with _a and _b */ - free(nde->_d); /* along with _rhs, _x12, and _x21 */ + free(nde->_d); /* along with _rhs */ nde->v = NULL; nde->_a = NULL; nde->_b = NULL; nde->_d = NULL; nde->_rhs = NULL; - nde->_x12 = NULL; - nde->_x21 = NULL; } } @@ -292,10 +290,8 @@ static void extnode_alloc_elements(Extnode* nde) { nde->_a = nde->v + nlayer; nde->_b = nde->_a + nlayer; - nde->_d = (double**) ecalloc(nlayer * 4, sizeof(double*)); + nde->_d = (double**) ecalloc(nlayer * 2, sizeof(double*)); nde->_rhs = nde->_d + nlayer; - nde->_x12 = nde->_rhs + nlayer; - nde->_x21 = nde->_x12 + nlayer; } } @@ -420,45 +416,45 @@ void nrn_rhs_ext(NrnThread* _nt) { } void nrn_setup_ext(NrnThread* _nt) { - int i, j, cnt; - Node *nd, *pnd, **ndlist; - double* pd; - double d, cfac, mfac; - Extnode *nde, *pnde; + double cfac, mfac; Memb_list* ml = _nt->_ecell_memb_list; if (!ml) { return; } /*printnode("begin setup");*/ - cnt = ml->nodecount; - ndlist = ml->nodelist; + int cnt = ml->nodecount; + Node** ndlist = ml->nodelist; cfac = .001 * _nt->cj; /* d contains all the membrane conductances (and capacitance) */ /* i.e. (cm/dt + di/dvm - dis/dvi)*[dvi] and (dis/dvi)*[dvx] */ - for (i = 0; i < cnt; ++i) { - nd = ndlist[i]; - nde = nd->extnode; - d = NODED(nd); + for (int i = 0; i < cnt; ++i) { + OcSparseMatrix& m = *_nt->_sp13mat; + Node* nd = ndlist[i]; + int index = nd->eqn_index_; + ExtNode* nde = nd->extnode; + double d = NODED(nd); /* nde->_d only has -ELECTRODE_CURRENT contribution */ d = (*nde->_d[0] += NODED(nd)); /* now d is only the membrane current contribution */ /* i.e. d = cm/dt + di/dvm */ - *nde->_x12[0] -= d; - *nde->_x21[0] -= d; + *m.mep(index - 1, index) -= d; + *m.mep(index, index - 1) -= d; #if I_MEMBRANE ml->data(i, sav_g_index) = d; #endif } /* series resistance, capacitance, and axial terms. */ - for (i = 0; i < cnt; ++i) { - nd = ndlist[i]; - nde = nd->extnode; - pnd = _nt->_v_parent[nd->v_node_index]; + for (int i = 0; i < cnt; ++i) { + OcSparseMatrix& m = *_nt->_sp13mat; + Node* nd = ndlist[i]; + int index = nd->eqn_index_; + ExtNode* nde = nd->extnode; + Node* pnd = _nt->_v_parent[nd->v_node_index]; if (pnd) { /* series resistance and capacitance to ground */ - j = 0; + int j = 0; for (;;) { /* between j and j+1 layer */ mfac = (*nde->param[xg_index_ext(j)] + *nde->param[xc_index_ext(j)] * cfac); *nde->_d[j] += mfac; @@ -467,10 +463,10 @@ void nrn_setup_ext(NrnThread* _nt) { break; } *nde->_d[j] += mfac; - *nde->_x12[j] -= mfac; - *nde->_x21[j] -= mfac; + *m.mep(index - 1 + j, index + j) -= mfac; + *m.mep(index + j, index - 1 + j) -= mfac; } - pnde = pnd->extnode; + ExtNode* pnde = pnd->extnode; /* axial connections */ if (pnde) { /* parent sec may not be extracellular */ OcSparseMatrix& m = *_nt->_sp13mat; diff --git a/src/nrnoc/section_fwd.hpp b/src/nrnoc/section_fwd.hpp index b39576575b..9e010014c1 100644 --- a/src/nrnoc/section_fwd.hpp +++ b/src/nrnoc/section_fwd.hpp @@ -42,8 +42,6 @@ struct Extnode { double* _b; double** _d; double** _rhs; /* d, rhs, a, and b are analogous to those in node */ - double** _x12; /* effect of v[layer] on eqn layer-1 (or internal)*/ - double** _x21; /* effect of v[layer-1 or internal] on eqn layer*/ }; #endif diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index 907066b21d..70570c519f 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -1954,18 +1954,8 @@ static void nrn_matrix_node_alloc(void) { int k = i + ie; nde->_d[ie] = nt->_sp13mat->mep(k, k); nde->_rhs[ie] = nt->_sp13_rhs + k + 1; - nde->_x21[ie] = nt->_sp13mat->mep(k, k - 1); - nde->_x12[ie] = nt->_sp13mat->mep(k - 1, k); } } - if (pnd) { - int j = pnd->eqn_index_; - if (nde && pnd->extnode) - for (int ie = 0; ie < nlayer; ++ie) { - int kp = j + ie; - int k = i + ie; - } - } } nrndae_alloc(); } else { From f2a1c37ae351d321fc7988af48c8eef6148d0b9a Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 28 Oct 2024 12:37:48 +0100 Subject: [PATCH 08/23] Replace _d by direct access --- src/nocmodl/noccout.cpp | 8 ++++++-- src/nrnoc/extcelln.cpp | 16 +++++++--------- src/nrnoc/section_fwd.hpp | 1 - src/nrnoc/treeset.cpp | 2 -- 4 files changed, 13 insertions(+), 14 deletions(-) diff --git a/src/nocmodl/noccout.cpp b/src/nocmodl/noccout.cpp index 31c83e2e73..f80f0ea553 100644 --- a/src/nocmodl/noccout.cpp +++ b/src/nocmodl/noccout.cpp @@ -260,7 +260,9 @@ void c_out() { P(" }\n"); P("#if EXTRACELLULAR\n"); P(" if (auto* const _extnode = _nrn_mechanism_access_extnode(_nd); _extnode) {\n"); - P(" *_extnode->_d[0] += _g;\n"); + P(" int index = _nd->eqn_index_;\n"); + P(" OcSparseMatrix& m = *_nt->_sp13mat;\n"); + P(" *m.mep(index, index) += _g;\n"); P(" }\n"); P("#endif\n"); } else { @@ -679,7 +681,9 @@ void c_out_vectorize() { P(" }\n"); P("#if EXTRACELLULAR\n"); P(" if (auto* const _extnode = _nrn_mechanism_access_extnode(_nd); _extnode) {\n"); - P(" *_extnode->_d[0] += _g;\n"); + P(" int index = _nd->eqn_index_;\n"); + P(" OcSparseMatrix& m = *_nt->_sp13mat;\n"); + P(" *m.mep(index, index) += _g;\n"); P(" }\n"); P("#endif\n"); } else { diff --git a/src/nrnoc/extcelln.cpp b/src/nrnoc/extcelln.cpp index 0292efd591..f73a86b42d 100644 --- a/src/nrnoc/extcelln.cpp +++ b/src/nrnoc/extcelln.cpp @@ -213,11 +213,9 @@ static void extcell_init(neuron::model_sorted_token const&, void extnode_free_elements(Extnode* nde) { if (nde->v) { free(nde->v); /* along with _a and _b */ - free(nde->_d); /* along with _rhs */ nde->v = NULL; nde->_a = NULL; nde->_b = NULL; - nde->_d = NULL; nde->_rhs = NULL; } } @@ -290,8 +288,7 @@ static void extnode_alloc_elements(Extnode* nde) { nde->_a = nde->v + nlayer; nde->_b = nde->_a + nlayer; - nde->_d = (double**) ecalloc(nlayer * 2, sizeof(double*)); - nde->_rhs = nde->_d + nlayer; + nde->_rhs = (double**) ecalloc(nlayer * 1, sizeof(double*)); } } @@ -436,7 +433,8 @@ void nrn_setup_ext(NrnThread* _nt) { ExtNode* nde = nd->extnode; double d = NODED(nd); /* nde->_d only has -ELECTRODE_CURRENT contribution */ - d = (*nde->_d[0] += NODED(nd)); + *m.mep(index, index) += NODED(nd); + d = m.getval(index, index); /* now d is only the membrane current contribution */ /* i.e. d = cm/dt + di/dvm */ *m.mep(index - 1, index) -= d; @@ -457,12 +455,12 @@ void nrn_setup_ext(NrnThread* _nt) { int j = 0; for (;;) { /* between j and j+1 layer */ mfac = (*nde->param[xg_index_ext(j)] + *nde->param[xc_index_ext(j)] * cfac); - *nde->_d[j] += mfac; + *m.mep(index + j, index + j) += mfac; ++j; if (j == nrn_nlayer_extracellular) { break; } - *nde->_d[j] += mfac; + *m.mep(index + j, index + j) += mfac; *m.mep(index - 1 + j, index + j) -= mfac; *m.mep(index + j, index - 1 + j) -= mfac; } @@ -473,8 +471,8 @@ void nrn_setup_ext(NrnThread* _nt) { int index = nd->eqn_index_; int parent_index = pnd->eqn_index_; for (j = 0; j < nrn_nlayer_extracellular; ++j) { - *nde->_d[j] -= nde->_b[j]; - *pnde->_d[j] -= nde->_a[j]; + *m.mep(index + j, index + j) -= nde->_b[j]; + *m.mep(parent_index + j, parent_index + j) -= nde->_a[j]; *m.mep(parent_index + j, index + j) += nde->_a[j]; *m.mep(index + j, parent_index + j) += nde->_b[j]; } diff --git a/src/nrnoc/section_fwd.hpp b/src/nrnoc/section_fwd.hpp index 9e010014c1..326acc2021 100644 --- a/src/nrnoc/section_fwd.hpp +++ b/src/nrnoc/section_fwd.hpp @@ -40,7 +40,6 @@ struct Extnode { double* v; /* v external. */ double* _a; double* _b; - double** _d; double** _rhs; /* d, rhs, a, and b are analogous to those in node */ }; #endif diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index 70570c519f..eecd1640d3 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -553,7 +553,6 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) { } } #if EXTRACELLULAR - /* nde->_d[0] contains the -ELECTRODE_CURRENT contribution to nd->_d */ nrn_setup_ext(_nt); #endif if (use_sparse13) { @@ -1952,7 +1951,6 @@ static void nrn_matrix_node_alloc(void) { if (nde) { for (int ie = 0; ie < nlayer; ++ie) { int k = i + ie; - nde->_d[ie] = nt->_sp13mat->mep(k, k); nde->_rhs[ie] = nt->_sp13_rhs + k + 1; } } From 29c4682aed09a3aba4035a20626172737e7e9091 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 28 Oct 2024 13:05:35 +0100 Subject: [PATCH 09/23] Fix compilation error --- src/nocmodl/noccout.cpp | 10 ++++------ src/nrnoc/extcelln.cpp | 7 ++++--- src/nrnoc/membfunc.cpp | 8 ++++++++ src/nrnoc/membfunc.h | 2 ++ src/nrnoc/multicore.cpp | 1 + 5 files changed, 19 insertions(+), 9 deletions(-) diff --git a/src/nocmodl/noccout.cpp b/src/nocmodl/noccout.cpp index f80f0ea553..0c6d26dd32 100644 --- a/src/nocmodl/noccout.cpp +++ b/src/nocmodl/noccout.cpp @@ -260,9 +260,8 @@ void c_out() { P(" }\n"); P("#if EXTRACELLULAR\n"); P(" if (auto* const _extnode = _nrn_mechanism_access_extnode(_nd); _extnode) {\n"); - P(" int index = _nd->eqn_index_;\n"); - P(" OcSparseMatrix& m = *_nt->_sp13mat;\n"); - P(" *m.mep(index, index) += _g;\n"); + P(" int ind = _nrn_mechanism_access_index(_nd);\n"); + P(" *_nrn_mechanism_get_matrix_elem(_nt, ind, ind) += _g;\n"); P(" }\n"); P("#endif\n"); } else { @@ -681,9 +680,8 @@ void c_out_vectorize() { P(" }\n"); P("#if EXTRACELLULAR\n"); P(" if (auto* const _extnode = _nrn_mechanism_access_extnode(_nd); _extnode) {\n"); - P(" int index = _nd->eqn_index_;\n"); - P(" OcSparseMatrix& m = *_nt->_sp13mat;\n"); - P(" *m.mep(index, index) += _g;\n"); + P(" int ind = _nrn_mechanism_access_index(_nd);\n"); + P(" *_nrn_mechanism_get_matrix_elem(_nt, ind, ind) += _g;\n"); P(" }\n"); P("#endif\n"); } else { diff --git a/src/nrnoc/extcelln.cpp b/src/nrnoc/extcelln.cpp index f73a86b42d..cea91a84be 100644 --- a/src/nrnoc/extcelln.cpp +++ b/src/nrnoc/extcelln.cpp @@ -7,6 +7,7 @@ #include "nrniv_mf.h" #include "hocassrt.h" #include "parse.hpp" +#include "ocmatrix.h" extern int nrn_use_daspk_; @@ -430,7 +431,7 @@ void nrn_setup_ext(NrnThread* _nt) { OcSparseMatrix& m = *_nt->_sp13mat; Node* nd = ndlist[i]; int index = nd->eqn_index_; - ExtNode* nde = nd->extnode; + Extnode* nde = nd->extnode; double d = NODED(nd); /* nde->_d only has -ELECTRODE_CURRENT contribution */ *m.mep(index, index) += NODED(nd); @@ -448,7 +449,7 @@ void nrn_setup_ext(NrnThread* _nt) { OcSparseMatrix& m = *_nt->_sp13mat; Node* nd = ndlist[i]; int index = nd->eqn_index_; - ExtNode* nde = nd->extnode; + Extnode* nde = nd->extnode; Node* pnd = _nt->_v_parent[nd->v_node_index]; if (pnd) { /* series resistance and capacitance to ground */ @@ -464,7 +465,7 @@ void nrn_setup_ext(NrnThread* _nt) { *m.mep(index - 1 + j, index + j) -= mfac; *m.mep(index + j, index - 1 + j) -= mfac; } - ExtNode* pnde = pnd->extnode; + Extnode* pnde = pnd->extnode; /* axial connections */ if (pnde) { /* parent sec may not be extracellular */ OcSparseMatrix& m = *_nt->_sp13mat; diff --git a/src/nrnoc/membfunc.cpp b/src/nrnoc/membfunc.cpp index e858536147..de3cc180f2 100644 --- a/src/nrnoc/membfunc.cpp +++ b/src/nrnoc/membfunc.cpp @@ -2,6 +2,7 @@ #include "multicore.h" #include "section.h" +#include "ocmatrix.h" #include @@ -49,6 +50,13 @@ double& _nrn_mechanism_access_rhs(Node* node) { double& _nrn_mechanism_access_voltage(Node* node) { return node->v(); } +int _nrn_mechanism_access_index(const Node* node) { + return node->eqn_index_; +} +double* _nrn_mechanism_get_matrix_elem(NrnThread* nt, int i, int j) { + OcSparseMatrix& m = *nt->_sp13mat; + return m.mep(i, j); +} neuron::container::data_handle _nrn_mechanism_get_area_handle(Node* node) { if (node) { return node->area_handle(); diff --git a/src/nrnoc/membfunc.h b/src/nrnoc/membfunc.h index 18647c19b7..a655828a8c 100644 --- a/src/nrnoc/membfunc.h +++ b/src/nrnoc/membfunc.h @@ -306,6 +306,8 @@ namespace _get { [[nodiscard]] double& _nrn_mechanism_access_param(Prop*, int field, int array_index = 0); [[nodiscard]] double& _nrn_mechanism_access_rhs(Node*); [[nodiscard]] double& _nrn_mechanism_access_voltage(Node*); +[[nodiscard]] int _nrn_mechanism_access_index(const Node*); +[[nodiscard]] double* _nrn_mechanism_get_matrix_elem(NrnThread* nt, int, int ); [[nodiscard]] neuron::container::data_handle _nrn_mechanism_get_area_handle(Node*); [[nodiscard]] Section* _nrn_mechanism_get_child(Section*); [[nodiscard]] int _nrn_mechanism_get_nnode(Section*); diff --git a/src/nrnoc/multicore.cpp b/src/nrnoc/multicore.cpp index a4a97c75e4..9c01395993 100644 --- a/src/nrnoc/multicore.cpp +++ b/src/nrnoc/multicore.cpp @@ -38,6 +38,7 @@ the handling of v_structure_change as long as possible. */ #include "nmodlmutex.h" +#include "ocmatrix.h" #include #include From a9fde15d28fa6e3d7b1f8660b7d5fc0068ec4bc5 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 28 Oct 2024 13:13:33 +0100 Subject: [PATCH 10/23] Stop pre-allocating _sp13mat --- src/nrnoc/treeset.cpp | 31 ------------------------------- 1 file changed, 31 deletions(-) diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index eecd1640d3..608ce70605 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -1912,40 +1912,9 @@ static void nrn_matrix_node_alloc(void) { i += nlayer; } } - // Creating for (int in = 0; in < nt->end; ++in) { Node *nd = nt->_v_node[in]; Extnode* nde = nd->extnode; - Node *pnd = nt->_v_parent[in]; - int i = nd->eqn_index_; - nt->_sp13mat->mep(i - 1, i - 1); - if (nde) { - for (int ie = 0; ie < nlayer; ++ie) { - int k = i + ie; - nt->_sp13mat->mep(k, k); - nt->_sp13mat->mep(k, k - 1); - nt->_sp13mat->mep(k - 1, k); - } - } - if (pnd) { - int j = pnd->eqn_index_; - nt->_sp13mat->mep(j - 1, i - 1); - nt->_sp13mat->mep(i - 1, j - 1); - if (nde && pnd->extnode) - for (int ie = 0; ie < nlayer; ++ie) { - int kp = j + ie; - int k = i + ie; - nt->_sp13mat->mep(kp, k); - nt->_sp13mat->mep(k, kp); - } - } - } - nt->_sp13mat->compress(); - // Collecting - for (int in = 0; in < nt->end; ++in) { - Node *nd = nt->_v_node[in]; - Extnode* nde = nd->extnode; - Node *pnd = nt->_v_parent[in]; int i = nd->eqn_index_; nt->_sp13_rhs[i] = nt->actual_rhs(in); if (nde) { From 2f964f1fe7a5a2432222e452a886d0d182e5d34c Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 29 Oct 2024 11:45:17 +0100 Subject: [PATCH 11/23] Fix MatrixMap --- src/nrniv/matrixmap.cpp | 55 +++++++++++++++++++---------------------- src/nrniv/matrixmap.h | 6 ++--- 2 files changed, 28 insertions(+), 33 deletions(-) diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp index 12b608e529..1bb2d149e4 100644 --- a/src/nrniv/matrixmap.cpp +++ b/src/nrniv/matrixmap.cpp @@ -18,48 +18,43 @@ void MatrixMap::mmfree() { } void MatrixMap::add(double fac) { - for (int i = 0; i < plen_; ++i) { - *ptree_[i] += fac * (*pm_[i]); + NrnThread* _nt = nrn_threads; + for (int i = 0; i < pm_.size(); ++i) { + auto [im, jm] = pm_[i]; + auto [it, jt] = ptree_[i]; + *_nt->_sp13mat->mep(it, jt) += fac * m_.getval(im, jm); + } +} + +int MatrixMap::compute_id(int i, int start, int nnode, Node** nodes, int* layer) const { + int it; + if (i < nnode) { + it = nodes[i]->eqn_index_ + layer[i]; + if (layer[i] > 0 && !nodes[i]->extnode) { + it = 0; + } + } else { + it = start + i - nnode; } + return it; } void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { - static double place_holder = 0; - NrnThread* _nt = nrn_threads; 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()); + pm_.reserve(nonzero_i.size()); + ptree_.reserve(nonzero_i.size()); for (int k = 0; k < nonzero_i.size(); k++) { const int i = nonzero_i[k]; const int j = nonzero_j[k]; - int it; - if (i < nnode) { - it = nodes[i]->eqn_index_ + layer[i]; - if (layer[i] > 0 && !nodes[i]->extnode) { - it = 0; - } - } else { - it = start + i - nnode; - } - int jt; - pm_[plen_] = m_.mep(i, j); - if (j < nnode) { - jt = nodes[j]->eqn_index_ + layer[j]; - if (layer[j] > 0 && !nodes[j]->extnode) { - jt = 0; - } - } else { - jt = start + j - nnode; - } + int it = compute_id(i, start, nnode, nodes, layer); + int jt = compute_id(j, start, nnode, nodes, layer); if (it == 0 || jt == 0) { - ptree_[plen_] = &place_holder; - } else { - ptree_[plen_] = _nt->_sp13mat->mep(it - 1, jt - 1); + continue; } - ++plen_; + pm_.emplace_back(i, j); + ptree_.emplace_back(it - 1, jt - 1); } } diff --git a/src/nrniv/matrixmap.h b/src/nrniv/matrixmap.h index f216dcc9c3..bbcb92466a 100644 --- a/src/nrniv/matrixmap.h +++ b/src/nrniv/matrixmap.h @@ -112,10 +112,10 @@ class MatrixMap { private: + int compute_id(int i, int start, int nnode, Node** nodes, int* layer) const; Matrix& m_; // the map - int plen_ = 0; - std::vector pm_{}; - std::vector ptree_{}; + std::vector> pm_{}; + std::vector> ptree_{}; }; From dc0525479872c67ee3259af90a6d015e1583e495 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 29 Oct 2024 12:32:52 +0100 Subject: [PATCH 12/23] Remove useless local variables --- src/nrnoc/extcelln.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/nrnoc/extcelln.cpp b/src/nrnoc/extcelln.cpp index cea91a84be..32b13ef618 100644 --- a/src/nrnoc/extcelln.cpp +++ b/src/nrnoc/extcelln.cpp @@ -468,8 +468,6 @@ void nrn_setup_ext(NrnThread* _nt) { Extnode* pnde = pnd->extnode; /* axial connections */ if (pnde) { /* parent sec may not be extracellular */ - OcSparseMatrix& m = *_nt->_sp13mat; - int index = nd->eqn_index_; int parent_index = pnd->eqn_index_; for (j = 0; j < nrn_nlayer_extracellular; ++j) { *m.mep(index + j, index + j) -= nde->_b[j]; From b2006eafca3d3d7884e846506c7180628a2ba926 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 30 Oct 2024 16:26:26 +0100 Subject: [PATCH 13/23] Modernize --- src/nrnoc/extcelln.cpp | 22 +++++++++++----------- src/nrnoc/membfunc.cpp | 2 +- src/nrnoc/treeset.cpp | 10 +++++----- 3 files changed, 17 insertions(+), 17 deletions(-) diff --git a/src/nrnoc/extcelln.cpp b/src/nrnoc/extcelln.cpp index 32b13ef618..404aef8b40 100644 --- a/src/nrnoc/extcelln.cpp +++ b/src/nrnoc/extcelln.cpp @@ -434,12 +434,12 @@ void nrn_setup_ext(NrnThread* _nt) { Extnode* nde = nd->extnode; double d = NODED(nd); /* nde->_d only has -ELECTRODE_CURRENT contribution */ - *m.mep(index, index) += NODED(nd); + m(index, index) += NODED(nd); d = m.getval(index, index); /* now d is only the membrane current contribution */ /* i.e. d = cm/dt + di/dvm */ - *m.mep(index - 1, index) -= d; - *m.mep(index, index - 1) -= d; + m(index - 1, index) -= d; + m(index, index - 1) -= d; #if I_MEMBRANE ml->data(i, sav_g_index) = d; #endif @@ -456,24 +456,24 @@ void nrn_setup_ext(NrnThread* _nt) { int j = 0; for (;;) { /* between j and j+1 layer */ mfac = (*nde->param[xg_index_ext(j)] + *nde->param[xc_index_ext(j)] * cfac); - *m.mep(index + j, index + j) += mfac; + m(index + j, index + j) += mfac; ++j; if (j == nrn_nlayer_extracellular) { break; } - *m.mep(index + j, index + j) += mfac; - *m.mep(index - 1 + j, index + j) -= mfac; - *m.mep(index + j, index - 1 + j) -= mfac; + m(index + j, index + j) += mfac; + m(index - 1 + j, index + j) -= mfac; + m(index + j, index - 1 + j) -= mfac; } Extnode* pnde = pnd->extnode; /* axial connections */ if (pnde) { /* parent sec may not be extracellular */ int parent_index = pnd->eqn_index_; for (j = 0; j < nrn_nlayer_extracellular; ++j) { - *m.mep(index + j, index + j) -= nde->_b[j]; - *m.mep(parent_index + j, parent_index + j) -= nde->_a[j]; - *m.mep(parent_index + j, index + j) += nde->_a[j]; - *m.mep(index + j, parent_index + j) += nde->_b[j]; + m(index + j, index + j) -= nde->_b[j]; + m(parent_index + j, parent_index + j) -= nde->_a[j]; + m(parent_index + j, index + j) += nde->_a[j]; + m(index + j, parent_index + j) += nde->_b[j]; } } } diff --git a/src/nrnoc/membfunc.cpp b/src/nrnoc/membfunc.cpp index de3cc180f2..a56680267a 100644 --- a/src/nrnoc/membfunc.cpp +++ b/src/nrnoc/membfunc.cpp @@ -55,7 +55,7 @@ int _nrn_mechanism_access_index(const Node* node) { } double* _nrn_mechanism_get_matrix_elem(NrnThread* nt, int i, int j) { OcSparseMatrix& m = *nt->_sp13mat; - return m.mep(i, j); + return m(i, j); } neuron::container::data_handle _nrn_mechanism_get_area_handle(Node* node) { if (node) { diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index 608ce70605..463f4e0626 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -365,7 +365,7 @@ void update_sp13_mat_based_on_actual_d(NrnThread* nt) { for (int i = 0; i < nt->end; ++i) { const int index = nt->_v_node[i]->eqn_index_; OcSparseMatrix& m = *nt->_sp13mat; - *m.mep(index - 1, index - 1) = nt->actual_d(i); + m(index - 1, index - 1) = nt->actual_d(i); } } @@ -580,15 +580,15 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) { // Update entries in sp13_mat { const int index = nd->eqn_index_; - *m.mep(index - 1, index - 1) -= nd_b; + m(index - 1, index - 1) -= nd_b; // used to update NODED (sparse13 matrix) using NODEA and NODEB ("SoA") Node* const parent_nd = _nt->_v_node[parent_i]; const int parent_index = parent_nd->eqn_index_; - *m.mep(parent_index - 1, parent_index - 1) -= nd_a; + m(parent_index - 1, parent_index - 1) -= nd_a; - *m.mep(parent_index - 1, index - 1) += nd_a; - *m.mep(index - 1, parent_index - 1) += nd_b; /* b may have value from lincir */ + m(parent_index - 1, index - 1) += nd_a; + m(index - 1, parent_index - 1) += nd_b; /* b may have value from lincir */ } // Also update the Node's d value in the SoA storage (is this needed?) vec_d[i] -= nd_b; From cfc1b3c6853dc435d8f1f9e0521e84a4b98b4d32 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 30 Oct 2024 16:40:41 +0100 Subject: [PATCH 14/23] Fix merge --- src/nrniv/matrixmap.cpp | 13 +++++-------- src/nrnoc/membfunc.cpp | 2 +- 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp index 9e1584fd6b..c3951ef80b 100644 --- a/src/nrniv/matrixmap.cpp +++ b/src/nrniv/matrixmap.cpp @@ -22,7 +22,7 @@ void MatrixMap::add(double fac) { for (int i = 0; i < pm_.size(); ++i) { auto [im, jm] = pm_[i]; auto [it, jt] = ptree_[i]; - _nt->_sp13mat(it, jt) += fac * m_(im, jm); + _nt->_sp13mat->coeff(it, jt) += fac * m_(im, jm); } } @@ -42,13 +42,10 @@ int MatrixMap::compute_id(int i, int start, int nnode, Node** nodes, int* layer) void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { mmfree(); - std::vector nonzero_i, nonzero_j; - m_.nonzeros(nonzero_i, nonzero_j); - pm_.reserve(nonzero_i.size()); - ptree_.reserve(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_.reserve(nzs.size()); + ptree_.reserve(nzs.size()); + for (const auto& [i, j]: nzs) { int it = compute_id(i, start, nnode, nodes, layer); int jt = compute_id(j, start, nnode, nodes, layer); if (it == 0 || jt == 0) { diff --git a/src/nrnoc/membfunc.cpp b/src/nrnoc/membfunc.cpp index a56680267a..de3cc180f2 100644 --- a/src/nrnoc/membfunc.cpp +++ b/src/nrnoc/membfunc.cpp @@ -55,7 +55,7 @@ int _nrn_mechanism_access_index(const Node* node) { } double* _nrn_mechanism_get_matrix_elem(NrnThread* nt, int i, int j) { OcSparseMatrix& m = *nt->_sp13mat; - return m(i, j); + return m.mep(i, j); } neuron::container::data_handle _nrn_mechanism_get_area_handle(Node* node) { if (node) { From 4667b0b89872f7e2ad395c43b0e322ce7ef1567d Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 31 Oct 2024 09:57:06 +0100 Subject: [PATCH 15/23] Remove useless function compress --- src/ivoc/ocmatrix.h | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/ivoc/ocmatrix.h b/src/ivoc/ocmatrix.h index 3a4ef4076c..9e874c86a3 100644 --- a/src/ivoc/ocmatrix.h +++ b/src/ivoc/ocmatrix.h @@ -232,10 +232,6 @@ class OcSparseMatrix final: public OcMatrix { // type 2 void zero() override; - void compress() { - m_.makeCompressed(); - } - private: Eigen::SparseMatrix m_{}; std::unique_ptr> lu_{}; From 6273e222dd1e0eef68a06b3141c6f5374c858c0a Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 31 Oct 2024 09:59:56 +0100 Subject: [PATCH 16/23] Remove useless include of sparse13 --- src/nrncvode/nrndaspk.cpp | 1 - src/nrniv/matrixmap.cpp | 2 -- src/nrnoc/fadvance.cpp | 1 - 3 files changed, 4 deletions(-) diff --git a/src/nrncvode/nrndaspk.cpp b/src/nrncvode/nrndaspk.cpp index 5434cdc031..255ca35733 100644 --- a/src/nrncvode/nrndaspk.cpp +++ b/src/nrncvode/nrndaspk.cpp @@ -6,7 +6,6 @@ #include #include -#include "spmatrix.h" #include "nrnoc2iv.h" #include "cvodeobj.h" #include "nrndaspk.h" diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp index c3951ef80b..0969b63a72 100644 --- a/src/nrniv/matrixmap.cpp +++ b/src/nrniv/matrixmap.cpp @@ -2,8 +2,6 @@ #include "matrixmap.h" #include -#include "spmatrix.h" - MatrixMap::MatrixMap(Matrix& mat) : m_(mat) {} diff --git a/src/nrnoc/fadvance.cpp b/src/nrnoc/fadvance.cpp index f51d0b9958..82c9e901c9 100644 --- a/src/nrnoc/fadvance.cpp +++ b/src/nrnoc/fadvance.cpp @@ -13,7 +13,6 @@ #include "utils/profile/profiler_interface.h" #include "nonvintblock.h" #include "nrncvode.h" -#include "spmatrix.h" #include From 2827fa028aafc39557b29f09f905c67b04b3b011 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 4 Nov 2024 16:05:03 +0100 Subject: [PATCH 17/23] Still free _rhs --- src/nrnoc/extcelln.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/nrnoc/extcelln.cpp b/src/nrnoc/extcelln.cpp index 404aef8b40..d93b008dcf 100644 --- a/src/nrnoc/extcelln.cpp +++ b/src/nrnoc/extcelln.cpp @@ -214,10 +214,11 @@ static void extcell_init(neuron::model_sorted_token const&, void extnode_free_elements(Extnode* nde) { if (nde->v) { free(nde->v); /* along with _a and _b */ - nde->v = NULL; - nde->_a = NULL; - nde->_b = NULL; - nde->_rhs = NULL; + free(nde->_rhs); + nde->v = nullptr; + nde->_a = nullptr; + nde->_b = nullptr; + nde->_rhs = nullptr; } } From 6fd562dc47eaf742c2d877b3f5e01c0f9afd1021 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 7 Nov 2024 14:58:48 +0100 Subject: [PATCH 18/23] Minimize change --- src/nrniv/matrixmap.cpp | 20 ++++++++++---------- src/nrnoc/extcelln.cpp | 29 ++++++++++++++++------------- 2 files changed, 26 insertions(+), 23 deletions(-) diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp index c9569f5d23..78239aa4f7 100644 --- a/src/nrniv/matrixmap.cpp +++ b/src/nrniv/matrixmap.cpp @@ -18,13 +18,13 @@ void MatrixMap::mmfree() { void MatrixMap::add(double fac) { NrnThread* _nt = nrn_threads; for (int i = 0; i < pm_.size(); ++i) { - auto [im, jm] = pm_[i]; - auto [it, jt] = ptree_[i]; - _nt->_sp13mat->coeff(it, jt) += fac * m_(im, jm); + auto [it, jt] = pm_[i]; + auto [im, jm] = ptree_[i]; + _nt->_sp13mat->coeff(im, jm) += fac * m_(it, jt); } } -int MatrixMap::compute_id(int i, int start, int nnode, Node** nodes, int* layer) const { +int MatrixMap::compute_index(int i, int start, int nnode, Node** nodes, int* layer) const { int it; if (i < nnode) { it = nodes[i]->eqn_index_ + layer[i]; @@ -44,11 +44,11 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { pm_.reserve(nzs.size()); ptree_.reserve(nzs.size()); for (const auto& [i, j]: nzs) { - int it = compute_id(i, start, nnode, nodes, layer); - int jt = compute_id(j, start, nnode, nodes, layer); - if (it == 0 || jt == 0) { - continue; + int it = compute_index(i, start, nnode, nodes, layer); + int jt = compute_index(j, start, nnode, nodes, layer); + if (it != 0 && jt != 0) { + pm_.emplace_back(i, j); + ptree_.emplace_back(it - 1, jt - 1); } - pm_.emplace_back(i, j); - ptree_.emplace_back(it - 1, jt - 1); + } } diff --git a/src/nrnoc/extcelln.cpp b/src/nrnoc/extcelln.cpp index d93b008dcf..9ef904fdd7 100644 --- a/src/nrnoc/extcelln.cpp +++ b/src/nrnoc/extcelln.cpp @@ -415,25 +415,28 @@ void nrn_rhs_ext(NrnThread* _nt) { } void nrn_setup_ext(NrnThread* _nt) { - double cfac, mfac; + int i, j, cnt; + Node *nd, *pnd, **ndlist; + double d, cfac, mfac; + Extnode *nde, *pnde; Memb_list* ml = _nt->_ecell_memb_list; if (!ml) { return; } /*printnode("begin setup");*/ - int cnt = ml->nodecount; - Node** ndlist = ml->nodelist; + cnt = ml->nodecount; + ndlist = ml->nodelist; cfac = .001 * _nt->cj; /* d contains all the membrane conductances (and capacitance) */ /* i.e. (cm/dt + di/dvm - dis/dvi)*[dvi] and (dis/dvi)*[dvx] */ - for (int i = 0; i < cnt; ++i) { + for (i = 0; i < cnt; ++i) { OcSparseMatrix& m = *_nt->_sp13mat; - Node* nd = ndlist[i]; + nd = ndlist[i]; int index = nd->eqn_index_; - Extnode* nde = nd->extnode; - double d = NODED(nd); + nde = nd->extnode; + d = NODED(nd); /* nde->_d only has -ELECTRODE_CURRENT contribution */ m(index, index) += NODED(nd); d = m.getval(index, index); @@ -446,15 +449,15 @@ void nrn_setup_ext(NrnThread* _nt) { #endif } /* series resistance, capacitance, and axial terms. */ - for (int i = 0; i < cnt; ++i) { + for (i = 0; i < cnt; ++i) { OcSparseMatrix& m = *_nt->_sp13mat; - Node* nd = ndlist[i]; + nd = ndlist[i]; int index = nd->eqn_index_; - Extnode* nde = nd->extnode; - Node* pnd = _nt->_v_parent[nd->v_node_index]; + nde = nd->extnode; + pnd = _nt->_v_parent[nd->v_node_index]; if (pnd) { /* series resistance and capacitance to ground */ - int j = 0; + j = 0; for (;;) { /* between j and j+1 layer */ mfac = (*nde->param[xg_index_ext(j)] + *nde->param[xc_index_ext(j)] * cfac); m(index + j, index + j) += mfac; @@ -466,7 +469,7 @@ void nrn_setup_ext(NrnThread* _nt) { m(index - 1 + j, index + j) -= mfac; m(index + j, index - 1 + j) -= mfac; } - Extnode* pnde = pnd->extnode; + pnde = pnd->extnode; /* axial connections */ if (pnde) { /* parent sec may not be extracellular */ int parent_index = pnd->eqn_index_; From 2399a433c4d23e4b91364c961e112347e087d836 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 7 Nov 2024 15:31:41 +0100 Subject: [PATCH 19/23] Reduce changes --- src/nrniv/matrixmap.cpp | 18 +++++++++--------- src/nrnoc/membfunc.cpp | 4 ++-- src/nrnoc/section_fwd.hpp | 2 +- src/nrnoc/treeset.cpp | 26 +++++++++++++++----------- 4 files changed, 27 insertions(+), 23 deletions(-) diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp index 78239aa4f7..69d3a85249 100644 --- a/src/nrniv/matrixmap.cpp +++ b/src/nrniv/matrixmap.cpp @@ -25,16 +25,16 @@ void MatrixMap::add(double fac) { } int MatrixMap::compute_index(int i, int start, int nnode, Node** nodes, int* layer) const { - int it; - if (i < nnode) { - it = nodes[i]->eqn_index_ + layer[i]; - if (layer[i] > 0 && !nodes[i]->extnode) { - it = 0; + int it; + if (i < nnode) { + it = nodes[i]->eqn_index_ + layer[i]; + if (layer[i] > 0 && !nodes[i]->extnode) { + it = 0; + } + } else { + it = start + i - nnode; } - } else { - it = start + i - nnode; - } - return it; + return it; } void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { diff --git a/src/nrnoc/membfunc.cpp b/src/nrnoc/membfunc.cpp index de3cc180f2..e6eeaf2970 100644 --- a/src/nrnoc/membfunc.cpp +++ b/src/nrnoc/membfunc.cpp @@ -53,9 +53,9 @@ double& _nrn_mechanism_access_voltage(Node* node) { int _nrn_mechanism_access_index(const Node* node) { return node->eqn_index_; } -double* _nrn_mechanism_get_matrix_elem(NrnThread* nt, int i, int j) { +double& _nrn_mechanism_get_matrix_elem(NrnThread* nt, int i, int j) { OcSparseMatrix& m = *nt->_sp13mat; - return m.mep(i, j); + return m.coeff(i, j); } neuron::container::data_handle _nrn_mechanism_get_area_handle(Node* node) { if (node) { diff --git a/src/nrnoc/section_fwd.hpp b/src/nrnoc/section_fwd.hpp index 326acc2021..54e7e367c1 100644 --- a/src/nrnoc/section_fwd.hpp +++ b/src/nrnoc/section_fwd.hpp @@ -40,7 +40,7 @@ struct Extnode { double* v; /* v external. */ double* _a; double* _b; - double** _rhs; /* d, rhs, a, and b are analogous to those in node */ + double** _rhs; /* rhs, a, and b are analogous to those in node */ }; #endif diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index 463f4e0626..eb7e6c25d2 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -1871,7 +1871,7 @@ and therefore is passed to spSolve as actual_rhs intead of actual_rhs-1. */ static void nrn_matrix_node_alloc(void) { - int b; + int i, b; Node* nd; NrnThread* nt; @@ -1896,9 +1896,10 @@ static void nrn_matrix_node_alloc(void) { } ++nrn_matrix_cnt_; if (use_sparse13) { + int in, extn, neqn; nt = nrn_threads; - int neqn = nt->end + nrndae_extra_eqn_count(); - int extn = 0; + neqn = nt->end + nrndae_extra_eqn_count(); + extn = 0; if (nt->_ecell_memb_list) { extn = nt->_ecell_memb_list->nodecount * nlayer; } @@ -1906,21 +1907,24 @@ static void nrn_matrix_node_alloc(void) { neqn += extn; nt->_sp13_rhs = (double*) ecalloc(neqn + 1, sizeof(double)); nt->_sp13mat = new OcSparseMatrix(neqn, neqn); - for (int in = 0, i = 1; in < nt->end; ++in, ++i) { + for (in = 0, i = 1; in < nt->end; ++in, ++i) { nt->_v_node[in]->eqn_index_ = i; if (nt->_v_node[in]->extnode) { i += nlayer; } } - for (int in = 0; in < nt->end; ++in) { - Node *nd = nt->_v_node[in]; - Extnode* nde = nd->extnode; - int i = nd->eqn_index_; + for (in = 0; in < nt->end; ++in) { + int ie, k; + Node *nd; + Extnode* nde; + nd = nt->_v_node[in]; + nde = nd->extnode; + i = nd->eqn_index_; nt->_sp13_rhs[i] = nt->actual_rhs(in); if (nde) { - for (int ie = 0; ie < nlayer; ++ie) { - int k = i + ie; - nde->_rhs[ie] = nt->_sp13_rhs + k + 1; + for (ie = 0; ie < nlayer; ++ie) { + k = i + ie + 1; + nde->_rhs[ie] = nt->_sp13_rhs + k; } } } From 2406e7c61ed775a44ead95c0061c8fdfec9cfac8 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 7 Nov 2024 15:37:36 +0100 Subject: [PATCH 20/23] Simplify --- src/nrnoc/treeset.cpp | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index eb7e6c25d2..d32f75b83c 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -575,20 +575,18 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) { Node* nd = _nt->_v_node[i]; OcSparseMatrix& m = *_nt->_sp13mat; auto const parent_i = _nt->_v_parent_index[i]; + auto* const parent_nd = _nt->_v_node[parent_i]; auto const nd_a = NODEA(nd); auto const nd_b = NODEB(nd); // Update entries in sp13_mat { const int index = nd->eqn_index_; - m(index - 1, index - 1) -= nd_b; - - // used to update NODED (sparse13 matrix) using NODEA and NODEB ("SoA") - Node* const parent_nd = _nt->_v_node[parent_i]; const int parent_index = parent_nd->eqn_index_; - m(parent_index - 1, parent_index - 1) -= nd_a; - m(parent_index - 1, index - 1) += nd_a; m(index - 1, parent_index - 1) += nd_b; /* b may have value from lincir */ + m(index - 1, index - 1) -= nd_b; + // used to update NODED (sparse13 matrix) using NODEA and NODEB ("SoA") + m(parent_index - 1, parent_index - 1) -= nd_a; } // Also update the Node's d value in the SoA storage (is this needed?) vec_d[i] -= nd_b; From 8fb9cb0eb7aa5934de493ed9dda0afab030c089e Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 7 Nov 2024 16:27:15 +0100 Subject: [PATCH 21/23] Fix compilation --- src/nocmodl/noccout.cpp | 4 ++-- src/nrnoc/membfunc.h | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/nocmodl/noccout.cpp b/src/nocmodl/noccout.cpp index 0c6d26dd32..9a7952b06d 100644 --- a/src/nocmodl/noccout.cpp +++ b/src/nocmodl/noccout.cpp @@ -261,7 +261,7 @@ void c_out() { P("#if EXTRACELLULAR\n"); P(" if (auto* const _extnode = _nrn_mechanism_access_extnode(_nd); _extnode) {\n"); P(" int ind = _nrn_mechanism_access_index(_nd);\n"); - P(" *_nrn_mechanism_get_matrix_elem(_nt, ind, ind) += _g;\n"); + P(" _nrn_mechanism_get_matrix_elem(_nt, ind, ind) += _g;\n"); P(" }\n"); P("#endif\n"); } else { @@ -681,7 +681,7 @@ void c_out_vectorize() { P("#if EXTRACELLULAR\n"); P(" if (auto* const _extnode = _nrn_mechanism_access_extnode(_nd); _extnode) {\n"); P(" int ind = _nrn_mechanism_access_index(_nd);\n"); - P(" *_nrn_mechanism_get_matrix_elem(_nt, ind, ind) += _g;\n"); + P(" _nrn_mechanism_get_matrix_elem(_nt, ind, ind) += _g;\n"); P(" }\n"); P("#endif\n"); } else { diff --git a/src/nrnoc/membfunc.h b/src/nrnoc/membfunc.h index a655828a8c..5d7d79b8be 100644 --- a/src/nrnoc/membfunc.h +++ b/src/nrnoc/membfunc.h @@ -307,7 +307,7 @@ namespace _get { [[nodiscard]] double& _nrn_mechanism_access_rhs(Node*); [[nodiscard]] double& _nrn_mechanism_access_voltage(Node*); [[nodiscard]] int _nrn_mechanism_access_index(const Node*); -[[nodiscard]] double* _nrn_mechanism_get_matrix_elem(NrnThread* nt, int, int ); +[[nodiscard]] double& _nrn_mechanism_get_matrix_elem(NrnThread* nt, int, int ); [[nodiscard]] neuron::container::data_handle _nrn_mechanism_get_area_handle(Node*); [[nodiscard]] Section* _nrn_mechanism_get_child(Section*); [[nodiscard]] int _nrn_mechanism_get_nnode(Section*); From 03d5c4bcdfc53502e358b4743ce03375cfd89d2b Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 7 Nov 2024 17:06:11 +0100 Subject: [PATCH 22/23] OcSparseMatrix => OcMatrix --- src/nrnoc/extcelln.cpp | 4 ++-- src/nrnoc/membfunc.cpp | 2 +- src/nrnoc/multicore.h | 4 ++-- src/nrnoc/treeset.cpp | 8 ++++---- 4 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/nrnoc/extcelln.cpp b/src/nrnoc/extcelln.cpp index 9ef904fdd7..9380924a1a 100644 --- a/src/nrnoc/extcelln.cpp +++ b/src/nrnoc/extcelln.cpp @@ -432,7 +432,7 @@ void nrn_setup_ext(NrnThread* _nt) { /* i.e. (cm/dt + di/dvm - dis/dvi)*[dvi] and (dis/dvi)*[dvx] */ for (i = 0; i < cnt; ++i) { - OcSparseMatrix& m = *_nt->_sp13mat; + OcMatrix& m = *_nt->_sp13mat; nd = ndlist[i]; int index = nd->eqn_index_; nde = nd->extnode; @@ -450,7 +450,7 @@ void nrn_setup_ext(NrnThread* _nt) { } /* series resistance, capacitance, and axial terms. */ for (i = 0; i < cnt; ++i) { - OcSparseMatrix& m = *_nt->_sp13mat; + OcMatrix& m = *_nt->_sp13mat; nd = ndlist[i]; int index = nd->eqn_index_; nde = nd->extnode; diff --git a/src/nrnoc/membfunc.cpp b/src/nrnoc/membfunc.cpp index e6eeaf2970..7d48a0ee0b 100644 --- a/src/nrnoc/membfunc.cpp +++ b/src/nrnoc/membfunc.cpp @@ -54,7 +54,7 @@ int _nrn_mechanism_access_index(const Node* node) { return node->eqn_index_; } double& _nrn_mechanism_get_matrix_elem(NrnThread* nt, int i, int j) { - OcSparseMatrix& m = *nt->_sp13mat; + OcMatrix& m = *nt->_sp13mat; return m.coeff(i, j); } neuron::container::data_handle _nrn_mechanism_get_area_handle(Node* node) { diff --git a/src/nrnoc/multicore.h b/src/nrnoc/multicore.h index ac8d27b718..0511fbb33b 100644 --- a/src/nrnoc/multicore.h +++ b/src/nrnoc/multicore.h @@ -29,7 +29,7 @@ actual_v, etc. #include "membfunc.h" #include -class OcSparseMatrix; +class OcMatrix; typedef struct NrnThreadMembList { /* patterned after CvMembList in cvodeobj.h */ struct NrnThreadMembList* next; @@ -92,7 +92,7 @@ struct NrnThread { Node** _v_parent; double* _sp13_rhs; /* rhs matrix for sparse13 solver. updates to and from this vector need to be transfered to actual_rhs */ - OcSparseMatrix* _sp13mat; /* handle to general sparse matrix */ + OcMatrix* _sp13mat; /* handle to general sparse matrix */ Memb_list* _ecell_memb_list; /* normally nullptr */ Node** _ecell_children; /* nodes with no extcell but parent has it */ void* _vcv; /* replaces old cvode_instance and nrn_cvode_ */ diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index d32f75b83c..3694e8a25a 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -353,7 +353,7 @@ void update_sp13_rhs_based_on_actual_rhs(NrnThread* nt) { void update_actual_d_based_on_sp13_mat(NrnThread* nt) { for (int i = 0; i < nt->end; ++i) { const int index = nt->_v_node[i]->eqn_index_; - OcSparseMatrix& m = *nt->_sp13mat; + OcMatrix& m = *nt->_sp13mat; nt->actual_d(i) = m.getval(index - 1, index - 1); } } @@ -364,7 +364,7 @@ void update_actual_d_based_on_sp13_mat(NrnThread* nt) { void update_sp13_mat_based_on_actual_d(NrnThread* nt) { for (int i = 0; i < nt->end; ++i) { const int index = nt->_v_node[i]->eqn_index_; - OcSparseMatrix& m = *nt->_sp13mat; + OcMatrix& m = *nt->_sp13mat; m(index - 1, index - 1) = nt->actual_d(i); } } @@ -573,7 +573,7 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) { update_sp13_mat_based_on_actual_d(_nt); // just because of activclamp_lhs for (i = i2; i < i3; ++i) { // note i2 Node* nd = _nt->_v_node[i]; - OcSparseMatrix& m = *_nt->_sp13mat; + OcMatrix& m = *_nt->_sp13mat; auto const parent_i = _nt->_v_parent_index[i]; auto* const parent_nd = _nt->_v_node[parent_i]; auto const nd_a = NODEA(nd); @@ -1904,7 +1904,7 @@ static void nrn_matrix_node_alloc(void) { /*printf(" %d extracellular nodes\n", extn);*/ neqn += extn; nt->_sp13_rhs = (double*) ecalloc(neqn + 1, sizeof(double)); - nt->_sp13mat = new OcSparseMatrix(neqn, neqn); + nt->_sp13mat = OcMatrix::instance(neqn, neqn, OcMatrix::MSPARSE); for (in = 0, i = 1; in < nt->end; ++in, ++i) { nt->_v_node[in]->eqn_index_ = i; if (nt->_v_node[in]->extnode) { From c47832bc62a41cf11cffebaaa0afc581c601ef76 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Fri, 8 Nov 2024 14:16:13 +0100 Subject: [PATCH 23/23] Add Extnode->eqn_index_ --- src/nrnoc/extcelln.cpp | 31 +++++++++++++++++-------------- src/nrnoc/section_fwd.hpp | 2 ++ src/nrnoc/treeset.cpp | 3 ++- 3 files changed, 21 insertions(+), 15 deletions(-) diff --git a/src/nrnoc/extcelln.cpp b/src/nrnoc/extcelln.cpp index 9380924a1a..552e9b58cb 100644 --- a/src/nrnoc/extcelln.cpp +++ b/src/nrnoc/extcelln.cpp @@ -431,53 +431,56 @@ void nrn_setup_ext(NrnThread* _nt) { /* d contains all the membrane conductances (and capacitance) */ /* i.e. (cm/dt + di/dvm - dis/dvi)*[dvi] and (dis/dvi)*[dvx] */ + // This loop handle conductances between the node and the first layer for (i = 0; i < cnt; ++i) { OcMatrix& m = *_nt->_sp13mat; nd = ndlist[i]; int index = nd->eqn_index_; nde = nd->extnode; + int ext_index = nde->eqn_index_; d = NODED(nd); /* nde->_d only has -ELECTRODE_CURRENT contribution */ - m(index, index) += NODED(nd); - d = m.getval(index, index); + m(ext_index + 0, ext_index + 0) += NODED(nd); + d = m.getval(ext_index + 0, ext_index + 0); /* now d is only the membrane current contribution */ /* i.e. d = cm/dt + di/dvm */ - m(index - 1, index) -= d; - m(index, index - 1) -= d; + m(index - 1, ext_index + 0) -= d; + m(ext_index + 0, index - 1) -= d; #if I_MEMBRANE ml->data(i, sav_g_index) = d; #endif } /* series resistance, capacitance, and axial terms. */ + // This look takes care of the conductances between layers for (i = 0; i < cnt; ++i) { OcMatrix& m = *_nt->_sp13mat; nd = ndlist[i]; - int index = nd->eqn_index_; nde = nd->extnode; + int ext_index = nde->eqn_index_; pnd = _nt->_v_parent[nd->v_node_index]; if (pnd) { /* series resistance and capacitance to ground */ j = 0; for (;;) { /* between j and j+1 layer */ mfac = (*nde->param[xg_index_ext(j)] + *nde->param[xc_index_ext(j)] * cfac); - m(index + j, index + j) += mfac; + m(ext_index + j, ext_index + j) += mfac; ++j; if (j == nrn_nlayer_extracellular) { break; } - m(index + j, index + j) += mfac; - m(index - 1 + j, index + j) -= mfac; - m(index + j, index - 1 + j) -= mfac; + m(ext_index + j, ext_index + j) += mfac; + m(ext_index + j - 1, ext_index + j) -= mfac; + m(ext_index + j, ext_index + j - 1) -= mfac; } pnde = pnd->extnode; /* axial connections */ if (pnde) { /* parent sec may not be extracellular */ - int parent_index = pnd->eqn_index_; + int parent_ext_index = pnde->eqn_index_; for (j = 0; j < nrn_nlayer_extracellular; ++j) { - m(index + j, index + j) -= nde->_b[j]; - m(parent_index + j, parent_index + j) -= nde->_a[j]; - m(parent_index + j, index + j) += nde->_a[j]; - m(index + j, parent_index + j) += nde->_b[j]; + m(ext_index + j, ext_index + j) -= nde->_b[j]; + m(parent_ext_index + j, parent_ext_index + j) -= nde->_a[j]; + m(parent_ext_index + j, ext_index + j) += nde->_a[j]; + m(ext_index + j, parent_ext_index + j) += nde->_b[j]; } } } diff --git a/src/nrnoc/section_fwd.hpp b/src/nrnoc/section_fwd.hpp index 54e7e367c1..ee86737e13 100644 --- a/src/nrnoc/section_fwd.hpp +++ b/src/nrnoc/section_fwd.hpp @@ -41,6 +41,8 @@ struct Extnode { double* _a; double* _b; double** _rhs; /* rhs, a, and b are analogous to those in node */ + + int eqn_index_; }; #endif diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index 3694e8a25a..cfb3f3c2d0 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -1906,7 +1906,7 @@ static void nrn_matrix_node_alloc(void) { nt->_sp13_rhs = (double*) ecalloc(neqn + 1, sizeof(double)); nt->_sp13mat = OcMatrix::instance(neqn, neqn, OcMatrix::MSPARSE); for (in = 0, i = 1; in < nt->end; ++in, ++i) { - nt->_v_node[in]->eqn_index_ = i; + nt->_v_node[in]->eqn_index_ = i; // 1-indexed if (nt->_v_node[in]->extnode) { i += nlayer; } @@ -1920,6 +1920,7 @@ static void nrn_matrix_node_alloc(void) { i = nd->eqn_index_; nt->_sp13_rhs[i] = nt->actual_rhs(in); if (nde) { + nde->eqn_index_ = i; // 0-indexed for (ie = 0; ie < nlayer; ++ie) { k = i + ie + 1; nde->_rhs[ie] = nt->_sp13_rhs + k;