Fix MatrixMap

Author: Nicolas Cornu

src/nrniv/matrixmap.cpp

@@ -18,48 +18,43 @@ void MatrixMap::mmfree() {
 }
 
 void MatrixMap::add(double fac) {
-    for (int i = 0; i < plen_; ++i) {
-        *ptree_[i] += fac * (*pm_[i]);
+    NrnThread* _nt = nrn_threads;
+    for (int i = 0; i < pm_.size(); ++i) {
+        auto [im, jm] = pm_[i];
+        auto [it, jt] = ptree_[i];
+        *_nt->_sp13mat->mep(it, jt) += fac * m_.getval(im, jm);
+    }
+}
+
+int MatrixMap::compute_id(int i, int start, int nnode, Node** nodes, int* layer) const {
+    int it;
+    if (i < nnode) {
+        it = nodes[i]->eqn_index_ + layer[i];
+        if (layer[i] > 0 && !nodes[i]->extnode) {
+            it = 0;
+        }
+    } else {
+        it = start + i - nnode;
     }
+    return it;
 }
 
 void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) {
-    static double place_holder = 0;
-    NrnThread* _nt = nrn_threads;
     mmfree();
 
-    plen_ = 0;
     std::vector<int> nonzero_i, nonzero_j;
     m_.nonzeros(nonzero_i, nonzero_j);
-    pm_.resize(nonzero_i.size());
-    ptree_.resize(nonzero_i.size());
+    pm_.reserve(nonzero_i.size());
+    ptree_.reserve(nonzero_i.size());
     for (int k = 0; k < nonzero_i.size(); k++) {
         const int i = nonzero_i[k];
         const int j = nonzero_j[k];
-        int it;
-        if (i < nnode) {
-            it = nodes[i]->eqn_index_ + layer[i];
-            if (layer[i] > 0 && !nodes[i]->extnode) {
-                it = 0;
-            }
-        } else {
-            it = start + i - nnode;
-        }
-        int jt;
-        pm_[plen_] = m_.mep(i, j);
-        if (j < nnode) {
-            jt = nodes[j]->eqn_index_ + layer[j];
-            if (layer[j] > 0 && !nodes[j]->extnode) {
-                jt = 0;
-            }
-        } else {
-            jt = start + j - nnode;
-        }
+        int it = compute_id(i, start, nnode, nodes, layer);
+        int jt = compute_id(j, start, nnode, nodes, layer);
         if (it == 0 || jt == 0) {
-            ptree_[plen_] = &place_holder;
-        } else {
-            ptree_[plen_] = _nt->_sp13mat->mep(it - 1, jt - 1);
+            continue;
         }
-        ++plen_;
+        pm_.emplace_back(i, j);
+        ptree_.emplace_back(it - 1, jt - 1);
     }
 }

src/nrniv/matrixmap.h

@@ -112,10 +112,10 @@ class MatrixMap {
 
 
   private:
+    int compute_id(int i, int start, int nnode, Node** nodes, int* layer) const;
     Matrix& m_;
 
     // the map
-    int plen_ = 0;
-    std::vector<double*> pm_{};
-    std::vector<double*> ptree_{};
+    std::vector<std::pair<int, int>> pm_{};
+    std::vector<std::pair<int, int>> ptree_{};
 };