From 181aef4849b915948b4d06b25114473e863094ce Mon Sep 17 00:00:00 2001
From: nicolas <nguidotti@nvidia.com>
Date: Tue, 9 Sep 2025 15:25:22 +0200
Subject: [PATCH 01/11] added a dual simplex method for calculating the phase 2
 from an existing basis_update_mpf_t.

---
 cpp/src/dual_simplex/basis_updates.cpp |  63 ++++++++++++++
 cpp/src/dual_simplex/basis_updates.hpp |  51 ++++++++++-
 cpp/src/dual_simplex/phase2.cpp        | 115 +++++++++++++------------
 cpp/src/dual_simplex/phase2.hpp        |  14 +++
 4 files changed, 183 insertions(+), 60 deletions(-)

diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp
index 18f74ccc2..0dba9910d 100644
--- 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
index b12950abb..3198ccd4b 100644
--- 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
index 6e468f915..f83dfe4a5 100644
--- 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
index 39311e060..0b5045311 100644
--- 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

From 0cb0dcf73c8fb3eb185080088929f3729d722f1a Mon Sep 17 00:00:00 2001
From: nicolas <nguidotti@nvidia.com>
Date: Thu, 9 Oct 2025 14:03:18 +0200
Subject: [PATCH 02/11] changed dual simplex to keep the `basis_update_mpf_t`
 for future use

---
 cpp/src/dual_simplex/basis_solves.cpp  |   8 +-
 cpp/src/dual_simplex/basis_updates.cpp |   2 +-
 cpp/src/dual_simplex/phase2.cpp        |  61 +++++++++--
 cpp/src/dual_simplex/phase2.hpp        |  10 ++
 cpp/src/dual_simplex/solve.cpp         | 140 +++++++++++++++++++++----
 cpp/src/dual_simplex/solve.hpp         |  15 +++
 6 files changed, 202 insertions(+), 34 deletions(-)

diff --git a/cpp/src/dual_simplex/basis_solves.cpp b/cpp/src/dual_simplex/basis_solves.cpp
index 2120152c5..6aac8a936 100644
--- a/cpp/src/dual_simplex/basis_solves.cpp
+++ b/cpp/src/dual_simplex/basis_solves.cpp
@@ -46,6 +46,13 @@ void get_basis_from_vstatus(i_t m,
   i_t n             = vstatus.size();
   i_t num_basic     = 0;
   i_t num_non_basic = 0;
+
+  assert(n >= m);
+  nonbasic_list.clear();
+  superbasic_list.clear();
+  nonbasic_list.reserve(n - m);
+  basis_list.resize(m);
+
   for (i_t j = 0; j < n; ++j) {
     if (vstatus[j] == variable_status_t::BASIC) {
       basis_list[num_basic++] = j;
@@ -61,7 +68,6 @@ void get_basis_from_vstatus(i_t m,
       superbasic_list.push_back(j);
     }
   }
-  i_t num_super_basic = superbasic_list.size();
   assert(num_basic == m);
 }
 
diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp
index 0dba9910d..f97413ff7 100644
--- a/cpp/src/dual_simplex/basis_updates.cpp
+++ b/cpp/src/dual_simplex/basis_updates.cpp
@@ -2100,7 +2100,7 @@ int basis_update_mpf_t<i_t, f_t>::factorize_basis(
                                       slacks_needed) == -1) {
       return deficient.size();
     }
-    settings.log.printf("Basis repaired\n");
+    settings.log.debug("Basis repaired\n");
   }
 
   assert(col_permutation_.size() == A.m);
diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp
index f83dfe4a5..2efac1246 100644
--- a/cpp/src/dual_simplex/phase2.cpp
+++ b/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];
@@ -2169,6 +2169,33 @@ class phase2_timers_t {
 };
 
 }  // namespace phase2
+template <typename i_t, typename f_t>
+dual::status_t factorize_basis(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,
+                               f_t start_time)
+{
+  const i_t m = lp.num_rows;
+  const i_t n = lp.num_cols;
+  assert(m <= n);
+  assert(vstatus.size() == n);
+
+  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.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::status_t::UNSET;
+}
 
 template <typename i_t, typename f_t>
 dual::status_t dual_phase2(i_t phase,
@@ -2183,14 +2210,6 @@ dual::status_t dual_phase2(i_t phase,
 {
   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<i_t> basic_list(m);
   std::vector<i_t> nonbasic_list;
   std::vector<i_t> superbasic_list;
@@ -2259,6 +2278,8 @@ dual::status_t dual_phase2_with_basis_update(i_t phase,
   std::vector<variable_status_t> vstatus_old = vstatus;
   std::vector<f_t> z_old                     = z;
 
+  phase2::bound_info(lp, settings);
+
   std::vector<f_t> c_basic(m);
   for (i_t k = 0; k < m; ++k) {
     const i_t j = basic_list[k];
@@ -3008,6 +3029,28 @@ template dual::status_t dual_phase2<int, double>(
   int& iter,
   std::vector<double>& steepest_edge_norms);
 
+template dual::status_t dual_phase2_with_basis_update<int, double>(
+  int phase,
+  int slack_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);
+
+template dual::status_t factorize_basis<int, double>(
+  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,
+  double start_time);
 #endif
 
 }  // namespace cuopt::linear_programming::dual_simplex
diff --git a/cpp/src/dual_simplex/phase2.hpp b/cpp/src/dual_simplex/phase2.hpp
index 0b5045311..bd1a80389 100644
--- a/cpp/src/dual_simplex/phase2.hpp
+++ b/cpp/src/dual_simplex/phase2.hpp
@@ -56,6 +56,15 @@ static std::string status_to_string(status_t status)
 }
 }  // namespace dual
 
+template <typename i_t, typename f_t>
+dual::status_t factorize_basis(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,
+                               f_t start_time);
+
 template <typename i_t, typename f_t>
 dual::status_t dual_phase2(i_t phase,
                            i_t slack_basis,
@@ -80,4 +89,5 @@ dual::status_t dual_phase2_with_basis_update(i_t phase,
                                              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
diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp
index dadb60cc6..4da964474 100644
--- a/cpp/src/dual_simplex/solve.cpp
+++ b/cpp/src/dual_simplex/solve.cpp
@@ -116,6 +116,35 @@ lp_status_t solve_linear_program_advanced(const lp_problem_t<i_t, f_t>& original
                                           lp_solution_t<i_t, f_t>& original_solution,
                                           std::vector<variable_status_t>& vstatus,
                                           std::vector<f_t>& edge_norms)
+{
+  const i_t m = original_lp.num_rows;
+  const i_t n = original_lp.num_cols;
+  assert(m <= n);
+  std::vector<i_t> basic_list(m);
+  std::vector<i_t> nonbasic_list;
+  basis_update_mpf_t<i_t, f_t> ft(m, settings.refactor_frequency);
+  return solve_linear_program_with_basis_update(original_lp,
+                                                start_time,
+                                                settings,
+                                                original_solution,
+                                                ft,
+                                                basic_list,
+                                                nonbasic_list,
+                                                vstatus,
+                                                edge_norms);
+}
+
+template <typename i_t, typename f_t>
+lp_status_t solve_linear_program_with_basis_update(
+  const lp_problem_t<i_t, f_t>& original_lp,
+  const f_t start_time,
+  const simplex_solver_settings_t<i_t, f_t>& settings,
+  lp_solution_t<i_t, f_t>& original_solution,
+  basis_update_mpf_t<i_t, f_t>& ft,
+  std::vector<i_t>& basic_list,
+  std::vector<i_t>& nonbasic_list,
+  std::vector<variable_status_t>& vstatus,
+  std::vector<f_t>& edge_norms)
 {
   lp_status_t lp_status = lp_status_t::UNSET;
   lp_problem_t<i_t, f_t> presolved_lp(original_lp.handle_ptr, 1, 1, 1);
@@ -138,21 +167,21 @@ lp_status_t solve_linear_program_advanced(const lp_problem_t<i_t, f_t>& original
   column_scaling(presolved_lp, settings, lp, column_scales);
   assert(presolved_lp.num_cols == lp.num_cols);
   lp_problem_t<i_t, f_t> phase1_problem(original_lp.handle_ptr, 1, 1, 1);
-  std::vector<variable_status_t> phase1_vstatus;
+
   f_t phase1_obj = -inf;
   create_phase1_problem(lp, phase1_problem);
   assert(phase1_problem.num_cols == presolved_lp.num_cols);
 
   // Set the vstatus for the phase1 problem based on a slack basis
-  phase1_vstatus.resize(phase1_problem.num_cols);
-  std::fill(phase1_vstatus.begin(), phase1_vstatus.end(), variable_status_t::NONBASIC_LOWER);
+  vstatus.resize(phase1_problem.num_cols);
+  std::fill(vstatus.begin(), vstatus.end(), variable_status_t::NONBASIC_LOWER);
   i_t num_basic = 0;
   for (i_t j = phase1_problem.num_cols - 1; j >= 0; --j) {
     const i_t col_start = phase1_problem.A.col_start[j];
     const i_t col_end   = phase1_problem.A.col_start[j + 1];
     const i_t nz        = col_end - col_start;
     if (nz == 1 && std::abs(phase1_problem.A.x[col_start]) == 1.0) {
-      phase1_vstatus[j] = variable_status_t::BASIC;
+      vstatus[j] = variable_status_t::BASIC;
       num_basic++;
     }
     if (num_basic == phase1_problem.num_rows) { break; }
@@ -160,9 +189,28 @@ lp_status_t solve_linear_program_advanced(const lp_problem_t<i_t, f_t>& original
   assert(num_basic == phase1_problem.num_rows);
   i_t iter = 0;
   lp_solution_t<i_t, f_t> phase1_solution(phase1_problem.num_rows, phase1_problem.num_cols);
-  std::vector<f_t> junk;
-  dual::status_t phase1_status = dual_phase2(
-    1, 1, start_time, phase1_problem, settings, phase1_vstatus, phase1_solution, iter, junk);
+
+  dual::status_t LU_status =
+    factorize_basis(lp, settings, vstatus, ft, basic_list, nonbasic_list, start_time);
+  if (LU_status == dual::status_t::NUMERICAL) {
+    settings.log.printf("Failed in factorizing the basis\n");
+    return lp_status_t::NUMERICAL_ISSUES;
+  }
+  if (LU_status == dual::status_t::TIME_LIMIT) { return lp_status_t::TIME_LIMIT; }
+
+  edge_norms.clear();
+  dual::status_t phase1_status = dual_phase2_with_basis_update(1,
+                                                               1,
+                                                               start_time,
+                                                               phase1_problem,
+                                                               settings,
+                                                               vstatus,
+                                                               ft,
+                                                               basic_list,
+                                                               nonbasic_list,
+                                                               phase1_solution,
+                                                               iter,
+                                                               edge_norms);
   if (phase1_status == dual::status_t::NUMERICAL ||
       phase1_status == dual::status_t::DUAL_UNBOUNDED) {
     settings.log.printf("Failed in Phase 1\n");
@@ -177,27 +225,62 @@ lp_status_t solve_linear_program_advanced(const lp_problem_t<i_t, f_t>& original
     lp_solution_t<i_t, f_t> solution(lp.num_rows, lp.num_cols);
     assert(lp.num_cols == phase1_problem.num_cols);
     assert(solution.x.size() == lp.num_cols);
-    vstatus = phase1_vstatus;
+
     edge_norms.clear();
-    dual::status_t status = dual_phase2(
-      2, iter == 0 ? 1 : 0, start_time, lp, settings, vstatus, solution, iter, edge_norms);
+    dual::status_t status = dual_phase2_with_basis_update(2,
+                                                          iter == 0 ? 1 : 0,
+                                                          start_time,
+                                                          lp,
+                                                          settings,
+                                                          vstatus,
+                                                          ft,
+                                                          basic_list,
+                                                          nonbasic_list,
+                                                          solution,
+                                                          iter,
+                                                          edge_norms);
+
     if (status == dual::status_t::NUMERICAL) {
       // Became dual infeasible. Try phase 1 again
-      phase1_vstatus = vstatus;
       settings.log.printf("Running Phase 1 again\n");
-      junk.clear();
-      dual_phase2(1,
-                  0,
-                  start_time,
-                  phase1_problem,
-                  settings,
-                  phase1_vstatus,
-                  phase1_solution,
-                  iter,
-                  edge_norms);
-      vstatus = phase1_vstatus;
+
+      dual::status_t LU_status =
+        factorize_basis(lp, settings, vstatus, ft, basic_list, nonbasic_list, start_time);
+
+      if (LU_status == dual::status_t::NUMERICAL) {
+        settings.log.printf("Failed in factorizing the basis\n");
+        return lp_status_t::NUMERICAL_ISSUES;
+      }
+
+      if (LU_status == dual::status_t::TIME_LIMIT) { return lp_status_t::TIME_LIMIT; }
+
       edge_norms.clear();
-      status = dual_phase2(2, 0, start_time, lp, settings, vstatus, solution, iter, edge_norms);
+      dual_phase2_with_basis_update(1,
+                                    0,
+                                    start_time,
+                                    phase1_problem,
+                                    settings,
+                                    vstatus,
+                                    ft,
+                                    basic_list,
+                                    nonbasic_list,
+                                    phase1_solution,
+                                    iter,
+                                    edge_norms);
+
+      edge_norms.clear();
+      status = dual_phase2_with_basis_update(2,
+                                             0,
+                                             start_time,
+                                             lp,
+                                             settings,
+                                             vstatus,
+                                             ft,
+                                             basic_list,
+                                             nonbasic_list,
+                                             solution,
+                                             iter,
+                                             edge_norms);
     }
     constexpr bool primal_cleanup = false;
     if (status == dual::status_t::OPTIMAL && primal_cleanup) {
@@ -602,6 +685,17 @@ template lp_status_t solve_linear_program_advanced(
   std::vector<variable_status_t>& vstatus,
   std::vector<double>& edge_norms);
 
+template lp_status_t solve_linear_program_with_basis_update(
+  const lp_problem_t<int, double>& original_lp,
+  const double start_time,
+  const simplex_solver_settings_t<int, double>& settings,
+  lp_solution_t<int, double>& original_solution,
+  basis_update_mpf_t<int, double>& ft,
+  std::vector<int>& basic_list,
+  std::vector<int>& nonbasic_list,
+  std::vector<variable_status_t>& vstatus,
+  std::vector<double>& edge_norms);
+
 template lp_status_t solve_linear_program_with_barrier(
   const user_problem_t<int, double>& user_problem,
   const simplex_solver_settings_t<int, double>& settings,
diff --git a/cpp/src/dual_simplex/solve.hpp b/cpp/src/dual_simplex/solve.hpp
index a54acf697..65418f5f2 100644
--- a/cpp/src/dual_simplex/solve.hpp
+++ b/cpp/src/dual_simplex/solve.hpp
@@ -17,6 +17,7 @@
 
 #pragma once
 
+#include <dual_simplex/basis_updates.hpp>
 #include <dual_simplex/initial_basis.hpp>
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/simplex_solver_settings.hpp>
@@ -56,6 +57,20 @@ lp_status_t solve_linear_program_advanced(const lp_problem_t<i_t, f_t>& original
                                           std::vector<variable_status_t>& vstatus,
                                           std::vector<f_t>& edge_norms);
 
+// Solve the LP using dual simplex and keep the `basis_update_mpf_t`
+// for future use.
+template <typename i_t, typename f_t>
+lp_status_t solve_linear_program_with_basis_update(
+  const lp_problem_t<i_t, f_t>& original_lp,
+  const f_t start_time,
+  const simplex_solver_settings_t<i_t, f_t>& settings,
+  lp_solution_t<i_t, f_t>& original_solution,
+  basis_update_mpf_t<i_t, f_t>& ft,
+  std::vector<i_t>& basic_list,
+  std::vector<i_t>& nonbasic_list,
+  std::vector<variable_status_t>& vstatus,
+  std::vector<f_t>& edge_norms);
+
 template <typename i_t, typename f_t>
 lp_status_t solve_linear_program_with_barrier(const user_problem_t<i_t, f_t>& user_problem,
                                               const simplex_solver_settings_t<i_t, f_t>& settings,

From 1cedfd6d8b9154030a29d0a086d0e4e58a93b3c8 Mon Sep 17 00:00:00 2001
From: nicolas <nguidotti@nvidia.com>
Date: Mon, 13 Oct 2025 17:31:02 +0200
Subject: [PATCH 03/11] small refactor

---
 cpp/src/dual_simplex/phase2.cpp | 14 ++------------
 1 file changed, 2 insertions(+), 12 deletions(-)

diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp
index b8becf3c8..37d4294b1 100644
--- a/cpp/src/dual_simplex/phase2.cpp
+++ b/cpp/src/dual_simplex/phase2.cpp
@@ -2209,23 +2209,13 @@ dual::status_t dual_phase2(i_t phase,
                            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;
-
-  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;
-  }
+  auto status = factorize_basis(lp, settings, vstatus, ft, basic_list, nonbasic_list, start_time);
+  if (status != dual::status_t::UNSET) { return status; }
 
-  if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; }
   return dual_phase2_with_basis_update(phase,
                                        slack_basis,
                                        start_time,

From 5768f50a681f85f803ce90a647eeb9b584094b7c Mon Sep 17 00:00:00 2001
From: nicolas <nguidotti@nvidia.com>
Date: Wed, 22 Oct 2025 14:32:01 +0200
Subject: [PATCH 04/11] fix incorrect matrix dimension in LU after LP presolve

---
 cpp/src/dual_simplex/basis_updates.cpp |  1 +
 cpp/src/dual_simplex/basis_updates.hpp | 17 +++++++++++++++++
 2 files changed, 18 insertions(+)

diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp
index 5efc86b58..e9cbfd345 100644
--- a/cpp/src/dual_simplex/basis_updates.cpp
+++ b/cpp/src/dual_simplex/basis_updates.cpp
@@ -2059,6 +2059,7 @@ int basis_update_mpf_t<i_t, f_t>::factorize_basis(
   std::vector<i_t> deficient;
   std::vector<i_t> slacks_needed;
 
+  if (L0_.m != A.m) { resize(A.m); }
   if (dual_simplex::factorize_basis(A,
                                     settings,
                                     basic_list,
diff --git a/cpp/src/dual_simplex/basis_updates.hpp b/cpp/src/dual_simplex/basis_updates.hpp
index 3198ccd4b..11f044d03 100644
--- a/cpp/src/dual_simplex/basis_updates.hpp
+++ b/cpp/src/dual_simplex/basis_updates.hpp
@@ -291,6 +291,23 @@ class basis_update_mpf_t {
     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++;

From 0ed38a907cb247d6879637342f56eb413a5b9c7e Mon Sep 17 00:00:00 2001
From: Christopher Maes <cmaes@nvidia.com>
Date: Wed, 22 Oct 2025 13:51:20 -0700
Subject: [PATCH 05/11] Refactor, simplify interface, and reduce number of code
 changes

---
 cpp/src/dual_simplex/basis_updates.cpp | 42 ++++++++--------
 cpp/src/dual_simplex/basis_updates.hpp | 10 ++--
 cpp/src/dual_simplex/phase2.cpp        | 63 +++++++++++++++++++-----
 cpp/src/dual_simplex/phase2.hpp        |  1 +
 cpp/src/dual_simplex/solve.cpp         | 68 +++++++-------------------
 5 files changed, 98 insertions(+), 86 deletions(-)

diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp
index e9cbfd345..b34ebcde8 100644
--- a/cpp/src/dual_simplex/basis_updates.cpp
+++ b/cpp/src/dual_simplex/basis_updates.cpp
@@ -2049,7 +2049,7 @@ void basis_update_mpf_t<i_t, f_t>::multiply_lu(csc_matrix_t<i_t, f_t>& out) cons
 }
 
 template <typename i_t, typename f_t>
-int basis_update_mpf_t<i_t, f_t>::factorize_basis(
+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,
@@ -2060,16 +2060,16 @@ int basis_update_mpf_t<i_t, f_t>::factorize_basis(
   std::vector<i_t> slacks_needed;
 
   if (L0_.m != A.m) { resize(A.m); }
-  if (dual_simplex::factorize_basis(A,
-                                    settings,
-                                    basic_list,
-                                    L0_,
-                                    U0_,
-                                    row_permutation_,
-                                    inverse_row_permutation_,
-                                    col_permutation_,
-                                    deficient,
-                                    slacks_needed) == -1) {
+  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);
 
@@ -2089,16 +2089,16 @@ int basis_update_mpf_t<i_t, f_t>::factorize_basis(
     }
 #endif
 
-    if (dual_simplex::factorize_basis(A,
-                                      settings,
-                                      basic_list,
-                                      L0_,
-                                      U0_,
-                                      row_permutation_,
-                                      inverse_row_permutation_,
-                                      col_permutation_,
-                                      deficient,
-                                      slacks_needed) == -1) {
+    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 (L.check_matrix() == -1) { settings.log.printf("Bad L after basis repair\n"); }
 #endif
diff --git a/cpp/src/dual_simplex/basis_updates.hpp b/cpp/src/dual_simplex/basis_updates.hpp
index 11f044d03..11e31732d 100644
--- a/cpp/src/dual_simplex/basis_updates.hpp
+++ b/cpp/src/dual_simplex/basis_updates.hpp
@@ -388,11 +388,11 @@ 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);
+  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()
diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp
index 37d4294b1..74caa44dd 100644
--- a/cpp/src/dual_simplex/phase2.cpp
+++ b/cpp/src/dual_simplex/phase2.cpp
@@ -2243,6 +2243,43 @@ dual::status_t dual_phase2_with_basis_update(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_basis_update(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_basis_update(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;
@@ -2269,6 +2306,18 @@ dual::status_t dual_phase2_with_basis_update(i_t phase,
   std::vector<f_t> z_old                     = z;
 
   phase2::bound_info(lp, settings);
+  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;
+    }
+
+    if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; }
+  }
 
   std::vector<f_t> c_basic(m);
   for (i_t k = 0; k < m; ++k) {
@@ -2905,14 +2954,14 @@ dual::status_t dual_phase2_with_basis_update(i_t phase,
 #endif
     if (should_refactor) {
       bool should_recompute_x = false;
-      if (ft.factorize_basis(lp.A, settings, basic_list, nonbasic_list, vstatus) > 0) {
+      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; }
         i_t count = 0;
         i_t deficient_size;
         while ((deficient_size =
-                  ft.factorize_basis(lp.A, settings, basic_list, nonbasic_list, vstatus)) > 0) {
+                  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));
@@ -3023,6 +3072,7 @@ template dual::status_t dual_phase2<int, double>(
 template dual::status_t dual_phase2_with_basis_update<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,
@@ -3033,15 +3083,6 @@ template dual::status_t dual_phase2_with_basis_update<int, double>(
   lp_solution_t<int, double>& sol,
   int& iter,
   std::vector<double>& steepest_edge_norms);
-
-template dual::status_t factorize_basis<int, double>(
-  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,
-  double start_time);
 #endif
 
 }  // namespace cuopt::linear_programming::dual_simplex
diff --git a/cpp/src/dual_simplex/phase2.hpp b/cpp/src/dual_simplex/phase2.hpp
index bd1a80389..4719c9d5b 100644
--- a/cpp/src/dual_simplex/phase2.hpp
+++ b/cpp/src/dual_simplex/phase2.hpp
@@ -79,6 +79,7 @@ dual::status_t dual_phase2(i_t phase,
 template <typename i_t, typename f_t>
 dual::status_t dual_phase2_with_basis_update(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,
diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp
index 4da964474..1b3e2d258 100644
--- a/cpp/src/dual_simplex/solve.cpp
+++ b/cpp/src/dual_simplex/solve.cpp
@@ -227,60 +227,30 @@ lp_status_t solve_linear_program_with_basis_update(
     assert(solution.x.size() == lp.num_cols);
 
     edge_norms.clear();
-    dual::status_t status = dual_phase2_with_basis_update(2,
-                                                          iter == 0 ? 1 : 0,
-                                                          start_time,
-                                                          lp,
-                                                          settings,
-                                                          vstatus,
-                                                          ft,
-                                                          basic_list,
-                                                          nonbasic_list,
-                                                          solution,
-                                                          iter,
-                                                          edge_norms);
-
+    bool initialize_basis_update = true;
+    dual::status_t status = dual_phase2_with_basis_update(
+      2, iter == 0 ? 1 : 0, initialize_basis_update, start_time, lp, settings, vstatus, ft, basic_list, nonbasic_list, solution, iter, edge_norms);
     if (status == dual::status_t::NUMERICAL) {
       // Became dual infeasible. Try phase 1 again
       settings.log.printf("Running Phase 1 again\n");
-
-      dual::status_t LU_status =
-        factorize_basis(lp, settings, vstatus, ft, basic_list, nonbasic_list, start_time);
-
-      if (LU_status == dual::status_t::NUMERICAL) {
-        settings.log.printf("Failed in factorizing the basis\n");
-        return lp_status_t::NUMERICAL_ISSUES;
-      }
-
-      if (LU_status == dual::status_t::TIME_LIMIT) { return lp_status_t::TIME_LIMIT; }
-
-      edge_norms.clear();
+      junk.clear();
+      initialize_basis_update = false;
       dual_phase2_with_basis_update(1,
-                                    0,
-                                    start_time,
-                                    phase1_problem,
-                                    settings,
-                                    vstatus,
-                                    ft,
-                                    basic_list,
-                                    nonbasic_list,
-                                    phase1_solution,
-                                    iter,
-                                    edge_norms);
-
+                  0,
+                  initialize_basis_update,
+                  start_time,
+                  phase1_problem,
+                  settings,
+                  phase1_vstatus,
+                  ft,
+                  basic_list,
+                  nonbasic_list,
+                  phase1_solution,
+                  iter,
+                  edge_norms);
+      vstatus = phase1_vstatus;
       edge_norms.clear();
-      status = dual_phase2_with_basis_update(2,
-                                             0,
-                                             start_time,
-                                             lp,
-                                             settings,
-                                             vstatus,
-                                             ft,
-                                             basic_list,
-                                             nonbasic_list,
-                                             solution,
-                                             iter,
-                                             edge_norms);
+      status = dual_phase2_with_basis_update(2, 0, initialize_basis_update, start_time, lp, settings, vstatus, ft, basic_list, nonbasic_list, solution, iter, edge_norms);
     }
     constexpr bool primal_cleanup = false;
     if (status == dual::status_t::OPTIMAL && primal_cleanup) {

From 166c44653ef4cc7943a82660c01be476403f9560 Mon Sep 17 00:00:00 2001
From: Christopher Maes <cmaes@nvidia.com>
Date: Wed, 22 Oct 2025 16:18:16 -0700
Subject: [PATCH 06/11] Further removal; missing when cherry-picking previous
 commit

---
 cpp/src/dual_simplex/basis_solves.cpp |  8 +---
 cpp/src/dual_simplex/phase2.cpp       | 63 ---------------------------
 cpp/src/dual_simplex/phase2.hpp       |  9 ----
 cpp/src/dual_simplex/solve.cpp        | 36 ++++-----------
 4 files changed, 10 insertions(+), 106 deletions(-)

diff --git a/cpp/src/dual_simplex/basis_solves.cpp b/cpp/src/dual_simplex/basis_solves.cpp
index 91a2566b9..c6f0ec486 100644
--- a/cpp/src/dual_simplex/basis_solves.cpp
+++ b/cpp/src/dual_simplex/basis_solves.cpp
@@ -46,13 +46,6 @@ void get_basis_from_vstatus(i_t m,
   i_t n             = vstatus.size();
   i_t num_basic     = 0;
   i_t num_non_basic = 0;
-
-  assert(n >= m);
-  nonbasic_list.clear();
-  superbasic_list.clear();
-  nonbasic_list.reserve(n - m);
-  basis_list.resize(m);
-
   for (i_t j = 0; j < n; ++j) {
     if (vstatus[j] == variable_status_t::BASIC) {
       basis_list[num_basic++] = j;
@@ -68,6 +61,7 @@ void get_basis_from_vstatus(i_t m,
       superbasic_list.push_back(j);
     }
   }
+  i_t num_super_basic = superbasic_list.size();
   assert(num_basic == m);
 }
 
diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp
index 74caa44dd..3509fa1be 100644
--- a/cpp/src/dual_simplex/phase2.cpp
+++ b/cpp/src/dual_simplex/phase2.cpp
@@ -2169,33 +2169,6 @@ class phase2_timers_t {
 };
 
 }  // namespace phase2
-template <typename i_t, typename f_t>
-dual::status_t factorize_basis(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,
-                               f_t start_time)
-{
-  const i_t m = lp.num_rows;
-  const i_t n = lp.num_cols;
-  assert(m <= n);
-  assert(vstatus.size() == n);
-
-  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.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::status_t::UNSET;
-}
 
 template <typename i_t, typename f_t>
 dual::status_t dual_phase2(i_t phase,
@@ -2207,42 +2180,6 @@ 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;
-  std::vector<i_t> basic_list(m);
-  std::vector<i_t> nonbasic_list;
-  basis_update_mpf_t<i_t, f_t> ft(m, settings.refactor_frequency);
-
-  auto status = factorize_basis(lp, settings, vstatus, ft, basic_list, nonbasic_list, start_time);
-  if (status != dual::status_t::UNSET) { return status; }
-
-  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;
diff --git a/cpp/src/dual_simplex/phase2.hpp b/cpp/src/dual_simplex/phase2.hpp
index 4719c9d5b..a4f8066d0 100644
--- a/cpp/src/dual_simplex/phase2.hpp
+++ b/cpp/src/dual_simplex/phase2.hpp
@@ -56,15 +56,6 @@ static std::string status_to_string(status_t status)
 }
 }  // namespace dual
 
-template <typename i_t, typename f_t>
-dual::status_t factorize_basis(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,
-                               f_t start_time);
-
 template <typename i_t, typename f_t>
 dual::status_t dual_phase2(i_t phase,
                            i_t slack_basis,
diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp
index 1b3e2d258..5850893bc 100644
--- a/cpp/src/dual_simplex/solve.cpp
+++ b/cpp/src/dual_simplex/solve.cpp
@@ -167,21 +167,21 @@ lp_status_t solve_linear_program_with_basis_update(
   column_scaling(presolved_lp, settings, lp, column_scales);
   assert(presolved_lp.num_cols == lp.num_cols);
   lp_problem_t<i_t, f_t> phase1_problem(original_lp.handle_ptr, 1, 1, 1);
-
+  std::vector<variable_status_t> phase1_vstatus;
   f_t phase1_obj = -inf;
   create_phase1_problem(lp, phase1_problem);
   assert(phase1_problem.num_cols == presolved_lp.num_cols);
 
   // Set the vstatus for the phase1 problem based on a slack basis
-  vstatus.resize(phase1_problem.num_cols);
-  std::fill(vstatus.begin(), vstatus.end(), variable_status_t::NONBASIC_LOWER);
+  phase1_vstatus.resize(phase1_problem.num_cols);
+  std::fill(phase1_vstatus.begin(), phase1_vstatus.end(), variable_status_t::NONBASIC_LOWER);
   i_t num_basic = 0;
   for (i_t j = phase1_problem.num_cols - 1; j >= 0; --j) {
     const i_t col_start = phase1_problem.A.col_start[j];
     const i_t col_end   = phase1_problem.A.col_start[j + 1];
     const i_t nz        = col_end - col_start;
     if (nz == 1 && std::abs(phase1_problem.A.x[col_start]) == 1.0) {
-      vstatus[j] = variable_status_t::BASIC;
+      phase1_vstatus[j] = variable_status_t::BASIC;
       num_basic++;
     }
     if (num_basic == phase1_problem.num_rows) { break; }
@@ -189,28 +189,9 @@ lp_status_t solve_linear_program_with_basis_update(
   assert(num_basic == phase1_problem.num_rows);
   i_t iter = 0;
   lp_solution_t<i_t, f_t> phase1_solution(phase1_problem.num_rows, phase1_problem.num_cols);
-
-  dual::status_t LU_status =
-    factorize_basis(lp, settings, vstatus, ft, basic_list, nonbasic_list, start_time);
-  if (LU_status == dual::status_t::NUMERICAL) {
-    settings.log.printf("Failed in factorizing the basis\n");
-    return lp_status_t::NUMERICAL_ISSUES;
-  }
-  if (LU_status == dual::status_t::TIME_LIMIT) { return lp_status_t::TIME_LIMIT; }
-
-  edge_norms.clear();
-  dual::status_t phase1_status = dual_phase2_with_basis_update(1,
-                                                               1,
-                                                               start_time,
-                                                               phase1_problem,
-                                                               settings,
-                                                               vstatus,
-                                                               ft,
-                                                               basic_list,
-                                                               nonbasic_list,
-                                                               phase1_solution,
-                                                               iter,
-                                                               edge_norms);
+  std::vector<f_t> junk;
+  dual::status_t phase1_status = dual_phase2(
+    1, 1, start_time, phase1_problem, settings, phase1_vstatus, phase1_solution, iter, junk);
   if (phase1_status == dual::status_t::NUMERICAL ||
       phase1_status == dual::status_t::DUAL_UNBOUNDED) {
     settings.log.printf("Failed in Phase 1\n");
@@ -225,13 +206,14 @@ lp_status_t solve_linear_program_with_basis_update(
     lp_solution_t<i_t, f_t> solution(lp.num_rows, lp.num_cols);
     assert(lp.num_cols == phase1_problem.num_cols);
     assert(solution.x.size() == lp.num_cols);
-
+    vstatus = phase1_vstatus;
     edge_norms.clear();
     bool initialize_basis_update = true;
     dual::status_t status = dual_phase2_with_basis_update(
       2, iter == 0 ? 1 : 0, initialize_basis_update, start_time, lp, settings, vstatus, ft, basic_list, nonbasic_list, solution, iter, edge_norms);
     if (status == dual::status_t::NUMERICAL) {
       // Became dual infeasible. Try phase 1 again
+      phase1_vstatus = vstatus;
       settings.log.printf("Running Phase 1 again\n");
       junk.clear();
       initialize_basis_update = false;

From 27884e46a1a1059cf8c9e146b9abf57f7eaded3b Mon Sep 17 00:00:00 2001
From: Christopher Maes <cmaes@nvidia.com>
Date: Wed, 22 Oct 2025 16:36:49 -0700
Subject: [PATCH 07/11] Rename functions. Address code-rabbit issues

---
 cpp/src/dual_simplex/basis_updates.cpp |  3 +-
 cpp/src/dual_simplex/basis_updates.hpp |  1 -
 cpp/src/dual_simplex/phase2.cpp        | 54 +++++++++++++-------------
 cpp/src/dual_simplex/phase2.hpp        | 26 ++++++-------
 cpp/src/dual_simplex/solve.cpp         | 28 ++++++-------
 5 files changed, 56 insertions(+), 56 deletions(-)

diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp
index b34ebcde8..2c1da6879 100644
--- a/cpp/src/dual_simplex/basis_updates.cpp
+++ b/cpp/src/dual_simplex/basis_updates.cpp
@@ -2074,6 +2074,7 @@ int basis_update_mpf_t<i_t, f_t>::refactor_basis(
     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) {
@@ -2100,7 +2101,7 @@ int basis_update_mpf_t<i_t, f_t>::refactor_basis(
                         deficient,
                         slacks_needed) == -1) {
 #ifdef CHECK_L_FACTOR
-      if (L.check_matrix() == -1) { settings.log.printf("Bad L after basis repair\n"); }
+      if (L0_.check_matrix() == -1) { settings.log.printf("Bad L after basis repair\n"); }
 #endif
       return deficient.size();
     }
diff --git a/cpp/src/dual_simplex/basis_updates.hpp b/cpp/src/dual_simplex/basis_updates.hpp
index 11e31732d..e13f4dddf 100644
--- a/cpp/src/dual_simplex/basis_updates.hpp
+++ b/cpp/src/dual_simplex/basis_updates.hpp
@@ -284,7 +284,6 @@ class basis_update_mpf_t {
 
   i_t reset()
   {
-    inverse_permutation(row_permutation_, inverse_row_permutation_);
     clear();
     compute_transposes();
     reset_stats();
diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp
index 3509fa1be..098c3b6e2 100644
--- a/cpp/src/dual_simplex/phase2.cpp
+++ b/cpp/src/dual_simplex/phase2.cpp
@@ -2188,35 +2188,35 @@ dual::status_t dual_phase2(i_t phase,
   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_basis_update(phase,
-                                       slack_basis,
-                                       initialize_basis,
-                                       start_time,
-                                       lp,
-                                       settings,
-                                       vstatus,
-                                       ft,
-                                       basic_list,
-                                       nonbasic_list,
-                                       sol,
-                                       iter,
-                                       delta_y_steepest_edge);
+  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_basis_update(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)
+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;
@@ -3006,7 +3006,7 @@ template dual::status_t dual_phase2<int, double>(
   int& iter,
   std::vector<double>& steepest_edge_norms);
 
-template dual::status_t dual_phase2_with_basis_update<int, double>(
+template dual::status_t dual_phase2_with_advanced_basis<int, double>(
   int phase,
   int slack_basis,
   bool initialize_basis,
diff --git a/cpp/src/dual_simplex/phase2.hpp b/cpp/src/dual_simplex/phase2.hpp
index a4f8066d0..e237ca372 100644
--- a/cpp/src/dual_simplex/phase2.hpp
+++ b/cpp/src/dual_simplex/phase2.hpp
@@ -68,18 +68,18 @@ dual::status_t dual_phase2(i_t phase,
                            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,
-                                             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);
+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);
 
 }  // namespace cuopt::linear_programming::dual_simplex
diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp
index 5850893bc..d286b0dbf 100644
--- a/cpp/src/dual_simplex/solve.cpp
+++ b/cpp/src/dual_simplex/solve.cpp
@@ -123,19 +123,19 @@ lp_status_t solve_linear_program_advanced(const lp_problem_t<i_t, f_t>& original
   std::vector<i_t> basic_list(m);
   std::vector<i_t> nonbasic_list;
   basis_update_mpf_t<i_t, f_t> ft(m, settings.refactor_frequency);
-  return solve_linear_program_with_basis_update(original_lp,
-                                                start_time,
-                                                settings,
-                                                original_solution,
-                                                ft,
-                                                basic_list,
-                                                nonbasic_list,
-                                                vstatus,
-                                                edge_norms);
+  return solve_linear_program_with_advanced_basis(original_lp,
+                                                  start_time,
+                                                  settings,
+                                                  original_solution,
+                                                  ft,
+                                                  basic_list,
+                                                  nonbasic_list,
+                                                  vstatus,
+                                                  edge_norms);
 }
 
 template <typename i_t, typename f_t>
-lp_status_t solve_linear_program_with_basis_update(
+lp_status_t solve_linear_program_with_advanced_basis(
   const lp_problem_t<i_t, f_t>& original_lp,
   const f_t start_time,
   const simplex_solver_settings_t<i_t, f_t>& settings,
@@ -209,7 +209,7 @@ lp_status_t solve_linear_program_with_basis_update(
     vstatus = phase1_vstatus;
     edge_norms.clear();
     bool initialize_basis_update = true;
-    dual::status_t status = dual_phase2_with_basis_update(
+    dual::status_t status = dual_phase2_with_advanced_basis(
       2, iter == 0 ? 1 : 0, initialize_basis_update, start_time, lp, settings, vstatus, ft, basic_list, nonbasic_list, solution, iter, edge_norms);
     if (status == dual::status_t::NUMERICAL) {
       // Became dual infeasible. Try phase 1 again
@@ -217,7 +217,7 @@ lp_status_t solve_linear_program_with_basis_update(
       settings.log.printf("Running Phase 1 again\n");
       junk.clear();
       initialize_basis_update = false;
-      dual_phase2_with_basis_update(1,
+      dual_phase2_with_advanced_basis(1,
                   0,
                   initialize_basis_update,
                   start_time,
@@ -232,7 +232,7 @@ lp_status_t solve_linear_program_with_basis_update(
                   edge_norms);
       vstatus = phase1_vstatus;
       edge_norms.clear();
-      status = dual_phase2_with_basis_update(2, 0, initialize_basis_update, start_time, lp, settings, vstatus, ft, basic_list, nonbasic_list, solution, iter, edge_norms);
+      status = dual_phase2_with_advanced_basis(2, 0, initialize_basis_update, start_time, lp, settings, vstatus, ft, basic_list, nonbasic_list, solution, iter, edge_norms);
     }
     constexpr bool primal_cleanup = false;
     if (status == dual::status_t::OPTIMAL && primal_cleanup) {
@@ -637,7 +637,7 @@ template lp_status_t solve_linear_program_advanced(
   std::vector<variable_status_t>& vstatus,
   std::vector<double>& edge_norms);
 
-template lp_status_t solve_linear_program_with_basis_update(
+template lp_status_t solve_linear_program_with_advanced_basis(
   const lp_problem_t<int, double>& original_lp,
   const double start_time,
   const simplex_solver_settings_t<int, double>& settings,

From deffe681242b3bfc5c44e9a146058bb776aa488f Mon Sep 17 00:00:00 2001
From: Christopher Maes <cmaes@nvidia.com>
Date: Wed, 22 Oct 2025 16:37:22 -0700
Subject: [PATCH 08/11] Style fixes

---
 cpp/src/dual_simplex/solve.cpp | 53 ++++++++++++++++++++++++----------
 1 file changed, 38 insertions(+), 15 deletions(-)

diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp
index d286b0dbf..ba6a05084 100644
--- a/cpp/src/dual_simplex/solve.cpp
+++ b/cpp/src/dual_simplex/solve.cpp
@@ -209,8 +209,19 @@ lp_status_t solve_linear_program_with_advanced_basis(
     vstatus = phase1_vstatus;
     edge_norms.clear();
     bool initialize_basis_update = true;
-    dual::status_t status = dual_phase2_with_advanced_basis(
-      2, iter == 0 ? 1 : 0, initialize_basis_update, start_time, lp, settings, vstatus, ft, basic_list, nonbasic_list, solution, iter, edge_norms);
+    dual::status_t status        = dual_phase2_with_advanced_basis(2,
+                                                            iter == 0 ? 1 : 0,
+                                                            initialize_basis_update,
+                                                            start_time,
+                                                            lp,
+                                                            settings,
+                                                            vstatus,
+                                                            ft,
+                                                            basic_list,
+                                                            nonbasic_list,
+                                                            solution,
+                                                            iter,
+                                                            edge_norms);
     if (status == dual::status_t::NUMERICAL) {
       // Became dual infeasible. Try phase 1 again
       phase1_vstatus = vstatus;
@@ -218,21 +229,33 @@ lp_status_t solve_linear_program_with_advanced_basis(
       junk.clear();
       initialize_basis_update = false;
       dual_phase2_with_advanced_basis(1,
-                  0,
-                  initialize_basis_update,
-                  start_time,
-                  phase1_problem,
-                  settings,
-                  phase1_vstatus,
-                  ft,
-                  basic_list,
-                  nonbasic_list,
-                  phase1_solution,
-                  iter,
-                  edge_norms);
+                                      0,
+                                      initialize_basis_update,
+                                      start_time,
+                                      phase1_problem,
+                                      settings,
+                                      phase1_vstatus,
+                                      ft,
+                                      basic_list,
+                                      nonbasic_list,
+                                      phase1_solution,
+                                      iter,
+                                      edge_norms);
       vstatus = phase1_vstatus;
       edge_norms.clear();
-      status = dual_phase2_with_advanced_basis(2, 0, initialize_basis_update, start_time, lp, settings, vstatus, ft, basic_list, nonbasic_list, solution, iter, edge_norms);
+      status = dual_phase2_with_advanced_basis(2,
+                                               0,
+                                               initialize_basis_update,
+                                               start_time,
+                                               lp,
+                                               settings,
+                                               vstatus,
+                                               ft,
+                                               basic_list,
+                                               nonbasic_list,
+                                               solution,
+                                               iter,
+                                               edge_norms);
     }
     constexpr bool primal_cleanup = false;
     if (status == dual::status_t::OPTIMAL && primal_cleanup) {

From f779cb71f01aea6096eb2100f891eb5f16316e18 Mon Sep 17 00:00:00 2001
From: Christopher Maes <cmaes@nvidia.com>
Date: Wed, 22 Oct 2025 17:23:36 -0700
Subject: [PATCH 09/11] Remove confusing col_permutation_ from middle-product
 update

---
 cpp/src/dual_simplex/basis_updates.cpp |  9 +++++----
 cpp/src/dual_simplex/basis_updates.hpp | 10 ----------
 2 files changed, 5 insertions(+), 14 deletions(-)

diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp
index 2c1da6879..d2127cc4d 100644
--- a/cpp/src/dual_simplex/basis_updates.cpp
+++ b/cpp/src/dual_simplex/basis_updates.cpp
@@ -2060,6 +2060,7 @@ int basis_update_mpf_t<i_t, f_t>::refactor_basis(
   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,
@@ -2067,7 +2068,7 @@ int basis_update_mpf_t<i_t, f_t>::refactor_basis(
                       U0_,
                       row_permutation_,
                       inverse_row_permutation_,
-                      col_permutation_,
+                      q,
                       deficient,
                       slacks_needed) == -1) {
     settings.log.debug("Initial factorization failed\n");
@@ -2097,7 +2098,7 @@ int basis_update_mpf_t<i_t, f_t>::refactor_basis(
                         U0_,
                         row_permutation_,
                         inverse_row_permutation_,
-                        col_permutation_,
+                        q,
                         deficient,
                         slacks_needed) == -1) {
 #ifdef CHECK_L_FACTOR
@@ -2108,8 +2109,8 @@ int basis_update_mpf_t<i_t, f_t>::refactor_basis(
     settings.log.debug("Basis repaired\n");
   }
 
-  assert(col_permutation_.size() == A.m);
-  reorder_basic_list(col_permutation_, basic_list);
+  assert(q.size() == A.m);
+  reorder_basic_list(q, basic_list);  // We no longer need q after reordering the basic list
   reset();
   return 0;
 }
diff --git a/cpp/src/dual_simplex/basis_updates.hpp b/cpp/src/dual_simplex/basis_updates.hpp
index e13f4dddf..d92a528a7 100644
--- a/cpp/src/dual_simplex/basis_updates.hpp
+++ b/cpp/src/dual_simplex/basis_updates.hpp
@@ -184,8 +184,6 @@ class basis_update_mpf_t {
       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),
@@ -214,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),
@@ -297,8 +293,6 @@ class basis_update_mpf_t {
     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);
@@ -398,8 +392,6 @@ class basis_update_mpf_t {
   {
     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;
@@ -452,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

From f669f642a900d28a5c9c439b91d922eb6d1ed8bc Mon Sep 17 00:00:00 2001
From: nicolas <nguidotti@nvidia.com>
Date: Thu, 23 Oct 2025 14:20:04 +0200
Subject: [PATCH 10/11] removed unused variable

---
 cpp/src/dual_simplex/solve.cpp | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp
index ba6a05084..b6f84a9dc 100644
--- a/cpp/src/dual_simplex/solve.cpp
+++ b/cpp/src/dual_simplex/solve.cpp
@@ -189,9 +189,9 @@ lp_status_t solve_linear_program_with_advanced_basis(
   assert(num_basic == phase1_problem.num_rows);
   i_t iter = 0;
   lp_solution_t<i_t, f_t> phase1_solution(phase1_problem.num_rows, phase1_problem.num_cols);
-  std::vector<f_t> junk;
+  edge_norms.clear();
   dual::status_t phase1_status = dual_phase2(
-    1, 1, start_time, phase1_problem, settings, phase1_vstatus, phase1_solution, iter, junk);
+    1, 1, start_time, phase1_problem, settings, phase1_vstatus, phase1_solution, iter, edge_norms);
   if (phase1_status == dual::status_t::NUMERICAL ||
       phase1_status == dual::status_t::DUAL_UNBOUNDED) {
     settings.log.printf("Failed in Phase 1\n");
@@ -226,7 +226,7 @@ lp_status_t solve_linear_program_with_advanced_basis(
       // Became dual infeasible. Try phase 1 again
       phase1_vstatus = vstatus;
       settings.log.printf("Running Phase 1 again\n");
-      junk.clear();
+      edge_norms.clear();
       initialize_basis_update = false;
       dual_phase2_with_advanced_basis(1,
                                       0,

From 8a06d3ea0fce2c443d37dc4c0cc02188bf757490 Mon Sep 17 00:00:00 2001
From: nicolas <nguidotti@nvidia.com>
Date: Thu, 23 Oct 2025 14:33:12 +0200
Subject: [PATCH 11/11] added missing function declaration

---
 cpp/src/dual_simplex/solve.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/cpp/src/dual_simplex/solve.hpp b/cpp/src/dual_simplex/solve.hpp
index 65418f5f2..53146eda9 100644
--- a/cpp/src/dual_simplex/solve.hpp
+++ b/cpp/src/dual_simplex/solve.hpp
@@ -60,7 +60,7 @@ lp_status_t solve_linear_program_advanced(const lp_problem_t<i_t, f_t>& original
 // Solve the LP using dual simplex and keep the `basis_update_mpf_t`
 // for future use.
 template <typename i_t, typename f_t>
-lp_status_t solve_linear_program_with_basis_update(
+lp_status_t solve_linear_program_with_advanced_basis(
   const lp_problem_t<i_t, f_t>& original_lp,
   const f_t start_time,
   const simplex_solver_settings_t<i_t, f_t>& settings,