Reverse order of children, to obtain same result as in serial

Author: filikat

highs/ipm/hipo/factorhighs/Analyse.cpp

@@ -1297,9 +1297,10 @@ void Analyse::findTreeSplitting() {
   // The tree is parallelised by creating a task for each single node and a task
   // for each group of subtrees.
 
-  // linked lists of children
-  std::vector<Int> head, next;
-  childrenLinkedList(sn_parent_, head, next);
+  // reverse linked lists of children
+  std::vector<Int> head_reverse, next_reverse;
+  childrenLinkedList(sn_parent_, head_reverse, next_reverse);
+  reverseLinkedList(head_reverse, next_reverse);
 
   // compute number of operations for each supernode
   std::vector<double> sn_ops(sn_count_);
@@ -1352,9 +1353,13 @@ void Analyse::findTreeSplitting() {
       // count is in the interval [small_thresh, 2*small_thresh]. Each group
       // corresponds to one task executed in parallel.
 
+      // The children are explored in reverse order, so that they can be synced
+      // in reverse order during the factorisation. In this way, the parallel
+      // and serial code obtain the exact same answer.
+
       double current_ops = 0.0;
       NodeData* current_nodedata = nullptr;
-      Int child = head[sn];
+      Int child = head_reverse[sn];
       while (child != -1) {
         bool is_small = subtree_ops[child] <= small_thresh;
 
@@ -1371,7 +1376,7 @@ void Analyse::findTreeSplitting() {
           if (current_ops > small_thresh) current_nodedata = nullptr;
         }
 
-        child = next[child];
+        child = next_reverse[child];
       }
 
     } else if (sn_parent_[sn] == -1) {

highs/ipm/hipo/factorhighs/Factorise.cpp

@@ -194,16 +194,16 @@ void Factorise::processSupernode(Int sn, const bool should_parallelise) {
       spawn(first_child_[sn], tg, false);
       do_parallelise = false;
     } else {
-      // spawn children of this supernode in reverse order
-      Int child_to_spawn = first_child_reverse_[sn];
+      // spawn children of this supernode in forward order
+      Int child_to_spawn = first_child_[sn];
       while (child_to_spawn != -1) {
         spawn(child_to_spawn, tg);
-        child_to_spawn = next_child_reverse_[child_to_spawn];
+        child_to_spawn = next_child_[child_to_spawn];
       }
 
       // wait for first child to finish, before starting the parent (if there is
       // a first child)
-      if (first_child_reverse_[sn] != -1) sync(first_child_[sn], tg);
+      if (first_child_reverse_[sn] != -1) sync(first_child_reverse_[sn], tg);
     }
   }
 
@@ -255,28 +255,20 @@ void Factorise::processSupernode(Int sn, const bool should_parallelise) {
   // ===================================================
   // Assemble frontal matrices of children
   // ===================================================
-  Int child_sn = first_child_[sn];
-  if (serial) child_sn = first_child_reverse_[sn];
-
+  Int child_sn = first_child_reverse_[sn];
   while (child_sn != -1) {
     // Child contribution is found:
     // - in cliquestack, if we are processing the tree in serial.
     // - in schur_contribution_ if we are processing the tree in parallel.
-    // Children are summed:
-    // - in reverse order in serial, so that the correct child is found on top
-    //   of the CliqueStack.
-    // - in forward order otherwise, so that if a small subtree is synced, and
-    //   it is not the first in its group, it will already have synced when it
-    //   is needed.
-
-    const double* child_clique;
+    // Children are summed always in reverse order.
 
     if (do_parallelise) {
       // sync with spawned child, apart from the first one
-      if (child_sn != first_child_[sn]) sync(child_sn, tg);
+      if (child_sn != first_child_reverse_[sn]) sync(child_sn, tg);
       if (flag_stop_.load(std::memory_order_relaxed)) return;
     }
 
+    const double* child_clique;
     if (!serial) {
       child_clique = schur_contribution_[child_sn].data();
       if (!child_clique) {
@@ -341,10 +333,7 @@ void Factorise::processSupernode(Int sn, const bool should_parallelise) {
     }
 
     // move on to the next child
-    if (!serial)
-      child_sn = next_child_[child_sn];
-    else
-      child_sn = next_child_reverse_[child_sn];
+    child_sn = next_child_reverse_[child_sn];
   }
 
   if (flag_stop_.load(std::memory_order_relaxed)) return;

highs/ipm/hipo/factorhighs/TreeSplitting.h

@@ -60,25 +60,29 @@ class TreeSplitting {
 // where Fj is the first descendant of node j.
 // belong_[j] is true for j=p,a,c,e,f and false for j=b,d,g.
 //
-// - When a is ran, a task is executed that executes the whole subtree of nodes
+// - When a is run, a task is executed that executes the whole subtree of nodes
 //   a, b and d.
-// - When b is ran, nothing happens, since it was already executed as part of
+// - When b is run, nothing happens, since it was already executed as part of
 //   the task that executed a.
-// - When c is ran, a task is created that executes only that supernode.
-// - When d is ran, nothing happens, since it was already executed as part of
+// - When c is run, a task is created that executes only that supernode.
+// - When d is run, nothing happens, since it was already executed as part of
 //   the task that executed a.
-// - When e is ran, a task is executed that executes the whole subtree of nodes
+// - When e is run, a task is executed that executes the whole subtree of nodes
 //   e and g.
-// - When f is ran, a task is created that executes only that supernode.
-// - When g is ran, nothing happens, since it was already executed as part of
+// - When f is run, a task is created that executes only that supernode.
+// - When g is run, nothing happens, since it was already executed as part of
 //   the task that executed e.
 //
 // It is important that node a is processed before b and d, so that when b or d
 // are synced, their operations are already performed. Otherwise, syncing node b
 // would complete even though the operations for node b have not been performed.
-// In other words, children should be synced in forward order, and thus spawned
-// in reverse order.
-//
+// When the factorisation is run in serial using a CliqueStack, the children are
+// processed in reverse order. To obtain the same result when running the code
+// in parallel, the children must be synced in reverse order, and so spawned in
+// forward order. To make the TreeSplitting work, the groups of supernodes
+// should be formed in reverse order, so that in the group {a,b,d}, d has the
+// lowest ordering and a the highest one. In this way, syncing in reverse order
+// produces the correct result.
 
 }  // namespace hipo
 #endif