Numeric is now only used within FHsolver, which now has a solve function. A bit cleaner

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_, sn_columns_);
-  return fact_obj.run(N);
+  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,6 +60,7 @@ class FHsolver {
   const Log* log_;
   DataCollector data_;
   Regul regul_;
+  Numeric N_;
 
   const Int nb_;  // block size
   static const Int default_nb_ = 128;
@@ -88,9 +85,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

@@ -388,6 +388,7 @@ bool Factorise::run(Numeric& num) {
   if (flag_stop_) return true;
 
   // move factorisation to numerical object
+  num.S_ = &S_;
   num.sn_columns_ = &sn_columns_;
   num.total_reg_ = std::move(total_reg_);
   num.swaps_ = std::move(swaps_);

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,15 +13,14 @@
 
 namespace hipo {
 
-Numeric::Numeric(const Symbolic& S) : S_{S} {}
-
-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
 
   assert(sn_columns_);
 
   // initialise solve handler
-  SH_.reset(new HybridSolveHandler(S_, *sn_columns_, swaps_, pivot_2x2_));
+  SH_.reset(new HybridSolveHandler(*S_, *sn_columns_, swaps_, pivot_2x2_));
 
   SH_->setData(data_);
 
@@ -33,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);
@@ -60,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());
@@ -70,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,
@@ -234,7 +237,7 @@ double Numeric::computeOmega(const std::vector<double>& b,
 void Numeric::conditionNumber() const {
   HighsRandom random;
 
-  const Int n = S_.size();
+  const Int n = S_->size();
 
   // estimate largest eigenvalue with power iteration:
   // x <- x / ||x||

highs/ipm/hipo/factorhighs/Numeric.h

@@ -22,7 +22,7 @@ class Numeric {
   std::vector<std::vector<double>> pivot_2x2_{};
 
   // symbolic object
-  const Symbolic& S_;
+  const Symbolic* S_;
 
   // object to handle solve phase in different formats
   mutable std::unique_ptr<SolveHandler> SH_;
@@ -43,11 +43,10 @@ class Numeric {
   // dynamic regularisation applied to the matrix
   std::vector<double> total_reg_{};
 
-  Numeric(const Symbolic& S);
-
   // Full solve with refinement
-  // Return pair (number of solves, final residual)
-  std::pair<Int, double> solve(std::vector<double>& x) const;
+  // Return also number of solves and final residual
+  Int solve(std::vector<double>& x, Int* solve_count = nullptr,
+            double* omega = nullptr) const;
 
   // Iterative refinement
   std::pair<Int, double> refine(const std::vector<double>& rhs,

highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp

@@ -13,7 +13,6 @@ FactorHiGHSSolver::FactorHiGHSSolver(Options& options, const Model& model,
                                      IpmData* record, const LogHighs& log)
     : FH_(&log, options.block_size),
       S_{},
-      N_(S_),
       regul_{regul},
       info_{info},
       data_{record},
@@ -298,7 +297,7 @@ Int FactorHiGHSSolver::factorAS(const HighsSparseMatrix& A,
 
   // factorise matrix
   clock.start();
-  if (FH_.factorise(N_, S_, rowsLower, ptrLower, valLower))
+  if (FH_.factorise(S_, rowsLower, ptrLower, valLower))
     return kStatusErrorFactorise;
   if (info_) {
     info_->factor_time += clock.stop();
@@ -333,8 +332,7 @@ Int FactorHiGHSSolver::factorNE(const HighsSparseMatrix& A,
   // modify them.
   std::vector<Int> ptrNE(ptrNE_);
   std::vector<Int> rowsNE(rowsNE_);
-  if (FH_.factorise(N_, S_, rowsNE, ptrNE, valNE_))
-    return kStatusErrorFactorise;
+  if (FH_.factorise(S_, rowsNE, ptrNE, valNE_)) return kStatusErrorFactorise;
   if (info_) {
     info_->factor_time += clock.stop();
     info_->factor_number++;
@@ -353,14 +351,16 @@ Int FactorHiGHSSolver::solveNE(const std::vector<double>& rhs,
   lhs = rhs;
 
   Clock clock;
-  auto solve_data = N_.solve(lhs);
+  Int solve_count;
+  double final_res;
+  Int solve_status = FH_.solve(lhs, &solve_count, &final_res);
   if (info_) {
     info_->solve_time += clock.stop();
-    info_->solve_number += solve_data.first;
+    info_->solve_number += solve_count;
   }
   if (data_) {
-    data_->back().num_solves += solve_data.first;
-    data_->back().omega = std::max(data_->back().omega, solve_data.second);
+    data_->back().num_solves += solve_count;
+    data_->back().omega = std::max(data_->back().omega, final_res);
   }
 
   return kStatusOk;
@@ -381,14 +381,16 @@ Int FactorHiGHSSolver::solveAS(const std::vector<double>& rhs_x,
   rhs.insert(rhs.end(), rhs_y.begin(), rhs_y.end());
 
   Clock clock;
-  auto solve_data = N_.solve(rhs);
+  Int solve_count;
+  double final_res;
+  Int solve_status = FH_.solve(rhs, &solve_count, &final_res);
   if (info_) {
     info_->solve_time += clock.stop();
-    info_->solve_number += solve_data.first;
+    info_->solve_number += solve_count;
   }
   if (data_) {
-    data_->back().num_solves += solve_data.first;
-    data_->back().omega = std::max(data_->back().omega, solve_data.second);
+    data_->back().num_solves += solve_count;
+    data_->back().omega = std::max(data_->back().omega, final_res);
   }
 
   // split lhs

highs/ipm/hipo/ipm/FactorHiGHSSolver.h

@@ -19,9 +19,6 @@ class FactorHiGHSSolver : public LinearSolver {
   // symbolic factorisation
   Symbolic S_;
 
-  // numeric factorisation
-  Numeric N_;
-
   // normal equations data
   std::vector<Int> ptrNE_, rowsNE_;
   std::vector<double> valNE_;