Set up tree splitting for solve

Author: filikat

highs/ipm/hipo/auxiliary/Auxiliary.cpp

@@ -223,6 +223,23 @@ double getDiagStart(Int n, Int k, Int nb, Int n_blocks, std::vector<Int>& start,
   return result;
 }
 
+void firstDescendant(const std::vector<Int>& parent, std::vector<Int>& first) {
+  // Given an elimination tree, find the first descendant of each node, i.e.
+  // given a node j, first[j] is the descendant of j with smallest number.
+  // Taken from Tim Davis "Direct Methods for Sparse Linear Systems".
+
+  const Int n = parent.size();
+  first.assign(n, -1);
+
+  for (Int i = 0; i < n; ++i) {
+    Int j = i;
+    while (j != -1 && first[j] == -1) {
+      first[j] = i;
+      j = parent[j];
+    }
+  }
+}
+
 Clock::Clock() { start(); }
 void Clock::start() { t0 = std::chrono::high_resolution_clock::now(); }
 double Clock::stop() const {

highs/ipm/hipo/auxiliary/Auxiliary.h

@@ -29,6 +29,7 @@ void processEdge(Int j, Int i, const std::vector<Int>& first,
                  std::vector<Int>& prevleaf, std::vector<Int>& ancestor);
 double getDiagStart(Int n, Int k, Int nb, Int n_blocks, std::vector<Int>& start,
                     bool triang = false);
+void firstDescendant(const std::vector<Int>& parent, std::vector<Int>& first);
 
 template <typename T>
 void permuteVector(std::vector<T>& v, const std::vector<Int>& perm) {

highs/ipm/hipo/factorhighs/Analyse.cpp

@@ -1226,6 +1226,118 @@ void Analyse::computeBlockStart() {
   }
 }
 
+void Analyse::findTreeSplittingSolve() {
+  // Split the tree into single nodes and subtrees for solve.
+  // The subtrees have at most 1% of total operations. They are grouped together
+  // so that each group of subtrees has enough operations.
+  // The tree is parallelised by creating a task for each single node and a task
+  // for each group of subtrees.
+
+  // compute number of operations for each supernode
+  std::vector<double> sn_ops(sn_count_);
+  double total_ops = 0;
+  for (Int sn = 0; sn < sn_count_; ++sn) {
+    // supernode size
+    const Int sz = sn_start_[sn + 1] - sn_start_[sn];
+
+    // frontal size
+    const Int fr = ptr_sn_[sn + 1] - ptr_sn_[sn];
+
+    // number of operations for this supernode for the solve
+    double ops_to_add = (double)sz * sz / 2 + (double)sz * fr;
+    sn_ops[sn] += ops_to_add;
+    total_ops += ops_to_add;
+
+    // add assembly operations
+    if (sn_parent_[sn] != -1) {
+      const Int ldc = fr - sz;
+      sn_ops[sn_parent_[sn]] += ldc;
+      total_ops += ldc;
+    }
+  }
+
+  // compute number of operations to process each subtree
+  std::vector<double> subtree_ops(sn_count_, 0.0);
+  for (Int sn = 0; sn < sn_count_; ++sn) {
+    subtree_ops[sn] += sn_ops[sn];
+    if (sn_parent_[sn] != -1) {
+      subtree_ops[sn_parent_[sn]] += subtree_ops[sn];
+    }
+  }
+
+  // Find first descendant of each supernode
+  std::vector<Int> first_desc;
+  firstDescendant(sn_parent_, first_desc);
+
+  // linked lists of children
+  std::vector<Int> head, next;
+  childrenLinkedList(sn_parent_, head, next);
+
+  node_data_ptr_.assign(sn_count_, nullptr);
+
+  // Divide the tree into single nodes and subtrees, such that each subtree has
+  // at most small_thresh operations overall. Group subtrees together, so that
+  // groups have enough operations.
+  const double small_thresh = 0.05 * total_ops;
+  for (Int sn = 0; sn < sn_count_; ++sn) {
+    if (subtree_ops[sn] > small_thresh) {
+      // sn is a single node
+      auto res_insert = tree_splitting_.insert({sn, {}});
+      node_data_ptr_[sn] = &res_insert.first->second;
+      res_insert.first->second.type = NodeType::single;
+      num_single_++;
+
+      // The children of this sn are either single nodes or head of subtrees.
+      // Divide the head of subtrees in groups, so that each group has enough
+      // operations. Each group corresponds to one task executed in parallel.
+
+      double current_ops = 0.0;
+      NodeData* current_nodedata = nullptr;
+      Int child = head[sn];
+      while (child != -1) {
+        bool is_small = subtree_ops[child] <= small_thresh;
+
+        if (is_small) {
+          num_subtrees_++;
+
+          if (!current_nodedata) {
+            auto res_insert = tree_splitting_.insert({child, {}});
+            current_nodedata = &res_insert.first->second;
+            node_data_ptr_[child] = current_nodedata;
+            current_nodedata->type = NodeType::subtree;
+            current_ops = 0.0;
+          }
+
+          current_ops += subtree_ops[child];
+          current_nodedata->group.push_back(child);
+          current_nodedata->firstdesc.push_back(first_desc[child]);
+
+          if (current_ops > small_thresh) current_nodedata = nullptr;
+        }
+
+        child = next[child];
+      }
+
+    } else if (sn_parent_[sn] == -1) {
+      // sn is small root: single task with whole subtree
+      auto res_insert = tree_splitting_.insert({sn, {}});
+      node_data_ptr_[sn] = &res_insert.first->second;
+      res_insert.first->second.type = NodeType::subtree;
+      res_insert.first->second.group.push_back(sn);
+      res_insert.first->second.firstdesc.push_back(first_desc[sn]);
+    }
+    /*
+    else if (subtree_ops[sn_parent_[sn]] > small_thresh) {
+      // sn is head of a subtree, processed as part of a group of subtrees
+      continue;
+    } else {
+      // sn is part of a subtree, but not the head
+      continue;
+    }
+    */
+  }
+}
+
 Int Analyse::run(Symbolic& S) {
   // Perform analyse phase and store the result into the symbolic object S.
   // After Run returns, the Analyse object is not valid.
@@ -1300,6 +1412,7 @@ Int Analyse::run(Symbolic& S) {
   computeStorage();
   computeBlockStart();
   computeCriticalPath();
+  findTreeSplittingSolve();
 
   // move relevant stuff into S
   S.n_ = n_;
@@ -1314,6 +1427,8 @@ Int Analyse::run(Symbolic& S) {
   S.serial_storage_ = serial_storage_;
   S.flops_ = dense_ops_;
   S.block_size_ = nb_;
+  S.num_single_ = num_single_;
+  S.num_subtrees_ = num_subtrees_;
 
   // compute largest supernode
   std::vector<Int> sn_size(sn_start_.begin() + 1, sn_start_.end());
@@ -1348,6 +1463,8 @@ Int Analyse::run(Symbolic& S) {
   S.relind_clique_ = std::move(relind_clique_);
   S.consecutive_sums_ = std::move(consecutive_sums_);
   S.clique_block_start_ = std::move(clique_block_start_);
+  S.tree_splitting_solve_ = std::move(tree_splitting_);
+  S.node_data_ptr_ = std::move(node_data_ptr_);
 
 #if HIPO_TIMING_LEVEL >= 1
   data_.sumTime(kTimeAnalyse, clock_total.stop());

highs/ipm/hipo/factorhighs/Analyse.h

@@ -3,6 +3,7 @@
 
 #include <algorithm>
 #include <vector>
+#include <map>
 
 #include "DataCollector.h"
 #include "Symbolic.h"
@@ -77,6 +78,11 @@ class Analyse {
 
   std::vector<std::vector<Int>> clique_block_start_{};
 
+  std::map<Int, NodeData> tree_splitting_;
+  std::vector<NodeData*> node_data_ptr_;
+  Int num_single_{};
+  Int num_subtrees_{};
+
   // block size
   Int nb_{};
 
@@ -102,6 +108,7 @@ class Analyse {
                       double& cl_entries) const;
   void computeCriticalPath();
   void computeBlockStart();
+  void findTreeSplittingSolve();
 
  public:
   // Constructor: matrix must be in lower triangular format

highs/ipm/hipo/factorhighs/Symbolic.cpp

@@ -44,6 +44,9 @@ Int Symbolic::cliqueSize(Int sn) const {
 bool Symbolic::parTree() const { return parallel_tree_; }
 bool Symbolic::parNode() const { return parallel_node_; }
 bool Symbolic::metisNo2hop() const { return metis_no2hop_; }
+const NodeData* Symbolic::nodeDataPtr(Int sn) const {
+  return node_data_ptr_[sn];
+}
 
 const std::vector<Int>& Symbolic::ptr() const { return ptr_; }
 const std::vector<Int>& Symbolic::iperm() const { return iperm_; }

highs/ipm/hipo/factorhighs/Symbolic.h

@@ -1,13 +1,21 @@
 #ifndef FACTORHIGHS_SYMBOLIC_H
 #define FACTORHIGHS_SYMBOLIC_H
 
+#include <map>
 #include <vector>
 
 #include "ipm/hipo/auxiliary/IntConfig.h"
 #include "ipm/hipo/auxiliary/Log.h"
 
 namespace hipo {
 
+enum NodeType { single, subtree };
+struct NodeData {
+  NodeType type;
+  std::vector<Int> firstdesc;
+  std::vector<Int> group;
+};
+
 // Symbolic factorisation object
 class Symbolic {
   // Options for parallelism
@@ -36,6 +44,7 @@ class Symbolic {
   Int sn_size_1_{};
   Int sn_size_10_{};
   Int sn_size_100_{};
+  Int num_single_, num_subtrees_;
 
   // Inverse permutation
   std::vector<Int> iperm_{};
@@ -97,6 +106,20 @@ class Symbolic {
   // Starting position of diagonal blocks for hybrid formats
   std::vector<std::vector<Int>> clique_block_start_{};
 
+  // Information to split the elimination tree for solve. Each entry in
+  // tree_splitting_solve_ corresponds to a task that is executed in parallel.
+  // tree_splitting_solve_ contains pairs (sn, data):
+  // - If data.type is single, then the task processes only the supernode sn.
+  // - If data.type is subtree, then the task processes each subtree rooted at
+  //    data.group[i]. Each subtree requires processing supernodes j,
+  //    data.firstdesc[i] <= j <= data.group[i].
+  std::map<Int, NodeData> tree_splitting_solve_;
+
+  // For each supernode, provide the pointer to the NodeData information, if the
+  // supernode is found in the tree splitting data structure. Otherwise,
+  // contains nullptr. Avoids too many lookups into the map.
+  std::vector<NodeData*> node_data_ptr_;
+
   friend class Analyse;
 
  public:
@@ -125,6 +148,7 @@ class Symbolic {
   bool parTree() const;
   bool parNode() const;
   bool metisNo2hop() const;
+  const NodeData* nodeDataPtr(Int sn) const;
   const std::vector<Int>& ptr() const;
   const std::vector<Int>& iperm() const;
   const std::vector<Int>& snParent() const;