From 5c5c730489bfd8e71e5a2bf5c73b44d266c31094 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 13 Dec 2023 17:09:08 +0100 Subject: [PATCH 01/26] Sparse13 --- src/nrncvode/occvode.cpp | 2 +- src/nrniv/matrixmap.cpp | 2 +- src/nrnoc/fadvance.cpp | 8 +++--- src/nrnoc/multicore.cpp | 6 ++-- src/nrnoc/multicore.h | 4 ++- src/nrnoc/solve.cpp | 18 ++++-------- src/nrnoc/treeset.cpp | 60 ++++++++++++++++++---------------------- 7 files changed, 45 insertions(+), 55 deletions(-) diff --git a/src/nrncvode/occvode.cpp b/src/nrncvode/occvode.cpp index 889f9fbc69..854e1972d4 100644 --- a/src/nrncvode/occvode.cpp +++ b/src/nrncvode/occvode.cpp @@ -356,7 +356,7 @@ void Cvode::daspk_init_eqn() { if (use_sparse13 == 0 || diam_changed != 0) { recalc_diam(); } - zneq = spGetSize(_nt->_sp13mat, 0); + zneq = _nt->_sp13mat->cols(); 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 02464c40a7..0fb76c2aff 100644 --- a/src/nrniv/matrixmap.cpp +++ b/src/nrniv/matrixmap.cpp @@ -73,7 +73,7 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { jt = start + j - nnode; } // printf("MatrixMap::alloc getElement(%d,%d)\n", it, jt); - ptree_[plen_] = spGetElement(_nt->_sp13mat, it, jt); + ptree_[plen_] = &_nt->_sp13mat->coeffRef(it-1, jt-1); ++plen_; } } diff --git a/src/nrnoc/fadvance.cpp b/src/nrnoc/fadvance.cpp index 11dc9b89df..a0a032ce1e 100644 --- a/src/nrnoc/fadvance.cpp +++ b/src/nrnoc/fadvance.cpp @@ -667,11 +667,11 @@ 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); + std::cout << &_nt->_sp13mat << std::endl; } else { - int i, n = spGetSize(_nt->_sp13mat, 0); - spPrint(_nt->_sp13mat, 1, 1, 1); - for (i = 1; i <= n; ++i) { + std::cout << &_nt->_sp13mat << std::endl; + int n = _nt->_sp13mat->cols(); + for (int i = 1; i <= n; ++i) { Printf("%d %g\n", i, _nt->actual_rhs(i)); } } diff --git a/src/nrnoc/multicore.cpp b/src/nrnoc/multicore.cpp index ac2dbdd70e..badf2f9b27 100644 --- a/src/nrnoc/multicore.cpp +++ b/src/nrnoc/multicore.cpp @@ -346,7 +346,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; @@ -442,8 +442,8 @@ void nrn_threads_free() { nt->_ecell_children = NULL; } if (nt->_sp13mat) { - spDestroy(nt->_sp13mat); - nt->_sp13mat = 0; + free(nt->_sp13mat); + nt->_sp13mat = nullptr; } nt->end = 0; nt->ncell = 0; diff --git a/src/nrnoc/multicore.h b/src/nrnoc/multicore.h index 398e9f43c9..4d23e0187d 100644 --- a/src/nrnoc/multicore.h +++ b/src/nrnoc/multicore.h @@ -30,6 +30,8 @@ actual_v, etc. #include +#include + typedef struct NrnThreadMembList { /* patterned after CvMembList in cvodeobj.h */ struct NrnThreadMembList* next; Memb_list* ml; @@ -91,7 +93,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 */ + Eigen::SparseMatrix* _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 08b104c1b5..b6b19ab331 100644 --- a/src/nrnoc/solve.cpp +++ b/src/nrnoc/solve.cpp @@ -367,19 +367,13 @@ 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"); - } - } + // Use a copy because makeCompressed can invalidate pointers + Eigen::SparseMatrix mat = *_nt->_sp13mat; + mat.makeCompressed(); update_sp13_rhs_based_on_actual_rhs(_nt); - spSolve(_nt->_sp13mat, _nt->_sp13_rhs, _nt->_sp13_rhs); + Eigen::SparseLU> lu{mat}; + Eigen::Map rhs(_nt->_sp13_rhs + 1, mat.cols()); + rhs = lu.solve(rhs); 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 b3b698851c..4ae0c09706 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -28,8 +28,6 @@ #include #include -extern spREAL* spGetElement(char*, int, int); - int nrn_shape_changed_; /* for notifying Shape class in nrniv */ double* nrn_mech_wtime_; @@ -391,13 +389,12 @@ void nrn_rhs(neuron::model_sorted_token const& cache_token, NrnThread& nt) { } auto* const vec_rhs = nt.node_rhs_storage(); if (use_sparse13) { - int i, neqn; nrn_thread_error("nrn_rhs use_sparse13"); - neqn = spGetSize(_nt->_sp13mat, 0); - for (i = 1; i <= neqn; ++i) { + int neqn = _nt->_sp13mat->cols(); + for (int i = 1; i <= neqn; ++i) { _nt->_sp13_rhs[i] = 0.; } - for (i = i1; i < i3; ++i) { + for (int i = i1; i < i3; ++i) { NODERHS(_nt->_v_node[i]) = 0.; } } else { @@ -499,7 +496,14 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) { if (use_sparse13) { // Zero the sparse13 matrix - spClear(_nt->_sp13mat); + // We cannot setZero() because it invalidates pointers + // _nt->_sp13mat->setZero(); + auto m = _nt->_sp13mat; + for (int k = 0; k < m->outerSize(); ++k) { + for (Eigen::SparseMatrix::InnerIterator it(*m, k); it; ++it) { + it.valueRef() = 0.; + } + } } // Make sure the SoA node diagonals are also zeroed (is this needed?) @@ -1871,8 +1875,9 @@ void nrn_matrix_node_free() { free(std::exchange(nt->_sp13_rhs, nullptr)); } if (nt->_sp13mat) { - spDestroy(nt->_sp13mat); - nt->_sp13mat = (char*) 0; + nt->_sp13mat->setZero(); + nt->_sp13mat->data().squeeze(); + nt->_sp13mat = nullptr; } } diam_changed = 1; @@ -1935,14 +1940,6 @@ int nrn_method_consistent(void) { } -/* -sparse13 uses the mathematical index scheme in which the rows and columns -range from 1...size instead of 0...size-1. Thus the calls to -spGetElement(i,j) have args that are 1 greater than a normal c style. -Also the actual_rhs_ uses this style, 1-neqn, when sparse13 is activated -and therefore is passed to spSolve as actual_rhs intead of actual_rhs-1. -*/ - static void nrn_matrix_node_alloc(void) { int i, b; Node* nd; @@ -1969,7 +1966,7 @@ static void nrn_matrix_node_alloc(void) { } ++nrn_matrix_cnt_; if (use_sparse13) { - int in, err, extn, neqn, j; + int in, extn, neqn, j; nt = nrn_threads; neqn = nt->end + nrndae_extra_eqn_count(); extn = 0; @@ -1979,10 +1976,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 = spCreate(neqn, 0, &err); - if (err != spOKAY) { - hoc_execerror("Couldn't create sparse matrix", (char*) 0); - } + nt->_sp13mat = new Eigen::SparseMatrix(neqn, neqn); for (in = 0, i = 1; in < nt->end; ++in, ++i) { nt->_v_node[in]->eqn_index_ = i; if (nt->_v_node[in]->extnode) { @@ -1998,26 +1992,26 @@ static void nrn_matrix_node_alloc(void) { pnd = nt->_v_parent[in]; 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->coeffRef(i-1, i-1); if (nde) { for (ie = 0; ie < nlayer; ++ie) { - k = i + ie + 1; - nde->_d[ie] = spGetElement(nt->_sp13mat, k, k); + k = i + ie; + nde->_d[ie] = &nt->_sp13mat->coeffRef(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); + nde->_x21[ie] = &nt->_sp13mat->coeffRef(k, k - 1); + nde->_x12[ie] = &nt->_sp13mat->coeffRef(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); + nd->_a_matelm = &nt->_sp13mat->coeffRef(j-1, i-1); + nd->_b_matelm = &nt->_sp13mat->coeffRef(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); + int kp = j + ie; + k = i + ie; + nde->_a_matelm[ie] = &nt->_sp13mat->coeffRef(kp, k); + nde->_b_matelm[ie] = &nt->_sp13mat->coeffRef(k, kp); } } else { /* not needed if index starts at 1 */ nd->_a_matelm = nullptr; From 42fa3537ddb380003063b003dc3a0beb95287f88 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 14 Dec 2023 11:48:05 +0100 Subject: [PATCH 02/26] Bump Eigen to stable version 3.4.0 --- external/eigen | 2 +- src/ivoc/ocmatrix.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/external/eigen b/external/eigen index 328b5f9085..3147391d94 160000 --- a/external/eigen +++ b/external/eigen @@ -1 +1 @@ -Subproject commit 328b5f90858f93344ebc1484df8cadfd2a5da6dd +Subproject commit 3147391d946bb4b6c68edd901f2add6ac1f31f8c diff --git a/src/ivoc/ocmatrix.cpp b/src/ivoc/ocmatrix.cpp index 2fd427fe11..7afa87498a 100644 --- a/src/ivoc/ocmatrix.cpp +++ b/src/ivoc/ocmatrix.cpp @@ -128,7 +128,7 @@ void OcFullMatrix::symmeigen(Matrix* mout, Vect* vout) { void OcFullMatrix::svd1(Matrix* u, Matrix* v, Vect* d) { auto v1 = Vect2VEC(d); - Eigen::JacobiSVD svd(m_); + Eigen::JacobiSVD svd(m_, Eigen::ComputeFullU | Eigen::ComputeFullV); v1 = svd.singularValues(); if (u) { u->full()->m_ = svd.matrixU().transpose(); From 7b93bc576dffa2518dcaa632f724a37e5080bc51 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 14 Dec 2023 14:18:52 +0100 Subject: [PATCH 03/26] Fix indices --- src/nrnoc/treeset.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index 4ae0c09706..b8eba01840 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -499,7 +499,7 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) { // We cannot setZero() because it invalidates pointers // _nt->_sp13mat->setZero(); auto m = _nt->_sp13mat; - for (int k = 0; k < m->outerSize(); ++k) { + for (int k = 0; k < m->outerSize(); ++k) { for (Eigen::SparseMatrix::InnerIterator it(*m, k); it; ++it) { it.valueRef() = 0.; } @@ -1995,11 +1995,11 @@ static void nrn_matrix_node_alloc(void) { nd->_d_matelm = &nt->_sp13mat->coeffRef(i-1, i-1); if (nde) { for (ie = 0; ie < nlayer; ++ie) { - k = i + ie; - nde->_d[ie] = &nt->_sp13mat->coeffRef(k, k); + k = i + ie + 1; + nde->_d[ie] = &nt->_sp13mat->coeffRef(k - 1, k - 1); nde->_rhs[ie] = nt->_sp13_rhs + k; - nde->_x21[ie] = &nt->_sp13mat->coeffRef(k, k - 1); - nde->_x12[ie] = &nt->_sp13mat->coeffRef(k - 1, k); + nde->_x21[ie] = &nt->_sp13mat->coeffRef(k - 1, k - 2); + nde->_x12[ie] = &nt->_sp13mat->coeffRef(k - 2, k - 1); } } if (pnd) { @@ -2008,10 +2008,10 @@ static void nrn_matrix_node_alloc(void) { nd->_b_matelm = &nt->_sp13mat->coeffRef(i-1, j-1); if (nde && pnd->extnode) for (ie = 0; ie < nlayer; ++ie) { - int kp = j + ie; - k = i + ie; - nde->_a_matelm[ie] = &nt->_sp13mat->coeffRef(kp, k); - nde->_b_matelm[ie] = &nt->_sp13mat->coeffRef(k, kp); + int kp = j + ie + 1; + k = i + ie + 1; + nde->_a_matelm[ie] = &nt->_sp13mat->coeffRef(kp - 1, k - 1); + nde->_b_matelm[ie] = &nt->_sp13mat->coeffRef(k - 1, kp - 1); } } else { /* not needed if index starts at 1 */ nd->_a_matelm = nullptr; From e85f753cbb01ae81a7ae9b2ef39d80eb854314a4 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 14 Dec 2023 14:31:00 +0100 Subject: [PATCH 04/26] Use not sparse matrix because of pointers --- src/nrniv/matrixmap.cpp | 2 +- src/nrnoc/multicore.h | 4 ++-- src/nrnoc/solve.cpp | 7 ++----- src/nrnoc/treeset.cpp | 21 ++++++--------------- 4 files changed, 11 insertions(+), 23 deletions(-) diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp index 0fb76c2aff..d26b70531e 100644 --- a/src/nrniv/matrixmap.cpp +++ b/src/nrniv/matrixmap.cpp @@ -73,7 +73,7 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { jt = start + j - nnode; } // printf("MatrixMap::alloc getElement(%d,%d)\n", it, jt); - ptree_[plen_] = &_nt->_sp13mat->coeffRef(it-1, jt-1); + ptree_[plen_] = &_nt->_sp13mat->coeffRef(it - 1, jt - 1); ++plen_; } } diff --git a/src/nrnoc/multicore.h b/src/nrnoc/multicore.h index 4d23e0187d..daa8cbdf2e 100644 --- a/src/nrnoc/multicore.h +++ b/src/nrnoc/multicore.h @@ -30,7 +30,7 @@ actual_v, etc. #include -#include +#include typedef struct NrnThreadMembList { /* patterned after CvMembList in cvodeobj.h */ struct NrnThreadMembList* next; @@ -93,7 +93,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 */ - Eigen::SparseMatrix* _sp13mat; /* handle to general sparse matrix */ + Eigen::MatrixXd* _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 b6b19ab331..627b51aecd 100644 --- a/src/nrnoc/solve.cpp +++ b/src/nrnoc/solve.cpp @@ -367,12 +367,9 @@ void nrn_solve(NrnThread* _nt) { int e; nrn_thread_error("solve use_sparse13"); update_sp13_mat_based_on_actual_d(_nt); - // Use a copy because makeCompressed can invalidate pointers - Eigen::SparseMatrix mat = *_nt->_sp13mat; - mat.makeCompressed(); update_sp13_rhs_based_on_actual_rhs(_nt); - Eigen::SparseLU> lu{mat}; - Eigen::Map rhs(_nt->_sp13_rhs + 1, mat.cols()); + Eigen::FullPivLU lu{*_nt->_sp13mat}; + Eigen::Map rhs(_nt->_sp13_rhs + 1, _nt->_sp13mat->cols()); rhs = lu.solve(rhs); update_actual_d_based_on_sp13_mat(_nt); update_actual_rhs_based_on_sp13_rhs(_nt); diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index b8eba01840..56e8fd6b76 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -495,15 +495,7 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) { } if (use_sparse13) { - // Zero the sparse13 matrix - // We cannot setZero() because it invalidates pointers - // _nt->_sp13mat->setZero(); - auto m = _nt->_sp13mat; - for (int k = 0; k < m->outerSize(); ++k) { - for (Eigen::SparseMatrix::InnerIterator it(*m, k); it; ++it) { - it.valueRef() = 0.; - } - } + _nt->_sp13mat->setZero(); } // Make sure the SoA node diagonals are also zeroed (is this needed?) @@ -1875,8 +1867,7 @@ void nrn_matrix_node_free() { free(std::exchange(nt->_sp13_rhs, nullptr)); } if (nt->_sp13mat) { - nt->_sp13mat->setZero(); - nt->_sp13mat->data().squeeze(); + delete nt->_sp13mat; nt->_sp13mat = nullptr; } } @@ -1976,7 +1967,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 Eigen::SparseMatrix(neqn, neqn); + nt->_sp13mat = new Eigen::MatrixXd(neqn, neqn); for (in = 0, i = 1; in < nt->end; ++in, ++i) { nt->_v_node[in]->eqn_index_ = i; if (nt->_v_node[in]->extnode) { @@ -1992,7 +1983,7 @@ static void nrn_matrix_node_alloc(void) { pnd = nt->_v_parent[in]; i = nd->eqn_index_; nt->_sp13_rhs[i] = nt->actual_rhs(in); - nd->_d_matelm = &nt->_sp13mat->coeffRef(i-1, i-1); + nd->_d_matelm = &nt->_sp13mat->coeffRef(i - 1, i - 1); if (nde) { for (ie = 0; ie < nlayer; ++ie) { k = i + ie + 1; @@ -2004,8 +1995,8 @@ static void nrn_matrix_node_alloc(void) { } if (pnd) { j = pnd->eqn_index_; - nd->_a_matelm = &nt->_sp13mat->coeffRef(j-1, i-1); - nd->_b_matelm = &nt->_sp13mat->coeffRef(i-1, j-1); + nd->_a_matelm = &nt->_sp13mat->coeffRef(j - 1, i - 1); + nd->_b_matelm = &nt->_sp13mat->coeffRef(i - 1, j - 1); if (nde && pnd->extnode) for (ie = 0; ie < nlayer; ++ie) { int kp = j + ie + 1; From 6ce1dcb884e5e6d306b2f2e0706cc5da36416434 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 14 Dec 2023 15:25:04 +0100 Subject: [PATCH 05/26] _sp13_rhs 0-indexed --- src/nrniv/nrndae.cpp | 6 +++--- src/nrnoc/multicore.cpp | 2 +- src/nrnoc/multicore.h | 2 +- src/nrnoc/solve.cpp | 2 +- src/nrnoc/treeset.cpp | 12 ++++++------ 5 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/nrniv/nrndae.cpp b/src/nrniv/nrndae.cpp index 4bbe44d8f6..962e13e08b 100644 --- a/src/nrniv/nrndae.cpp +++ b/src/nrniv/nrndae.cpp @@ -221,7 +221,7 @@ void NrnDAE::dkmap(std::vector>& pv, y_.data() + i}; pvdot[bmap_[i] - 1] = neuron::container::data_handle{neuron::container::do_not_search, - _nt->_sp13_rhs + bmap_[i]}; + _nt->_sp13_rhs + bmap_[i] - 1}; } } @@ -231,7 +231,7 @@ void NrnDAE::update() { // note that the following is correct also for states that refer // to the internal potential of a segment. i.e rhs is v + vext[0] for (int i = 0; i < size_; ++i) { - y_[i] += _nt->_sp13_rhs[bmap_[i]]; + y_[i] += _nt->_sp13_rhs[bmap_[i] - 1]; } // for (int i=0; i < size_; ++i) printf(" i=%d bmap_[i]=%d y_[i]=%g\n", i, bmap_[i], // y_->elem(i)); @@ -311,7 +311,7 @@ void NrnDAE::rhs() { v2y(); f_(y_, yptmp_, size_); for (int i = 0; i < size_; ++i) { - _nt->_sp13_rhs[bmap_[i]] += yptmp_[i]; + _nt->_sp13_rhs[bmap_[i] - 1] += yptmp_[i]; } } diff --git a/src/nrnoc/multicore.cpp b/src/nrnoc/multicore.cpp index badf2f9b27..2092db836d 100644 --- a/src/nrnoc/multicore.cpp +++ b/src/nrnoc/multicore.cpp @@ -442,7 +442,7 @@ void nrn_threads_free() { nt->_ecell_children = NULL; } if (nt->_sp13mat) { - free(nt->_sp13mat); + delete nt->_sp13mat; nt->_sp13mat = nullptr; } nt->end = 0; diff --git a/src/nrnoc/multicore.h b/src/nrnoc/multicore.h index daa8cbdf2e..bfa100d967 100644 --- a/src/nrnoc/multicore.h +++ b/src/nrnoc/multicore.h @@ -93,7 +93,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 */ - Eigen::MatrixXd* _sp13mat; /* handle to general sparse matrix */ + Eigen::MatrixXd* _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 627b51aecd..72128c0915 100644 --- a/src/nrnoc/solve.cpp +++ b/src/nrnoc/solve.cpp @@ -369,7 +369,7 @@ void nrn_solve(NrnThread* _nt) { update_sp13_mat_based_on_actual_d(_nt); update_sp13_rhs_based_on_actual_rhs(_nt); Eigen::FullPivLU lu{*_nt->_sp13mat}; - Eigen::Map rhs(_nt->_sp13_rhs + 1, _nt->_sp13mat->cols()); + Eigen::Map rhs(_nt->_sp13_rhs, _nt->_sp13mat->cols()); rhs = lu.solve(rhs); update_actual_d_based_on_sp13_mat(_nt); update_actual_rhs_based_on_sp13_rhs(_nt); diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index 56e8fd6b76..ccbddbe169 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -332,7 +332,7 @@ vm += dvi-dvx */ void update_actual_rhs_based_on_sp13_rhs(NrnThread* nt) { for (int i = 0; i < nt->end; i++) { - nt->actual_rhs(i) = nt->_sp13_rhs[nt->_v_node[i]->eqn_index_]; + nt->actual_rhs(i) = nt->_sp13_rhs[nt->_v_node[i]->eqn_index_ - 1]; } } @@ -341,7 +341,7 @@ void update_actual_rhs_based_on_sp13_rhs(NrnThread* nt) { */ void update_sp13_rhs_based_on_actual_rhs(NrnThread* nt) { for (int i = 0; i < nt->end; i++) { - nt->_sp13_rhs[nt->_v_node[i]->eqn_index_] = nt->actual_rhs(i); + nt->_sp13_rhs[nt->_v_node[i]->eqn_index_ - 1] = nt->actual_rhs(i); } } @@ -391,7 +391,7 @@ void nrn_rhs(neuron::model_sorted_token const& cache_token, NrnThread& nt) { if (use_sparse13) { nrn_thread_error("nrn_rhs use_sparse13"); int neqn = _nt->_sp13mat->cols(); - for (int i = 1; i <= neqn; ++i) { + for (int i = 0; i < neqn; ++i) { _nt->_sp13_rhs[i] = 0.; } for (int i = i1; i < i3; ++i) { @@ -1966,7 +1966,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->_sp13_rhs = (double*) ecalloc(neqn, sizeof(double)); nt->_sp13mat = new Eigen::MatrixXd(neqn, neqn); for (in = 0, i = 1; in < nt->end; ++in, ++i) { nt->_v_node[in]->eqn_index_ = i; @@ -1982,13 +1982,13 @@ static void nrn_matrix_node_alloc(void) { nde = nd->extnode; pnd = nt->_v_parent[in]; i = nd->eqn_index_; - nt->_sp13_rhs[i] = nt->actual_rhs(in); + nt->_sp13_rhs[i - 1] = nt->actual_rhs(in); nd->_d_matelm = &nt->_sp13mat->coeffRef(i - 1, i - 1); if (nde) { for (ie = 0; ie < nlayer; ++ie) { k = i + ie + 1; nde->_d[ie] = &nt->_sp13mat->coeffRef(k - 1, k - 1); - nde->_rhs[ie] = nt->_sp13_rhs + k; + nde->_rhs[ie] = nt->_sp13_rhs + k - 1; nde->_x21[ie] = &nt->_sp13mat->coeffRef(k - 1, k - 2); nde->_x12[ie] = &nt->_sp13mat->coeffRef(k - 2, k - 1); } From 20e595bc9394140af3a0ff54e9d9ea2c1f8c7909 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 14 Dec 2023 17:47:04 +0100 Subject: [PATCH 06/26] Let avoid include headers --- src/nrncvode/occvode.cpp | 2 ++ src/nrniv/matrixmap.cpp | 2 +- src/nrniv/nrndae.cpp | 2 ++ src/nrnoc/fadvance.cpp | 2 +- src/nrnoc/multicore.h | 6 +++++- src/nrnoc/solve.cpp | 2 +- src/nrnoc/treeset.cpp | 2 +- 7 files changed, 13 insertions(+), 5 deletions(-) diff --git a/src/nrncvode/occvode.cpp b/src/nrncvode/occvode.cpp index 854e1972d4..c748b895b5 100644 --- a/src/nrncvode/occvode.cpp +++ b/src/nrncvode/occvode.cpp @@ -14,6 +14,8 @@ #include +#include + extern void setup_topology(), v_setup_vectors(); extern void recalc_diam(); extern int nrn_errno_check(int); diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp index d26b70531e..a0d426fb26 100644 --- a/src/nrniv/matrixmap.cpp +++ b/src/nrniv/matrixmap.cpp @@ -3,7 +3,7 @@ #include using std::vector; -#include "spmatrix.h" +#include MatrixMap::MatrixMap(Matrix& mat) : m_(mat) diff --git a/src/nrniv/nrndae.cpp b/src/nrniv/nrndae.cpp index 962e13e08b..58f8df170a 100644 --- a/src/nrniv/nrndae.cpp +++ b/src/nrniv/nrndae.cpp @@ -6,6 +6,8 @@ #include "treeset.h" #include "utils/enumerate.h" +#include + extern int secondorder; static NrnDAEPtrList nrndae_list; diff --git a/src/nrnoc/fadvance.cpp b/src/nrnoc/fadvance.cpp index a0a032ce1e..64c59ab4bc 100644 --- a/src/nrnoc/fadvance.cpp +++ b/src/nrnoc/fadvance.cpp @@ -12,7 +12,7 @@ #include "utils/profile/profiler_interface.h" #include "nonvintblock.h" #include "nrncvode.h" -#include "spmatrix.h" +#include #include diff --git a/src/nrnoc/multicore.h b/src/nrnoc/multicore.h index bfa100d967..58c50a66fb 100644 --- a/src/nrnoc/multicore.h +++ b/src/nrnoc/multicore.h @@ -30,7 +30,11 @@ actual_v, etc. #include -#include +namespace Eigen { + template + class Matrix; + using MatrixXd = Matrix; +} typedef struct NrnThreadMembList { /* patterned after CvMembList in cvodeobj.h */ struct NrnThreadMembList* next; diff --git a/src/nrnoc/solve.cpp b/src/nrnoc/solve.cpp index 72128c0915..d5862f9375 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 #include "treeset.h" #include diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index ccbddbe169..2f03da2839 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -16,7 +16,7 @@ #include "ocnotify.h" #include "partrans.h" #include "section.h" -#include "spmatrix.h" +#include #include "utils/profile/profiler_interface.h" #include "multicore.h" From 2abf5003cd8f942838985ca22b0fa29a8994a031 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 14 Dec 2023 17:49:27 +0100 Subject: [PATCH 07/26] Format --- src/nrnoc/multicore.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/nrnoc/multicore.h b/src/nrnoc/multicore.h index 58c50a66fb..6e12c51a65 100644 --- a/src/nrnoc/multicore.h +++ b/src/nrnoc/multicore.h @@ -31,10 +31,10 @@ actual_v, etc. #include namespace Eigen { - template - class Matrix; - using MatrixXd = Matrix; -} +template +class Matrix; +using MatrixXd = Matrix; +} // namespace Eigen typedef struct NrnThreadMembList { /* patterned after CvMembList in cvodeobj.h */ struct NrnThreadMembList* next; From 6e03f48a9d288017f811a297623855668552ea97 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 14 Dec 2023 17:58:21 +0100 Subject: [PATCH 08/26] Format --- src/nrnoc/multicore.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/nrnoc/multicore.cpp b/src/nrnoc/multicore.cpp index 2092db836d..0efab3460e 100644 --- a/src/nrnoc/multicore.cpp +++ b/src/nrnoc/multicore.cpp @@ -47,6 +47,8 @@ the handling of v_structure_change as long as possible. #include #include +#include + #define CACHELINE_ALLOC(name, type, size) \ name = (type*) nrn_cacheline_alloc((void**) &name, size * sizeof(type)) #define CACHELINE_CALLOC(name, type, size) \ From 135fc529175853ab6930e20dffa7dcf208740a5f Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 4 Jun 2024 11:02:15 +0200 Subject: [PATCH 09/26] more -1+1 --- src/nrncvode/occvode.cpp | 3 --- src/nrniv/matrixmap.cpp | 2 +- src/nrniv/nrndae.cpp | 16 ++++++------- src/nrnoc/section.h | 2 +- src/nrnoc/treeset.cpp | 49 +++++++++++++++++++++------------------- 5 files changed, 36 insertions(+), 36 deletions(-) diff --git a/src/nrncvode/occvode.cpp b/src/nrncvode/occvode.cpp index c748b895b5..1b7b705cfe 100644 --- a/src/nrncvode/occvode.cpp +++ b/src/nrncvode/occvode.cpp @@ -25,9 +25,6 @@ extern int nrn_errno_check(int); extern Symlist* hoc_built_in_symlist; -#include "spmatrix.h" -extern double* sp13mat; - #if 1 || NRNMPI extern void (*nrnthread_v_transfer_)(NrnThread*); extern void (*nrnmpi_v_transfer_)(); diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp index a0d426fb26..0739291587 100644 --- a/src/nrniv/matrixmap.cpp +++ b/src/nrniv/matrixmap.cpp @@ -73,7 +73,7 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { jt = start + j - nnode; } // printf("MatrixMap::alloc getElement(%d,%d)\n", it, jt); - ptree_[plen_] = &_nt->_sp13mat->coeffRef(it - 1, jt - 1); + ptree_[plen_] = &_nt->_sp13mat->coeffRef(it, jt); ++plen_; } } diff --git a/src/nrniv/nrndae.cpp b/src/nrniv/nrndae.cpp index 58f8df170a..3bf18f7b5a 100644 --- a/src/nrniv/nrndae.cpp +++ b/src/nrniv/nrndae.cpp @@ -49,7 +49,7 @@ void nrndae_alloc() { neqn += _nt->_ecell_memb_list->nodecount * nlayer; } for (NrnDAE* item: nrndae_list) { - item->alloc(neqn + 1); + item->alloc(neqn); neqn += item->extra_eqn_count(); } } @@ -219,11 +219,11 @@ void NrnDAE::dkmap(std::vector>& pv, NrnThread* _nt = nrn_threads; for (int i = nnode_; i < size_; ++i) { // printf("bmap_[%d] = %d\n", i, bmap_[i]); - pv[bmap_[i] - 1] = neuron::container::data_handle{neuron::container::do_not_search, + pv[bmap_[i]] = neuron::container::data_handle{neuron::container::do_not_search, y_.data() + i}; - pvdot[bmap_[i] - 1] = + pvdot[bmap_[i]] = neuron::container::data_handle{neuron::container::do_not_search, - _nt->_sp13_rhs + bmap_[i] - 1}; + _nt->_sp13_rhs + bmap_[i]}; } } @@ -233,7 +233,7 @@ void NrnDAE::update() { // note that the following is correct also for states that refer // to the internal potential of a segment. i.e rhs is v + vext[0] for (int i = 0; i < size_; ++i) { - y_[i] += _nt->_sp13_rhs[bmap_[i] - 1]; + y_[i] += _nt->_sp13_rhs[bmap_[i]]; } // for (int i=0; i < size_; ++i) printf(" i=%d bmap_[i]=%d y_[i]=%g\n", i, bmap_[i], // y_->elem(i)); @@ -289,7 +289,7 @@ void NrnDAE::dkres(double* y, double* yprime, double* delta) { for (int i = 0; i < size_; ++i) { // printf("%d %d %g %g\n", i, bmap_[i]-1, y[bmap_[i]-1], yprime[bmap_[i]-1]); - yptmp_[i] = yprime[bmap_[i] - 1]; + yptmp_[i] = yprime[bmap_[i]]; } Vect* cyp; if (assumed_identity_) { @@ -303,7 +303,7 @@ void NrnDAE::dkres(double* y, double* yprime, double* delta) { cyp = &cyp_; } for (int i = 0; i < size_; ++i) { - delta[bmap_[i] - 1] -= (*cyp)[i]; + delta[bmap_[i]] -= (*cyp)[i]; } } @@ -313,7 +313,7 @@ void NrnDAE::rhs() { v2y(); f_(y_, yptmp_, size_); for (int i = 0; i < size_; ++i) { - _nt->_sp13_rhs[bmap_[i] - 1] += yptmp_[i]; + _nt->_sp13_rhs[bmap_[i]] += yptmp_[i]; } } diff --git a/src/nrnoc/section.h b/src/nrnoc/section.h index 0f496d1f46..27620c8ccd 100644 --- a/src/nrnoc/section.h +++ b/src/nrnoc/section.h @@ -182,7 +182,7 @@ struct Node { double _rinv{}; /* conductance uS from node to parent */ double* _a_matelm; double* _b_matelm; - double* _d_matelm; + int diag_value; 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 6e44d6c7bb..564f381a68 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -333,7 +333,7 @@ vm += dvi-dvx */ void update_actual_rhs_based_on_sp13_rhs(NrnThread* nt) { for (int i = 0; i < nt->end; i++) { - nt->actual_rhs(i) = nt->_sp13_rhs[nt->_v_node[i]->eqn_index_ - 1]; + nt->actual_rhs(i) = nt->_sp13_rhs[nt->_v_node[i]->eqn_index_]; } } @@ -342,7 +342,7 @@ void update_actual_rhs_based_on_sp13_rhs(NrnThread* nt) { */ void update_sp13_rhs_based_on_actual_rhs(NrnThread* nt) { for (int i = 0; i < nt->end; i++) { - nt->_sp13_rhs[nt->_v_node[i]->eqn_index_ - 1] = nt->actual_rhs(i); + nt->_sp13_rhs[nt->_v_node[i]->eqn_index_] = nt->actual_rhs(i); } } @@ -351,7 +351,8 @@ 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; + int index = nt->_v_node[i]->diag_value; + nt->actual_d(i) = nt->_sp13mat->coeffRef(index, index); } } @@ -360,7 +361,8 @@ 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); + int index = nt->_v_node[i]->diag_value; + nt->_sp13mat->coeffRef(index, index) = nt->actual_d(i); } } @@ -574,9 +576,11 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) { // 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; + int index = nd->diag_value; + _nt->_sp13mat->coeffRef(index, index) -= nd_b; // used to update NODED (sparse13 matrix) using NODEA and NODEB ("SoA") - *parent_nd->_d_matelm -= nd_a; + index = parent_nd->diag_value; + _nt->_sp13mat->coeffRef(index, index) -= 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; @@ -1987,41 +1991,40 @@ static void nrn_matrix_node_alloc(void) { neqn += extn; nt->_sp13_rhs = (double*) ecalloc(neqn, sizeof(double)); nt->_sp13mat = new Eigen::MatrixXd(neqn, neqn); - for (in = 0, i = 1; in < nt->end; ++in, ++i) { + for (in = 0, i = 0; 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_; - nt->_sp13_rhs[i - 1] = nt->actual_rhs(in); - nd->_d_matelm = &nt->_sp13mat->coeffRef(i - 1, i - 1); + nt->_sp13_rhs[i] = nt->actual_rhs(in); + nd->diag_value = i; if (nde) { - for (ie = 0; ie < nlayer; ++ie) { - k = i + ie + 1; - nde->_d[ie] = &nt->_sp13mat->coeffRef(k - 1, k - 1); - nde->_rhs[ie] = nt->_sp13_rhs + k - 1; - nde->_x21[ie] = &nt->_sp13mat->coeffRef(k - 1, k - 2); - nde->_x12[ie] = &nt->_sp13mat->coeffRef(k - 2, k - 1); + for (int ie = 0; ie < nlayer; ++ie) { + int k = i + ie; + nde->_d[ie] = &nt->_sp13mat->coeffRef(k, k); + nde->_rhs[ie] = nt->_sp13_rhs + k; + nde->_x21[ie] = &nt->_sp13mat->coeffRef(k, k - 1); + nde->_x12[ie] = &nt->_sp13mat->coeffRef(k - 1, k); } } if (pnd) { j = pnd->eqn_index_; - nd->_a_matelm = &nt->_sp13mat->coeffRef(j - 1, i - 1); - nd->_b_matelm = &nt->_sp13mat->coeffRef(i - 1, j - 1); + nd->_a_matelm = &nt->_sp13mat->coeffRef(j, i); + nd->_b_matelm = &nt->_sp13mat->coeffRef(i, j); if (nde && pnd->extnode) - for (ie = 0; ie < nlayer; ++ie) { - int kp = j + ie + 1; - k = i + ie + 1; - nde->_a_matelm[ie] = &nt->_sp13mat->coeffRef(kp - 1, k - 1); - nde->_b_matelm[ie] = &nt->_sp13mat->coeffRef(k - 1, kp - 1); + for (int ie = 0; ie < nlayer; ++ie) { + int kp = j + ie; + int k = i + ie; + nde->_a_matelm[ie] = &nt->_sp13mat->coeffRef(kp, k); + nde->_b_matelm[ie] = &nt->_sp13mat->coeffRef(k, kp); } } else { /* not needed if index starts at 1 */ nd->_a_matelm = nullptr; From d52b91e14bdc3cf9bbf93d94f435c2d66f17cd87 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Tue, 4 Jun 2024 11:58:11 +0000 Subject: [PATCH 10/26] Fix formatting --- src/nrniv/nrndae.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/nrniv/nrndae.cpp b/src/nrniv/nrndae.cpp index 3bf18f7b5a..4b0db8ad5f 100644 --- a/src/nrniv/nrndae.cpp +++ b/src/nrniv/nrndae.cpp @@ -220,10 +220,9 @@ void NrnDAE::dkmap(std::vector>& pv, for (int i = nnode_; i < size_; ++i) { // printf("bmap_[%d] = %d\n", i, bmap_[i]); pv[bmap_[i]] = neuron::container::data_handle{neuron::container::do_not_search, - y_.data() + i}; - pvdot[bmap_[i]] = - neuron::container::data_handle{neuron::container::do_not_search, - _nt->_sp13_rhs + bmap_[i]}; + y_.data() + i}; + pvdot[bmap_[i]] = neuron::container::data_handle{neuron::container::do_not_search, + _nt->_sp13_rhs + bmap_[i]}; } } From c4ae1a066adb77496fd857f670087bff9c5c2179 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 4 Jun 2024 17:09:29 +0200 Subject: [PATCH 11/26] Simplify MatrixMap --- src/nrniv/matrixmap.cpp | 37 +++++++++---------------------------- src/nrniv/matrixmap.h | 7 +++---- 2 files changed, 12 insertions(+), 32 deletions(-) diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp index 02464c40a7..f07cfbed2b 100644 --- a/src/nrniv/matrixmap.cpp +++ b/src/nrniv/matrixmap.cpp @@ -1,37 +1,24 @@ #include <../../nrnconf.h> #include "matrixmap.h" #include -using std::vector; #include "spmatrix.h" MatrixMap::MatrixMap(Matrix& mat) - : m_(mat) - , plen_(0) - , ptree_(NULL) - , pm_(NULL) {} + : m_(mat) {} MatrixMap::MatrixMap(Matrix* mat) - : m_(*mat) - , plen_(0) - , ptree_(NULL) - , pm_(NULL) {} - -MatrixMap::~MatrixMap() { - mmfree(); -} + : m_(*mat) {} void MatrixMap::mmfree() { - // safe to delete NULL pointers - delete[] ptree_; - delete[] pm_; - ptree_ = NULL; - pm_ = NULL; + pm_.clear(); + pm_.shrink_to_fit(); + ptree_.clear(); + ptree_.shrink_to_fit(); } void MatrixMap::add(double fac) { for (int i = 0; i < plen_; ++i) { - // printf("i=%d %g += %g * %g\n", i, *ptree_[i], fac, *pm_[i]); *ptree_[i] += fac * (*pm_[i]); } } @@ -39,23 +26,18 @@ void MatrixMap::add(double fac) { void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { NrnThread* _nt = nrn_threads; mmfree(); - // how many elements - int nrow = m_.nrow(); - int ncol = m_.ncol(); - // printf("MatrixMap::alloc nrow=%d ncol=%d\n", nrow, ncol); plen_ = 0; - vector nonzero_i, nonzero_j; + std::vector nonzero_i, nonzero_j; m_.nonzeros(nonzero_i, nonzero_j); - pm_ = new double*[nonzero_i.size()]; - ptree_ = new double*[nonzero_i.size()]; + 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]; int it; if (i < nnode) { it = nodes[i]->eqn_index_ + layer[i]; - // printf("i=%d it=%d area=%g\n", i, it, nodes[i]->area); if (layer[i] > 0 && !nodes[i]->extnode) { it = 0; } @@ -72,7 +54,6 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { } else { jt = start + j - nnode; } - // printf("MatrixMap::alloc getElement(%d,%d)\n", it, jt); ptree_[plen_] = spGetElement(_nt->_sp13mat, it, jt); ++plen_; } diff --git a/src/nrniv/matrixmap.h b/src/nrniv/matrixmap.h index eba8b892d7..f216dcc9c3 100644 --- a/src/nrniv/matrixmap.h +++ b/src/nrniv/matrixmap.h @@ -7,7 +7,6 @@ class MatrixMap { public: MatrixMap(Matrix*); MatrixMap(Matrix&); - ~MatrixMap(); void alloc(int, int, Node**, int*); void mmfree(); @@ -116,7 +115,7 @@ class MatrixMap { Matrix& m_; // the map - int plen_; - double** pm_; - double** ptree_; + int plen_ = 0; + std::vector pm_{}; + std::vector ptree_{}; }; From 94276eb30cf331e65b29fa1c989a50ace8f78cec Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 5 Jun 2024 13:00:43 +0200 Subject: [PATCH 12/26] Revert change --- src/nrniv/matrixmap.cpp | 36 +++++++++++++++++++++++++++--------- src/nrniv/matrixmap.h | 7 ++++--- src/nrniv/multisplit.cpp | 5 +++++ 3 files changed, 36 insertions(+), 12 deletions(-) diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp index 3696f482c5..918ecd72a7 100644 --- a/src/nrniv/matrixmap.cpp +++ b/src/nrniv/matrixmap.cpp @@ -1,24 +1,37 @@ #include <../../nrnconf.h> #include "matrixmap.h" #include +using std::vector; #include MatrixMap::MatrixMap(Matrix& mat) - : m_(mat) {} + : m_(mat) + , plen_(0) + , ptree_(NULL) + , pm_(NULL) {} MatrixMap::MatrixMap(Matrix* mat) - : m_(*mat) {} + : m_(*mat) + , plen_(0) + , ptree_(NULL) + , pm_(NULL) {} + +MatrixMap::~MatrixMap() { + mmfree(); +} void MatrixMap::mmfree() { - pm_.clear(); - pm_.shrink_to_fit(); - ptree_.clear(); - ptree_.shrink_to_fit(); + // safe to delete NULL pointers + delete[] ptree_; + delete[] pm_; + ptree_ = NULL; + pm_ = NULL; } void MatrixMap::add(double fac) { for (int i = 0; i < plen_; ++i) { + // printf("i=%d %g += %g * %g\n", i, *ptree_[i], fac, *pm_[i]); *ptree_[i] += fac * (*pm_[i]); } } @@ -26,18 +39,23 @@ void MatrixMap::add(double fac) { void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { NrnThread* _nt = nrn_threads; mmfree(); + // how many elements + int nrow = m_.nrow(); + int ncol = m_.ncol(); + // printf("MatrixMap::alloc nrow=%d ncol=%d\n", nrow, ncol); plen_ = 0; - std::vector nonzero_i, nonzero_j; + vector nonzero_i, nonzero_j; m_.nonzeros(nonzero_i, nonzero_j); - pm_.resize(nonzero_i.size()); - ptree_.resize(nonzero_i.size()); + pm_ = new double*[nonzero_i.size()]; + ptree_ = new double*[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]; + // printf("i=%d it=%d area=%g\n", i, it, nodes[i]->area); if (layer[i] > 0 && !nodes[i]->extnode) { it = 0; } diff --git a/src/nrniv/matrixmap.h b/src/nrniv/matrixmap.h index f216dcc9c3..eba8b892d7 100644 --- a/src/nrniv/matrixmap.h +++ b/src/nrniv/matrixmap.h @@ -7,6 +7,7 @@ class MatrixMap { public: MatrixMap(Matrix*); MatrixMap(Matrix&); + ~MatrixMap(); void alloc(int, int, Node**, int*); void mmfree(); @@ -115,7 +116,7 @@ class MatrixMap { Matrix& m_; // the map - int plen_ = 0; - std::vector pm_{}; - std::vector ptree_{}; + int plen_; + double** pm_; + double** ptree_; }; diff --git a/src/nrniv/multisplit.cpp b/src/nrniv/multisplit.cpp index 11ea02b357..3ca23e4dcc 100644 --- a/src/nrniv/multisplit.cpp +++ b/src/nrniv/multisplit.cpp @@ -3367,6 +3367,11 @@ for (i=i1; i < i3; ++i) { nt->_v_node[i]->eqn_index_ = -1; } for (i = i1; i < backbone_end; ++i) { +#if 0 +printf("backbone i=%d %d %s %d", i, node[i]->v_node_index, secname(node[i]->sec), node[i]->sec_node_index_); +printf(" -> %s %d\n", parent[i]?secname(parent[i]->sec):"root", +parent[i]?parent[i]->sec_node_index_:-1); +#endif node[i - i1]->eqn_index_ = i; } // the rest are in order From 92b08570f8744feb654a8f5466ed742fc0ceac6e Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 5 Jun 2024 13:19:18 +0200 Subject: [PATCH 13/26] Simplify & fix weird things --- src/nrncvode/nrndaspk.cpp | 2 +- src/nrncvode/occvode.cpp | 2 +- src/nrniv/matrixmap.cpp | 2 +- src/nrniv/nrndae.cpp | 30 ++++++++++++++++++++--------- src/nrnoc/fadvance.cpp | 2 +- src/nrnoc/section.h | 2 +- src/nrnoc/solve.cpp | 2 +- src/nrnoc/treeset.cpp | 40 ++++++++++++++++++--------------------- 8 files changed, 45 insertions(+), 37 deletions(-) diff --git a/src/nrncvode/nrndaspk.cpp b/src/nrncvode/nrndaspk.cpp index 7087d82351..5434cdc031 100644 --- a/src/nrncvode/nrndaspk.cpp +++ b/src/nrncvode/nrndaspk.cpp @@ -500,7 +500,7 @@ for (i=0; i < z.nvsize_; ++i) { auto const vec_sav_rhs = nt->node_sav_rhs_storage(); for (i = 0; i < n; ++i) { Node* nd = ml->nodelist[i]; - int j = nd->eqn_index_; + int j = nd->eqn_index_ - 1; Extnode* nde = nd->extnode; double cdvm; if (nde) { diff --git a/src/nrncvode/occvode.cpp b/src/nrncvode/occvode.cpp index 5a1d686f61..1b7b705cfe 100644 --- a/src/nrncvode/occvode.cpp +++ b/src/nrncvode/occvode.cpp @@ -400,7 +400,7 @@ void Cvode::daspk_init_eqn() { Extnode* nde; nd = _nt->_v_node[in]; nde = nd->extnode; - i = nd->eqn_index_; // the sparse matrix index starts at 1 + i = nd->eqn_index_ - 1; // the sparse matrix index starts at 1 z.pv_[i] = nd->v_handle(); z.pvdot_[i] = nd->rhs_handle(); if (nde) { diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp index 918ecd72a7..dd175f31c3 100644 --- a/src/nrniv/matrixmap.cpp +++ b/src/nrniv/matrixmap.cpp @@ -72,7 +72,7 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { } else { jt = start + j - nnode; } - ptree_[plen_] = &_nt->_sp13mat->coeffRef(it, jt); + ptree_[plen_] = &_nt->_sp13mat->coeffRef(it - 1, jt - 1); ++plen_; } } diff --git a/src/nrniv/nrndae.cpp b/src/nrniv/nrndae.cpp index ce4ea9e08b..4bbe44d8f6 100644 --- a/src/nrniv/nrndae.cpp +++ b/src/nrniv/nrndae.cpp @@ -6,8 +6,6 @@ #include "treeset.h" #include "utils/enumerate.h" -#include - extern int secondorder; static NrnDAEPtrList nrndae_list; @@ -49,7 +47,7 @@ void nrndae_alloc() { neqn += _nt->_ecell_memb_list->nodecount * nlayer; } for (NrnDAE* item: nrndae_list) { - item->alloc(neqn); + item->alloc(neqn + 1); neqn += item->extra_eqn_count(); } } @@ -200,35 +198,48 @@ NrnDAE::~NrnDAE() { if (elayer_) { delete[] elayer_; } + // if (nrndae_list->count() == 0) { + // use_sparse13 = 0; + // } nrn_matrix_node_free(); } int NrnDAE::extra_eqn_count() { + // printf("NrnDAE::extra_eqn_count %lx\n", (long)this); + // printf(" nnode_=%d g_->nrow()=%d\n", nnode_, g_->nrow()); return c_->nrow() - nnode_; } void NrnDAE::dkmap(std::vector>& pv, std::vector>& pvdot) { + // printf("NrnDAE::dkmap\n"); NrnThread* _nt = nrn_threads; for (int i = nnode_; i < size_; ++i) { - pv[bmap_[i]] = neuron::container::data_handle{neuron::container::do_not_search, - y_.data() + i}; - pvdot[bmap_[i]] = neuron::container::data_handle{neuron::container::do_not_search, - _nt->_sp13_rhs + bmap_[i]}; + // printf("bmap_[%d] = %d\n", i, bmap_[i]); + pv[bmap_[i] - 1] = neuron::container::data_handle{neuron::container::do_not_search, + y_.data() + i}; + pvdot[bmap_[i] - 1] = + neuron::container::data_handle{neuron::container::do_not_search, + _nt->_sp13_rhs + bmap_[i]}; } } void NrnDAE::update() { + // printf("NrnDAE::update %lx\n", (long)this); NrnThread* _nt = nrn_threads; // note that the following is correct also for states that refer // to the internal potential of a segment. i.e rhs is v + vext[0] for (int i = 0; i < size_; ++i) { y_[i] += _nt->_sp13_rhs[bmap_[i]]; } + // for (int i=0; i < size_; ++i) printf(" i=%d bmap_[i]=%d y_[i]=%g\n", i, bmap_[i], + // y_->elem(i)); } void NrnDAE::init() { + // printf("NrnDAE::init %lx\n", (long)this); + // printf("init size_=%d %d %d %d\n", size_, y_->size(), y0_->size(), b_->size()); Vect& y0 = *y0_; v2y(); @@ -275,7 +286,8 @@ void NrnDAE::dkres(double* y, double* yprime, double* delta) { // representation of the y and yprime for (int i = 0; i < size_; ++i) { - yptmp_[i] = yprime[bmap_[i]]; + // printf("%d %d %g %g\n", i, bmap_[i]-1, y[bmap_[i]-1], yprime[bmap_[i]-1]); + yptmp_[i] = yprime[bmap_[i] - 1]; } Vect* cyp; if (assumed_identity_) { @@ -289,7 +301,7 @@ void NrnDAE::dkres(double* y, double* yprime, double* delta) { cyp = &cyp_; } for (int i = 0; i < size_; ++i) { - delta[bmap_[i]] -= (*cyp)[i]; + delta[bmap_[i] - 1] -= (*cyp)[i]; } } diff --git a/src/nrnoc/fadvance.cpp b/src/nrnoc/fadvance.cpp index 91f1f08ff0..e064e1cd88 100644 --- a/src/nrnoc/fadvance.cpp +++ b/src/nrnoc/fadvance.cpp @@ -676,7 +676,7 @@ void nrn_print_matrix(NrnThread* _nt) { } else { std::cout << &_nt->_sp13mat << std::endl; int n = _nt->_sp13mat->cols(); - for (int i = 0; i < n; ++i) { + for (int i = 1; i <= n; ++i) { Printf("%d %g\n", i, _nt->actual_rhs(i)); } } diff --git a/src/nrnoc/section.h b/src/nrnoc/section.h index 27620c8ccd..0f496d1f46 100644 --- a/src/nrnoc/section.h +++ b/src/nrnoc/section.h @@ -182,7 +182,7 @@ struct Node { double _rinv{}; /* conductance uS from node to parent */ double* _a_matelm; double* _b_matelm; - int diag_value; + 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/solve.cpp b/src/nrnoc/solve.cpp index fbe8e6e0df..508f7763c6 100644 --- a/src/nrnoc/solve.cpp +++ b/src/nrnoc/solve.cpp @@ -369,7 +369,7 @@ void nrn_solve(NrnThread* _nt) { update_sp13_mat_based_on_actual_d(_nt); update_sp13_rhs_based_on_actual_rhs(_nt); Eigen::FullPivLU lu{*_nt->_sp13mat}; - Eigen::Map rhs(_nt->_sp13_rhs, _nt->_sp13mat->cols()); + Eigen::Map rhs(_nt->_sp13_rhs + 1, _nt->_sp13mat->cols()); rhs = lu.solve(rhs); update_actual_d_based_on_sp13_mat(_nt); update_actual_rhs_based_on_sp13_rhs(_nt); diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index 564f381a68..83d3384591 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -351,8 +351,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) { - int index = nt->_v_node[i]->diag_value; - nt->actual_d(i) = nt->_sp13mat->coeffRef(index, index); + nt->actual_d(i) = *nt->_v_node[i]->_d_matelm; } } @@ -361,8 +360,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) { - int index = nt->_v_node[i]->diag_value; - nt->_sp13mat->coeffRef(index, index) = nt->actual_d(i); + *nt->_v_node[i]->_d_matelm = nt->actual_d(i); } } @@ -394,7 +392,7 @@ void nrn_rhs(neuron::model_sorted_token const& cache_token, NrnThread& nt) { if (use_sparse13) { nrn_thread_error("nrn_rhs use_sparse13"); int neqn = _nt->_sp13mat->cols(); - for (int i = 0; i < neqn; ++i) { + for (int i = 1; i <= neqn; ++i) { _nt->_sp13_rhs[i] = 0.; } for (int i = i1; i < i3; ++i) { @@ -576,11 +574,9 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) { // Update entries in sp13_mat *nd->_a_matelm += nd_a; *nd->_b_matelm += nd_b; /* b may have value from lincir */ - int index = nd->diag_value; - _nt->_sp13mat->coeffRef(index, index) -= nd_b; + *nd->_d_matelm -= nd_b; // used to update NODED (sparse13 matrix) using NODEA and NODEB ("SoA") - index = parent_nd->diag_value; - _nt->_sp13mat->coeffRef(index, index) -= nd_a; + *parent_nd->_d_matelm -= 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; @@ -1989,9 +1985,9 @@ static void nrn_matrix_node_alloc(void) { } /*printf(" %d extracellular nodes\n", extn);*/ neqn += extn; - nt->_sp13_rhs = (double*) ecalloc(neqn, sizeof(double)); + nt->_sp13_rhs = (double*) ecalloc(neqn + 1, sizeof(double)); nt->_sp13mat = new Eigen::MatrixXd(neqn, neqn); - for (in = 0, i = 0; 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; @@ -2005,26 +2001,26 @@ static void nrn_matrix_node_alloc(void) { pnd = nt->_v_parent[in]; i = nd->eqn_index_; nt->_sp13_rhs[i] = nt->actual_rhs(in); - nd->diag_value = i; + nd->_d_matelm = &nt->_sp13mat->coeffRef(i - 1, i - 1); if (nde) { for (int ie = 0; ie < nlayer; ++ie) { - int k = i + ie; - nde->_d[ie] = &nt->_sp13mat->coeffRef(k, k); + int k = i + ie + 1; + nde->_d[ie] = &nt->_sp13mat->coeffRef(k - 1, k - 1); nde->_rhs[ie] = nt->_sp13_rhs + k; - nde->_x21[ie] = &nt->_sp13mat->coeffRef(k, k - 1); - nde->_x12[ie] = &nt->_sp13mat->coeffRef(k - 1, k); + nde->_x21[ie] = &nt->_sp13mat->coeffRef(k - 1, k - 2); + nde->_x12[ie] = &nt->_sp13mat->coeffRef(k - 2, k - 1); } } if (pnd) { j = pnd->eqn_index_; - nd->_a_matelm = &nt->_sp13mat->coeffRef(j, i); - nd->_b_matelm = &nt->_sp13mat->coeffRef(i, j); + nd->_a_matelm = &nt->_sp13mat->coeffRef(j - 1, i - 1); + nd->_b_matelm = &nt->_sp13mat->coeffRef(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->coeffRef(kp, k); - nde->_b_matelm[ie] = &nt->_sp13mat->coeffRef(k, kp); + int kp = j + ie + 1; + int k = i + ie + 1; + nde->_a_matelm[ie] = &nt->_sp13mat->coeffRef(kp - 1, k - 1); + nde->_b_matelm[ie] = &nt->_sp13mat->coeffRef(k - 1, kp - 1); } } else { /* not needed if index starts at 1 */ nd->_a_matelm = nullptr; From 3e64192ad5324fd02b7103f10a4d3d2c4f682b33 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 5 Jun 2024 14:45:42 +0200 Subject: [PATCH 14/26] Place holder & renaming --- src/nrncvode/occvode.cpp | 2 +- src/nrniv/matrixmap.cpp | 7 +++++- src/nrniv/nrndae.cpp | 6 ++--- src/nrnoc/fadvance.cpp | 6 ++--- src/nrnoc/multicore.cpp | 10 ++++---- src/nrnoc/multicore.h | 4 +-- src/nrnoc/solve.cpp | 4 +-- src/nrnoc/treeset.cpp | 54 ++++++++++++++++++++-------------------- 8 files changed, 49 insertions(+), 44 deletions(-) diff --git a/src/nrncvode/occvode.cpp b/src/nrncvode/occvode.cpp index 1b7b705cfe..15181d0d25 100644 --- a/src/nrncvode/occvode.cpp +++ b/src/nrncvode/occvode.cpp @@ -355,7 +355,7 @@ void Cvode::daspk_init_eqn() { if (use_sparse13 == 0 || diam_changed != 0) { recalc_diam(); } - zneq = _nt->_sp13mat->cols(); + zneq = _nt->_sparseMat->cols(); 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 e20420a3a8..8921c02b25 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_] = &_nt->_sp13mat->coeffRef(it - 1, jt - 1); + if (it == 0 || jt == 0) { + ptree_[plen_] = &_nt->_sparseMat->coeffRef(it - 1, jt - 1); + } else { + ptree_[plen_] = &place_holder; + } ++plen_; } } diff --git a/src/nrniv/nrndae.cpp b/src/nrniv/nrndae.cpp index 4bbe44d8f6..b323321bed 100644 --- a/src/nrniv/nrndae.cpp +++ b/src/nrniv/nrndae.cpp @@ -221,7 +221,7 @@ void NrnDAE::dkmap(std::vector>& pv, y_.data() + i}; pvdot[bmap_[i] - 1] = neuron::container::data_handle{neuron::container::do_not_search, - _nt->_sp13_rhs + bmap_[i]}; + _nt->_sparse_rhs + bmap_[i]}; } } @@ -231,7 +231,7 @@ void NrnDAE::update() { // note that the following is correct also for states that refer // to the internal potential of a segment. i.e rhs is v + vext[0] for (int i = 0; i < size_; ++i) { - y_[i] += _nt->_sp13_rhs[bmap_[i]]; + y_[i] += _nt->_sparse_rhs[bmap_[i]]; } // for (int i=0; i < size_; ++i) printf(" i=%d bmap_[i]=%d y_[i]=%g\n", i, bmap_[i], // y_->elem(i)); @@ -311,7 +311,7 @@ void NrnDAE::rhs() { v2y(); f_(y_, yptmp_, size_); for (int i = 0; i < size_; ++i) { - _nt->_sp13_rhs[bmap_[i]] += yptmp_[i]; + _nt->_sparse_rhs[bmap_[i]] += yptmp_[i]; } } diff --git a/src/nrnoc/fadvance.cpp b/src/nrnoc/fadvance.cpp index e064e1cd88..305f0ec748 100644 --- a/src/nrnoc/fadvance.cpp +++ b/src/nrnoc/fadvance.cpp @@ -672,10 +672,10 @@ void nrn_print_matrix(NrnThread* _nt) { Node* nd; if (use_sparse13) { if (ifarg(1) && chkarg(1, 0., 1.) == 0.) { - std::cout << &_nt->_sp13mat << std::endl; + std::cout << &_nt->_sparseMat << std::endl; } else { - std::cout << &_nt->_sp13mat << std::endl; - int n = _nt->_sp13mat->cols(); + std::cout << &_nt->_sparseMat << std::endl; + int n = _nt->_sparseMat->cols(); for (int i = 1; i <= n; ++i) { Printf("%d %g\n", i, _nt->actual_rhs(i)); } diff --git a/src/nrnoc/multicore.cpp b/src/nrnoc/multicore.cpp index 4aa34130f3..c4964a8c50 100644 --- a/src/nrnoc/multicore.cpp +++ b/src/nrnoc/multicore.cpp @@ -339,14 +339,14 @@ void nrn_threads_create(int n, bool parallel) { for (j = 0; j < BEFORE_AFTER_SIZE; ++j) { nt->tbl[j] = (NrnThreadBAList*) 0; } - nt->_sp13_rhs = nullptr; + nt->_sparse_rhs = nullptr; nt->_v_parent_index = 0; nt->_v_node = 0; nt->_v_parent = 0; nt->_ecell_memb_list = 0; nt->_ecell_child_cnt = 0; nt->_ecell_children = NULL; - nt->_sp13mat = nullptr; + nt->_sparseMat = nullptr; nt->_ctime = 0.0; nt->_vcv = 0; nt->_node_data_offset = 0; @@ -440,9 +440,9 @@ void nrn_threads_free() { free(nt->_ecell_children); nt->_ecell_children = NULL; } - if (nt->_sp13mat) { - delete nt->_sp13mat; - nt->_sp13mat = nullptr; + if (nt->_sparseMat) { + delete nt->_sparseMat; + nt->_sparseMat = nullptr; } nt->end = 0; nt->ncell = 0; diff --git a/src/nrnoc/multicore.h b/src/nrnoc/multicore.h index 6e12c51a65..5ec151fba4 100644 --- a/src/nrnoc/multicore.h +++ b/src/nrnoc/multicore.h @@ -95,9 +95,9 @@ struct NrnThread { int* _v_parent_index; Node** _v_node; Node** _v_parent; - double* _sp13_rhs; /* rhs matrix for sparse13 solver. updates to and from this vector + double* _sparse_rhs; /* rhs matrix for sparse13 solver. updates to and from this vector need to be transfered to actual_rhs */ - Eigen::MatrixXd* _sp13mat{}; /* handle to general sparse matrix */ + Eigen::MatrixXd* _sparseMat{}; /* 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 508f7763c6..b435e4c90d 100644 --- a/src/nrnoc/solve.cpp +++ b/src/nrnoc/solve.cpp @@ -368,8 +368,8 @@ void nrn_solve(NrnThread* _nt) { nrn_thread_error("solve use_sparse13"); update_sp13_mat_based_on_actual_d(_nt); update_sp13_rhs_based_on_actual_rhs(_nt); - Eigen::FullPivLU lu{*_nt->_sp13mat}; - Eigen::Map rhs(_nt->_sp13_rhs + 1, _nt->_sp13mat->cols()); + Eigen::FullPivLU lu{*_nt->_sparseMat}; + Eigen::Map rhs(_nt->_sparse_rhs + 1, _nt->_sparseMat->cols()); rhs = lu.solve(rhs); update_actual_d_based_on_sp13_mat(_nt); update_actual_rhs_based_on_sp13_rhs(_nt); diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index 83d3384591..27e9998d6b 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -329,20 +329,20 @@ vm += dvi-dvx */ /* - * Update actual_rhs based on _sp13_rhs used for sparse13 solver + * Update actual_rhs based on _sparse_rhs used for sparse13 solver */ void update_actual_rhs_based_on_sp13_rhs(NrnThread* nt) { for (int i = 0; i < nt->end; i++) { - nt->actual_rhs(i) = nt->_sp13_rhs[nt->_v_node[i]->eqn_index_]; + nt->actual_rhs(i) = nt->_sparse_rhs[nt->_v_node[i]->eqn_index_]; } } /* - * Update _sp13_rhs used for sparse13 solver based on changes on actual_rhs + * Update _sparse_rhs used for sparse13 solver based on changes on actual_rhs */ void update_sp13_rhs_based_on_actual_rhs(NrnThread* nt) { for (int i = 0; i < nt->end; i++) { - nt->_sp13_rhs[nt->_v_node[i]->eqn_index_] = nt->actual_rhs(i); + nt->_sparse_rhs[nt->_v_node[i]->eqn_index_] = nt->actual_rhs(i); } } @@ -391,9 +391,9 @@ void nrn_rhs(neuron::model_sorted_token const& cache_token, NrnThread& nt) { auto* const vec_rhs = nt.node_rhs_storage(); if (use_sparse13) { nrn_thread_error("nrn_rhs use_sparse13"); - int neqn = _nt->_sp13mat->cols(); + int neqn = _nt->_sparseMat->cols(); for (int i = 1; i <= neqn; ++i) { - _nt->_sp13_rhs[i] = 0.; + _nt->_sparse_rhs[i] = 0.; } for (int i = i1; i < i3; ++i) { NODERHS(_nt->_v_node[i]) = 0.; @@ -496,7 +496,7 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) { } if (use_sparse13) { - _nt->_sp13mat->setZero(); + _nt->_sparseMat->setZero(); } // Make sure the SoA node diagonals are also zeroed (is this needed?) @@ -1882,12 +1882,12 @@ void node_data(void) { void nrn_matrix_node_free() { NrnThread* nt; FOR_THREADS(nt) { - if (nt->_sp13_rhs) { - free(std::exchange(nt->_sp13_rhs, nullptr)); + if (nt->_sparse_rhs) { + free(std::exchange(nt->_sparse_rhs, nullptr)); } - if (nt->_sp13mat) { - delete nt->_sp13mat; - nt->_sp13mat = nullptr; + if (nt->_sparseMat) { + delete nt->_sparseMat; + nt->_sparseMat = nullptr; } } diam_changed = 1; @@ -1957,16 +1957,16 @@ static void nrn_matrix_node_alloc(void) { nrn_method_consistent(); nt = nrn_threads; - /*printf("use_sparse13=%d sp13mat=%lx rhs=%lx\n", use_sparse13, (long)nt->_sp13mat, + /*printf("use_sparse13=%d sp13mat=%lx rhs=%lx\n", use_sparse13, (long)nt->_sparseMat, * (long)nt->_actual_rhs);*/ if (use_sparse13) { - if (nt->_sp13mat) { + if (nt->_sparseMat) { return; } else { nrn_matrix_node_free(); } } else { - if (nt->_sp13mat) { + if (nt->_sparseMat) { v_structure_change = 1; v_setup_vectors(); return; @@ -1985,8 +1985,8 @@ 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 Eigen::MatrixXd(neqn, neqn); + nt->_sparse_rhs = (double*) ecalloc(neqn + 1, sizeof(double)); + nt->_sparseMat = new Eigen::MatrixXd(neqn, neqn); for (in = 0, i = 1; in < nt->end; ++in, ++i) { nt->_v_node[in]->eqn_index_ = i; if (nt->_v_node[in]->extnode) { @@ -2000,27 +2000,27 @@ static void nrn_matrix_node_alloc(void) { nde = nd->extnode; pnd = nt->_v_parent[in]; i = nd->eqn_index_; - nt->_sp13_rhs[i] = nt->actual_rhs(in); - nd->_d_matelm = &nt->_sp13mat->coeffRef(i - 1, i - 1); + nt->_sparse_rhs[i] = nt->actual_rhs(in); + nd->_d_matelm = &nt->_sparseMat->coeffRef(i - 1, i - 1); if (nde) { for (int ie = 0; ie < nlayer; ++ie) { int k = i + ie + 1; - nde->_d[ie] = &nt->_sp13mat->coeffRef(k - 1, k - 1); - nde->_rhs[ie] = nt->_sp13_rhs + k; - nde->_x21[ie] = &nt->_sp13mat->coeffRef(k - 1, k - 2); - nde->_x12[ie] = &nt->_sp13mat->coeffRef(k - 2, k - 1); + nde->_d[ie] = &nt->_sparseMat->coeffRef(k - 1, k - 1); + nde->_rhs[ie] = nt->_sparse_rhs + k; + nde->_x21[ie] = &nt->_sparseMat->coeffRef(k - 1, k - 2); + nde->_x12[ie] = &nt->_sparseMat->coeffRef(k - 2, k - 1); } } if (pnd) { j = pnd->eqn_index_; - nd->_a_matelm = &nt->_sp13mat->coeffRef(j - 1, i - 1); - nd->_b_matelm = &nt->_sp13mat->coeffRef(i - 1, j - 1); + nd->_a_matelm = &nt->_sparseMat->coeffRef(j - 1, i - 1); + nd->_b_matelm = &nt->_sparseMat->coeffRef(i - 1, j - 1); if (nde && pnd->extnode) for (int ie = 0; ie < nlayer; ++ie) { int kp = j + ie + 1; int k = i + ie + 1; - nde->_a_matelm[ie] = &nt->_sp13mat->coeffRef(kp - 1, k - 1); - nde->_b_matelm[ie] = &nt->_sp13mat->coeffRef(k - 1, kp - 1); + nde->_a_matelm[ie] = &nt->_sparseMat->coeffRef(kp - 1, k - 1); + nde->_b_matelm[ie] = &nt->_sparseMat->coeffRef(k - 1, kp - 1); } } else { /* not needed if index starts at 1 */ nd->_a_matelm = nullptr; From e08cef16ecd8b87412d581b064682f5cb0724a52 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 5 Jun 2024 18:45:49 +0200 Subject: [PATCH 15/26] Let's go for green! --- src/nrncvode/occvode.cpp | 3 ++- src/nrniv/matrixmap.cpp | 4 ++-- src/nrnoc/multicore.h | 9 ++------- src/nrnoc/solve.cpp | 4 ++-- src/nrnoc/treeset.cpp | 37 +++++++++++++++++++++++++++++++++++-- 5 files changed, 43 insertions(+), 14 deletions(-) diff --git a/src/nrncvode/occvode.cpp b/src/nrncvode/occvode.cpp index 15181d0d25..608784e604 100644 --- a/src/nrncvode/occvode.cpp +++ b/src/nrncvode/occvode.cpp @@ -1,3 +1,5 @@ +#include + #include <../../nrnconf.h> #include #include "nrn_ansi.h" @@ -14,7 +16,6 @@ #include -#include extern void setup_topology(), v_setup_vectors(); extern void recalc_diam(); diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp index 8921c02b25..050685ce18 100644 --- a/src/nrniv/matrixmap.cpp +++ b/src/nrniv/matrixmap.cpp @@ -56,9 +56,9 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { jt = start + j - nnode; } if (it == 0 || jt == 0) { - ptree_[plen_] = &_nt->_sparseMat->coeffRef(it - 1, jt - 1); - } else { ptree_[plen_] = &place_holder; + } else { + ptree_[plen_] = &_nt->_sparseMat->coeffRef(it - 1, jt - 1); } ++plen_; } diff --git a/src/nrnoc/multicore.h b/src/nrnoc/multicore.h index 5ec151fba4..bdbe73c6ed 100644 --- a/src/nrnoc/multicore.h +++ b/src/nrnoc/multicore.h @@ -27,15 +27,10 @@ actual_v, etc. */ #include "membfunc.h" +#include #include -namespace Eigen { -template -class Matrix; -using MatrixXd = Matrix; -} // namespace Eigen - typedef struct NrnThreadMembList { /* patterned after CvMembList in cvodeobj.h */ struct NrnThreadMembList* next; Memb_list* ml; @@ -97,7 +92,7 @@ struct NrnThread { Node** _v_parent; double* _sparse_rhs; /* rhs matrix for sparse13 solver. updates to and from this vector need to be transfered to actual_rhs */ - Eigen::MatrixXd* _sparseMat{}; /* handle to general sparse matrix */ + Eigen::SparseMatrix* _sparseMat{}; /* 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 b435e4c90d..10fb8b2506 100644 --- a/src/nrnoc/solve.cpp +++ b/src/nrnoc/solve.cpp @@ -58,6 +58,7 @@ node.v + extnode.v[0] #include "ocnotify.h" #include "section.h" #include +#include #include "treeset.h" #include @@ -364,11 +365,10 @@ void nrn_solve(NrnThread* _nt) { } #else if (use_sparse13) { - int e; nrn_thread_error("solve use_sparse13"); update_sp13_mat_based_on_actual_d(_nt); update_sp13_rhs_based_on_actual_rhs(_nt); - Eigen::FullPivLU lu{*_nt->_sparseMat}; + Eigen::SparseLU> lu(*_nt->_sparseMat); Eigen::Map rhs(_nt->_sparse_rhs + 1, _nt->_sparseMat->cols()); rhs = lu.solve(rhs); update_actual_d_based_on_sp13_mat(_nt); diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index 27e9998d6b..15fd57a54c 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -496,7 +496,12 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) { } if (use_sparse13) { - _nt->_sparseMat->setZero(); + Eigen::SparseMatrix& m_ = *_nt->_sparseMat; + for (int k = 0; k < m_.outerSize(); ++k) { + for (Eigen::SparseMatrix::InnerIterator it(m_, k); it; ++it) { + it.valueRef() = 0.; + } + } } // Make sure the SoA node diagonals are also zeroed (is this needed?) @@ -1986,13 +1991,41 @@ static void nrn_matrix_node_alloc(void) { /*printf(" %d extracellular nodes\n", extn);*/ neqn += extn; nt->_sparse_rhs = (double*) ecalloc(neqn + 1, sizeof(double)); - nt->_sparseMat = new Eigen::MatrixXd(neqn, neqn); + nt->_sparseMat = new Eigen::SparseMatrix(neqn, neqn); 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 (in = 0; in < nt->end; ++in) { + const Node *nd = nt->_v_node[in]; + const Extnode* nde = nd->extnode; + const Node *pnd = nt->_v_parent[in]; + int i = nd->eqn_index_; + nt->_sparseMat->coeffRef(i - 1, i - 1) = 0.; + if (nde) { + for (int ie = 0; ie < nlayer; ++ie) { + int k = i + ie + 1; + nt->_sparseMat->coeffRef(k - 1, k - 1) = 0.; + nt->_sparseMat->coeffRef(k - 1, k - 2) = 0; + nt->_sparseMat->coeffRef(k - 2, k - 1) = 0; + } + } + if (pnd) { + int j = pnd->eqn_index_; + nt->_sparseMat->coeffRef(j - 1, i - 1); + nt->_sparseMat->coeffRef(i - 1, j - 1); + if (nde && pnd->extnode) + for (int ie = 0; ie < nlayer; ++ie) { + int kp = j + ie + 1; + int k = i + ie + 1; + nt->_sparseMat->coeffRef(kp - 1, k - 1) = 0; + nt->_sparseMat->coeffRef(k - 1, kp - 1) = 0; + } + } + } + // nt->_sparseMat->makeCompressed(); for (in = 0; in < nt->end; ++in) { Node *nd, *pnd; Extnode* nde; From f17b925996e7d6e2e861db076ce933edd7edd410 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Wed, 5 Jun 2024 17:25:18 +0000 Subject: [PATCH 16/26] Fix formatting --- src/nrnoc/multicore.h | 13 +++++++------ src/nrnoc/treeset.cpp | 4 ++-- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/nrnoc/multicore.h b/src/nrnoc/multicore.h index bdbe73c6ed..5934ba925a 100644 --- a/src/nrnoc/multicore.h +++ b/src/nrnoc/multicore.h @@ -90,12 +90,13 @@ struct NrnThread { int* _v_parent_index; Node** _v_node; Node** _v_parent; - double* _sparse_rhs; /* rhs matrix for sparse13 solver. updates to and from this vector - need to be transfered to actual_rhs */ - Eigen::SparseMatrix* _sparseMat{}; /* 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_ */ + double* _sparse_rhs; /* rhs matrix for sparse13 solver. updates to and from this vector + need to be transfered to actual_rhs */ + Eigen::SparseMatrix* _sparseMat{}; /* 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_ */ #if 1 double _ctime; /* computation time in seconds (using nrnmpi_wtime) */ diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index 15fd57a54c..d0fc16ca61 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -1999,9 +1999,9 @@ static void nrn_matrix_node_alloc(void) { } } for (in = 0; in < nt->end; ++in) { - const Node *nd = nt->_v_node[in]; + const Node* nd = nt->_v_node[in]; const Extnode* nde = nd->extnode; - const Node *pnd = nt->_v_parent[in]; + const Node* pnd = nt->_v_parent[in]; int i = nd->eqn_index_; nt->_sparseMat->coeffRef(i - 1, i - 1) = 0.; if (nde) { From cbd699507393110608b53d4303935a10682c88e7 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 6 Jun 2024 01:01:39 +0200 Subject: [PATCH 17/26] Avoid including eigen there --- src/nrnoc/multicore.h | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/nrnoc/multicore.h b/src/nrnoc/multicore.h index bdbe73c6ed..a87f72929f 100644 --- a/src/nrnoc/multicore.h +++ b/src/nrnoc/multicore.h @@ -27,7 +27,10 @@ actual_v, etc. */ #include "membfunc.h" -#include +namespace Eigen { + template + class SparseMatrix; +} #include @@ -92,7 +95,7 @@ struct NrnThread { Node** _v_parent; double* _sparse_rhs; /* rhs matrix for sparse13 solver. updates to and from this vector need to be transfered to actual_rhs */ - Eigen::SparseMatrix* _sparseMat{}; /* handle to general sparse matrix */ + Eigen::SparseMatrix* _sparseMat{}; /* 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_ */ From 3706ba01f1204669ec4459be262808564d7bb1d6 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 6 Jun 2024 01:04:01 +0200 Subject: [PATCH 18/26] format --- src/nrnoc/multicore.h | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/nrnoc/multicore.h b/src/nrnoc/multicore.h index a87f72929f..8272dfca9e 100644 --- a/src/nrnoc/multicore.h +++ b/src/nrnoc/multicore.h @@ -28,8 +28,8 @@ actual_v, etc. #include "membfunc.h" namespace Eigen { - template - class SparseMatrix; +template +class SparseMatrix; } #include @@ -93,12 +93,12 @@ struct NrnThread { int* _v_parent_index; Node** _v_node; Node** _v_parent; - double* _sparse_rhs; /* rhs matrix for sparse13 solver. updates to and from this vector - need to be transfered to actual_rhs */ + double* _sparse_rhs; /* rhs matrix for sparse13 solver. updates to and from this vector + need to be transfered to actual_rhs */ Eigen::SparseMatrix* _sparseMat{}; /* 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_ */ + 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_ */ #if 1 double _ctime; /* computation time in seconds (using nrnmpi_wtime) */ From 61aca57fe2aa197737fd46da4a4b515f551b3d47 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 6 Jun 2024 10:04:20 +0200 Subject: [PATCH 19/26] Prepare the matrix for solving with SparseLU --- src/nrnoc/solve.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/nrnoc/solve.cpp b/src/nrnoc/solve.cpp index 10fb8b2506..0a8029a531 100644 --- a/src/nrnoc/solve.cpp +++ b/src/nrnoc/solve.cpp @@ -368,6 +368,7 @@ void nrn_solve(NrnThread* _nt) { nrn_thread_error("solve use_sparse13"); update_sp13_mat_based_on_actual_d(_nt); update_sp13_rhs_based_on_actual_rhs(_nt); + _nt->_sparseMat->makeCompressed(); Eigen::SparseLU> lu(*_nt->_sparseMat); Eigen::Map rhs(_nt->_sparse_rhs + 1, _nt->_sparseMat->cols()); rhs = lu.solve(rhs); From 846364d63f1c9d0d0e545b645523079582ef26ed Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 6 Jun 2024 11:18:05 +0200 Subject: [PATCH 20/26] Compress as soon as possible --- src/nrnoc/treeset.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index d0fc16ca61..1428b411fc 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -2025,7 +2025,7 @@ static void nrn_matrix_node_alloc(void) { } } } - // nt->_sparseMat->makeCompressed(); + nt->_sparseMat->makeCompressed(); for (in = 0; in < nt->end; ++in) { Node *nd, *pnd; Extnode* nde; From f4bc02d2774eb531d7561bada4289241e9556a9a Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 6 Jun 2024 12:56:29 +0200 Subject: [PATCH 21/26] Use triplets to set the matrix --- src/nrniv/nrndae.cpp | 2 +- src/nrnoc/treeset.cpp | 119 ++++++++++++++++++++++++------------------ 2 files changed, 68 insertions(+), 53 deletions(-) diff --git a/src/nrniv/nrndae.cpp b/src/nrniv/nrndae.cpp index b323321bed..2ab296ffab 100644 --- a/src/nrniv/nrndae.cpp +++ b/src/nrniv/nrndae.cpp @@ -8,7 +8,7 @@ extern int secondorder; -static NrnDAEPtrList nrndae_list; +static std::list nrndae_list; int nrndae_list_is_empty() { return nrndae_list.empty() ? 1 : 0; diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index 1428b411fc..ef58e3130c 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -1998,69 +1998,84 @@ static void nrn_matrix_node_alloc(void) { i += nlayer; } } - for (in = 0; in < nt->end; ++in) { - const Node* nd = nt->_v_node[in]; - const Extnode* nde = nd->extnode; - const Node* pnd = nt->_v_parent[in]; - int i = nd->eqn_index_; - nt->_sparseMat->coeffRef(i - 1, i - 1) = 0.; - if (nde) { - for (int ie = 0; ie < nlayer; ++ie) { - int k = i + ie + 1; - nt->_sparseMat->coeffRef(k - 1, k - 1) = 0.; - nt->_sparseMat->coeffRef(k - 1, k - 2) = 0; - nt->_sparseMat->coeffRef(k - 2, k - 1) = 0; - } - } - if (pnd) { - int j = pnd->eqn_index_; - nt->_sparseMat->coeffRef(j - 1, i - 1); - nt->_sparseMat->coeffRef(i - 1, j - 1); - if (nde && pnd->extnode) + { + // sparse13 was a linked list of single values and so we were able to have + // stable pointers on internal values. + // Eigen::SparseMatrix is more like a vector, that is reallocated when needed + // and so pointers will be invalidate. + // More over, when solving with SparseLU, the matrix is compressed and the + // internal state is fully changed. + // setFromTriplets is an efficient way of having an already compressed matrix. + // We can safely get stable pointers as long as we don't insert new value with + // coeff() and coeffRef(). + std::vector> triplets{}; + for (in = 0; in < nt->end; ++in) { + const Node* nd = nt->_v_node[in]; + const Extnode* nde = nd->extnode; + const Node* pnd = nt->_v_parent[in]; + int i = nd->eqn_index_; + triplets.emplace_back(i - 1, i - 1, 0.); + if (nde) { for (int ie = 0; ie < nlayer; ++ie) { - int kp = j + ie + 1; int k = i + ie + 1; - nt->_sparseMat->coeffRef(kp - 1, k - 1) = 0; - nt->_sparseMat->coeffRef(k - 1, kp - 1) = 0; + triplets.emplace_back(k - 1, k - 1, 0.); + triplets.emplace_back(k - 1, k - 2, 0.); + triplets.emplace_back(k - 2, k - 1, 0.); } + } + if (pnd) { + int j = pnd->eqn_index_; + triplets.emplace_back(j - 1, i - 1, 0.); + triplets.emplace_back(i - 1, j - 1, 0.); + if (nde && pnd->extnode) + for (int ie = 0; ie < nlayer; ++ie) { + triplets.emplace_back(j + ie, i + ie, 0.); + triplets.emplace_back(i + ie, j + ie, 0.); + } + } } + nt->_sparseMat->setFromTriplets(triplets.begin(), triplets.end()); } + // Add some new nodes + nrndae_alloc(); nt->_sparseMat->makeCompressed(); - for (in = 0; in < nt->end; ++in) { - Node *nd, *pnd; - Extnode* nde; - nd = nt->_v_node[in]; - nde = nd->extnode; - pnd = nt->_v_parent[in]; - i = nd->eqn_index_; - nt->_sparse_rhs[i] = nt->actual_rhs(in); - nd->_d_matelm = &nt->_sparseMat->coeffRef(i - 1, i - 1); - if (nde) { - for (int ie = 0; ie < nlayer; ++ie) { - int k = i + ie + 1; - nde->_d[ie] = &nt->_sparseMat->coeffRef(k - 1, k - 1); - nde->_rhs[ie] = nt->_sparse_rhs + k; - nde->_x21[ie] = &nt->_sparseMat->coeffRef(k - 1, k - 2); - nde->_x12[ie] = &nt->_sparseMat->coeffRef(k - 2, k - 1); - } - } - if (pnd) { - j = pnd->eqn_index_; - nd->_a_matelm = &nt->_sparseMat->coeffRef(j - 1, i - 1); - nd->_b_matelm = &nt->_sparseMat->coeffRef(i - 1, j - 1); - if (nde && pnd->extnode) + { + // Now that we have our stable matrix, let's extract pointers. + for (in = 0; in < nt->end; ++in) { + Node *nd, *pnd; + Extnode* nde; + nd = nt->_v_node[in]; + nde = nd->extnode; + pnd = nt->_v_parent[in]; + i = nd->eqn_index_; + nt->_sparse_rhs[i] = nt->actual_rhs(in); + nd->_d_matelm = &nt->_sparseMat->coeffRef(i - 1, i - 1); + if (nde) { for (int ie = 0; ie < nlayer; ++ie) { - int kp = j + ie + 1; int k = i + ie + 1; - nde->_a_matelm[ie] = &nt->_sparseMat->coeffRef(kp - 1, k - 1); - nde->_b_matelm[ie] = &nt->_sparseMat->coeffRef(k - 1, kp - 1); + nde->_d[ie] = &nt->_sparseMat->coeffRef(k - 1, k - 1); + nde->_rhs[ie] = nt->_sparse_rhs + k; + nde->_x21[ie] = &nt->_sparseMat->coeffRef(k - 1, k - 2); + nde->_x12[ie] = &nt->_sparseMat->coeffRef(k - 2, k - 1); } - } else { /* not needed if index starts at 1 */ - nd->_a_matelm = nullptr; - nd->_b_matelm = nullptr; + } + if (pnd) { + j = pnd->eqn_index_; + nd->_a_matelm = &nt->_sparseMat->coeffRef(j - 1, i - 1); + nd->_b_matelm = &nt->_sparseMat->coeffRef(i - 1, j - 1); + if (nde && pnd->extnode) + for (int ie = 0; ie < nlayer; ++ie) { + int kp = j + ie + 1; + int k = i + ie + 1; + nde->_a_matelm[ie] = &nt->_sparseMat->coeffRef(kp - 1, k - 1); + nde->_b_matelm[ie] = &nt->_sparseMat->coeffRef(k - 1, kp - 1); + } + } else { /* not needed if index starts at 1 */ + nd->_a_matelm = nullptr; + nd->_b_matelm = nullptr; + } } } - nrndae_alloc(); } else { FOR_THREADS(nt) { assert(nrndae_extra_eqn_count() == 0); From 40d6844f4afb6c3ac18642d772b95de9c7da0c92 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 6 Jun 2024 15:58:38 +0200 Subject: [PATCH 22/26] Be expressive --- src/nrnoc/treeset.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index ef58e3130c..e7219ea9b6 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -2029,8 +2029,10 @@ static void nrn_matrix_node_alloc(void) { triplets.emplace_back(i - 1, j - 1, 0.); if (nde && pnd->extnode) for (int ie = 0; ie < nlayer; ++ie) { - triplets.emplace_back(j + ie, i + ie, 0.); - triplets.emplace_back(i + ie, j + ie, 0.); + int kp = j + ie + 1; + int k = i + ie + 1; + triplets.emplace_back(kp - 1, k - 1, 0.); + triplets.emplace_back(k - 1, kp - 1, 0.); } } } From bd7248da27240fb08ccece541ba8bef00e5c3223 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 6 Jun 2024 15:58:21 +0200 Subject: [PATCH 23/26] Any fixed typo is good --- share/lib/nrn.defaults.in | 12 ++++++------ src/ivoc/idraw.cpp | 2 +- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/share/lib/nrn.defaults.in b/share/lib/nrn.defaults.in index 332c3e0a83..22049eb11a 100755 --- a/share/lib/nrn.defaults.in +++ b/share/lib/nrn.defaults.in @@ -78,13 +78,13 @@ *xvalue_field_size_increase: 10. // following is historical standard for unix. -@nrndef_unix@*font: *helvetica-medium-r-normal*--14* -@nrndef_unix@*MenuBar*font: *-helvetica-bold-o-normal--14* -@nrndef_unix@*MenuItem*font: *-helvetica-bold-o-normal--14* +@nrndef_unix@*font: *fixed-medium-r-normal*--14* +@nrndef_unix@*MenuBar*font: *-fixed-bold-o-normal--14* +@nrndef_unix@*MenuItem*font: *-fixed-bold-o-normal--14* // Mac and mswin internal default is -@nrndef_mac@*font: *helvetica-medium-r-normal*--14* -@nrndef_mac@*MenuBar*font: *-helvetica-bold-o-normal--14* -@nrndef_mac@*MenuItem*font: *-helvetica-bold-o-normal--14* +@nrndef_mac@*font: *fixed-medium-r-normal*--14* +@nrndef_mac@*MenuBar*font: *-fixed-bold-o-normal--14* +@nrndef_mac@*MenuItem*font: *-fixed-bold-o-normal--14* @nrndef_mswin@*font: *Arial*bold*--12* //*MenuBar*font: *Arial*bold*--12* //*MenuItem*font: *Arial*bold*--12* diff --git a/src/ivoc/idraw.cpp b/src/ivoc/idraw.cpp index 1216b226d1..f2798032eb 100644 --- a/src/ivoc/idraw.cpp +++ b/src/ivoc/idraw.cpp @@ -152,7 +152,7 @@ void OcIdraw::text(Canvas*, out << font->name() << font->size() << "SetF\n"; } else { out << "\ -%I f -*-helvetica-medium-r-normal-*-12-*-*-*-*-*-*-*\n\ +%I f -*-fixed-medium-r-normal-*-12-*-*-*-*-*-*-*\n\ Helvetica 12 SetF\n\ "; }; From 299d2edf6817a8c1c999eaa3637540064504806a Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 6 Jun 2024 18:33:17 +0200 Subject: [PATCH 24/26] Fix plen --- src/nrniv/matrixmap.cpp | 22 +++++++++++++--------- src/nrniv/matrixmap.h | 2 +- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp index 050685ce18..98a12eab79 100644 --- a/src/nrniv/matrixmap.cpp +++ b/src/nrniv/matrixmap.cpp @@ -19,20 +19,23 @@ void MatrixMap::mmfree() { void MatrixMap::add(double fac) { for (int i = 0; i < plen_; ++i) { - *ptree_[i] += fac * (*pm_[i]); + auto [j, k] = ptree_[i]; + if (j != -1) { + // std::cerr << "modifying '" << nrn_threads->_sparseMat->coeff(j -1, k - 1) << "' (" << j - 1 << ", " << k - 1 << ") with " << fac * (*pm_[i]) << std::endl; + nrn_threads->_sparseMat->coeffRef(j - 1, k - 1) += fac * (*pm_[i]); + } } } 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()); + plen_ = nonzero_i.size(); for (int k = 0; k < nonzero_i.size(); k++) { const int i = nonzero_i[k]; const int j = nonzero_j[k]; @@ -46,7 +49,7 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { it = start + i - nnode; } int jt; - pm_[plen_] = m_.mep(i, j); + pm_.push_back(m_.mep(i, j)); if (j < nnode) { jt = nodes[j]->eqn_index_ + layer[j]; if (layer[j] > 0 && !nodes[j]->extnode) { @@ -56,10 +59,11 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { jt = start + j - nnode; } if (it == 0 || jt == 0) { - ptree_[plen_] = &place_holder; + ptree_.emplace_back(-1, -1); } else { - ptree_[plen_] = &_nt->_sparseMat->coeffRef(it - 1, jt - 1); + std::cerr << "someone touch me in the heart!: (" << it << ", " << jt << ")" << std::endl; + ptree_.emplace_back(it, jt); + _nt->_sparseMat->coeffRef(it - 1, jt - 1) = 0.; } - ++plen_; } } diff --git a/src/nrniv/matrixmap.h b/src/nrniv/matrixmap.h index f216dcc9c3..0a21e5ac0b 100644 --- a/src/nrniv/matrixmap.h +++ b/src/nrniv/matrixmap.h @@ -117,5 +117,5 @@ class MatrixMap { // the map int plen_ = 0; std::vector pm_{}; - std::vector ptree_{}; + std::vector> ptree_{}; }; From e1df5f231ac6511e1671c8d6278c5055e815d1a8 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Thu, 6 Jun 2024 16:35:03 +0000 Subject: [PATCH 25/26] Fix formatting --- src/nrniv/matrixmap.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp index 98a12eab79..96478b57db 100644 --- a/src/nrniv/matrixmap.cpp +++ b/src/nrniv/matrixmap.cpp @@ -21,7 +21,8 @@ void MatrixMap::add(double fac) { for (int i = 0; i < plen_; ++i) { auto [j, k] = ptree_[i]; if (j != -1) { - // std::cerr << "modifying '" << nrn_threads->_sparseMat->coeff(j -1, k - 1) << "' (" << j - 1 << ", " << k - 1 << ") with " << fac * (*pm_[i]) << std::endl; + // std::cerr << "modifying '" << nrn_threads->_sparseMat->coeff(j -1, k - 1) << "' (" << + // j - 1 << ", " << k - 1 << ") with " << fac * (*pm_[i]) << std::endl; nrn_threads->_sparseMat->coeffRef(j - 1, k - 1) += fac * (*pm_[i]); } } @@ -61,7 +62,8 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { if (it == 0 || jt == 0) { ptree_.emplace_back(-1, -1); } else { - std::cerr << "someone touch me in the heart!: (" << it << ", " << jt << ")" << std::endl; + std::cerr << "someone touch me in the heart!: (" << it << ", " << jt << ")" + << std::endl; ptree_.emplace_back(it, jt); _nt->_sparseMat->coeffRef(it - 1, jt - 1) = 0.; } From 6b411b67e7fa907785f6f289b8703333205bc817 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 22 Oct 2024 11:18:10 +0200 Subject: [PATCH 26/26] Remove problem during merge --- share/lib/nrn.defaults.in | 12 ++++++------ src/ivoc/idraw.cpp | 2 +- src/nrnoc/treeset.cpp | 2 -- 3 files changed, 7 insertions(+), 9 deletions(-) diff --git a/share/lib/nrn.defaults.in b/share/lib/nrn.defaults.in index 22049eb11a..332c3e0a83 100755 --- a/share/lib/nrn.defaults.in +++ b/share/lib/nrn.defaults.in @@ -78,13 +78,13 @@ *xvalue_field_size_increase: 10. // following is historical standard for unix. -@nrndef_unix@*font: *fixed-medium-r-normal*--14* -@nrndef_unix@*MenuBar*font: *-fixed-bold-o-normal--14* -@nrndef_unix@*MenuItem*font: *-fixed-bold-o-normal--14* +@nrndef_unix@*font: *helvetica-medium-r-normal*--14* +@nrndef_unix@*MenuBar*font: *-helvetica-bold-o-normal--14* +@nrndef_unix@*MenuItem*font: *-helvetica-bold-o-normal--14* // Mac and mswin internal default is -@nrndef_mac@*font: *fixed-medium-r-normal*--14* -@nrndef_mac@*MenuBar*font: *-fixed-bold-o-normal--14* -@nrndef_mac@*MenuItem*font: *-fixed-bold-o-normal--14* +@nrndef_mac@*font: *helvetica-medium-r-normal*--14* +@nrndef_mac@*MenuBar*font: *-helvetica-bold-o-normal--14* +@nrndef_mac@*MenuItem*font: *-helvetica-bold-o-normal--14* @nrndef_mswin@*font: *Arial*bold*--12* //*MenuBar*font: *Arial*bold*--12* //*MenuItem*font: *Arial*bold*--12* diff --git a/src/ivoc/idraw.cpp b/src/ivoc/idraw.cpp index f2798032eb..1216b226d1 100644 --- a/src/ivoc/idraw.cpp +++ b/src/ivoc/idraw.cpp @@ -152,7 +152,7 @@ void OcIdraw::text(Canvas*, out << font->name() << font->size() << "SetF\n"; } else { out << "\ -%I f -*-fixed-medium-r-normal-*-12-*-*-*-*-*-*-*\n\ +%I f -*-helvetica-medium-r-normal-*-12-*-*-*-*-*-*-*\n\ Helvetica 12 SetF\n\ "; }; diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index be25013db2..f3637addb0 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -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_;