Lvardt cell CvMembList.ml.size() should be the number of contiguous regions of the fixed step thread Memb_list

src/nrncvode/cvodeobj.h

@@ -30,7 +30,8 @@ struct model_sorted_token;
  *   contiguous
  * - with ml.size() >= 1 and ml[i].nodecount == 1 when non-contiguous instances need to be processed
  *
- * generic configurations with ml.size() and ml[i].nodecount both larger than one are not supported.
+ * generic configurations with ml.size() and ml[i].nodecount both larger than one are only
+ * supported for the local variable time step method.
  */
 struct CvMembList {
     CvMembList(int type)

src/nrncvode/netcvode.cpp

@@ -1615,6 +1615,20 @@
                 }
             }
 
+            // Modified to also count the nodes and set the offsets for
+            // each CvMembList.ml[contig_region].
+            // The sum of the ml[i].nodecount must equal the mechanism
+            // nodecount for the cell and each ml[i] data must be contiguous.
+            // Ideally the node permutation would be such that each cell
+            // is contiguous. So only needing a ml[0]. That is sadly not
+            // the case with the default permutation. The cell root nodes are
+            // all at the beginning, and thereafter only Section nodes are
+            // contiguous. It would be easy to permute nodes so that each cell
+            // is contiguous (except root node). This would result in a
+            // CvMembList.ml.size() == 1 almost always with an exception of
+            // size() == 2 only for extracellular and for POINT_PROCESSes
+            // located both in the root node and other cell nodes.
+
             for (NrnThreadMembList* tml = _nt->tml; tml; tml = tml->next) {
                 i = tml->index;
                 const Memb_func& mf = memb_func[i];
@@ -1625,42 +1639,57 @@
                     // singly linked list built below
                     int j;
                     for (j = 0; j < ml->nodecount; ++j) {
+                        auto offset = ml->get_storage_offset() + j;
+                        // for each Memb_list instance constructed, keep
+                        // track of its initial storage offset (i.e. offset)
+
                         int inode = ml->nodelist[j]->v_node_index;
+                        int icell = cellnum[inode];
                         Cvode& cv = d.lcv_[cellnum[inode]];
                         CvodeThreadData& z = cv.ctd_[0];
-                        if (!z.cv_memb_list_) {
+
+                        // Circumstances for creating a new CvMembList
+                        // or (due to non-contiguity of a cell),
+                        // appending a Memb_list instance to cml->ml
+                        if (!z.cv_memb_list_) {  // initialize the first
                             cml = new CvMembList{i};
-                            cml->next = nullptr;
-                            assert(cml->ml.size() == 1);
-                            cml->ml[0].nodecount = 0;
                             z.cv_memb_list_ = cml;
+                            cml->next = nullptr;
                             last[cellnum[inode]] = cml;
-                        }
-                        if (last[cellnum[inode]]->index == i) {
-                            assert(last[cellnum[inode]]->ml.size() == 1);
-                            ++last[cellnum[inode]]->ml[0].nodecount;
-                        } else {
+                            assert(cml->ml.size() == 1);
+                            assert(cml->ml[0].nodecount == 0);
+                        } else if (last[cellnum[inode]]->index != i) {  // initialize next
                             cml = new CvMembList{i};
                             last[cellnum[inode]]->next = cml;
                             cml->next = nullptr;
                             last[cellnum[inode]] = cml;
                             assert(cml->ml.size() == 1);
-                            cml->ml[0].nodecount = 1;
+                            assert(cml->ml[0].nodecount == 0);
+                        } else {  // if non-contiguous, append Memb_list
+                            cml = last[cellnum[inode]];
+                            auto& cvml = cml->ml.back();
+                            auto cvml_offset = cvml.get_storage_offset() + cvml.nodecount;
+                            if (cvml_offset != offset) {
+                                // not contiguous, add another Memb_list
+                                // instance to cml->ml
+                                cml->ml.emplace_back(cml->index);
+                                assert(cml->ml.back().nodecount == 0);
+                            }
+                        }
+
+                        auto& cvml = cml->ml.back();
+                        if (cvml.nodecount == 0) {  // first time for this Memb_List
+                            cvml.set_storage_offset(offset);
                         }
+                        // Increment count of last Memb_list in cml->ml.
+                        ++cvml.nodecount;
                     }
                 }
             }
-            // allocate and re-initialize count
+
             std::vector<CvMembList*> cvml(d.nlcv_);
             for (i = 0; i < d.nlcv_; ++i) {
-                cvml[i] = d.lcv_[i].ctd_[0].cv_memb_list_;
-                for (cml = cvml[i]; cml; cml = cml->next) {
-                    // non-contiguous mode, so we're going to create a lot of 1-element Memb_list
-                    // inside cml->ml
-                    cml->ml.reserve(cml->ml[0].nodecount);
-                    // remove the single entry from contiguous mode
-                    cml->ml.clear();
-                }
+                cvml[i] = d.lcv_[i].ctd_[0].cv_memb_list_;  // whole cell in thread
             }
             // fill pointers (and nodecount)
             // now list order is from 0 to n_memb_func
@@ -1674,20 +1703,39 @@
                         int icell = cellnum[ml->nodelist[j]->v_node_index];
                         if (cvml[icell]->index != i) {
                             cvml[icell] = cvml[icell]->next;
-                            assert(cvml[icell] && cvml[icell]->index);
+                            assert(cvml[icell] && cvml[icell]->index == i);
                         }
-                        cml = cvml[icell];
-                        auto& newml = cml->ml.emplace_back(cml->index /* mechanism type */);
-                        newml.nodecount = 1;
-                        newml.nodelist = new Node*[1];
-                        newml.nodelist[0] = ml->nodelist[j];
-                        newml.nodeindices = new int[1]{ml->nodeindices[j]};
-                        newml.prop = new Prop* [1] { ml->prop[j] };
-                        if (!mf.hoc_mech) {
-                            newml.set_storage_offset(ml->get_storage_offset() + j);
-                            newml.pdata = new Datum* [1] { ml->pdata[j] };
+                        auto& cml = cvml[icell];
+                        int kk = j;
+                        for (auto& newml: cml->ml) {
+                            auto nodecount = newml.nodecount;
+                            if (!newml.nodelist) {
+                                // do nodecount of these for ml and then
+                                // skip forward by nodecount in the outer
+                                // ml->nodecount j loop (i.e. a contiguity
+                                // region)
+                                newml.nodelist = new Node*[nodecount];
+                                newml.nodeindices = new int[nodecount];
+                                newml.prop = new Prop*[nodecount];
+                                if (!mf.hoc_mech) {
+                                    newml.pdata = new Datum*[nodecount];
+                                }
+                                for (int k = 0; k < nodecount; ++k) {
+                                    newml.nodelist[k] = ml->nodelist[kk + k];
+                                    newml.nodeindices[k] = ml->nodeindices[kk + k];
+                                    assert(cellnum[newml.nodeindices[k]] ==
+                                           cellnum[ml->nodeindices[j]]);
+                                    newml.prop[k] = ml->prop[kk + k];
+                                    if (!mf.hoc_mech) {
+                                        newml.pdata[k] = ml->pdata[kk + k];
+                                    }
+                                }
+                                kk += nodecount;
+                                newml._thread = ml->_thread;
+                                j += nodecount - 1;
+                                break;
+                            }
                         }
-                        newml._thread = ml->_thread;
                     }
                 }
             }

src/nrncvode/occvode.cpp

@@ -176,7 +176,8 @@ printf("%d Cvode::init_eqn id=%d neq_v_=%d #nonvint=%d #nonvint_extra=%d nvsize=
         if (z.cmlcap_) {
             for (auto& ml: z.cmlcap_->ml) {
                 // support `1 x n` and `n x 1` but not `n x m`
-                assert(z.cmlcap_->ml.size() == 1 || ml.nodecount == 1);
+                // why not? (I'm trying to improve lvardt performance)
+                // assert(z.cmlcap_->ml.size() == 1 || ml.nodecount == 1);
                 zneq_cap_v += ml.nodecount;
             }
         }
@@ -207,13 +208,17 @@ printf("%d Cvode::init_eqn id=%d neq_v_=%d #nonvint=%d #nonvint_extra=%d nvsize=
             // sentinal values for determining no_cap
             NODERHS(z.v_node_[i]) = 1.;
         }
-        for (i = 0; i < zneq_cap_v; ++i) {
-            auto* const node = z.cmlcap_->ml.size() == 1 ? z.cmlcap_->ml[0].nodelist[i]
-                                                         : z.cmlcap_->ml[i].nodelist[0];
-            z.pv_[i] = node->v_handle();
-            z.pvdot_[i] = node->rhs_handle();
-            *z.pvdot_[i] = 0.;  // only ones = 1 are no_cap
-        }
+        i = 0;
+        if (zneq_cap_v)
+            for (auto& ml: z.cmlcap_->ml) {
+                for (int j = 0; j < ml.nodecount; ++j) {
+                    auto* const node = ml.nodelist[j];
+                    z.pv_[i] = node->v_handle();
+                    z.pvdot_[i] = node->rhs_handle();
+                    *z.pvdot_[i] = 0.;  // only ones = 1 are no_cap
+                    ++i;
+                }
+            }
 
         // the remainder are no_cap nodes
         if (z.no_cap_node_) {