Eigen is better with ColMajor

neuronsimulator · Jun 6, 2024 · 3a12308 · 3a12308
1 parent f4bc02d
commit 3a12308
Show file tree

Hide file tree

Showing 4 changed files with 15 additions and 15 deletions.
diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp
@@ -58,7 +58,7 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) {
         if (it == 0 || jt == 0) {
             ptree_[plen_] = &place_holder;
         } else {
-            ptree_[plen_] = &_nt->_sparseMat->coeffRef(it - 1, jt - 1);
+            ptree_[plen_] = &_nt->_sparseMat->coeffRef(jt - 1, it - 1);
         }
         ++plen_;
     }

diff --git a/src/nrnoc/multicore.h b/src/nrnoc/multicore.h
@@ -95,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<double, 1, int>* _sparseMat{}; /* handle to general sparse matrix */
+    Eigen::SparseMatrix<double, 0, 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_ */

diff --git 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);
         _nt->_sparseMat->makeCompressed();
-        Eigen::SparseLU<Eigen::SparseMatrix<double, Eigen::RowMajor>> lu(*_nt->_sparseMat);
+        Eigen::SparseLU<Eigen::SparseMatrix<double>> 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);

diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp
@@ -496,9 +496,9 @@ void nrn_lhs(neuron::model_sorted_token const& sorted_token, NrnThread& nt) {
     }
 
     if (use_sparse13) {
-        Eigen::SparseMatrix<double, Eigen::RowMajor>& m_ = *_nt->_sparseMat;
+        Eigen::SparseMatrix<double>& m_ = *_nt->_sparseMat;
         for (int k = 0; k < m_.outerSize(); ++k) {
-            for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(m_, k); it; ++it) {
+            for (Eigen::SparseMatrix<double>::InnerIterator it(m_, k); it; ++it) {
                 it.valueRef() = 0.;
             }
         }
@@ -1991,7 +1991,7 @@ 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::SparseMatrix<double, Eigen::RowMajor>(neqn, neqn);
+        nt->_sparseMat = new Eigen::SparseMatrix<double>(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) {
@@ -2019,18 +2019,18 @@ static void nrn_matrix_node_alloc(void) {
                     for (int ie = 0; ie < nlayer; ++ie) {
                         int k = i + ie + 1;
                         triplets.emplace_back(k - 1, k - 1, 0.);
-                        triplets.emplace_back(k - 1, k - 2, 0.);
                         triplets.emplace_back(k - 2, k - 1, 0.);
+                        triplets.emplace_back(k - 1, k - 2, 0.);
                     }
                 }
                 if (pnd) {
                     int j = pnd->eqn_index_;
-                    triplets.emplace_back(j - 1, i - 1, 0.);
                     triplets.emplace_back(i - 1, j - 1, 0.);
+                    triplets.emplace_back(j - 1, i - 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.);
+                            triplets.emplace_back(j + ie, i + ie, 0.);
                         }
                 }
             }
@@ -2055,20 +2055,20 @@ static void nrn_matrix_node_alloc(void) {
                         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);
+                        nde->_x21[ie] = &nt->_sparseMat->coeffRef(k - 2, k - 1);
+                        nde->_x12[ie] = &nt->_sparseMat->coeffRef(k - 1, k - 2);
                     }
                 }
                 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);
+                    nd->_a_matelm = &nt->_sparseMat->coeffRef(i - 1, j - 1);
+                    nd->_b_matelm = &nt->_sparseMat->coeffRef(j - 1, i - 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);
+                            nde->_a_matelm[ie] = &nt->_sparseMat->coeffRef(k - 1, kp - 1);
+                            nde->_b_matelm[ie] = &nt->_sparseMat->coeffRef(kp - 1, k - 1);
                         }
                 } else { /* not needed if index starts at 1 */
                     nd->_a_matelm = nullptr;