Added parallel forward solve

Author: filikat

cmake/sources.cmake

@@ -214,7 +214,6 @@ set(factor_highs_sources
     ipm/hipo/factorhighs/FormatHandler.cpp
     ipm/hipo/factorhighs/HybridHybridFormatHandler.cpp
     ipm/hipo/factorhighs/HybridSolveHandler.cpp
-    ipm/hipo/factorhighs/KrylovMethodsIpm.cpp
     ipm/hipo/factorhighs/Numeric.cpp
     ipm/hipo/factorhighs/SolveHandler.cpp
     ipm/hipo/factorhighs/Swaps.cpp
@@ -233,7 +232,6 @@ set(factor_highs_headers
     ipm/hipo/factorhighs/FormatHandler.h
     ipm/hipo/factorhighs/HybridHybridFormatHandler.h
     ipm/hipo/factorhighs/HybridSolveHandler.h
-    ipm/hipo/factorhighs/KrylovMethodsIpm.h
     ipm/hipo/factorhighs/Numeric.h
     ipm/hipo/factorhighs/ReturnValues.h
     ipm/hipo/factorhighs/SolveHandler.h

highs/ipm/hipo/auxiliary/Auxiliary.cpp

@@ -248,4 +248,24 @@ double Clock::stop() const {
   return d.count();
 }
 
+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.
+
+  // No virtual destructor in TaskGroup. Do not call this class via pointer to
+  // the base!
+
+  cancel();
+
+  // re-call taskWait if it throws, until it succeeds
+  while (true) {
+    try {
+      taskWait();
+      break;
+    } catch (HighsTask::Interrupt) {
+      continue;
+    }
+  }
+}
+
 }  // namespace hipo

highs/ipm/hipo/auxiliary/Auxiliary.h

@@ -8,6 +8,7 @@
 #include <vector>
 
 #include "ipm/hipo/auxiliary/IntConfig.h"
+#include "parallel/HighsParallel.h"
 
 namespace hipo {
 
@@ -66,6 +67,11 @@ class Clock {
   double stop() const;
 };
 
+class TaskGroupSpecial : public highs::parallel::TaskGroup {
+ public:
+  ~TaskGroupSpecial();
+};
+
 }  // namespace hipo
 
 #endif

highs/ipm/hipo/factorhighs/FactorHiGHS.cpp

@@ -48,7 +48,10 @@ Int FHsolver::factorise(const Symbolic& S, const std::vector<Int>& rows,
   return fact_obj.run(N_);
 }
 
-Int FHsolver::solve(std::vector<double>& x) { return N_.solve(x); }
+Int FHsolver::solve(std::vector<double>& x) {
+  N_.setup();
+  return N_.solve(x);
+}
 
 void FHsolver::getRegularisation(std::vector<double>& reg) { N_.getReg(reg); }
 

highs/ipm/hipo/factorhighs/Factorise.cpp

@@ -156,29 +156,6 @@ void Factorise::permute(const std::vector<Int>& iperm) {
   valA_ = std::move(new_val);
 }
 
-class TaskGroupSpecial : public highs::parallel::TaskGroup {
-  // 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!
-
-    cancel();
-
-    // re-call taskWait if it throws, until it succeeds
-    while (true) {
-      try {
-        taskWait();
-        break;
-      } catch (HighsTask::Interrupt) {
-        continue;
-      }
-    }
-  }
-};
-
 void Factorise::processSupernode(Int sn) {
   // Assemble frontal matrix for supernode sn, perform partial factorisation and
   // store the result.

highs/ipm/hipo/factorhighs/HybridSolveHandler.cpp

@@ -12,8 +12,18 @@ namespace hipo {
 HybridSolveHandler::HybridSolveHandler(
     const Symbolic& S, const std::vector<std::vector<double>>& sn_columns,
     const std::vector<std::vector<Int>>& swaps,
-    const std::vector<std::vector<double>>& pivot_2x2)
-    : SolveHandler(S, sn_columns), swaps_{swaps}, pivot_2x2_{pivot_2x2} {}
+    const std::vector<std::vector<double>>& pivot_2x2,
+    const std::vector<Int>& fc, const std::vector<Int>& nc,
+    const std::vector<Int>& fcr, const std::vector<Int>& ncr,
+    std::vector<std::vector<double>>& local)
+    : SolveHandler(S, sn_columns),
+      swaps_{swaps},
+      pivot_2x2_{pivot_2x2},
+      first_child_{fc},
+      next_child_{nc},
+      first_child_reverse_{fcr},
+      next_child_reverse_{ncr},
+      local_{local} {}
 
 void HybridSolveHandler::forwardSolve(std::vector<double>& x) const {
   // Forward solve.
@@ -314,4 +324,212 @@ void HybridSolveHandler::diagSolve(std::vector<double>& x) const {
   }
 }
 
+void HybridSolveHandler::parForwardSolve(std::vector<double>& x) {
+  if (S_.parTree()) {
+    TaskGroupSpecial tg;
+
+    for (Int sn = 0; sn < S_.sn(); ++sn) {
+      if (S_.snParent(sn) == -1) spawnNode(sn, x, tg);
+    }
+
+    tg.taskWait();
+
+  } else {
+    for (Int sn = 0; sn < S_.sn(); ++sn) {
+      processSupernode(sn, x, false);
+    }
+  }
+
+  for (int sn = 0; sn < S_.sn(); ++sn) {
+    const int sn_size = S_.snStart(sn + 1) - S_.snStart(sn);
+    std::memcpy(&x[S_.snStart(sn)], local_[sn].data(),
+                sn_size * sizeof(double));
+  }
+}
+
+void HybridSolveHandler::spawnNode(Int sn, const std::vector<double>& x,
+                                   const TaskGroupSpecial& tg, bool do_spawn) {
+  // if do_spawn is true, a task is actually spawned, otherwise, it is executed
+  // immediately. This avoids the overhead of spawning a task if a supernode has
+  // a single child.
+
+  const NodeData* ptr = S_.nodeDataPtr(sn);
+  if (!ptr) return;
+
+  if (ptr->type == NodeType::single) {
+    // sn is single node; spawn only that
+
+    auto f = [this, &x, sn]() { processSupernode(sn, x, true); };
+
+    if (do_spawn)
+      tg.spawn(std::move(f));
+    else
+      f();
+
+  } else {
+    // sn is head of the first subtree in a group of small subtrees; spawn all
+    // of them
+
+    auto f = [this, &x, ptr]() {
+      for (Int i = 0; i < ptr->group.size(); ++i) {
+        Int st_head = ptr->group[i];
+        Int start = ptr->firstdesc[i];
+        Int end = st_head + 1;
+        for (Int sn = start; sn < end; ++sn) {
+          processSupernode(sn, x, false);
+        }
+      }
+    };
+
+    if (do_spawn)
+      tg.spawn(std::move(f));
+    else
+      f();
+  }
+}
+
+void HybridSolveHandler::syncNode(Int sn, const TaskGroupSpecial& tg) {
+  // If spawnNode(sn,tg) created a task, then sync it.
+  // This happens only if sn is found in the treeSplitting data structure.
+
+  if (S_.nodeDataPtr(sn)) tg.sync();
+}
+
+void HybridSolveHandler::processSupernode(Int sn, const std::vector<double>& x,
+                                          bool parallelise) {
+  // Parallel forward solve.
+  // Blas calls: dtrsv, dgemv
+
+  // supernode columns in format FH
+
+#if HIPO_TIMING_LEVEL >= 2
+  Clock clock;
+#endif
+
+  TaskGroupSpecial tg;
+
+  if (parallelise) {
+    // if there is only one child, do not parallelise
+    if (first_child_[sn] != -1 && next_child_[first_child_[sn]] == -1) {
+      spawnNode(first_child_[sn], x, tg, false);
+      parallelise = false;
+    } else {
+      // spawn children of this supernode in reverse order
+      int child_to_spawn = first_child_reverse_[sn];
+      while (child_to_spawn != -1) {
+        spawnNode(child_to_spawn, x, tg);
+        child_to_spawn = next_child_reverse_[child_to_spawn];
+      }
+    }
+  }
+
+  const Int nb = S_.blockSize();
+
+  // leading size of supernode
+  const Int ldSn = S_.ptr(sn + 1) - S_.ptr(sn);
+
+  // number of columns in the supernode
+  const Int sn_size = S_.snStart(sn + 1) - S_.snStart(sn);
+
+  // first colums of the supernode
+  const Int sn_start = S_.snStart(sn);
+
+  // index to access S->rows for this supernode
+  const Int start_row = S_.ptr(sn);
+
+  // number of blocks of columns
+  const Int n_blocks = (sn_size - 1) / nb + 1;
+
+  // index to access snColumns[sn]
+  Int SnCol_ind{};
+
+  // initialize local storage for this supernode
+  double* local = local_[sn].data();
+  std::memset(local, 0, local_[sn].size() * sizeof(double));
+
+  // contribution from original vector
+  std::memcpy(local, &x[sn_start], sn_size * sizeof(double));
+
+  // contributions from children
+  int child = first_child_[sn];
+  while (child != -1) {
+    if (parallelise) {
+      // wait for child to be ready
+      syncNode(child, tg);
+    }
+
+    std::vector<double>& child_x = local_[child];
+    const int child_size = S_.snStart(child + 1) - S_.snStart(child);
+
+    // number of entries to assemble into local
+    const int nc = child_x.size() - child_size;
+
+#if HIPO_TIMING_LEVEL >= 2
+    clock.start();
+#endif
+    // assemble each contribution of this child
+    for (int i = 0; i < nc; ++i) {
+      const int j = S_.relindClique(child, i);
+      local[j] += child_x[child_size + i];
+    }
+#if HIPO_TIMING_LEVEL >= 2
+    if (data_) data_->sumTime(kTimeSolveSolve_sparse, clock.stop());
+#endif
+
+    child = next_child_[child];
+  }
+
+  // go through blocks of columns for this supernode
+  for (int j = 0; j < n_blocks; ++j) {
+    // number of columns in the block
+    const int jb = std::min(nb, sn_size - nb * j);
+
+    // number of entries in diagonal part
+    const int diag_entries = jb * jb;
+
+    // index to access vector x
+    const int x_start = sn_start + nb * j;
+
+#ifdef HIPO_PIVOTING
+#if HIPO_TIMING_LEVEL >= 2
+    clock.start();
+#endif
+    // apply swaps to portion of rhs that is affected
+    const int* current_swaps = &swaps_[sn][nb * j];
+    permuteWithSwaps(&local[nb * j], current_swaps, jb);
+#if HIPO_TIMING_LEVEL >= 2
+    if (data_) data_->sumTime(kTimeSolveSolve_swap, clock.stop());
+#endif
+#endif
+
+#if HIPO_TIMING_LEVEL >= 2
+    clock.start();
+#endif
+    callAndTime_dtrsv('U', 'T', 'U', jb, &sn_columns_[sn][SnCol_ind], jb,
+                      &local[nb * j], 1, *data_);
+
+    SnCol_ind += diag_entries;
+
+    const int gemv_size = ldSn - nb * j - jb;
+
+    callAndTime_dgemv('T', jb, gemv_size, -1.0, &sn_columns_[sn][SnCol_ind], jb,
+                      &local[nb * j], 1, 1.0, &local[nb * j + jb], 1, *data_);
+    SnCol_ind += jb * gemv_size;
+#if HIPO_TIMING_LEVEL >= 2
+    if (data_) data_->sumTime(kTimeSolveSolve_dense, clock.stop());
+#endif
+
+#ifdef HIPO_PIVOTING
+#if HIPO_TIMING_LEVEL >= 2
+    clock.start();
+#endif
+    // apply inverse swaps
+    permuteWithSwaps(&local[nb * j], current_swaps, jb, true);
+#if HIPO_TIMING_LEVEL >= 2
+    if (data_) data_->sumTime(kTimeSolveSolve_swap, clock.stop());
+#endif
+#endif
+  }
+}
+
 }  // namespace hipo

highs/ipm/hipo/factorhighs/HybridSolveHandler.h

@@ -2,22 +2,40 @@
 #define FACTORHIGHS_HYBRID_SOLVE_HANDLER_H
 
 #include "SolveHandler.h"
+#include "ipm/hipo/auxiliary/Auxiliary.h"
 
 namespace hipo {
 
 class HybridSolveHandler : public SolveHandler {
   const std::vector<std::vector<Int>>& swaps_;
   const std::vector<std::vector<double>>& pivot_2x2_;
 
+  const std::vector<Int>& first_child_;
+  const std::vector<Int>& next_child_;
+  const std::vector<Int>& first_child_reverse_;
+  const std::vector<Int>& next_child_reverse_;
+
+  std::vector<std::vector<double>>& local_;
+
   void forwardSolve(std::vector<double>& x) const override;
   void backwardSolve(std::vector<double>& x) const override;
   void diagSolve(std::vector<double>& x) const override;
 
+  void parForwardSolve(std::vector<double>& x) override;
+
+  void processSupernode(Int sn, const std::vector<double>& x, bool parallelise);
+  void spawnNode(Int sn, const std::vector<double>& x,
+                 const TaskGroupSpecial& tg, bool do_spawn = true);
+  void syncNode(Int sn, const TaskGroupSpecial& tg);
+
  public:
   HybridSolveHandler(const Symbolic& S,
                      const std::vector<std::vector<double>>& sn_columns,
                      const std::vector<std::vector<Int>>& swaps,
-                     const std::vector<std::vector<double>>& pivot_2x2);
+                     const std::vector<std::vector<double>>& pivot_2x2,
+                     const std::vector<Int>& fc, const std::vector<Int>& nc,
+                     const std::vector<Int>& fcr, const std::vector<Int>& ncr,
+                     std::vector<std::vector<double>>& local);
 };
 
 }  // namespace hipo