treeset: use eigen instead of sparse13

src/nrncvode/occvode.cpp

@@ -1,3 +1,5 @@
+#include <Eigen/Eigen>
+
 #include <../../nrnconf.h>
 #include <errno.h>
 #include "nrn_ansi.h"
@@ -14,6 +16,7 @@
 
 #include <utility>
 
+
 extern void setup_topology(), v_setup_vectors();
 extern void recalc_diam();
 extern int nrn_errno_check(int);
@@ -23,9 +26,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_)();
@@ -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->_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) {

src/nrniv/matrixmap.cpp

@@ -2,7 +2,7 @@
 #include "matrixmap.h"
 #include <vector>
 
-#include "spmatrix.h"
+#include <Eigen/Eigen>
 
 MatrixMap::MatrixMap(Matrix& mat)
     : m_(mat) {}
@@ -24,6 +24,7 @@ void MatrixMap::add(double fac) {
 }
 
 void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) {
+    static double place_holder = 0.;
     NrnThread* _nt = nrn_threads;
     mmfree();
 
@@ -54,7 +55,11 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) {
         } else {
             jt = start + j - nnode;
         }
-        ptree_[plen_] = spGetElement(_nt->_sp13mat, it, jt);
+        if (it == 0 || jt == 0) {
+            ptree_[plen_] = &place_holder;
+        } else {
+            ptree_[plen_] = &_nt->_sparseMat->coeffRef(it - 1, jt - 1);
+        }
         ++plen_;
     }
 }

src/nrniv/nrndae.cpp

@@ -221,7 +221,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->_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];
     }
 }
 

src/nrnoc/fadvance.cpp

@@ -13,7 +13,7 @@
 #include "utils/profile/profiler_interface.h"
 #include "nonvintblock.h"
 #include "nrncvode.h"
-#include "spmatrix.h"
+#include <Eigen/Eigen>
 
 #include <vector>
 
@@ -672,11 +672,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->_sparseMat << 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->_sparseMat << std::endl;
+            int n = _nt->_sparseMat->cols();
+            for (int i = 1; i <= n; ++i) {
                 Printf("%d %g\n", i, _nt->actual_rhs(i));
             }
         }

src/nrnoc/multicore.cpp

@@ -46,6 +46,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) \
@@ -337,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 = 0;
+                nt->_sparseMat = nullptr;
                 nt->_ctime = 0.0;
                 nt->_vcv = 0;
                 nt->_node_data_offset = 0;
@@ -438,9 +440,9 @@ void nrn_threads_free() {
             free(nt->_ecell_children);
             nt->_ecell_children = NULL;
         }
-        if (nt->_sp13mat) {
-            spDestroy(nt->_sp13mat);
-            nt->_sp13mat = 0;
+        if (nt->_sparseMat) {
+            delete nt->_sparseMat;
+            nt->_sparseMat = nullptr;
         }
         nt->end = 0;
         nt->ncell = 0;

src/nrnoc/multicore.h

@@ -27,6 +27,10 @@ actual_v, etc.
 */
 
 #include "membfunc.h"
+namespace Eigen {
+template <typename _Scalar, int _Options, typename _StorageIndex>
+class SparseMatrix;
+}
 
 #include <cstddef>
 
@@ -89,12 +93,12 @@ 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
-                                    need to be transfered to actual_rhs */
-    char* _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_ */
+    double* _sparse_rhs; /* rhs matrix for sparse13 solver. updates to and from this vector
+                          need to be transfered to actual_rhs */
+    Eigen::SparseMatrix<double, 1, int>* _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) */

src/nrnoc/solve.cpp

@@ -57,7 +57,8 @@ node.v + extnode.v[0]
 #include "nrnmpiuse.h"
 #include "ocnotify.h"
 #include "section.h"
-#include "spmatrix.h"
+#include <Eigen/Eigen>
+#include <Eigen/LU>
 #include "treeset.h"
 
 #include <stdio.h>
@@ -364,22 +365,13 @@ 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);
-        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);
+        _nt->_sparseMat->makeCompressed();
+        Eigen::SparseLU<Eigen::SparseMatrix<double, Eigen::RowMajor>> lu(*_nt->_sparseMat);
+        Eigen::Map<Eigen::VectorXd> 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);
     } else {