Improved comments and removed unused function in Numeric

Author: filikat

highs/ipm/hipo/factorhighs/FactorHiGHS.h

@@ -65,7 +65,9 @@ class FHsolver {
   const Int nb_;  // block size
   static const Int default_nb_ = 128;
 
-  // columns of factorisation, stored by supernode
+  // 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:
@@ -107,34 +109,6 @@ class FHsolver {
   void setRegularisation(double reg_p, double reg_d);
 };
 
-/* To do
-
-- At the moment, the format handler allocates new space for frontal each time,
-  and then moves the values into sn_columns_, with the previous space being
-  deallocated. This is very inefficient, because it causes many memory
-  allocation/deallocation. It may also cause fragmentation.
-- Have sn_columns_ outside of factorise, into FHsolver. Then, pass it into
-  Factorise and assemble the contributions directly in place. Numeric then
-  receives a pointer to the data, so that it can access the data in the same way
-  as now. Format handler needs to assign zeros during initialization each time.
-- In this way, the space for the columns of the factor is persistent throughout
-  the ipm iterations, no extra allocation/deallocation is needed, no copies are
-  needed. Potentially, there may be more contention of the cache lines, but this
-  is to be seen.
-
-- For the moment, the same can be done with cliques. No need to
-  allocate/deallocate them. They can stay allocated and be reused. This would
-  probably considerably increase the memory usage though and the stack approach
-  is to be preferred. However, the latter requires a different parallelisation.
-
-
-- at the moment, I managed to keep the same interface to Numeric, but it's a bit
-  ugly. Maybe Numeric should be an object that is crested by Factorise, which
-  only contains pointers to data. Then I can just copy into the general numeric
-  object within FactorHiGHSSolver and include the pointers to the most up to
-  date data.
-*/
-
 }  // namespace hipo
 
 #endif

highs/ipm/hipo/factorhighs/Factorise.cpp

@@ -358,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

highs/ipm/hipo/factorhighs/Factorise.h

@@ -34,6 +34,8 @@ class Factorise {
   std::vector<std::vector<double>> schur_contribution_{};
 
   // columns of L, stored as dense supernodes
+  // 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

highs/ipm/hipo/factorhighs/HybridHybridFormatHandler.cpp

@@ -22,6 +22,9 @@ void HybridHybridFormatHandler::initFrontal() {
   Int frontal_size = getDiagStart(ldf_, sn_size_, nb_, n_blocks, diag_start_);
   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/Numeric.cpp

@@ -234,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

highs/ipm/hipo/factorhighs/Numeric.h

@@ -9,6 +9,10 @@
 #include "Symbolic.h"
 #include "ipm/hipo/auxiliary/IntConfig.h"
 
+// Numeric allows to perform solve, though a pointer to the numerical factor,
+// that is stored in FHsolver. It also holds auxiliary data about swaps,
+// pivots...
+
 namespace hipo {
 
 class Numeric {
@@ -58,8 +62,6 @@ class Numeric {
   double computeOmega(const std::vector<double>& b,
                       const std::vector<double>& x,
                       const std::vector<double>& res) const;
-
-  void conditionNumber() const;
 };
 
 }  // namespace hipo