Added stack of cliques for subtrees

Author: filikat

cmake/sources.cmake

@@ -204,6 +204,7 @@ set(hipo_headers
 set(factor_highs_sources
     ipm/hipo/factorhighs/Analyse.cpp
     ipm/hipo/factorhighs/CallAndTimeBlas.cpp
+    ipm/hipo/factorhighs/CliqueStack.cpp
     ipm/hipo/factorhighs/DataCollector.cpp
     ipm/hipo/factorhighs/DenseFactHybrid.cpp
     ipm/hipo/factorhighs/DenseFactKernel.cpp
@@ -223,6 +224,7 @@ set(factor_highs_sources
 set(factor_highs_headers
     ipm/hipo/factorhighs/Analyse.h
     ipm/hipo/factorhighs/CallAndTimeBlas.h
+    ipm/hipo/factorhighs/CliqueStack.h
     ipm/hipo/factorhighs/DataCollector.h
     ipm/hipo/factorhighs/DenseFact.h
     ipm/hipo/factorhighs/DgemmParallel.h

highs/ipm/hipo/auxiliary/Auxiliary.cpp

@@ -250,8 +250,8 @@ void processEdge(Int j, Int i, const std::vector<Int>& first,
   prevleaf[i] = j;
 }
 
-double getDiagStart(Int n, Int k, Int nb, Int n_blocks, std::vector<Int>& start,
-                    bool triang) {
+int64_t getDiagStart(Int n, Int k, Int nb, Int n_blocks,
+                     std::vector<Int>& start, bool triang) {
   // start position of diagonal blocks for blocked dense formats
   start.assign(n_blocks, 0);
   for (Int i = 1; i < n_blocks; ++i) {
@@ -260,8 +260,9 @@ double getDiagStart(Int n, Int k, Int nb, Int n_blocks, std::vector<Int>& start,
   }
 
   Int jb = std::min(nb, k - (n_blocks - 1) * nb);
-  double result = (double)start.back() + (double)(n - (n_blocks - 1) * nb) * jb;
-  if (triang) result -= (double)jb * (jb - 1) / 2;
+  int64_t result =
+      (int64_t)start.back() + (int64_t)(n - (n_blocks - 1) * nb) * jb;
+  if (triang) result -= (int64_t)jb * (jb - 1) / 2;
   return result;
 }
 

highs/ipm/hipo/auxiliary/Auxiliary.h

@@ -35,8 +35,8 @@ void dfsPostorder(Int node, Int& start, std::vector<Int>& head,
 void processEdge(Int j, Int i, const std::vector<Int>& first,
                  std::vector<Int>& maxfirst, std::vector<Int>& delta,
                  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);
+int64_t 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>

highs/ipm/hipo/factorhighs/Analyse.cpp

@@ -908,8 +908,8 @@ void Analyse::relativeIndClique() {
   }
 }
 
-void Analyse::computeStorage(Int fr, Int sz, double& fr_entries,
-                             double& cl_entries) const {
+void Analyse::computeStorage(Int fr, Int sz, int64_t& fr_entries,
+                             int64_t& cl_entries) const {
   // compute storage required by frontal and clique, based on the format used
 
   const Int cl = fr - sz;
@@ -920,17 +920,17 @@ void Analyse::computeStorage(Int fr, Int sz, double& fr_entries,
 
   // clique is stored as a collection of rectangles
   n_blocks = (cl - 1) / nb_ + 1;
-  double schur_size{};
+  int64_t schur_size{};
   for (Int j = 0; j < n_blocks; ++j) {
     const Int jb = std::min(nb_, cl - j * nb_);
-    schur_size += (double)(cl - j * nb_) * jb;
+    schur_size += (ino64_t)(cl - j * nb_) * jb;
   }
   cl_entries = schur_size;
 }
 
 void Analyse::computeStorage() {
-  std::vector<double> clique_entries(sn_count_);
-  std::vector<double> frontal_entries(sn_count_);
+  std::vector<int64_t> clique_entries(sn_count_);
+  std::vector<int64_t> frontal_entries(sn_count_);
   std::vector<double> storage(sn_count_);
   std::vector<double> storage_factors(sn_count_);
 
@@ -1050,8 +1050,8 @@ void Analyse::computeCriticalPath() {
 }
 
 void Analyse::reorderChildren() {
-  std::vector<double> clique_entries(sn_count_);
-  std::vector<double> frontal_entries(sn_count_);
+  std::vector<int64_t> clique_entries(sn_count_);
+  std::vector<int64_t> frontal_entries(sn_count_);
   std::vector<double> storage(sn_count_);
   std::vector<double> storage_factors(sn_count_);
 
@@ -1304,12 +1304,15 @@ void Analyse::findTreeSplitting() {
             is_in_tree_splitting_[child] = true;
             current_nodedata = &res_insert.first->second;
             current_nodedata->type = NodeType::subtree;
+            current_nodedata->stack_size = 0;
             current_ops = 0.0;
           }
 
           current_ops += subtree_ops[child];
           current_nodedata->group.push_back(child);
           current_nodedata->firstdesc.push_back(first_desc[child]);
+          current_nodedata->stack_size =
+              std::max(current_nodedata->stack_size, stack_subtrees_[child]);
 
           if (current_ops > small_thresh) current_nodedata = nullptr;
         }
@@ -1324,6 +1327,7 @@ void Analyse::findTreeSplitting() {
       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]);
+      res_insert.first->second.stack_size = stack_subtrees_[sn];
     }
     /*
     else if (subtree_ops[sn_parent_[sn]] > small_thresh) {
@@ -1337,6 +1341,62 @@ void Analyse::findTreeSplitting() {
   }
 }
 
+void Analyse::computeStackSize() {
+  // Compute the minimum size of the stack to process each subtree.
+
+  std::vector<int64_t> clique_entries(sn_count_);
+  std::vector<int64_t> frontal_entries(sn_count_);
+  stack_subtrees_.assign(sn_count_, 0);
+
+  // initialise data of supernodes
+  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];
+
+    // compute storage based on format used
+    computeStorage(fr, sz, frontal_entries[sn], clique_entries[sn]);
+  }
+
+  // linked lists of children
+  std::vector<Int> head, next;
+  childrenLinkedList(sn_parent_, head, next);
+
+  // go through the supernodes
+  for (Int sn = 0; sn < sn_count_; ++sn) {
+    // leaf node
+    if (head[sn] == -1) {
+      stack_subtrees_[sn] = clique_entries[sn];
+      continue;
+    }
+
+    // Compute storage
+    // storage is found as max(storage_1,storage_2), where
+    // storage_1 = max_j stack_size[j] + \sum_{k up to j-1} clique_entries[k]
+    // storage_2 = clique_total_entries (including node itself)
+
+    int64_t clique_partial_entries{};
+    int64_t storage_1{};
+
+    Int child = head[sn];
+    while (child != -1) {
+      int64_t current = stack_subtrees_[child] + clique_partial_entries;
+
+      clique_partial_entries += clique_entries[child];
+      storage_1 = std::max(storage_1, current);
+
+      child = next[child];
+    }
+
+    int64_t storage_2 = clique_partial_entries + clique_entries[sn];
+
+    stack_subtrees_[sn] = std::max(storage_1, storage_2);
+    max_stack_size_ = std::max(max_stack_size_, stack_subtrees_[sn]);
+  }
+}
+
 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.
@@ -1411,6 +1471,7 @@ Int Analyse::run(Symbolic& S) {
   computeStorage();
   computeBlockStart();
   computeCriticalPath();
+  computeStackSize();
 
   findTreeSplitting();
 
@@ -1427,6 +1488,7 @@ Int Analyse::run(Symbolic& S) {
   S.serial_storage_ = serial_storage_;
   S.flops_ = dense_ops_;
   S.block_size_ = nb_;
+  S.max_stack_size_ = max_stack_size_;
 
   // compute largest supernode
   std::vector<Int> sn_size(sn_start_.begin() + 1, sn_start_.end());

highs/ipm/hipo/factorhighs/Analyse.h

@@ -80,6 +80,9 @@ class Analyse {
   std::map<Int, NodeData> tree_splitting_;
   std::vector<bool> is_in_tree_splitting_;
 
+  std::vector<int64_t> stack_subtrees_{};
+  int64_t max_stack_size_{};
+
   // block size
   Int nb_{};
 
@@ -101,10 +104,11 @@ class Analyse {
   void relativeIndClique();
   void reorderChildren();
   void computeStorage();
-  void computeStorage(Int fr, Int sz, double& fr_entries,
-                      double& cl_entries) const;
+  void computeStorage(Int fr, Int sz, int64_t& fr_entries,
+                      int64_t& cl_entries) const;
   void computeCriticalPath();
   void computeBlockStart();
+  void computeStackSize();
 
   void findTreeSplitting();
 

highs/ipm/hipo/factorhighs/CliqueStack.cpp

@@ -0,0 +1,75 @@
+#include "CliqueStack.h"
+
+#include <cassert>
+#include <cstring>
+
+namespace hipo {
+
+CliqueStack::CliqueStack(Int stack_size) {
+  stack_.assign(stack_size, 0.0);
+  top_ = 0;
+  workspace_ = nullptr;
+  worksize_ = 0;
+}
+
+double* CliqueStack::setup(Int clique_size) {
+  // Clear workspace
+
+  assert(!workspace_ && !worksize_);
+
+  // This should not trigger reallocation, because the reserve in init is done
+  // with the maximum possible size of a clique.
+  if (top_ + clique_size > stack_.size()) {
+    printf("=== Warning, reallocation of workspace\n");
+    stack_.resize(top_ + clique_size, 0.0);
+  }
+
+  workspace_ = &stack_[top_];
+  worksize_ = clique_size;
+
+  // initialize workspace to zero
+  std::memset(workspace_, 0, worksize_ * sizeof(double));
+
+  return workspace_;
+}
+
+const double* CliqueStack::get() const { return stack_.data(); }
+Int CliqueStack::getTop() const { return top_; }
+
+const double* CliqueStack::getChild(Int& child_sn) const {
+  // Get the top of the stack, in terms of supernode ID of the child and pointer
+  // to its data.
+
+  child_sn = sn_pushed_.top().first;
+  Int child_size = sn_pushed_.top().second;
+  const double* child = &stack_[top_ - child_size];
+
+  return child;
+}
+
+void CliqueStack::popChild() {
+  // Remove top child from the stack
+
+  Int child_size = sn_pushed_.top().second;
+  sn_pushed_.pop();
+
+  top_ -= child_size;
+}
+
+void CliqueStack::pushWork(Int sn) {
+  // Put the content of the workspace at the top of the stack
+
+  // stack_[top_] has lower address than workspace, so no need to resize.
+  // workspace_ and stack_[top_] do not overlap, so use memcpy
+  std::memcpy(&stack_[top_], workspace_, worksize_ * sizeof(double));
+
+  top_ += worksize_;
+
+  // keep track of supernodes pushed
+  sn_pushed_.push({sn, worksize_});
+
+  worksize_ = 0;
+  workspace_ = nullptr;
+}
+
+}  // namespace hipo

highs/ipm/hipo/factorhighs/CliqueStack.h

@@ -0,0 +1,82 @@
+#ifndef FACTORHIGHS_CLIQUE_STACK_H
+#define FACTORHIGHS_CLIQUE_STACK_H
+
+#include <stack>
+#include <vector>
+
+#include "ipm/hipo/auxiliary/IntConfig.h"
+
+// Class to manage the stack of cliques.
+//
+// The stack is used as follows:
+// - use init to initialize the stack, providing the max stack size. If the
+//   parameter is correct, there will be no reallocation of stack during its
+//   operation.
+// - use setup, with the size of the clique that is being computed, to
+//   initialize enough elements in the workspace. This also returns a pointer to
+//   write the new clique. The workspace is simply a chunk of space at the top
+//   of the stack where the new clique is being computed before being pushed in
+//   the right position.
+// - use getChild to obtain information about the child clique that is at the
+//   top of the stack. This also returns a pointer to read the data of the child
+//   clique.
+// - use popChild when the top child clique is no longer needed, to move on to
+//   the next one.
+// - use pushWork, with the ID of supernode that is being processed, when all
+//   children clique have been assembled, the dense factorization is completed,
+//   and the clique in the workspace is ready to be pushed onto the stack.
+
+namespace hipo {
+
+class CliqueStack {
+  std::vector<double> stack_;
+  double* workspace_;
+  Int worksize_;
+
+  // pairs (sn, size) of supernodes that got pushed
+  std::stack<std::pair<Int, Int>> sn_pushed_{};
+
+  Int top_{};
+
+ public:
+  CliqueStack(Int size);
+  double* setup(Int clique_size);
+  const double* getChild(Int& child_sn) const;
+  void popChild();
+  void pushWork(Int sn);
+  const double* get() const;
+  Int getTop() const;
+};
+
+}  // namespace hipo
+
+// Example: sn 7 has children 5,6 and the stack looks like this:
+//
+// | - - - | - - - - - - - | - | - - - - - - - ...
+// |   4   |       5       | 6 | t
+// | - - - | - - - - - - - | - | - - - - - - - ...
+//
+// where t points to the top of the stack.
+// The information about the supernodes pushed (in sn_pushed_) is:
+// (sn 4, size 7), (sn 5, size 15), (sn 6, size 3)
+//
+// The workspace starts at t and fits in the stack.
+// The clique of sn 7 is constructed in the workspace. After a child clique is
+// assembled, it is popped from the stack, by moving t back. After all children
+// are assembled and the dense factorization is done, the stack looks like this:
+//
+// | - - - | - - - - - - - - - | - - - - - - - - - - - | - - ...
+// |   4   | t      free       |           7           |
+// | - - - | - - - - - - - - - | - - - - - - - - - - - | - - ...
+//
+// The clique 7 can then be pushed onto the stack by copying it at position t.
+// At the end, the stack looks like this:
+//
+// | - - - | - - - - - - - - - - - | - - - - - ...
+// |   4   |           7           | t
+// | - - - | - - - - - - - - - - - | - - - - - ...
+//
+// The information in sn_pushed_ is:
+// (sn 4, size 7), (sn 7, size 23)
+
+#endif