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,9 +16,11 @@ void MatrixMap::mmfree() {
 }
 
 void MatrixMap::add(double fac) {
+    NrnThread* _nt = nrn_threads;
     for (int i = 0; i < pm_.size(); ++i) {
         auto [it, jt] = pm_[i];
-        *ptree_[i] += fac * m_(it, jt);
+        auto [im, jm] = ptree_[i];
+        _nt->_sp13mat->coeff(im, jm) += fac * m_(it, jt);
     }
 }
 
@@ -38,7 +38,6 @@ int MatrixMap::compute_index(int i, int start, int nnode, Node** nodes, int* lay
 }
 
 void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) {
-    NrnThread* _nt = nrn_threads;
     mmfree();
 
     std::vector<std::pair<int, int>> nzs = m_.nonzeros();
@@ -49,7 +48,7 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) {
         int jt = compute_index(j, start, nnode, nodes, layer);
         if (it != 0 && jt != 0) {
             pm_.emplace_back(i, j);
-            ptree_.emplace_back(spGetElement(_nt->_sp13mat, it, jt));
+            ptree_.emplace_back(it - 1, jt - 1);
         }
     }
 }

src/nrniv/matrixmap.h

@@ -116,7 +116,7 @@ class MatrixMap {
 
     // the map
     std::vector<std::pair<int, int>> pm_{};
-    std::vector<double*> ptree_{};
+    std::vector<std::pair<int, int>> ptree_{};
 
     int compute_index(int, int, int, Node**, int*) const;
 };

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*));
     }
 }
 
@@ -426,7 +417,6 @@ 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;
     Memb_list* ml = _nt->_ecell_memb_list;
@@ -441,47 +431,56 @@ void nrn_setup_ext(NrnThread* _nt) {
     /* d contains all the membrane conductances (and capacitance) */
     /* i.e. (cm/dt + di/dvm - dis/dvi)*[dvi] and
         (dis/dvi)*[dvx] */
+    // This loop handle conductances between the node and the first layer
     for (i = 0; i < cnt; ++i) {
+        OcMatrix& m = *_nt->_sp13mat;
         nd = ndlist[i];
+        int index = nd->eqn_index_;
         nde = nd->extnode;
+        int ext_index = nde->eqn_index_;
         d = NODED(nd);
         /* nde->_d only has -ELECTRODE_CURRENT contribution */
-        d = (*nde->_d[0] += NODED(nd));
+        m(ext_index + 0, ext_index + 0) += NODED(nd);
+        d = m.getval(ext_index + 0, ext_index + 0);
         /* now d is only the membrane current contribution */
         /* i.e. d =  cm/dt + di/dvm */
-        *nde->_x12[0] -= d;
-        *nde->_x21[0] -= d;
+        m(index - 1, ext_index + 0) -= d;
+        m(ext_index + 0, index - 1) -= d;
 #if I_MEMBRANE
         ml->data(i, sav_g_index) = d;
 #endif
     }
     /* series resistance, capacitance, and axial terms. */
+    // This look takes care of the conductances between layers
     for (i = 0; i < cnt; ++i) {
+        OcMatrix& m = *_nt->_sp13mat;
         nd = ndlist[i];
         nde = nd->extnode;
+        int ext_index = nde->eqn_index_;
         pnd = _nt->_v_parent[nd->v_node_index];
         if (pnd) {
             /* series resistance and capacitance to ground */
             j = 0;
             for (;;) { /* between j and j+1 layer */
                 mfac = (*nde->param[xg_index_ext(j)] + *nde->param[xc_index_ext(j)] * cfac);
-                *nde->_d[j] += mfac;
+                m(ext_index + j, ext_index + j) += mfac;
                 ++j;
                 if (j == nrn_nlayer_extracellular) {
                     break;
                 }
-                *nde->_d[j] += mfac;
-                *nde->_x12[j] -= mfac;
-                *nde->_x21[j] -= mfac;
+                m(ext_index + j, ext_index + j) += mfac;
+                m(ext_index + j - 1, ext_index + j) -= mfac;
+                m(ext_index + j, ext_index + j - 1) -= mfac;
             }
             pnde = pnd->extnode;
             /* axial connections */
             if (pnde) { /* parent sec may not be extracellular */
+                int parent_ext_index = pnde->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(ext_index + j, ext_index + j) -= nde->_b[j];
+                    m(parent_ext_index + j, parent_ext_index + j) -= nde->_a[j];
+                    m(parent_ext_index + j, ext_index + j) += nde->_a[j];
+                    m(ext_index + j, parent_ext_index + j) += nde->_b[j];
                 }
             }
         }

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) {
+    OcMatrix& m = *nt->_sp13mat;
+    return m.coeff(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 OcMatrix;
 
 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 */
+    OcMatrix* _sp13mat;              /* handle to general sparse matrix */
     Memb_list* _ecell_memb_list; /* normally nullptr */
     Node** _ecell_children;      /* nodes with no extcell but parent has it */
     void* _vcv;                  /* replaces old cvode_instance and nrn_cvode_ */

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 */

src/nrnoc/section_fwd.hpp

@@ -40,12 +40,9 @@ struct Extnode {
     double* v; /* v external. */
     double* _a;
     double* _b;
-    double** _d;
-    double** _rhs; /* d, rhs, a, and b are analogous to those in node */
-    double** _a_matelm;
-    double** _b_matelm;
-    double** _x12; /* effect of v[layer] on eqn layer-1 (or internal)*/
-    double** _x21; /* effect of v[layer-1 or internal] on eqn layer*/
+    double** _rhs; /* rhs, a, and b are analogous to those in node */
+
+    int eqn_index_;
 };
 #endif
 

src/nrnoc/solve.cpp

@@ -57,8 +57,9 @@ node.v + extnode.v[0]
 #include "nrnmpiuse.h"
 #include "ocnotify.h"
 #include "section.h"
-#include "spmatrix.h"
+#include "ocmatrix.h"
 #include "treeset.h"
+#include "ivocvect.h"
 
 #include <stdio.h>
 #include <stdlib.h>
@@ -367,19 +368,17 @@ void nrn_solve(NrnThread* _nt) {
         int e;
         nrn_thread_error("solve use_sparse13");
         update_sp13_mat_based_on_actual_d(_nt);
-        e = spFactor(_nt->_sp13mat);
-        if (e != spOKAY) {
-            switch (e) {
-            case spZERO_DIAG:
-                hoc_execerror("spFactor error:", "Zero Diagonal");
-            case spNO_MEMORY:
-                hoc_execerror("spFactor error:", "No Memory");
-            case spSINGULAR:
-                hoc_execerror("spFactor error:", "Singular");
-            }
-        }
+
         update_sp13_rhs_based_on_actual_rhs(_nt);
-        spSolve(_nt->_sp13mat, _nt->_sp13_rhs, _nt->_sp13_rhs);
+        int size = _nt->_sp13mat->ncol();
+        std::vector<double> vector_in(_nt->_sp13_rhs + 1, _nt->_sp13_rhs + size + 1);
+        IvocVect in; in.vec() = vector_in;
+        IvocVect out(size);
+        _nt->_sp13mat->solv(&in, &out, false);
+        for (int i = 0; i < size; ++i) {
+            _nt->_sp13_rhs[i + 1] = out[i];
+
+        }
         update_actual_d_based_on_sp13_mat(_nt);
         update_actual_rhs_based_on_sp13_rhs(_nt);
     } else {