Split tree in single nodes and subtrees for parallelisation

Author: filikat

highs/ipm/hipo/auxiliary/Auxiliary.cpp

@@ -265,6 +265,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

@@ -37,6 +37,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

@@ -2,6 +2,7 @@
 
 #include <fstream>
 #include <iostream>
+#include <map>
 #include <random>
 #include <stack>
 
@@ -1284,6 +1285,72 @@ void Analyse::computeBlockStart() {
   }
 }
 
+void Analyse::findTreeSplitting() {
+  // Split the tree into single nodes and subtrees.
+  // The subtrees have at most 1% of total operations.
+  // The tree is parallelised by creating a task for each single node and a task
+  // for each subtree.
+
+  // linked lists of children
+  std::vector<Int> head, next;
+  childrenLinkedList(sn_parent_, head, next);
+
+  // compute number of operations for each supernode
+  std::vector<double> sn_ops(sn_count_);
+  double total_ops = dense_ops_;
+  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 (Int i = 0; i < sz; ++i) {
+      sn_ops[sn] += (double)(fr - i - 1) * (fr - i - 1);
+    }
+
+    // add assembly operations times spops_weight to the parent
+    if (sn_parent_[sn] != -1) {
+      const Int ldc = fr - sz;
+      sn_ops[sn_parent_[sn]] += ldc * (ldc + 1) / 2 * spops_weight;
+      total_ops += ldc * (ldc + 1) / 2 * spops_weight;
+    }
+  }
+
+  // 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);
+
+  // Divide the tree into single nodes and subtrees, such that each subtree has
+  // at most small_thresh operations overall.
+  const double small_thresh = small_thresh_coeff * total_ops;
+  for (Int sn = 0; sn < sn_count_; ++sn) {
+    if (subtree_ops[sn] > small_thresh) {
+      // sn is a single node
+      tree_splitting_.insert({sn, {NodeType::single, -1}});
+
+    } else if (sn_parent_[sn] == -1 ||
+               subtree_ops[sn_parent_[sn]] > small_thresh) {
+      // sn is head of a subtree
+      tree_splitting_.insert({sn, {NodeType::subtree, first_desc[sn]}});
+
+    } 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.
@@ -1359,6 +1426,8 @@ Int Analyse::run(Symbolic& S) {
   computeBlockStart();
   computeCriticalPath();
 
+  findTreeSplitting();
+
   // move relevant stuff into S
   S.n_ = n_;
   S.sn_ = sn_count_;
@@ -1406,6 +1475,7 @@ 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_ = std::move(tree_splitting_);
 
 #if HIPO_TIMING_LEVEL >= 1
   data_.sumTime(kTimeAnalyse, clock_total.stop());

highs/ipm/hipo/factorhighs/Analyse.h

@@ -77,6 +77,8 @@ class Analyse {
 
   std::vector<std::vector<Int>> clique_block_start_{};
 
+  std::map<Int, NodeData> tree_splitting_;
+
   // block size
   Int nb_{};
 
@@ -103,6 +105,8 @@ class Analyse {
   void computeCriticalPath();
   void computeBlockStart();
 
+  void findTreeSplitting();
+
  public:
   // Constructor: matrix must be in lower triangular format
   Analyse(const std::vector<Int>& rows, const std::vector<Int>& ptr,

highs/ipm/hipo/factorhighs/FactorHiGHSSettings.h

@@ -51,6 +51,10 @@ const double kDynamicDiagCoeff = 1e-24;
 const Int kMaxRefinementIter = 3;
 const double kRefinementTolerance = 1e-12;
 
+// tree splitting
+const double small_thresh_coeff = 0.01;
+const double spops_weight = 30.0;
+
 struct Regul {
   double primal{};
   double dual{};

highs/ipm/hipo/factorhighs/Factorise.cpp

@@ -11,7 +11,6 @@
 #include "SymScaling.h"
 #include "ipm/hipo/auxiliary/Auxiliary.h"
 #include "ipm/hipo/auxiliary/Log.h"
-#include "parallel/HighsParallel.h"
 
 namespace hipo {
 
@@ -168,42 +167,40 @@ void Factorise::permute(const std::vector<Int>& iperm) {
   valA_ = std::move(new_val);
 }
 
-class TaskGroupSpecial : public highs::parallel::TaskGroup {
+TaskGroupSpecial::~TaskGroupSpecial() {
   // Using TaskGroup may throw an exception when tasks are cancelled. Not sure
   // exactly why this happens, but for now this fix seems to work.
 
- public:
-  ~TaskGroupSpecial() {
-    // No virtual destructor in TaskGroup. Do not call this class via pointer to
-    // the base!
+  // No virtual destructor in TaskGroup. Do not call this class via pointer to
+  // the base!
 
-    cancel();
+  cancel();
 
-    // re-call taskWait if it throws, until it succeeds
-    while (true) {
-      try {
-        taskWait();
-        break;
-      } catch (HighsTask::Interrupt) {
-        continue;
-      }
+  // re-call taskWait if it throws, until it succeeds
+  while (true) {
+    try {
+      taskWait();
+      break;
+    } catch (HighsTask::Interrupt) {
+      continue;
     }
   }
-};
+}
 
-void Factorise::processSupernode(Int sn) {
+void Factorise::processSupernode(Int sn, bool parallelise) {
   // Assemble frontal matrix for supernode sn, perform partial factorisation and
   // store the result.
 
   TaskGroupSpecial tg;
 
   if (flag_stop_) return;
 
-  if (S_.parTree()) {
+  if (parallelise) {
     // spawn children of this supernode in reverse order
     Int child_to_spawn = first_child_reverse_[sn];
     while (child_to_spawn != -1) {
-      tg.spawn([=]() { processSupernode(child_to_spawn); });
+      if (spawnNode(child_to_spawn, tg)) return;
+
       child_to_spawn = next_child_reverse_[child_to_spawn];
     }
 
@@ -263,7 +260,7 @@ void Factorise::processSupernode(Int sn) {
     // Schur contribution of the current child
     std::vector<double>& child_clique = schur_contribution_[child_sn];
 
-    if (S_.parTree()) {
+    if (parallelise) {
       // sync with spawned child, apart from the first one
       if (child_sn != first_child_[sn]) tg.sync();
 
@@ -376,6 +373,35 @@ void Factorise::processSupernode(Int sn) {
 #endif
 }
 
+void Factorise::processSubtree(Int start, Int end) {
+  for (Int sn = start; sn < end; ++sn) {
+    processSupernode(sn, false);
+  }
+}
+
+bool Factorise::spawnNode(Int sn, const TaskGroupSpecial& tg) {
+  auto it = S_.treeSplitting().find(sn);
+
+  if (it == S_.treeSplitting().end()) {
+    log_->printDevInfo("Missing supernode from tree splitting\n");
+    flag_stop_ = true;
+    return true;
+  }
+
+  if (it->second.type == NodeType::single) {
+    // sn is single node, spawn only that
+    tg.spawn([=]() { processSupernode(sn, true); });
+
+  } else {
+    // sn is subtree, spawn the whole subtree
+    Int start = it->second.first;
+    Int end = sn + 1;
+    tg.spawn([=]() { processSubtree(start, end); });
+  }
+
+  return false;
+}
+
 bool Factorise::run(Numeric& num) {
 #if HIPO_TIMING_LEVEL >= 1
   Clock clock;
@@ -395,12 +421,10 @@ bool Factorise::run(Numeric& num) {
   sn_columns_.resize(S_.sn());
 
   if (S_.parTree()) {
-    Int spawned_roots{};
     // spawn tasks for root supernodes
     for (Int sn = 0; sn < S_.sn(); ++sn) {
       if (S_.snParent(sn) == -1) {
-        tg.spawn([=]() { processSupernode(sn); });
-        ++spawned_roots;
+        if (spawnNode(sn, tg)) return true;
       }
     }
 
@@ -409,7 +433,7 @@ bool Factorise::run(Numeric& num) {
   } else {
     // go through each supernode serially
     for (Int sn = 0; sn < S_.sn(); ++sn) {
-      processSupernode(sn);
+      processSupernode(sn, false);
     }
   }
 

highs/ipm/hipo/factorhighs/Factorise.h

@@ -7,9 +7,15 @@
 #include "Symbolic.h"
 #include "ipm/hipo/auxiliary/IntConfig.h"
 #include "ipm/hipo/auxiliary/Log.h"
+#include "parallel/HighsParallel.h"
 
 namespace hipo {
 
+class TaskGroupSpecial : public highs::parallel::TaskGroup {
+ public:
+  ~TaskGroupSpecial();
+};
+
 class Factorise {
  public:
   // matrix to factorise
@@ -68,7 +74,9 @@ class Factorise {
 
  public:
   void permute(const std::vector<Int>& iperm);
-  void processSupernode(Int sn);
+  void processSupernode(Int sn, bool parallelise);
+  void processSubtree(Int start, Int end);
+  bool spawnNode(Int sn, const TaskGroupSpecial& tg);
 
  public:
   Factorise(const Symbolic& S, const std::vector<Int>& rowsA,

highs/ipm/hipo/factorhighs/Symbolic.cpp

@@ -44,6 +44,9 @@ const std::vector<Int>& Symbolic::iperm() const { return iperm_; }
 const std::vector<Int>& Symbolic::snParent() const { return sn_parent_; }
 const std::vector<Int>& Symbolic::snStart() const { return sn_start_; }
 const std::vector<Int>& Symbolic::pivotSign() const { return pivot_sign_; }
+const std::map<Int, NodeData>& Symbolic::treeSplitting() const {
+  return tree_splitting_;
+}
 
 std::string memoryString(double mem) {
   std::stringstream ss;
@@ -74,6 +77,8 @@ void Symbolic::print(const Log& log, bool verbose) const {
     log_stream << textline("Critical ops:") << sci(critops_, 0, 1) << '\n';
     log_stream << textline("Max tree speedup:") << fix(flops_ / critops_, 0, 2)
                << '\n';
+    log_stream << textline("Number of tasks:") << integer(tree_splitting_.size())
+               << '\n';
     log_stream << textline("Artificial nz:") << sci(artificial_nz_, 0, 1)
                << '\n';
     log_stream << textline("Artificial ops:") << sci(artificial_ops_, 0, 1)

highs/ipm/hipo/factorhighs/Symbolic.h

@@ -1,13 +1,20 @@
 #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;
+  Int first;
+};
+
 // Symbolic factorisation object
 class Symbolic {
   // Options for parallelism
@@ -94,6 +101,8 @@ class Symbolic {
   // Starting position of diagonal blocks for hybrid formats
   std::vector<std::vector<Int>> clique_block_start_{};
 
+  std::map<Int, NodeData> tree_splitting_;
+
   friend class Analyse;
 
  public:
@@ -124,6 +133,7 @@ class Symbolic {
   const std::vector<Int>& snParent() const;
   const std::vector<Int>& snStart() const;
   const std::vector<Int>& pivotSign() const;
+  const std::map<Int, NodeData>& treeSplitting() const;
 
   void print(const Log& log, bool verbose = false) const;
 };