cvode no_cap uses direct backing store instead of Node*

src/nrncvode/cvodeobj.cpp

@@ -1050,6 +1050,11 @@ void NetCvode::set_CVRhsFn() {
     }
 }
 
+static Section* cv_rootsec(const Cvode* cv) {
+    NrnThread* nt = cv->nth_ ? cv->nth_ : nrn_threads;
+    return nt->_v_node[cv->ctd_[0].vnode_begin_index_]->sec;
+}
+
 int Cvode::cvode_init(double) {
     int err = SUCCESS;
     // note, a change in stiff_ due to call of stiff() destroys mem_
@@ -1062,7 +1067,7 @@ int Cvode::cvode_init(double) {
         if (err != SUCCESS) {
             Printf("Cvode %p %s CVReInit error %d\n",
                    fmt::ptr(this),
-                   secname(ctd_[0].v_node_[ctd_[0].rootnodecount_]->sec),
+                   secname(cv_rootsec(this)),
                    err);
             return err;
         }
@@ -1078,7 +1083,7 @@ int Cvode::cvode_init(double) {
         if (err != SUCCESS) {
             Printf("Cvode %p %s CVodeMalloc error %d\n",
                    fmt::ptr(this),
-                   secname(ctd_[0].v_node_[ctd_[0].rootnodecount_]->sec),
+                   secname(cv_rootsec(this)),
                    err);
             return err;
         }
@@ -1317,7 +1322,7 @@ int Cvode::cvode_advance_tn(neuron::model_sorted_token const& sorted_token) {
     if (err < 0) {
         Printf("CVode %p %s advance_tn failed, err=%d.\n",
                fmt::ptr(this),
-               secname(ctd_[0].v_node_[ctd_[0].rootnodecount_]->sec),
+               secname(cv_rootsec(this)),
                err);
         pf_(t_, y_, nullptr, &opaque);
         return err;
@@ -1360,7 +1365,7 @@ int Cvode::cvode_interpolate(double tout) {
     if (err < 0) {
         Printf("CVode %p %s interpolate failed, err=%d.\n",
                fmt::ptr(this),
-               secname(ctd_[0].v_node_[ctd_[0].rootnodecount_]->sec),
+               secname(cv_rootsec(this)),
                err);
         return err;
     }
@@ -1404,7 +1409,7 @@ void Cvode::statistics() {
 #if 1
     Printf("\nCvode instance %p %s statistics : %d %s states\n",
            fmt::ptr(this),
-           secname(ctd_[0].v_node_[ctd_[0].rootnodecount_]->sec),
+           secname(cv_rootsec(this)),
            neq_,
            (use_daspk_ ? "IDA" : "CVode"));
     Printf("   %d advance_tn, %d interpolate, %d init (%d due to at_time)\n",

src/nrncvode/cvodeobj.h

@@ -58,21 +58,29 @@ class CvodeThreadData {
     virtual ~CvodeThreadData();
     void delete_memb_list(CvMembList*);
 
-    int no_cap_count_;  // number of nodes with no capacitance
-    int no_cap_child_count_;
-    Node** no_cap_node_;
-    Node** no_cap_child_;  // connected to nodes that have no capacitance
+    std::vector<int> no_cap_indices_;
+    std::vector<int> no_cap_child_indices_;
     CvMembList* cv_memb_list_;
     CvMembList* cmlcap_;
     CvMembList* cmlext_;       // used only by daspk
     CvMembList* no_cap_memb_;  // used only by cvode, point processes in the no cap nodes
     BAMechList* before_breakpoint_;
     BAMechList* after_solve_;
     BAMechList* before_step_;
-    int rootnodecount_;
-    int v_node_count_;
-    Node** v_node_;
-    Node** v_parent_;
+
+    // Analogous to NrnThread ncell and end to allow similar SoA container
+    // indexing as in fixed step triang and bksub in nrnoc/solve.cpp.
+    // The v, rhs, d, a, b, parent storage containers are assumed to be
+    // sorted for cell contiguity (when lvardt is used) though roots are
+    // all at the beginning of each thread. The following 4 indices are
+    // used to directly access regions of those containers for tree setup
+    // and solving. This eliminates the use of low performance
+    // Node** v_node_, v_parent_ to access the containers.
+    int rootnode_begin_index_;
+    int rootnode_end_index_;
+    int vnode_begin_index_;
+    int vnode_end_index_;
+
     PreSynList* psl_th_;  // with a threshold
     std::list<WatchCondition*>* watch_list_;
     // since scatter/gather are hot loops, don't want to use data_handle

src/nrncvode/cvtrset.cpp

@@ -14,28 +14,35 @@ void Cvode::rhs(neuron::model_sorted_token const& sorted_token, NrnThread* _nt)
     if (diam_changed) {
         recalc_diam();
     }
-    if (z.v_node_count_ == 0) {
+    if (z.rootnode_end_index_ == 0) {
         return;
     }
-    for (int i = 0; i < z.v_node_count_; ++i) {
-        NODERHS(z.v_node_[i]) = 0.;
+    auto* const vec_rhs = _nt->node_rhs_storage();
+    for (int i = z.rootnode_begin_index_; i < z.rootnode_end_index_; ++i) {
+        vec_rhs[i] = 0.;
+    }
+    for (int i = z.vnode_begin_index_; i < z.vnode_end_index_; ++i) {
+        vec_rhs[i] = 0.;
     }
     auto const vec_sav_rhs = _nt->node_sav_rhs_storage();
     if (vec_sav_rhs) {
-        for (int i = 0; i < z.v_node_count_; ++i) {
-            Node* nd = z.v_node_[i];
-            vec_sav_rhs[nd->v_node_index] = 0;
+        for (int i = z.rootnode_begin_index_; i < z.rootnode_end_index_; ++i) {
+            vec_sav_rhs[i] = 0.;
+        }
+        for (int i = z.vnode_begin_index_; i < z.vnode_end_index_; ++i) {
+            vec_sav_rhs[i] = 0.;
         }
     }
 
     rhs_memb(sorted_token, z.cv_memb_list_, _nt);
-    auto const vec_rhs = _nt->node_rhs_storage();
     nrn_nonvint_block_current(_nt->end, vec_rhs, _nt->id);
 
     if (vec_sav_rhs) {
-        for (int i = 0; i < z.v_node_count_; ++i) {
-            auto const node_index = z.v_node_[i]->v_node_index;
-            vec_sav_rhs[node_index] -= vec_rhs[node_index];
+        for (int i = z.rootnode_begin_index_; i < z.rootnode_end_index_; ++i) {
+            vec_sav_rhs[i] -= vec_rhs[i];
+        }
+        for (int i = z.vnode_begin_index_; i < z.vnode_end_index_; ++i) {
+            vec_sav_rhs[i] -= vec_rhs[i];
         }
     }
 
@@ -46,13 +53,13 @@ void Cvode::rhs(neuron::model_sorted_token const& sorted_token, NrnThread* _nt)
     auto const vec_a = _nt->node_a_storage();
     auto const vec_b = _nt->node_b_storage();
     auto const vec_v = _nt->node_voltage_storage();
-    for (auto i = z.rootnodecount_; i < z.v_node_count_; ++i) {
-        auto const node_i = z.v_node_[i]->v_node_index;
-        auto const parent_i = z.v_parent_[i]->v_node_index;
-        auto const dv = vec_v[parent_i] - vec_v[node_i];
-        // our connection coefficients are negative
-        vec_rhs[node_i] -= vec_b[node_i] * dv;
-        vec_rhs[parent_i] += vec_a[node_i] * dv;
+    auto* const parents = _nt->_v_parent_index;
+    for (int i = z.vnode_begin_index_; i < z.vnode_end_index_; ++i) {
+        auto const parent_i = parents[i];
+        auto const dv = vec_v[parent_i] - vec_v[i];
+        // our connection coefficients are negative so
+        vec_rhs[i] -= vec_b[i] * dv;
+        vec_rhs[parent_i] += vec_a[i] * dv;
     }
 }
 
@@ -80,11 +87,15 @@ void Cvode::lhs(neuron::model_sorted_token const& sorted_token, NrnThread* _nt)
     int i;
 
     CvodeThreadData& z = CTD(_nt->id);
-    if (z.v_node_count_ == 0) {
+    if (z.rootnode_end_index_ == 0) {
         return;
     }
-    for (i = 0; i < z.v_node_count_; ++i) {
-        NODED(z.v_node_[i]) = 0.;
+    auto* const vec_d = _nt->node_d_storage();
+    for (int i = z.rootnode_begin_index_; i < z.rootnode_end_index_; ++i) {
+        vec_d[i] = 0.;
+    }
+    for (int i = z.vnode_begin_index_; i < z.vnode_end_index_; ++i) {
+        vec_d[i] = 0.;
     }
 
     lhs_memb(sorted_token, z.cv_memb_list_, _nt);
@@ -96,11 +107,12 @@ void Cvode::lhs(neuron::model_sorted_token const& sorted_token, NrnThread* _nt)
     // fast_imem not needed since exact icap added in nrn_div_capacity
 
     /* now add the axial currents */
-    for (i = 0; i < z.v_node_count_; ++i) {
-        NODED(z.v_node_[i]) -= NODEB(z.v_node_[i]);
-    }
-    for (i = z.rootnodecount_; i < z.v_node_count_; ++i) {
-        NODED(z.v_parent_[i]) -= NODEA(z.v_node_[i]);
+    auto* const vec_a = _nt->node_a_storage();
+    auto* const vec_b = _nt->node_b_storage();
+    auto* const parent_i = _nt->_v_parent_index;
+    for (int i = z.vnode_begin_index_; i < z.vnode_end_index_; ++i) {
+        vec_d[i] -= vec_b[i];
+        vec_d[parent_i[i]] -= vec_a[i];
     }
 }
 
@@ -125,33 +137,33 @@ void Cvode::lhs_memb(neuron::model_sorted_token const& sorted_token,
 
 /* triangularization of the matrix equations */
 void Cvode::triang(NrnThread* _nt) {
-    Node *nd, *pnd;
-    double p;
-    int i;
     CvodeThreadData& z = CTD(_nt->id);
 
-    for (i = z.v_node_count_ - 1; i >= z.rootnodecount_; --i) {
-        nd = z.v_node_[i];
-        pnd = z.v_parent_[i];
-        p = NODEA(nd) / NODED(nd);
-        NODED(pnd) -= p * NODEB(nd);
-        NODERHS(pnd) -= p * NODERHS(nd);
+    auto* const vec_a = _nt->node_a_storage();
+    auto* const vec_b = _nt->node_b_storage();
+    auto* const vec_d = _nt->node_d_storage();
+    auto* const vec_rhs = _nt->node_rhs_storage();
+    for (int i = z.vnode_end_index_ - 1; i >= z.vnode_begin_index_; --i) {
+        double const p = vec_a[i] / vec_d[i];
+        auto const parent_i = _nt->_v_parent_index[i];
+        vec_d[parent_i] -= p * vec_b[i];
+        vec_rhs[parent_i] -= p * vec_rhs[i];
     }
 }
 
 /* back substitution to finish solving the matrix equations */
 void Cvode::bksub(NrnThread* _nt) {
-    Node *nd, *cnd;
-    int i;
     CvodeThreadData& z = CTD(_nt->id);
 
-    for (i = 0; i < z.rootnodecount_; ++i) {
-        NODERHS(z.v_node_[i]) /= NODED(z.v_node_[i]);
+    auto* const vec_b = _nt->node_b_storage();
+    auto* const vec_d = _nt->node_d_storage();
+    auto* const vec_rhs = _nt->node_rhs_storage();
+    auto* const parent_i = _nt->_v_parent_index;
+    for (int i = z.rootnode_begin_index_; i < z.rootnode_end_index_; ++i) {
+        vec_rhs[i] /= vec_d[i];
     }
-    for (i = z.rootnodecount_; i < z.v_node_count_; ++i) {
-        cnd = z.v_node_[i];
-        nd = z.v_parent_[i];
-        NODERHS(cnd) -= NODEB(cnd) * NODERHS(nd);
-        NODERHS(cnd) /= NODED(cnd);
+    for (int i = z.vnode_begin_index_; i < z.vnode_end_index_; ++i) {
+        vec_rhs[i] -= vec_b[i] * vec_rhs[parent_i[i]];
+        vec_rhs[i] /= vec_d[i];
     }
 }