Persistent basis_update_mpf_t in the Dual Simplex Phase 2

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,72 @@ 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); }
+  if (factorize_basis(A,
+                      settings,
+                      basic_list,
+                      L0_,
+                      U0_,
+                      row_permutation_,
+                      inverse_row_permutation_,
+                      col_permutation_,
+                      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_,
+                        col_permutation_,
+                        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(col_permutation_.size() == A.m);
+  reorder_basic_list(col_permutation_, 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,33 @@ 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),
+      col_permutation_(n),
+      inverse_col_permutation_(n),
+      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,
@@ -205,7 +234,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 +255,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 +278,35 @@ 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);
+    col_permutation_.resize(n);
+    inverse_col_permutation_.resize(n);
+    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,6 +386,13 @@ 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()
   {

cpp/src/dual_simplex/phase2.cpp

@@ -778,7 +778,7 @@ i_t steepest_edge_pricing_with_infeasibilities(const lp_problem_t<i_t, f_t>& lp,
     const i_t j              = infeasibility_indices[k];
     const f_t squared_infeas = squared_infeasibilities[j];
     const f_t val            = squared_infeas / dy_steepest_edge[j];
-    if (val > max_val || val == max_val && j > leaving_index) {
+    if (val > max_val || (val == max_val && j > leaving_index)) {
       max_val                = val;
       leaving_index          = j;
       const f_t lower_infeas = lp.lower[j] - x[j];
@@ -2180,6 +2180,43 @@ dual::status_t dual_phase2(i_t phase,
                            lp_solution_t<i_t, f_t>& sol,
                            i_t& iter,
                            std::vector<f_t>& delta_y_steepest_edge)
+{
+  const i_t m = lp.num_rows;
+  const i_t n = lp.num_cols;
+  std::vector<i_t> basic_list(m);
+  std::vector<i_t> nonbasic_list;
+  std::vector<i_t> superbasic_list;
+  basis_update_mpf_t<i_t, f_t> ft(m, settings.refactor_frequency);
+  const bool initialize_basis = true;
+  return dual_phase2_with_advanced_basis(phase,
+                                         slack_basis,
+                                         initialize_basis,
+                                         start_time,
+                                         lp,
+                                         settings,
+                                         vstatus,
+                                         ft,
+                                         basic_list,
+                                         nonbasic_list,
+                                         sol,
+                                         iter,
+                                         delta_y_steepest_edge);
+}
+
+template <typename i_t, typename f_t>
+dual::status_t dual_phase2_with_advanced_basis(i_t phase,
+                                               i_t slack_basis,
+                                               bool initialize_basis,
+                                               f_t start_time,
+                                               const lp_problem_t<i_t, f_t>& lp,
+                                               const simplex_solver_settings_t<i_t, f_t>& settings,
+                                               std::vector<variable_status_t>& vstatus,
+                                               basis_update_mpf_t<i_t, f_t>& ft,
+                                               std::vector<i_t>& basic_list,
+                                               std::vector<i_t>& nonbasic_list,
+                                               lp_solution_t<i_t, f_t>& sol,
+                                               i_t& iter,
+                                               std::vector<f_t>& delta_y_steepest_edge)
 {
   const i_t m = lp.num_rows;
   const i_t n = lp.num_cols;
@@ -2191,9 +2228,6 @@ dual::status_t dual_phase2(i_t phase,
   assert(lp.lower.size() == n);
   assert(lp.upper.size() == n);
   assert(lp.rhs.size() == m);
-  std::vector<i_t> basic_list(m);
-  std::vector<i_t> nonbasic_list;
-  std::vector<i_t> superbasic_list;
 
   std::vector<f_t>& x = sol.x;
   std::vector<f_t>& y = sol.y;
@@ -2209,33 +2243,18 @@ dual::status_t dual_phase2(i_t phase,
   std::vector<f_t> z_old                     = z;
 
   phase2::bound_info(lp, settings);
-  get_basis_from_vstatus(m, vstatus, basic_list, nonbasic_list, superbasic_list);
-  assert(superbasic_list.size() == 0);
-  assert(nonbasic_list.size() == n - m);
-
-  // Compute L*U = A(p, basic_list)
-  csc_matrix_t<i_t, f_t> L(m, m, 1);
-  csc_matrix_t<i_t, f_t> U(m, m, 1);
-  std::vector<i_t> pinv(m);
-  std::vector<i_t> p;
-  std::vector<i_t> q;
-  std::vector<i_t> deficient;
-  std::vector<i_t> slacks_needed;
-
-  if (factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) ==
-      -1) {
-    settings.log.debug("Initial factorization failed\n");
-    basis_repair(lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
-    if (factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) ==
-        -1) {
+  if (initialize_basis) {
+    std::vector<i_t> superbasic_list;
+    get_basis_from_vstatus(m, vstatus, basic_list, nonbasic_list, superbasic_list);
+    assert(superbasic_list.size() == 0);
+    assert(nonbasic_list.size() == n - m);
+
+    if (ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus) > 0) {
       return dual::status_t::NUMERICAL;
     }
-    settings.log.printf("Basis repaired\n");
+
+    if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; }
   }
-  if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; }
-  assert(q.size() == m);
-  reorder_basic_list(q, basic_list);
-  basis_update_mpf_t<i_t, f_t> ft(L, U, p, settings.refactor_frequency);
 
   std::vector<f_t> c_basic(m);
   for (i_t k = 0; k < m; ++k) {
@@ -2872,54 +2891,28 @@ dual::status_t dual_phase2(i_t phase,
 #endif
     if (should_refactor) {
       bool should_recompute_x = false;
-      if (factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) ==
-          -1) {
+      if (ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus) > 0) {
         should_recompute_x = true;
         settings.log.printf("Failed to factorize basis. Iteration %d\n", iter);
         if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; }
-        basis_repair(lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
         i_t count = 0;
-        while (factorize_basis(
-                 lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) == -1) {
+        i_t deficient_size;
+        while ((deficient_size =
+                  ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus)) > 0) {
           settings.log.printf("Failed to repair basis. Iteration %d. %d deficient columns.\n",
                               iter,
-                              static_cast<int>(deficient.size()));
-#ifdef CHECK_L_FACTOR
-          if (L.check_matrix() == -1) { settings.log.printf("Bad L after basis repair\n"); }
-#endif
+                              static_cast<int>(deficient_size));
+
           if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; }
           settings.threshold_partial_pivoting_tol = 1.0;
+
           count++;
           if (count > 10) { return dual::status_t::NUMERICAL; }
-          basis_repair(
-            lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
-
-#ifdef CHECK_BASIS_REPAIR
-          csc_matrix_t<i_t, f_t> B(m, m, 0);
-          form_b(lp.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
         }
 
         settings.log.printf("Successfully repaired basis. Iteration %d\n", iter);
       }
-      reorder_basic_list(q, basic_list);
-#ifdef CHECK_L_FACTOR
-      if (L.check_matrix() == -1) { settings.log.printf("Bad L factor\n"); }
-#endif
-      ft.reset(L, U, p);
+
       phase2::reset_basis_mark(basic_list, nonbasic_list, basic_mark, nonbasic_mark);
       if (should_recompute_x) {
         std::vector<f_t> unperturbed_x(n);
@@ -3013,6 +3006,20 @@ template dual::status_t dual_phase2<int, double>(
   int& iter,
   std::vector<double>& steepest_edge_norms);
 
+template dual::status_t dual_phase2_with_advanced_basis<int, double>(
+  int phase,
+  int slack_basis,
+  bool initialize_basis,
+  double start_time,
+  const lp_problem_t<int, double>& lp,
+  const simplex_solver_settings_t<int, double>& settings,
+  std::vector<variable_status_t>& vstatus,
+  basis_update_mpf_t<int, double>& ft,
+  std::vector<int>& basic_list,
+  std::vector<int>& nonbasic_list,
+  lp_solution_t<int, double>& sol,
+  int& iter,
+  std::vector<double>& steepest_edge_norms);
 #endif
 
 }  // namespace cuopt::linear_programming::dual_simplex