Merge branch 'hipo-dev' into hipo

Author: filikat

highs/ipm/hipo/factorhighs/FactorHiGHS.cpp

@@ -41,12 +41,15 @@ Int FHsolver::analyse(Symbolic& S, const std::vector<Int>& rows,
   return an_obj.run(S);
 }
 
-Int FHsolver::factorise(Numeric& N, const Symbolic& S,
-                        const std::vector<Int>& rows,
+Int FHsolver::factorise(const Symbolic& S, const std::vector<Int>& rows,
                         const std::vector<Int>& ptr,
                         const std::vector<double>& vals) {
-  Factorise fact_obj(S, rows, ptr, vals, regul_, log_, data_);
-  return fact_obj.run(N);
+  Factorise fact_obj(S, rows, ptr, vals, regul_, log_, data_, sn_columns_);
+  return fact_obj.run(N_);
+}
+
+Int FHsolver::solve(std::vector<double>& x, Int* solve_count, double* omega) {
+  return N_.solve(x, solve_count, omega);
 }
 
 }  // namespace hipo

highs/ipm/hipo/factorhighs/FactorHiGHS.h

@@ -22,7 +22,6 @@ It is stored using three arrays:
 
 The direct solver uses the following objects:
 - Symbolic, to store the symbolic factorization;
-- Numeric, to store the numeric factorization;
 - FHsolver, to perform analyse and factorise phases.
 
 Define a vector signs that contains the expected sign of each pivot (1 or -1).
@@ -32,13 +31,10 @@ M^{-1} * rhs.
 Then, the factorization is performed as follows.
 
     Symbolic S;
-    Numeric N(S);
-
     FHsolver FH;
     FH.analyse(S, rows, ptr, signs);
-    FH.factorise(N, S, rows, ptr, val);
-
-    N.solve(rhs);
+    FH.factorise(S, rows, ptr, val);
+    FH.solve(x);
 
 Printing to screen is achieved using the interface in auxiliary/Log.h. Pass an
 object of type Log for normal printing:
@@ -64,10 +60,16 @@ class FHsolver {
   const Log* log_;
   DataCollector data_;
   Regul regul_;
+  Numeric N_;
 
   const Int nb_;  // block size
   static const Int default_nb_ = 128;
 
+  // Columns of factorisation, stored by supernode.
+  // This memory is allocated the first time that it is used. Subsequent
+  // factorisations reuse the same memory.
+  std::vector<std::vector<double>> sn_columns_;
+
  public:
   // Create object and initialise DataCollector
   FHsolver(const Log* log = nullptr, Int block_size = default_nb_);
@@ -85,9 +87,18 @@ class FHsolver {
   // numerical factorisation in object N. Matrix is moved into the object, so
   // rows, ptr, vals are invalid afterwards.
   // See ReturnValues.h for errors.
-  Int factorise(Numeric& N, const Symbolic& S, const std::vector<Int>& rows,
+  Int factorise(const Symbolic& S, const std::vector<Int>& rows,
                 const std::vector<Int>& ptr, const std::vector<double>& vals);
 
+  // Perform solve phase with rhs given by x, which is overwritten with the
+  // solution. solve_count returns the number of solves performed during the
+  // phase (including refinement). omega returns the final residual after
+  // refinement.
+  // For now refinement is performed automatically as part of solve,
+  // this will change in the future.
+  Int solve(std::vector<double>& x, Int* solve_count = nullptr,
+            double* omega = nullptr);
+
   // If multiple factorisation are performed, call newIter() before each
   // factorisation. This is used only to collect data for debugging, if
   // expensive data collection is turned on at compile time.

highs/ipm/hipo/factorhighs/Factorise.cpp

@@ -18,8 +18,9 @@ namespace hipo {
 Factorise::Factorise(const Symbolic& S, const std::vector<Int>& rowsA,
                      const std::vector<Int>& ptrA,
                      const std::vector<double>& valA, const Regul& regul,
-                     const Log* log, DataCollector& data)
-    : S_{S}, regul_{regul}, log_{log}, data_{data} {
+                     const Log* log, DataCollector& data,
+                     std::vector<std::vector<double>>& sn_columns)
+    : S_{S}, regul_{regul}, log_{log}, data_{data}, sn_columns_{sn_columns} {
   // Input the symmetric matrix to be factorised in CSC format and the symbolic
   // factorisation coming from Analyse.
   // Only the lower triangular part of the matrix is used.
@@ -200,7 +201,7 @@ void Factorise::processSupernode(Int sn) {
   // initialise the format handler
   // this also allocates space for the frontal matrix and schur complement
   std::unique_ptr<FormatHandler> FH(
-      new HybridHybridFormatHandler(S_, sn, regul_, data_));
+      new HybridHybridFormatHandler(S_, sn, regul_, data_, sn_columns_[sn]));
 
 #if HIPO_TIMING_LEVEL >= 2
   data_.sumTime(kTimeFactorisePrepare, clock.stop());
@@ -343,8 +344,8 @@ void Factorise::processSupernode(Int sn) {
   FH->extremeEntries();
 
   // terminate the format handler
-  FH->terminate(sn_columns_[sn], schur_contribution_[sn], total_reg_,
-                swaps_[sn], pivot_2x2_[sn]);
+  FH->terminate(schur_contribution_[sn], total_reg_, swaps_[sn],
+                pivot_2x2_[sn]);
 #if HIPO_TIMING_LEVEL >= 2
   data_.sumTime(kTimeFactoriseTerminate, clock.stop());
 #endif
@@ -357,12 +358,15 @@ bool Factorise::run(Numeric& num) {
 
   total_reg_.assign(n_, 0.0);
 
-  // allocate space for list of generated elements and columns of L
+  // allocate space
   schur_contribution_.resize(S_.sn());
-  sn_columns_.resize(S_.sn());
   swaps_.resize(S_.sn());
   pivot_2x2_.resize(S_.sn());
 
+  // This should actually allocate only the first time, then sn_columns_ reuses
+  // the memory of previous factorisations.
+  sn_columns_.resize(S_.sn());
+
   if (S_.parTree()) {
     Int spawned_roots{};
     // spawn tasks for root supernodes
@@ -387,7 +391,8 @@ bool Factorise::run(Numeric& num) {
   if (flag_stop_) return true;
 
   // move factorisation to numerical object
-  num.sn_columns_ = std::move(sn_columns_);
+  num.S_ = &S_;
+  num.sn_columns_ = &sn_columns_;
   num.total_reg_ = std::move(total_reg_);
   num.swaps_ = std::move(swaps_);
   num.pivot_2x2_ = std::move(pivot_2x2_);

highs/ipm/hipo/factorhighs/Factorise.h

@@ -34,7 +34,9 @@ class Factorise {
   std::vector<std::vector<double>> schur_contribution_{};
 
   // columns of L, stored as dense supernodes
-  std::vector<std::vector<double>> sn_columns_{};
+  // This memory is managed outside of Factorise, so that it can be reused for
+  // all ipm iterations.
+  std::vector<std::vector<double>>& sn_columns_;
 
   // swaps of columns for each supernode, ordered locally within a block
   std::vector<std::vector<Int>> swaps_{};
@@ -71,7 +73,8 @@ class Factorise {
  public:
   Factorise(const Symbolic& S, const std::vector<Int>& rowsA,
             const std::vector<Int>& ptrA, const std::vector<double>& valA,
-            const Regul& regul, const Log* log, DataCollector& data);
+            const Regul& regul, const Log* log, DataCollector& data,
+            std::vector<std::vector<double>>& sn_columns);
 
   bool run(Numeric& num);
 };

highs/ipm/hipo/factorhighs/FormatHandler.cpp

@@ -6,21 +6,22 @@
 
 namespace hipo {
 
-FormatHandler::FormatHandler(const Symbolic& S, Int sn, const Regul& regul)
+FormatHandler::FormatHandler(const Symbolic& S, Int sn, const Regul& regul,
+                             std::vector<double>& frontal)
     : S_{&S},
       regul_{regul},
       sn_{sn},
       nb_{S_->blockSize()},
       sn_size_{S_->snStart(sn_ + 1) - S_->snStart(sn_)},
       ldf_{S_->ptr(sn_ + 1) - S_->ptr(sn_)},
-      ldc_{ldf_ - sn_size_} {
+      ldc_{ldf_ - sn_size_},
+      frontal_{frontal} {
   local_reg_.resize(sn_size_);
   swaps_.resize(sn_size_);
   pivot_2x2_.resize(sn_size_);
 }
 
-void FormatHandler::terminate(std::vector<double>& frontal,
-                              std::vector<double>& clique,
+void FormatHandler::terminate(std::vector<double>& clique,
                               std::vector<double>& total_reg,
                               std::vector<Int>& swaps,
                               std::vector<double>& pivot_2x2) {
@@ -29,7 +30,6 @@ void FormatHandler::terminate(std::vector<double>& frontal,
   // accessed only here, while a local copy is used for the assembly and dense
   // factorisation. This should avoid the problem of false sharing.
 
-  frontal = std::move(frontal_);
   clique = std::move(clique_);
   swaps = std::move(swaps_);
   pivot_2x2 = std::move(pivot_2x2_);

highs/ipm/hipo/factorhighs/FormatHandler.h

@@ -48,17 +48,17 @@ class FormatHandler {
   const Int ldc_{};
 
   // local copies to be moved at the end
-  std::vector<double> frontal_{};
+  std::vector<double>& frontal_;
   std::vector<double> clique_{};
   std::vector<double> local_reg_{};
   std::vector<Int> swaps_{};
   std::vector<double> pivot_2x2_{};
 
  public:
-  FormatHandler(const Symbolic& S, Int sn, const Regul& regul);
-  void terminate(std::vector<double>& frontal, std::vector<double>& clique,
-                 std::vector<double>& total_reg, std::vector<Int>& swaps,
-                 std::vector<double>& pivot_2x2);
+  FormatHandler(const Symbolic& S, Int sn, const Regul& regul,
+                std::vector<double>& frontal);
+  void terminate(std::vector<double>& clique, std::vector<double>& total_reg,
+                 std::vector<Int>& swaps, std::vector<double>& pivot_2x2);
 
   // avoid copies
   FormatHandler(const FormatHandler&) = delete;

highs/ipm/hipo/factorhighs/HybridHybridFormatHandler.cpp

@@ -7,10 +7,10 @@
 
 namespace hipo {
 
-HybridHybridFormatHandler::HybridHybridFormatHandler(const Symbolic& S, Int sn,
-                                                     const Regul& regul,
-                                                     DataCollector& data)
-    : FormatHandler(S, sn, regul), data_{data} {
+HybridHybridFormatHandler::HybridHybridFormatHandler(
+    const Symbolic& S, Int sn, const Regul& regul, DataCollector& data,
+    std::vector<double>& frontal)
+    : FormatHandler(S, sn, regul, frontal), data_{data} {
   // initialise frontal and clique
   initFrontal();
   initClique();
@@ -20,8 +20,11 @@ void HybridHybridFormatHandler::initFrontal() {
   const Int n_blocks = (sn_size_ - 1) / nb_ + 1;
   diag_start_.resize(n_blocks);
   Int frontal_size = getDiagStart(ldf_, sn_size_, nb_, n_blocks, diag_start_);
-  frontal_.resize(frontal_size + extra_space);
+  frontal_.assign(frontal_size + extra_space, 0.0);
   // NB: the plus 10 is not needed, but it avoids weird problems later on.
+
+  // frontal_ is actually allocated just the first time, then the memory is
+  // reused from the previous factorisations and just initialised.
 }
 
 void HybridHybridFormatHandler::initClique() {

highs/ipm/hipo/factorhighs/HybridHybridFormatHandler.h

@@ -23,7 +23,7 @@ class HybridHybridFormatHandler : public FormatHandler {
 
  public:
   HybridHybridFormatHandler(const Symbolic& S, Int sn, const Regul& regul,
-                            DataCollector& data);
+                            DataCollector& data, std::vector<double>& frontal);
 };
 
 }  // namespace hipo

highs/ipm/hipo/factorhighs/Numeric.cpp

@@ -3,6 +3,7 @@
 #include "DataCollector.h"
 #include "FactorHiGHSSettings.h"
 #include "HybridSolveHandler.h"
+#include "ReturnValues.h"
 #include "Timing.h"
 #include "ipm/hipo/auxiliary/Auxiliary.h"
 #include "ipm/hipo/auxiliary/Log.h"
@@ -12,14 +13,15 @@
 
 namespace hipo {
 
-Numeric::Numeric(const Symbolic& S) : S_{S} {
-  // initialise solve handler
-  SH_.reset(new HybridSolveHandler(S_, sn_columns_, swaps_, pivot_2x2_));
-}
-
-std::pair<Int, double> Numeric::solve(std::vector<double>& x) const {
+Int Numeric::solve(std::vector<double>& x, Int* solve_count,
+                   double* omega) const {
   // Return the number of solves performed
 
+  if (!sn_columns_ || !S_) return kRetInvalidPointer;
+
+  // initialise solve handler
+  SH_.reset(new HybridSolveHandler(*S_, *sn_columns_, swaps_, pivot_2x2_));
+
   SH_->setData(data_);
 
 #if HIPO_TIMING_LEVEL >= 1
@@ -31,7 +33,7 @@ std::pair<Int, double> Numeric::solve(std::vector<double>& x) const {
 #endif
 
   // permute rhs
-  permuteVectorInverse(x, S_.iperm());
+  permuteVectorInverse(x, S_->iperm());
 
   // make a copy of permuted rhs, for refinement
   const std::vector<double> rhs(x);
@@ -58,7 +60,7 @@ std::pair<Int, double> Numeric::solve(std::vector<double>& x) const {
 #endif
 
   // unpermute solution
-  permuteVector(x, S_.iperm());
+  permuteVector(x, S_->iperm());
 
 #if HIPO_TIMING_LEVEL >= 2
   if (data_) data_->sumTime(kTimeSolvePrepare, clock_fine.stop());
@@ -68,7 +70,10 @@ std::pair<Int, double> Numeric::solve(std::vector<double>& x) const {
   if (data_) data_->sumTime(kTimeSolve, clock.stop());
 #endif
 
-  return {refine_data.first + 1, refine_data.second};
+  if (solve_count) *solve_count = refine_data.first + 1;
+  if (omega) *omega = refine_data.second;
+
+  return kRetOk;
 }
 
 std::vector<double> Numeric::residual(const std::vector<double>& rhs,
@@ -229,73 +234,4 @@ double Numeric::computeOmega(const std::vector<double>& b,
   return omega_1 + omega_2;
 }
 
-void Numeric::conditionNumber() const {
-  HighsRandom random;
-
-  const Int n = S_.size();
-
-  // estimate largest eigenvalue with power iteration:
-  // x <- x / ||x||
-  // y <- M * x
-  // lambda = ||y||
-
-  double lambda_large{};
-  std::vector<double> x(n);
-  for (Int i = 0; i < n; ++i) x[i] = random.fraction();
-
-  for (Int iter = 0; iter < 10; ++iter) {
-    // normalize x
-    double x_norm = norm2(x);
-    vectorScale(x, 1.0 / x_norm);
-
-    // multiply by matrix
-    std::vector<double> y(n);
-    symProduct(ptrA_, rowsA_, valA_, x, y);
-
-    double norm_y = norm2(y);
-
-    if (std::abs(lambda_large - norm_y) / std::max(1.0, lambda_large) < 1e-2) {
-      // converged
-      break;
-    }
-
-    lambda_large = norm_y;
-    x = std::move(y);
-  }
-
-  // estimate inverse of smallest eigenvalue with inverse power iteration:
-  // x <- x / ||x||
-  // y <- M^-1 * x
-  // lambda_inv = ||y||
-
-  double lambda_small_inv{};
-  for (Int i = 0; i < n; ++i) x[i] = random.fraction();
-
-  for (Int iter = 0; iter < 10; ++iter) {
-    // normalize x
-    double x_norm = norm2(x);
-    vectorScale(x, 1.0 / x_norm);
-
-    // solve with matrix
-    std::vector<double> y(x);
-    SH_->forwardSolve(y);
-    SH_->diagSolve(y);
-    SH_->backwardSolve(y);
-
-    double norm_y = norm2(y);
-
-    if (std::abs(lambda_small_inv - norm_y) / std::max(1.0, lambda_small_inv) <
-        1e-2) {
-      // converged
-      break;
-    }
-
-    lambda_small_inv = norm_y;
-    x = std::move(y);
-  }
-
-  // condition number: lambda_large / lambda_small
-  double cond = lambda_large * lambda_small_inv;
-}
-
 }  // namespace hipo