Merge pull request #2642 from ERGO-Code/hipo

Better refinement for HiPO

Author: jajhall

cmake/sources.cmake

@@ -184,6 +184,7 @@ set(hipo_sources
     ipm/hipo/ipm/Iterate.cpp 
     ipm/hipo/ipm/LogHighs.cpp
     ipm/hipo/ipm/Model.cpp
+    ipm/hipo/ipm/Refine.cpp
     ipm/hipo/ipm/Solver.cpp)
 
 set(hipo_headers

docs/src/installation.md

@@ -39,8 +39,8 @@ git checkout 521-ts
 
 Then build with
 ```
-cmake -S. -B build
--DGKLIB_PATH=/path_to_METIS_repo/GKlib
+cmake -S. -B build \
+-DGKLIB_PATH=/path_to_METIS_repo/GKlib \
 -DCMAKE_INSTALL_PREFIX=path_to_installs_dir
 cmake --build build
 cmake --install build

highs/ipm/hipo/auxiliary/Auxiliary.cpp

@@ -97,48 +97,6 @@ void transpose(const std::vector<Int>& ptr, const std::vector<Int>& rows,
   }
 }
 
-void symProduct(const std::vector<Int>& ptr, const std::vector<Int>& rows,
-                const std::vector<double>& vals, const std::vector<double>& x,
-                std::vector<double>& y, double alpha) {
-  // Matrix-vector product in CSC format, for symmetric matrix which stores only
-  // the lower triangle.
-  // Compute y = y + alpha * M * x
-
-  const Int n = ptr.size() - 1;
-
-  for (Int col = 0; col < n; ++col) {
-    for (Int el = ptr[col]; el < ptr[col + 1]; ++el) {
-      Int row = rows[el];
-      double val = vals[el];
-
-      y[row] += alpha * val * x[col];
-      if (row != col) y[col] += alpha * val * x[row];
-    }
-  }
-}
-
-void symProductQuad(const std::vector<Int>& ptr, const std::vector<Int>& rows,
-                    const std::vector<double>& vals,
-                    const std::vector<double>& x, std::vector<HighsCDouble>& y,
-                    double alpha) {
-  // Matrix-vector product in CSC format, for symmetric matrix which stores only
-  // the lower triangle.
-  // Compute y = y + alpha * M * x
-
-  const Int n = ptr.size() - 1;
-
-  for (Int col = 0; col < n; ++col) {
-    for (Int el = ptr[col]; el < ptr[col + 1]; ++el) {
-      Int row = rows[el];
-      HighsCDouble val = vals[el];
-
-      y[row] += val * (HighsCDouble)x[col] * (HighsCDouble)alpha;
-      if (row != col)
-        y[col] += val * (HighsCDouble)x[row] * (HighsCDouble)alpha;
-    }
-  }
-}
-
 void childrenLinkedList(const std::vector<Int>& parent, std::vector<Int>& head,
                         std::vector<Int>& next) {
   // Create linked lists of children in elimination tree.

highs/ipm/hipo/auxiliary/Auxiliary.h

@@ -8,7 +8,6 @@
 #include <vector>
 
 #include "ipm/hipo/auxiliary/IntConfig.h"
-#include "util/HighsCDouble.h"
 
 namespace hipo {
 
@@ -20,13 +19,6 @@ void transpose(const std::vector<Int>& ptr, const std::vector<Int>& rows,
 void transpose(const std::vector<Int>& ptr, const std::vector<Int>& rows,
                const std::vector<double>& val, std::vector<Int>& ptrT,
                std::vector<Int>& rowsT, std::vector<double>& valT);
-void symProduct(const std::vector<Int>& ptr, const std::vector<Int>& rows,
-                const std::vector<double>& vals, const std::vector<double>& x,
-                std::vector<double>& y, double alpha = 1.0);
-void symProductQuad(const std::vector<Int>& ptr, const std::vector<Int>& rows,
-                    const std::vector<double>& vals,
-                    const std::vector<double>& x, std::vector<HighsCDouble>& y,
-                    double alpha);
 void childrenLinkedList(const std::vector<Int>& parent, std::vector<Int>& head,
                         std::vector<Int>& next);
 void reverseLinkedList(std::vector<Int>& head, std::vector<Int>& next);

highs/ipm/hipo/factorhighs/DataCollector.cpp

@@ -141,7 +141,7 @@ void DataCollector::printTimes(const Log& log) const {
              << "\n";
 
 #if HIPO_TIMING_LEVEL >= 2
-  log_stream << "\tPrepare:                "
+  log_stream << "\tPrepare fact:           "
              << fix(times[kTimeFactorisePrepare], 8, 4) << " ("
              << fix(times[kTimeFactorisePrepare] / times[kTimeFactorise] * 100,
                     4, 1)
@@ -189,7 +189,7 @@ void DataCollector::printTimes(const Log& log) const {
              << "\n";
 
 #if HIPO_TIMING_LEVEL >= 2
-  log_stream << "\tPrepare:                "
+  log_stream << "\tPrepare solve:          "
              << fix(times[kTimeSolvePrepare], 8, 4) << " ("
              << fix(times[kTimeSolvePrepare] / times[kTimeSolve] * 100, 4, 1)
              << "%)\n";
@@ -203,14 +203,6 @@ void DataCollector::printTimes(const Log& log) const {
              << fix(times[kTimeSolveSolve_sparse], 8, 4) << "\n";
   log_stream << "\t\tswap:           " << fix(times[kTimeSolveSolve_swap], 8, 4)
              << "\n";
-  log_stream << "\tResidual:               "
-             << fix(times[kTimeSolveResidual], 8, 4) << " ("
-             << fix(times[kTimeSolveResidual] / times[kTimeSolve] * 100, 4, 1)
-             << "%)\n";
-  log_stream << "\tOmega:                  "
-             << fix(times[kTimeSolveOmega], 8, 4) << " ("
-             << fix(times[kTimeSolveOmega] / times[kTimeSolve] * 100, 4, 1)
-             << "%)\n";
 #endif
   log_stream << "----------------------------------------------------\n";
 
@@ -225,68 +217,57 @@ void DataCollector::printTimes(const Log& log) const {
   log_stream << "BLAS time               \t" << fix(total_blas_time, 8, 4)
              << '\n';
   log_stream << "\tcopy:           \t" << fix(times[kTimeBlas_copy], 8, 4)
-             << " ("
-             << fix(times[kTimeBlas_copy] / times[kTimeSolve] * 100, 4, 1)
+             << " (" << fix(times[kTimeBlas_copy] / total_blas_time * 100, 4, 1)
              << "%) in "
              << integer(blas_calls[kTimeBlas_copy - kTimeBlasStart], 10)
              << " calls\n";
   log_stream << "\taxpy:           \t" << fix(times[kTimeBlas_axpy], 8, 4)
-             << " ("
-             << fix(times[kTimeBlas_axpy] / times[kTimeSolve] * 100, 4, 1)
+             << " (" << fix(times[kTimeBlas_axpy] / total_blas_time * 100, 4, 1)
              << "%) in "
              << integer(blas_calls[kTimeBlas_axpy - kTimeBlasStart], 10)
              << " calls\n";
   log_stream << "\tscal:           \t" << fix(times[kTimeBlas_scal], 8, 4)
-             << " ("
-             << fix(times[kTimeBlas_scal] / times[kTimeSolve] * 100, 4, 1)
+             << " (" << fix(times[kTimeBlas_scal] / total_blas_time * 100, 4, 1)
              << "%) in "
              << integer(blas_calls[kTimeBlas_scal - kTimeBlasStart], 10)
              << " calls\n";
   log_stream << "\tswap:           \t" << fix(times[kTimeBlas_swap], 8, 4)
-             << " ("
-             << fix(times[kTimeBlas_swap] / times[kTimeSolve] * 100, 4, 1)
+             << " (" << fix(times[kTimeBlas_swap] / total_blas_time * 100, 4, 1)
              << "%) in "
              << integer(blas_calls[kTimeBlas_swap - kTimeBlasStart], 10)
              << " calls\n";
   log_stream << "\tgemv:           \t" << fix(times[kTimeBlas_gemv], 8, 4)
-             << " ("
-             << fix(times[kTimeBlas_gemv] / times[kTimeSolve] * 100, 4, 1)
+             << " (" << fix(times[kTimeBlas_gemv] / total_blas_time * 100, 4, 1)
              << "%) in "
              << integer(blas_calls[kTimeBlas_gemv - kTimeBlasStart], 10)
              << " calls\n";
   log_stream << "\ttrsv:           \t" << fix(times[kTimeBlas_trsv], 8, 4)
-             << " ("
-             << fix(times[kTimeBlas_trsv] / times[kTimeSolve] * 100, 4, 1)
+             << " (" << fix(times[kTimeBlas_trsv] / total_blas_time * 100, 4, 1)
              << "%) in "
              << integer(blas_calls[kTimeBlas_trsv - kTimeBlasStart], 10)
              << " calls\n";
   log_stream << "\ttpsv:           \t" << fix(times[kTimeBlas_tpsv], 8, 4)
-             << " ("
-             << fix(times[kTimeBlas_tpsv] / times[kTimeSolve] * 100, 4, 1)
+             << " (" << fix(times[kTimeBlas_tpsv] / total_blas_time * 100, 4, 1)
              << "%) in "
              << integer(blas_calls[kTimeBlas_tpsv - kTimeBlasStart], 10)
              << " calls\n";
   log_stream << "\tger:            \t" << fix(times[kTimeBlas_ger], 8, 4)
-             << " ("
-             << fix(times[kTimeBlas_ger] / times[kTimeSolve] * 100, 4, 1)
+             << " (" << fix(times[kTimeBlas_ger] / total_blas_time * 100, 4, 1)
              << "%) in "
              << integer(blas_calls[kTimeBlas_ger - kTimeBlasStart], 10)
              << " calls\n";
   log_stream << "\ttrsm:           \t" << fix(times[kTimeBlas_trsm], 8, 4)
-             << " ("
-             << fix(times[kTimeBlas_trsm] / times[kTimeSolve] * 100, 4, 1)
+             << " (" << fix(times[kTimeBlas_trsm] / total_blas_time * 100, 4, 1)
              << "%) in "
              << integer(blas_calls[kTimeBlas_trsm - kTimeBlasStart], 10)
              << " calls\n";
   log_stream << "\tsyrk:           \t" << fix(times[kTimeBlas_syrk], 8, 4)
-             << " ("
-             << fix(times[kTimeBlas_syrk] / times[kTimeSolve] * 100, 4, 1)
+             << " (" << fix(times[kTimeBlas_syrk] / total_blas_time * 100, 4, 1)
              << "%) in "
              << integer(blas_calls[kTimeBlas_syrk - kTimeBlasStart], 10)
              << " calls\n";
   log_stream << "\tgemm:           \t" << fix(times[kTimeBlas_gemm], 8, 4)
-             << " ("
-             << fix(times[kTimeBlas_gemm] / times[kTimeSolve] * 100, 4, 1)
+             << " (" << fix(times[kTimeBlas_gemm] / total_blas_time * 100, 4, 1)
              << "%) in "
              << integer(blas_calls[kTimeBlas_gemm - kTimeBlasStart], 10)
              << " calls\n";

highs/ipm/hipo/factorhighs/FactorHiGHS.cpp

@@ -48,8 +48,8 @@ Int FHsolver::factorise(const Symbolic& S, const std::vector<Int>& rows,
   return fact_obj.run(N_);
 }
 
-Int FHsolver::solve(std::vector<double>& x, Int* solve_count, double* omega) {
-  return N_.solve(x, solve_count, omega);
-}
+Int FHsolver::solve(std::vector<double>& x) { return N_.solve(x); }
+
+void FHsolver::getRegularisation(std::vector<double>& reg) { N_.getReg(reg); }
 
 }  // namespace hipo

highs/ipm/hipo/factorhighs/FactorHiGHS.h

@@ -91,13 +91,8 @@ class FHsolver {
                 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);
+  // solution.
+  Int solve(std::vector<double>& x);
 
   // If multiple factorisation are performed, call newIter() before each
   // factorisation. This is used only to collect data for debugging, if
@@ -107,6 +102,8 @@ class FHsolver {
   // Set values for static regularisation to be added when a pivot is selected.
   // If regularisation is already added to the matrix, ignore.
   void setRegularisation(double reg_p, double reg_d);
+
+  void getRegularisation(std::vector<double>& reg);
 };
 
 }  // namespace hipo

highs/ipm/hipo/factorhighs/FactorHiGHSSettings.h

@@ -47,10 +47,6 @@ const Int kBlockParallelThreshold = 5;
 // regularisation
 const double kDynamicDiagCoeff = 1e-24;
 
-// refinement
-const Int kMaxRefinementIter = 3;
-const double kRefinementTolerance = 1e-12;
-
 // metis
 const Int kMetisSeed = 42;
 

highs/ipm/hipo/factorhighs/Factorise.cpp

@@ -74,29 +74,17 @@ Factorise::Factorise(const Symbolic& S, const std::vector<Int>& rowsA,
     min_diag_ = std::min(min_diag_, val);
   }
 
-  // infinity norm of columns of A
-  inf_norm_cols_.assign(n_, 0.0);
-  for (Int col = 0; col < n_; ++col) {
-    for (Int el = ptrA_[col]; el < ptrA_[col + 1]; ++el) {
-      Int row = rowsA_[el];
-      double val = valA_[el];
-      inf_norm_cols_[col] = std::max(inf_norm_cols_[col], std::abs(val));
-      if (row != col)
-        inf_norm_cols_[row] = std::max(inf_norm_cols_[row], std::abs(val));
-    }
-  }
-
   // one norm of columns of A
-  one_norm_cols_.assign(n_, 0.0);
+  std::vector<double> one_norm_cols(n_, 0.0);
   for (Int col = 0; col < n_; ++col) {
     for (Int el = ptrA_[col]; el < ptrA_[col + 1]; ++el) {
       Int row = rowsA_[el];
       double val = valA_[el];
-      one_norm_cols_[col] += std::abs(val);
-      if (row != col) one_norm_cols_[row] += std::abs(val);
+      one_norm_cols[col] += std::abs(val);
+      if (row != col) one_norm_cols[row] += std::abs(val);
     }
   }
-  A_norm1_ = *std::max_element(one_norm_cols_.begin(), one_norm_cols_.end());
+  A_norm1_ = *std::max_element(one_norm_cols.begin(), one_norm_cols.end());
 
   data_.setNorms(A_norm1_, max_diag_);
 }
@@ -421,11 +409,6 @@ bool Factorise::run(Numeric& num) {
   num.total_reg_ = std::move(total_reg_);
   num.swaps_ = std::move(swaps_);
   num.pivot_2x2_ = std::move(pivot_2x2_);
-  num.ptrA_ = std::move(ptrA_);
-  num.rowsA_ = std::move(rowsA_);
-  num.valA_ = std::move(valA_);
-  num.one_norm_cols_ = std::move(one_norm_cols_);
-  num.inf_norm_cols_ = std::move(inf_norm_cols_);
   num.data_ = &data_;
 
 #if HIPO_TIMING_LEVEL >= 1

highs/ipm/hipo/factorhighs/Factorise.h

@@ -51,8 +51,6 @@ class Factorise {
   double max_diag_{};
   double min_diag_{};
   double A_norm1_{};
-  std::vector<double> one_norm_cols_{};
-  std::vector<double> inf_norm_cols_{};
 
   // regularisation
   std::vector<double> total_reg_{};