From 3a1230871b94b6ec353249ebcc8401b40756e033 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 6 Jun 2024 14:00:47 +0200 Subject: [PATCH] Eigen is better with ColMajor --- src/nrniv/matrixmap.cpp | 2 +- src/nrnoc/multicore.h | 2 +- src/nrnoc/solve.cpp | 2 +- src/nrnoc/treeset.cpp | 24 ++++++++++++------------ 4 files changed, 15 insertions(+), 15 deletions(-) diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp index 050685ce18..c8a66e0a30 100644 --- 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 index 8272dfca9e..39b6173540 100644 --- 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* _sparseMat{}; /* handle to general sparse matrix */ + Eigen::SparseMatrix* _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 index 0a8029a531..83d27b7011 100644 --- 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> lu(*_nt->_sparseMat); + Eigen::SparseLU> lu(*_nt->_sparseMat); Eigen::Map 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 index ef58e3130c..b5984e2f94 100644 --- 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& m_ = *_nt->_sparseMat; + Eigen::SparseMatrix& m_ = *_nt->_sparseMat; for (int k = 0; k < m_.outerSize(); ++k) { - for (Eigen::SparseMatrix::InnerIterator it(m_, k); it; ++it) { + for (Eigen::SparseMatrix::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(neqn, neqn); + nt->_sparseMat = new Eigen::SparseMatrix(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;