Hines Matrix solver and cell permute docs (#726)

* Initial comments added for the cell permutation algorithms
* Additional comments on permutation strategies

Co-authored-by: Christos Kotsalos <christos.kotsalos@epfl.ch>

Author: iomaganaris

coreneuron/network/tnode.hpp

@@ -17,20 +17,42 @@ class TNode;
 
 using VecTNode = std::vector<TNode*>;
 
+/**
+ * \class TNode
+ * \brief TNode is the tree node that represents the tree of the compartments
+ */
 class TNode {
   public:
     TNode(int ix);
     virtual ~TNode();
     TNode* parent;
     VecTNode children;
-    size_t mkhash();
-    size_t hash;
-    size_t treesize;
-    size_t nodevec_index;
-    size_t treenode_order;
-    size_t level;
-    size_t cellindex;
-    size_t groupindex;
+    size_t mkhash();  /// Hash algorith that generates a hash based on the hash of the children and
+                      /// the number of compartments of the children
+    size_t hash;      /// Hash value generated by mkhash
+    size_t treesize;  /// Total number of compartments from the current node and below
+    size_t nodevec_index;   /// index in nodevec that is set in check()
+                            /// In cell permute 2 this is set as Breadth First traversal
+    size_t treenode_order;  /// For cell permute 1 (Interleaved):
+                            /// - This is the id given to the compartments based on a Breadth First
+                            /// access on the tree that is created in the original circuit
+                            /// - This is what makes the cell ordering interleaved
+                            /// For cell permute 2 (Constant Depth):
+                            /// VVVTN: Vector (groups of cells) of vector (levels of this group of
+                            /// cells. Maxsize = maxlevel) of vector of TNodes This changes 3 times
+                            /// during cell permute 2:
+                            /// 1. According to the sorting of the nodes of each level
+                            /// 2. According to the sorting of the parents' treenode_order of the
+                            /// previous ordering
+                            /// 3. According to children and parents data races. Parents and
+                            /// children of the tree are moved by question2() so that threads that
+                            /// exist on the same warp don't have data races when updating the
+                            /// children and parent variables, so that threads have to wait in
+                            /// atomic instructions. If there are any races then those are solved by
+                            /// atomic instructions.
+    size_t level;           /// level of of this compartment in the tree
+    size_t cellindex;       /// Cell ID that this compartment belongs to
+    size_t groupindex;      /// Initialized index / groupsize
     int nodeindex;
 };
 

coreneuron/permute/cellorder.cpp

@@ -616,7 +616,7 @@ void solve_interleaved2(int ith) {
             int ic = icore & (warpsize - 1);  // figure out the & mask
             int ncycle = ncycles[iwarp];
             int* stride = strides + stridedispl[iwarp];
-            int root = rootbegin[iwarp];
+            int root = rootbegin[iwarp];  // cell ID -> [0, ncell)
             int lastroot = rootbegin[iwarp + 1];
             int firstnode = nodebegin[iwarp];
             int lastnode = nodebegin[iwarp + 1];

coreneuron/permute/cellorder1.cpp

@@ -396,6 +396,7 @@ void check(VecTNode& nodevec) {
             ncell++;
         }
     }
+    ///  Check that the first compartments of nodevec are the root nodes (cells)
     for (size_t i = 0; i < ncell; ++i) {
         nrn_assert(nodevec[i]->parent == nullptr);
     }
@@ -543,12 +544,16 @@ static void admin1(int ncell,
         }
     }
 
+    // this vector is used to move from one compartment to the other (per cell)
+    // its length is equal to the cell with the highest number of compartments
     stride = (int*) ecalloc_align(nstride + 1, sizeof(int));
     for (int i = 0; i <= nstride; ++i) {
         stride[i] = 0;
     }
     for (size_t i = ncell; i < nodevec.size(); ++i) {
         TNode& nd = *nodevec[i];
+        // compute how many compartments with the same order
+        // treenode_order : defined in breadth first fashion (for each cell separately)
         stride[nd.treenode_order - 1] += 1;  // -1 because treenode order includes root
     }
 }
@@ -590,6 +595,26 @@ static size_t stride_length(size_t begin, size_t end, VecTNode& nodevec) {
     return end - begin;
 }
 
+/**
+ * \brief Prepare for solve_interleaved2
+ *
+ * One group of cells per warp.
+ *
+ * warp[i] has a number of compute cycles (ncycle[i])
+ * the index of its first root (rootbegin[i], last rootbegin[nwarp] = ncell)
+ * the index of its first node (nodebegin[i], last nodebegin[nwarp] = nnode)
+ *
+ * Each compute cycle has a stride
+ * A stride is how many nodes are processed by a warp in one compute cycle
+ * There are nstride strides. nstride is the sum of ncycles of all warps.
+ * warp[i] has ncycle[i] strides
+ * same as sum of ncycle
+ * warp[i] has a stridedispl[i] which is stridedispl[i-1] + ncycle[i].
+ * ie. The zeroth cycle of warp[j] works on stride[stridedispl[j]]
+ * The value of a stride beginning at node i (node i is computed by core 0 of
+ * some warp for some cycle) is determined by stride_length(i, j, nodevec)
+ *
+ */
 static void admin2(int ncell,
                    VecTNode& nodevec,
                    int& nwarp,
@@ -628,6 +653,7 @@ static void admin2(int ncell,
         size_t j = size_t(nodebegin[iwarp + 1]);
         int nc = 0;
         size_t i = nodebegin[iwarp];
+        // in this loop we traverse all the children of all the cells in the current warp (iwarp)
         while (i < j) {
             i += stride_length(i, j, nodevec);
             ++nc;  // how many times the warp should loop in order to finish with all the tree

coreneuron/permute/cellorder2.cpp

@@ -69,6 +69,7 @@ static void sortlevel(VTN& level) {
     }
 }
 
+// TODO: refactor since sortlevel() is traversing the nodes in same order
 static void set_treenode_order(VVTN& levels) {
     size_t order = 0;
     for (auto& level: levels) {
@@ -344,8 +345,8 @@ static void question2(VVTN& levels) {
     // belong to separate warps
 
     // work backward and check the distance from parent to children.
-    // if parent in different group then there is no vitiating race.
-    // if children in different group then ther is no race (satisfied by
+    // if parent in different group (warp?) then there is no vitiating race.
+    // if children in different group (warp?) then ther is no race (satisfied by
     // atomic).
     // If there is a vitiating race, then figure out how many nodes
     // need to be inserted just before the parent to avoid the race.