Persistent basis_update_mpf_t in the Dual Simplex Phase 2 (#383)

This PR extract the LU factorisation from the `dual_phase2` to a separated method. It also added a new methods in the dual simplex to allow the `basis_update_mpf_t` to persist between calls. 
This allows the branch-and-bound to pass the `basis_update` from parent to child in the future.

Authors:
  - Nicolas L. Guidotti (https://github.com/nguidotti)
  - Chris Maes (https://github.com/chris-maes)

Approvers:
  - Chris Maes (https://github.com/chris-maes)
  - Akif ÇÖRDÜK (https://github.com/akifcorduk)

URL: #383

Author: nguidotti

cpp/src/dual_simplex/basis_updates.cpp

@@ -15,7 +15,9 @@
  * limitations under the License.
  */
 
+#include <dual_simplex/basis_solves.hpp>
 #include <dual_simplex/basis_updates.hpp>
+#include <dual_simplex/initial_basis.hpp>
 #include <dual_simplex/triangle_solve.hpp>
 
 #include <cmath>
@@ -2046,6 +2048,73 @@ void basis_update_mpf_t<i_t, f_t>::multiply_lu(csc_matrix_t<i_t, f_t>& out) cons
   out.nz_max = B_nz;
 }
 
+template <typename i_t, typename f_t>
+int basis_update_mpf_t<i_t, f_t>::refactor_basis(
+  const csc_matrix_t<i_t, f_t>& A,
+  const simplex_solver_settings_t<i_t, f_t>& settings,
+  std::vector<i_t>& basic_list,
+  std::vector<i_t>& nonbasic_list,
+  std::vector<variable_status_t>& vstatus)
+{
+  std::vector<i_t> deficient;
+  std::vector<i_t> slacks_needed;
+
+  if (L0_.m != A.m) { resize(A.m); }
+  std::vector<i_t> q;
+  if (factorize_basis(A,
+                      settings,
+                      basic_list,
+                      L0_,
+                      U0_,
+                      row_permutation_,
+                      inverse_row_permutation_,
+                      q,
+                      deficient,
+                      slacks_needed) == -1) {
+    settings.log.debug("Initial factorization failed\n");
+    basis_repair(A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
+
+#ifdef CHECK_BASIS_REPAIR
+    const i_t m = A.m;
+    csc_matrix_t<i_t, f_t> B(m, m, 0);
+    form_b(A, basic_list, B);
+    for (i_t k = 0; k < deficient.size(); ++k) {
+      const i_t j         = deficient[k];
+      const i_t col_start = B.col_start[j];
+      const i_t col_end   = B.col_start[j + 1];
+      const i_t col_nz    = col_end - col_start;
+      if (col_nz != 1) { settings.log.printf("Deficient column %d has %d nonzeros\n", j, col_nz); }
+      const i_t i = B.i[col_start];
+      if (i != slacks_needed[k]) {
+        settings.log.printf("Slack %d needed but found %d instead\n", slacks_needed[k], i);
+      }
+    }
+#endif
+
+    if (factorize_basis(A,
+                        settings,
+                        basic_list,
+                        L0_,
+                        U0_,
+                        row_permutation_,
+                        inverse_row_permutation_,
+                        q,
+                        deficient,
+                        slacks_needed) == -1) {
+#ifdef CHECK_L_FACTOR
+      if (L0_.check_matrix() == -1) { settings.log.printf("Bad L after basis repair\n"); }
+#endif
+      return deficient.size();
+    }
+    settings.log.debug("Basis repaired\n");
+  }
+
+  assert(q.size() == A.m);
+  reorder_basic_list(q, basic_list);  // We no longer need q after reordering the basic list
+  reset();
+  return 0;
+}
+
 #ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE
 template class basis_update_t<int, double>;
 template class basis_update_mpf_t<int, double>;

cpp/src/dual_simplex/basis_updates.hpp

@@ -17,6 +17,8 @@
 
 #pragma once
 
+#include <dual_simplex/initial_basis.hpp>
+#include <dual_simplex/simplex_solver_settings.hpp>
 #include <dual_simplex/sparse_matrix.hpp>
 #include <dual_simplex/sparse_vector.hpp>
 #include <dual_simplex/types.hpp>
@@ -176,6 +178,31 @@ class basis_update_t {
 template <typename i_t, typename f_t>
 class basis_update_mpf_t {
  public:
+  basis_update_mpf_t(i_t n, const i_t refactor_frequency)
+    : L0_(n, n, 1),
+      U0_(n, n, 1),
+      row_permutation_(n),
+      inverse_row_permutation_(n),
+      S_(n, 0, 0),
+      xi_workspace_(2 * n, 0),
+      x_workspace_(n, 0.0),
+      U0_transpose_(1, 1, 1),
+      L0_transpose_(1, 1, 1),
+      refactor_frequency_(refactor_frequency),
+      total_sparse_L_transpose_(0),
+      total_dense_L_transpose_(0),
+      total_sparse_L_(0),
+      total_dense_L_(0),
+      total_sparse_U_transpose_(0),
+      total_dense_U_transpose_(0),
+      total_sparse_U_(0),
+      total_dense_U_(0),
+      hypersparse_threshold_(0.05)
+  {
+    clear();
+    reset_stats();
+  }
+
   basis_update_mpf_t(const csc_matrix_t<i_t, f_t>& Linit,
                      const csc_matrix_t<i_t, f_t>& Uinit,
                      const std::vector<i_t>& p,
@@ -185,8 +212,6 @@ class basis_update_mpf_t {
       row_permutation_(p),
       inverse_row_permutation_(p.size()),
       S_(Linit.m, 0, 0),
-      col_permutation_(Linit.m),
-      inverse_col_permutation_(Linit.m),
       xi_workspace_(2 * Linit.m, 0),
       x_workspace_(Linit.m, 0.0),
       U0_transpose_(1, 1, 1),
@@ -205,7 +230,7 @@ class basis_update_mpf_t {
     inverse_permutation(row_permutation_, inverse_row_permutation_);
     clear();
     compute_transposes();
-    reset_stas();
+    reset_stats();
   }
 
   void print_stats() const
@@ -226,7 +251,7 @@ class basis_update_mpf_t {
     // clang-format on
   }
 
-  void reset_stas()
+  void reset_stats()
   {
     num_calls_L_           = 0;
     num_calls_U_           = 0;
@@ -249,10 +274,33 @@ class basis_update_mpf_t {
     inverse_permutation(row_permutation_, inverse_row_permutation_);
     clear();
     compute_transposes();
-    reset_stas();
+    reset_stats();
+    return 0;
+  }
+
+  i_t reset()
+  {
+    clear();
+    compute_transposes();
+    reset_stats();
     return 0;
   }
 
+  void resize(i_t n)
+  {
+    L0_.resize(n, n, 1);
+    U0_.resize(n, n, 1);
+    row_permutation_.resize(n);
+    inverse_row_permutation_.resize(n);
+    S_.resize(n, 0, 0);
+    xi_workspace_.resize(2 * n, 0);
+    x_workspace_.resize(n, 0.0);
+    U0_transpose_.resize(1, 1, 1);
+    L0_transpose_.resize(1, 1, 1);
+    clear();
+    reset_stats();
+  }
+
   f_t estimate_solution_density(f_t rhs_nz, f_t sum, i_t& num_calls, bool& use_hypersparse) const
   {
     num_calls++;
@@ -332,13 +380,18 @@ class basis_update_mpf_t {
 
   void multiply_lu(csc_matrix_t<i_t, f_t>& out) const;
 
+  // Compute L*U = A(p, basic_list)
+  int refactor_basis(const csc_matrix_t<i_t, f_t>& A,
+                     const simplex_solver_settings_t<i_t, f_t>& settings,
+                     std::vector<i_t>& basic_list,
+                     std::vector<i_t>& nonbasic_list,
+                     std::vector<variable_status_t>& vstatus);
+
  private:
   void clear()
   {
     pivot_indices_.clear();
     pivot_indices_.reserve(L0_.m);
-    std::iota(col_permutation_.begin(), col_permutation_.end(), 0);
-    std::iota(inverse_col_permutation_.begin(), inverse_col_permutation_.end(), 0);
     S_.col_start.resize(refactor_frequency_ + 1);
     S_.col_start[0] = 0;
     S_.col_start[1] = 0;
@@ -391,8 +444,6 @@ class basis_update_mpf_t {
   std::vector<i_t> pivot_indices_;  // indicies for rank-1 updates to L
   csc_matrix_t<i_t, f_t> S_;        // stores information about the rank-1 updates to L
   std::vector<f_t> mu_values_;      // stores information about the rank-1 updates to L
-  std::vector<i_t> col_permutation_;          // symmetric permuation q used in U(q, q) represents Q
-  std::vector<i_t> inverse_col_permutation_;  // inverse permutation represents Q'
   mutable std::vector<i_t> xi_workspace_;
   mutable std::vector<f_t> x_workspace_;
   mutable csc_matrix_t<i_t, f_t> U0_transpose_;  // Needed for sparse solves