treeset: replace sparse13 by eigen

src/nocmodl/noccout.cpp

@@ -260,7 +260,8 @@ 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 ind = _nrn_mechanism_access_index(_nd);\n");
+            P("   *_nrn_mechanism_get_matrix_elem(_nt, ind, ind) += _g;\n");
             P(" }\n");
             P("#endif\n");
         } else {
@@ -679,7 +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("   *_extnode->_d[0] += _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 {

src/nrncvode/nrndaspk.cpp

@@ -6,7 +6,6 @@
 
 #include <stdio.h>
 #include <math.h>
-#include "spmatrix.h"
 #include "nrnoc2iv.h"
 #include "cvodeobj.h"
 #include "nrndaspk.h"

src/nrncvode/occvode.cpp

@@ -18,7 +18,7 @@
 #include <numeric>
 
 
-#include "spmatrix.h"
+#include "ocmatrix.h"
 extern double* sp13mat;
 
 #if 1 || NRNMPI
@@ -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) {

src/nrniv/matrixmap.cpp

@@ -2,8 +2,6 @@
 #include "matrixmap.h"
 #include <vector>
 
-#include "spmatrix.h"
-
 MatrixMap::MatrixMap(Matrix& mat)
     : m_(mat) {}
 
@@ -18,41 +16,40 @@ void MatrixMap::mmfree() {
 }
 
 void MatrixMap::add(double fac) {
-    for (int i = 0; i < plen_; ++i) {
-        auto [it, jt] = pm_[i];
-        *ptree_[i] += fac * m_(it, jt);
+    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);
+    }
+}
+
+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) {
-    NrnThread* _nt = nrn_threads;
     mmfree();
 
-    plen_ = 0;
     std::vector<std::pair<int, int>> nzs = m_.nonzeros();
-    pm_.resize(nzs.size());
-    ptree_.resize(nzs.size());
-    for (const auto [i, j]: nzs) {
-        int it;
-        if (i < nnode) {
-            it = nodes[i]->eqn_index_ + layer[i];
-            if (layer[i] > 0 && !nodes[i]->extnode) {
-                it = 0;
-            }
-        } else {
-            it = start + i - nnode;
-        }
-        int jt;
-        pm_[plen_] = std::make_pair(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;
+    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;
         }
-        ptree_[plen_] = spGetElement(_nt->_sp13mat, it, jt);
-        ++plen_;
+        pm_.emplace_back(i, j);
+        ptree_.emplace_back(it - 1, jt - 1);
     }
 }

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<std::pair<int, int>> pm_{};
-    std::vector<double*> ptree_{};
+    std::vector<std::pair<int, int>> ptree_{};
 };

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_;
@@ -213,16 +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 */
-        free(nde->_d); /* along with _rhs, _a_matelm, _b_matelm, _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;
+        free(nde->_rhs);
+        nde->v = nullptr;
+        nde->_a = nullptr;
+        nde->_b = nullptr;
+        nde->_rhs = nullptr;
     }
 }
 
@@ -294,12 +290,7 @@ 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->_rhs = nde->_d + nlayer;
-        nde->_a_matelm = nde->_rhs + nlayer;
-        nde->_b_matelm = nde->_a_matelm + nlayer;
-        nde->_x12 = nde->_b_matelm + nlayer;
-        nde->_x21 = nde->_x12 + nlayer;
+        nde->_rhs = (double**) ecalloc(nlayer * 1, sizeof(double*));
     }
 }
 
@@ -424,64 +415,66 @@ 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));
+        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 */
-        *nde->_x12[0] -= d;
-        *nde->_x21[0] -= d;
+        m(index - 1, index) -= d;
+        m(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;
+                m(index + j, index + j) += mfac;
                 ++j;
                 if (j == nrn_nlayer_extracellular) {
                     break;
                 }
-                *nde->_d[j] += mfac;
-                *nde->_x12[j] -= mfac;
-                *nde->_x21[j] -= mfac;
+                m(index + j, index + j) += mfac;
+                m(index - 1 + j, index + j) -= mfac;
+                m(index + j, index - 1 + j) -= mfac;
             }
-            pnde = pnd->extnode;
+            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) {
-                    *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(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];
                 }
             }
         }

src/nrnoc/fadvance.cpp

@@ -13,7 +13,6 @@
 #include "utils/profile/profiler_interface.h"
 #include "nonvintblock.h"
 #include "nrncvode.h"
-#include "spmatrix.h"
 
 #include <vector>
 
@@ -672,13 +671,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) {

src/nrnoc/membfunc.cpp

@@ -2,6 +2,7 @@
 
 #include "multicore.h"
 #include "section.h"
+#include "ocmatrix.h"
 
 #include <cassert>
 
@@ -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<double> _nrn_mechanism_get_area_handle(Node* node) {
     if (node) {
         return node->area_handle();

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<double> _nrn_mechanism_get_area_handle(Node*);
 [[nodiscard]] Section* _nrn_mechanism_get_child(Section*);
 [[nodiscard]] int _nrn_mechanism_get_nnode(Section*);

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 <cstdint>
 #include <condition_variable>
@@ -345,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;
@@ -440,8 +441,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;

src/nrnoc/multicore.h

@@ -29,6 +29,7 @@ actual_v, etc.
 #include "membfunc.h"
 
 #include <cstddef>
+class OcSparseMatrix;
 
 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_ */

src/nrnoc/section.h

@@ -181,9 +181,6 @@ struct Node {
         return _node_handle.non_owning_handle();
     }
     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 */