added a dual simplex method for calculating the phase 2 from an existing basis_update_mpf_t.

nguidotti · nguidotti · commit 181aef4849b9 · 2025-10-07T16:14:13.000+02:00
diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/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,67 @@ 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>::factorize_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 (dual_simplex::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
+    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 (dual_simplex::factorize_basis(A,
+                                      settings,
+                                      basic_list,
+                                      L0_,
+                                      U0_,
+                                      row_permutation_,
+                                      inverse_row_permutation_,
+                                      col_permutation_,
+                                      deficient,
+                                      slacks_needed) == -1) {
+      return deficient.size();
+    }
+    settings.log.printf("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>;
diff --git a/cpp/src/dual_simplex/basis_updates.hpp b/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,7 +278,16 @@ class basis_update_mpf_t {
     inverse_permutation(row_permutation_, inverse_row_permutation_);
     clear();
     compute_transposes();
-    reset_stas();
+    reset_stats();
+    return 0;
+  }
+
+  i_t reset()
+  {
+    inverse_permutation(row_permutation_, inverse_row_permutation_);
+    clear();
+    compute_transposes();
+    reset_stats();
     return 0;
   }
 
@@ -332,6 +370,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 factorize_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()
   {
diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp
@@ -2195,6 +2195,57 @@ dual::status_t dual_phase2(i_t phase,
   std::vector<i_t> nonbasic_list;
   std::vector<i_t> superbasic_list;
 
+  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);
+
+  basis_update_mpf_t<i_t, f_t> ft(m, settings.refactor_frequency);
+
+  if (ft.factorize_basis(lp.A, settings, basic_list, nonbasic_list, vstatus) > 0) {
+    return dual::status_t::NUMERICAL;
+  }
+
+  if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; }
+  return dual_phase2_with_basis_update(phase,
+                                       slack_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_basis_update(i_t phase,
+                                             i_t slack_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;
+  assert(m <= n);
+  assert(vstatus.size() == n);
+  assert(lp.A.m == m);
+  assert(lp.A.n == n);
+  assert(lp.objective.size() == n);
+  assert(lp.lower.size() == n);
+  assert(lp.upper.size() == n);
+  assert(lp.rhs.size() == m);
+
   std::vector<f_t>& x = sol.x;
   std::vector<f_t>& y = sol.y;
   std::vector<f_t>& z = sol.z;
@@ -2208,35 +2259,6 @@ dual::status_t dual_phase2(i_t phase,
   std::vector<variable_status_t> vstatus_old = vstatus;
   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) {
-      return dual::status_t::NUMERICAL;
-    }
-    settings.log.printf("Basis repaired\n");
-  }
-  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) {
     const i_t j = basic_list[k];
@@ -2872,48 +2894,27 @@ 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.factorize_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.factorize_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()));
+                              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);
-      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);
diff --git a/cpp/src/dual_simplex/phase2.hpp b/cpp/src/dual_simplex/phase2.hpp
@@ -17,6 +17,7 @@
 
 #pragma once
 
+#include <dual_simplex/basis_updates.hpp>
 #include <dual_simplex/initial_basis.hpp>
 #include <dual_simplex/logger.hpp>
 #include <dual_simplex/presolve.hpp>
@@ -66,4 +67,17 @@ dual::status_t dual_phase2(i_t phase,
                            i_t& iter,
                            std::vector<f_t>& steepest_edge_norms);
 
+template <typename i_t, typename f_t>
+dual::status_t dual_phase2_with_basis_update(i_t phase,
+                                             i_t slack_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);
 }  // namespace cuopt::linear_programming::dual_simplex
