diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp
index 18f74ccc2..d2127cc4d 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,73 @@ void basis_update_mpf_t<i_t, f_t>::multiply_lu(csc_matrix_t<i_t, f_t>& out) cons
   out.nz_max = B_nz;
 }
 
+template <typename i_t, typename f_t>
+int basis_update_mpf_t<i_t, f_t>::refactor_basis(
+  const csc_matrix_t<i_t, f_t>& A,
+  const simplex_solver_settings_t<i_t, f_t>& settings,
+  std::vector<i_t>& basic_list,
+  std::vector<i_t>& nonbasic_list,
+  std::vector<variable_status_t>& vstatus)
+{
+  std::vector<i_t> deficient;
+  std::vector<i_t> slacks_needed;
+
+  if (L0_.m != A.m) { resize(A.m); }
+  std::vector<i_t> q;
+  if (factorize_basis(A,
+                      settings,
+                      basic_list,
+                      L0_,
+                      U0_,
+                      row_permutation_,
+                      inverse_row_permutation_,
+                      q,
+                      deficient,
+                      slacks_needed) == -1) {
+    settings.log.debug("Initial factorization failed\n");
+    basis_repair(A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
+
+#ifdef CHECK_BASIS_REPAIR
+    const i_t m = A.m;
+    csc_matrix_t<i_t, f_t> B(m, m, 0);
+    form_b(A, basic_list, B);
+    for (i_t k = 0; k < deficient.size(); ++k) {
+      const i_t j         = deficient[k];
+      const i_t col_start = B.col_start[j];
+      const i_t col_end   = B.col_start[j + 1];
+      const i_t col_nz    = col_end - col_start;
+      if (col_nz != 1) { settings.log.printf("Deficient column %d has %d nonzeros\n", j, col_nz); }
+      const i_t i = B.i[col_start];
+      if (i != slacks_needed[k]) {
+        settings.log.printf("Slack %d needed but found %d instead\n", slacks_needed[k], i);
+      }
+    }
+#endif
+
+    if (factorize_basis(A,
+                        settings,
+                        basic_list,
+                        L0_,
+                        U0_,
+                        row_permutation_,
+                        inverse_row_permutation_,
+                        q,
+                        deficient,
+                        slacks_needed) == -1) {
+#ifdef CHECK_L_FACTOR
+      if (L0_.check_matrix() == -1) { settings.log.printf("Bad L after basis repair\n"); }
+#endif
+      return deficient.size();
+    }
+    settings.log.debug("Basis repaired\n");
+  }
+
+  assert(q.size() == A.m);
+  reorder_basic_list(q, basic_list);  // We no longer need q after reordering the basic list
+  reset();
+  return 0;
+}
+
 #ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE
 template class basis_update_t<int, double>;
 template class basis_update_mpf_t<int, double>;
diff --git a/cpp/src/dual_simplex/basis_updates.hpp b/cpp/src/dual_simplex/basis_updates.hpp
index b12950abb..d92a528a7 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,31 @@ class basis_update_t {
 template <typename i_t, typename f_t>
 class basis_update_mpf_t {
  public:
+  basis_update_mpf_t(i_t n, const i_t refactor_frequency)
+    : L0_(n, n, 1),
+      U0_(n, n, 1),
+      row_permutation_(n),
+      inverse_row_permutation_(n),
+      S_(n, 0, 0),
+      xi_workspace_(2 * n, 0),
+      x_workspace_(n, 0.0),
+      U0_transpose_(1, 1, 1),
+      L0_transpose_(1, 1, 1),
+      refactor_frequency_(refactor_frequency),
+      total_sparse_L_transpose_(0),
+      total_dense_L_transpose_(0),
+      total_sparse_L_(0),
+      total_dense_L_(0),
+      total_sparse_U_transpose_(0),
+      total_dense_U_transpose_(0),
+      total_sparse_U_(0),
+      total_dense_U_(0),
+      hypersparse_threshold_(0.05)
+  {
+    clear();
+    reset_stats();
+  }
+
   basis_update_mpf_t(const csc_matrix_t<i_t, f_t>& Linit,
                      const csc_matrix_t<i_t, f_t>& Uinit,
                      const std::vector<i_t>& p,
@@ -185,8 +212,6 @@ class basis_update_mpf_t {
       row_permutation_(p),
       inverse_row_permutation_(p.size()),
       S_(Linit.m, 0, 0),
-      col_permutation_(Linit.m),
-      inverse_col_permutation_(Linit.m),
       xi_workspace_(2 * Linit.m, 0),
       x_workspace_(Linit.m, 0.0),
       U0_transpose_(1, 1, 1),
@@ -205,7 +230,7 @@ class basis_update_mpf_t {
     inverse_permutation(row_permutation_, inverse_row_permutation_);
     clear();
     compute_transposes();
-    reset_stas();
+    reset_stats();
   }
 
   void print_stats() const
@@ -226,7 +251,7 @@ class basis_update_mpf_t {
     // clang-format on
   }
 
-  void reset_stas()
+  void reset_stats()
   {
     num_calls_L_           = 0;
     num_calls_U_           = 0;
@@ -249,10 +274,33 @@ class basis_update_mpf_t {
     inverse_permutation(row_permutation_, inverse_row_permutation_);
     clear();
     compute_transposes();
-    reset_stas();
+    reset_stats();
+    return 0;
+  }
+
+  i_t reset()
+  {
+    clear();
+    compute_transposes();
+    reset_stats();
     return 0;
   }
 
+  void resize(i_t n)
+  {
+    L0_.resize(n, n, 1);
+    U0_.resize(n, n, 1);
+    row_permutation_.resize(n);
+    inverse_row_permutation_.resize(n);
+    S_.resize(n, 0, 0);
+    xi_workspace_.resize(2 * n, 0);
+    x_workspace_.resize(n, 0.0);
+    U0_transpose_.resize(1, 1, 1);
+    L0_transpose_.resize(1, 1, 1);
+    clear();
+    reset_stats();
+  }
+
   f_t estimate_solution_density(f_t rhs_nz, f_t sum, i_t& num_calls, bool& use_hypersparse) const
   {
     num_calls++;
@@ -332,13 +380,18 @@ class basis_update_mpf_t {
 
   void multiply_lu(csc_matrix_t<i_t, f_t>& out) const;
 
+  // Compute L*U = A(p, basic_list)
+  int refactor_basis(const csc_matrix_t<i_t, f_t>& A,
+                     const simplex_solver_settings_t<i_t, f_t>& settings,
+                     std::vector<i_t>& basic_list,
+                     std::vector<i_t>& nonbasic_list,
+                     std::vector<variable_status_t>& vstatus);
+
  private:
   void clear()
   {
     pivot_indices_.clear();
     pivot_indices_.reserve(L0_.m);
-    std::iota(col_permutation_.begin(), col_permutation_.end(), 0);
-    std::iota(inverse_col_permutation_.begin(), inverse_col_permutation_.end(), 0);
     S_.col_start.resize(refactor_frequency_ + 1);
     S_.col_start[0] = 0;
     S_.col_start[1] = 0;
@@ -391,8 +444,6 @@ class basis_update_mpf_t {
   std::vector<i_t> pivot_indices_;  // indicies for rank-1 updates to L
   csc_matrix_t<i_t, f_t> S_;        // stores information about the rank-1 updates to L
   std::vector<f_t> mu_values_;      // stores information about the rank-1 updates to L
-  std::vector<i_t> col_permutation_;          // symmetric permuation q used in U(q, q) represents Q
-  std::vector<i_t> inverse_col_permutation_;  // inverse permutation represents Q'
   mutable std::vector<i_t> xi_workspace_;
   mutable std::vector<f_t> x_workspace_;
   mutable csc_matrix_t<i_t, f_t> U0_transpose_;  // Needed for sparse solves
diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp
index c5a5dfe71..098c3b6e2 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];
@@ -2180,6 +2180,43 @@ dual::status_t dual_phase2(i_t phase,
                            lp_solution_t<i_t, f_t>& sol,
                            i_t& iter,
                            std::vector<f_t>& delta_y_steepest_edge)
+{
+  const i_t m = lp.num_rows;
+  const i_t n = lp.num_cols;
+  std::vector<i_t> basic_list(m);
+  std::vector<i_t> nonbasic_list;
+  std::vector<i_t> superbasic_list;
+  basis_update_mpf_t<i_t, f_t> ft(m, settings.refactor_frequency);
+  const bool initialize_basis = true;
+  return dual_phase2_with_advanced_basis(phase,
+                                         slack_basis,
+                                         initialize_basis,
+                                         start_time,
+                                         lp,
+                                         settings,
+                                         vstatus,
+                                         ft,
+                                         basic_list,
+                                         nonbasic_list,
+                                         sol,
+                                         iter,
+                                         delta_y_steepest_edge);
+}
+
+template <typename i_t, typename f_t>
+dual::status_t dual_phase2_with_advanced_basis(i_t phase,
+                                               i_t slack_basis,
+                                               bool initialize_basis,
+                                               f_t start_time,
+                                               const lp_problem_t<i_t, f_t>& lp,
+                                               const simplex_solver_settings_t<i_t, f_t>& settings,
+                                               std::vector<variable_status_t>& vstatus,
+                                               basis_update_mpf_t<i_t, f_t>& ft,
+                                               std::vector<i_t>& basic_list,
+                                               std::vector<i_t>& nonbasic_list,
+                                               lp_solution_t<i_t, f_t>& sol,
+                                               i_t& iter,
+                                               std::vector<f_t>& delta_y_steepest_edge)
 {
   const i_t m = lp.num_rows;
   const i_t n = lp.num_cols;
@@ -2191,9 +2228,6 @@ dual::status_t dual_phase2(i_t phase,
   assert(lp.lower.size() == n);
   assert(lp.upper.size() == n);
   assert(lp.rhs.size() == m);
-  std::vector<i_t> basic_list(m);
-  std::vector<i_t> nonbasic_list;
-  std::vector<i_t> superbasic_list;
 
   std::vector<f_t>& x = sol.x;
   std::vector<f_t>& y = sol.y;
@@ -2209,33 +2243,18 @@ dual::status_t dual_phase2(i_t phase,
   std::vector<f_t> z_old                     = z;
 
   phase2::bound_info(lp, settings);
-  get_basis_from_vstatus(m, vstatus, basic_list, nonbasic_list, superbasic_list);
-  assert(superbasic_list.size() == 0);
-  assert(nonbasic_list.size() == n - m);
-
-  // Compute L*U = A(p, basic_list)
-  csc_matrix_t<i_t, f_t> L(m, m, 1);
-  csc_matrix_t<i_t, f_t> U(m, m, 1);
-  std::vector<i_t> pinv(m);
-  std::vector<i_t> p;
-  std::vector<i_t> q;
-  std::vector<i_t> deficient;
-  std::vector<i_t> slacks_needed;
-
-  if (factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) ==
-      -1) {
-    settings.log.debug("Initial factorization failed\n");
-    basis_repair(lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
-    if (factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) ==
-        -1) {
+  if (initialize_basis) {
+    std::vector<i_t> superbasic_list;
+    get_basis_from_vstatus(m, vstatus, basic_list, nonbasic_list, superbasic_list);
+    assert(superbasic_list.size() == 0);
+    assert(nonbasic_list.size() == n - m);
+
+    if (ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus) > 0) {
       return dual::status_t::NUMERICAL;
     }
-    settings.log.printf("Basis repaired\n");
+
+    if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; }
   }
-  if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; }
-  assert(q.size() == m);
-  reorder_basic_list(q, basic_list);
-  basis_update_mpf_t<i_t, f_t> ft(L, U, p, settings.refactor_frequency);
 
   std::vector<f_t> c_basic(m);
   for (i_t k = 0; k < m; ++k) {
@@ -2872,54 +2891,28 @@ dual::status_t dual_phase2(i_t phase,
 #endif
     if (should_refactor) {
       bool should_recompute_x = false;
-      if (factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) ==
-          -1) {
+      if (ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus) > 0) {
         should_recompute_x = true;
         settings.log.printf("Failed to factorize basis. Iteration %d\n", iter);
         if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; }
-        basis_repair(lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
         i_t count = 0;
-        while (factorize_basis(
-                 lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) == -1) {
+        i_t deficient_size;
+        while ((deficient_size =
+                  ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus)) > 0) {
           settings.log.printf("Failed to repair basis. Iteration %d. %d deficient columns.\n",
                               iter,
-                              static_cast<int>(deficient.size()));
-#ifdef CHECK_L_FACTOR
-          if (L.check_matrix() == -1) { settings.log.printf("Bad L after basis repair\n"); }
-#endif
+                              static_cast<int>(deficient_size));
+
           if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; }
           settings.threshold_partial_pivoting_tol = 1.0;
+
           count++;
           if (count > 10) { return dual::status_t::NUMERICAL; }
-          basis_repair(
-            lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
-
-#ifdef CHECK_BASIS_REPAIR
-          csc_matrix_t<i_t, f_t> B(m, m, 0);
-          form_b(lp.A, basic_list, B);
-          for (i_t k = 0; k < deficient.size(); ++k) {
-            const i_t j         = deficient[k];
-            const i_t col_start = B.col_start[j];
-            const i_t col_end   = B.col_start[j + 1];
-            const i_t col_nz    = col_end - col_start;
-            if (col_nz != 1) {
-              settings.log.printf("Deficient column %d has %d nonzeros\n", j, col_nz);
-            }
-            const i_t i = B.i[col_start];
-            if (i != slacks_needed[k]) {
-              settings.log.printf("Slack %d needed but found %d instead\n", slacks_needed[k], i);
-            }
-          }
-#endif
         }
 
         settings.log.printf("Successfully repaired basis. Iteration %d\n", iter);
       }
-      reorder_basic_list(q, basic_list);
-#ifdef CHECK_L_FACTOR
-      if (L.check_matrix() == -1) { settings.log.printf("Bad L factor\n"); }
-#endif
-      ft.reset(L, U, p);
+
       phase2::reset_basis_mark(basic_list, nonbasic_list, basic_mark, nonbasic_mark);
       if (should_recompute_x) {
         std::vector<f_t> unperturbed_x(n);
@@ -3013,6 +3006,20 @@ template dual::status_t dual_phase2<int, double>(
   int& iter,
   std::vector<double>& steepest_edge_norms);
 
+template dual::status_t dual_phase2_with_advanced_basis<int, double>(
+  int phase,
+  int slack_basis,
+  bool initialize_basis,
+  double start_time,
+  const lp_problem_t<int, double>& lp,
+  const simplex_solver_settings_t<int, double>& settings,
+  std::vector<variable_status_t>& vstatus,
+  basis_update_mpf_t<int, double>& ft,
+  std::vector<int>& basic_list,
+  std::vector<int>& nonbasic_list,
+  lp_solution_t<int, double>& sol,
+  int& iter,
+  std::vector<double>& steepest_edge_norms);
 #endif
 
 }  // namespace cuopt::linear_programming::dual_simplex
diff --git a/cpp/src/dual_simplex/phase2.hpp b/cpp/src/dual_simplex/phase2.hpp
index 39311e060..e237ca372 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,19 @@ 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_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 dadb60cc6..b6f84a9dc 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_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_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,
+  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);
@@ -160,9 +189,9 @@ 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;
+  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");
@@ -179,25 +208,54 @@ lp_status_t solve_linear_program_advanced(const lp_problem_t<i_t, f_t>& original
     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);
+    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);
     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);
+      edge_norms.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);
       vstatus = phase1_vstatus;
       edge_norms.clear();
-      status = dual_phase2(2, 0, start_time, lp, settings, vstatus, 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) {
@@ -602,6 +660,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_advanced_basis(
+  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..53146eda9 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_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,
+  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,