treeset: use eigen instead of sparse13

src/nrncvode/occvode.cpp

@@ -14,6 +14,8 @@
 
 #include <utility>
 
+#include <Eigen/Eigen>
+
 extern void setup_topology(), v_setup_vectors();
 extern void recalc_diam();
 extern int nrn_errno_check(int);
@@ -356,7 +358,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) {

src/nrniv/matrixmap.cpp

@@ -3,7 +3,7 @@
 #include <vector>
 using std::vector;
 
-#include "spmatrix.h"
+#include <Eigen/Eigen>
 
 MatrixMap::MatrixMap(Matrix& mat)
     : m_(mat)
@@ -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_;
     }
 }

src/nrniv/nrndae.cpp

@@ -6,6 +6,8 @@
 #include "treeset.h"
 #include "utils/enumerate.h"
 
+#include <Eigen/Eigen>
+
 extern int secondorder;
 
 static NrnDAEPtrList nrndae_list;
@@ -221,7 +223,7 @@ void NrnDAE::dkmap(std::vector<neuron::container::data_handle<double>>& pv,
                                                                   y_.data() + i};
         pvdot[bmap_[i] - 1] =
             neuron::container::data_handle<double>{neuron::container::do_not_search,
-                                                   _nt->_sp13_rhs + bmap_[i]};
+                                                   _nt->_sp13_rhs + bmap_[i] - 1};
     }
 }
 
@@ -231,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]];
+        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 +313,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];
     }
 }
 

src/nrnoc/fadvance.cpp

@@ -12,7 +12,7 @@
 #include "utils/profile/profiler_interface.h"
 #include "nonvintblock.h"
 #include "nrncvode.h"
-#include "spmatrix.h"
+#include <Eigen/Eigen>
 
 #include <vector>
 
@@ -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));
             }
         }

src/nrnoc/multicore.cpp

@@ -47,6 +47,8 @@ the handling of v_structure_change as long as possible.
 #include <vector>
 #include <iostream>
 
+#include <Eigen/Eigen>
+
 #define CACHELINE_ALLOC(name, type, size) \
     name = (type*) nrn_cacheline_alloc((void**) &name, size * sizeof(type))
 #define CACHELINE_CALLOC(name, type, size) \
@@ -346,7 +348,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 +444,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

@@ -30,6 +30,12 @@ actual_v, etc.
 
 #include <cstddef>
 
+namespace Eigen {
+template <typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
+class Matrix;
+using MatrixXd = Matrix<double, -1, -1, 0, -1, -1>;
+}  // namespace Eigen
+
 typedef struct NrnThreadMembList { /* patterned after CvMembList in cvodeobj.h */
     struct NrnThreadMembList* next;
     Memb_list* ml;
@@ -91,7 +97,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::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_ */

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 <Eigen/Eigen>
 #include "treeset.h"
 
 #include <stdio.h>
@@ -367,19 +367,10 @@ 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);
+        Eigen::FullPivLU<Eigen::MatrixXd> lu{*_nt->_sp13mat};
+        Eigen::Map<Eigen::VectorXd> 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);
     } else {

src/nrnoc/treeset.cpp

@@ -16,7 +16,7 @@
 #include "ocnotify.h"
 #include "partrans.h"
 #include "section.h"
-#include "spmatrix.h"
+#include <Eigen/Eigen>
 #include "utils/profile/profiler_interface.h"
 #include "multicore.h"
 
@@ -28,8 +28,6 @@
 #include <algorithm>
 #include <string>
 
-extern spREAL* spGetElement(char*, int, int);
-
 int nrn_shape_changed_; /* for notifying Shape class in nrniv */
 double* nrn_mech_wtime_;
 
@@ -334,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];
     }
 }
 
@@ -343,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,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 = 0; 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 {
@@ -498,8 +495,7 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) {
     }
 
     if (use_sparse13) {
-        // Zero the sparse13 matrix
-        spClear(_nt->_sp13mat);
+        _nt->_sp13mat->setZero();
     }
 
     // Make sure the SoA node diagonals are also zeroed (is this needed?)
@@ -1871,8 +1867,8 @@ void nrn_matrix_node_free() {
             free(std::exchange(nt->_sp13_rhs, nullptr));
         }
         if (nt->_sp13mat) {
-            spDestroy(nt->_sp13mat);
-            nt->_sp13mat = (char*) 0;
+            delete nt->_sp13mat;
+            nt->_sp13mat = nullptr;
         }
     }
     diam_changed = 1;
@@ -1935,14 +1931,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 +1957,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;
@@ -1978,11 +1966,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 = spCreate(neqn, 0, &err);
-        if (err != spOKAY) {
-            hoc_execerror("Couldn't create sparse matrix", (char*) 0);
-        }
+        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;
             if (nt->_v_node[in]->extnode) {
@@ -1997,27 +1982,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 = spGetElement(nt->_sp13mat, i, i);
+            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] = spGetElement(nt->_sp13mat, k, k);
-                    nde->_rhs[ie] = nt->_sp13_rhs + k;
-                    nde->_x21[ie] = spGetElement(nt->_sp13mat, k, k - 1);
-                    nde->_x12[ie] = spGetElement(nt->_sp13mat, k - 1, k);
+                    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);
                 }
             }
             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);
+                        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;