diff --git a/benchmarks/linear_programming/run_mps_files.sh b/benchmarks/linear_programming/run_mps_files.sh
index 3f911ff806..e409cacaad 100755
--- a/benchmarks/linear_programming/run_mps_files.sh
+++ b/benchmarks/linear_programming/run_mps_files.sh
@@ -209,7 +209,7 @@ OUTPUT_DIR=${OUTPUT_DIR:-.}
 RELAXATION=${RELAXATION:-false}
 MIP_HEURISTICS_ONLY=${MIP_HEURISTICS_ONLY:-false}
 WRITE_LOG_FILE=${WRITE_LOG_FILE:-false}
-NUM_CPU_THREADS=${NUM_CPU_THREADS:-1}
+NUM_CPU_THREADS=${NUM_CPU_THREADS:--1}
 BATCH_NUM=${BATCH_NUM:-0}
 N_BATCHES=${N_BATCHES:-1}
 LOG_TO_CONSOLE=${LOG_TO_CONSOLE:-true}
diff --git a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp
index 2c62f1b443..4f6320752a 100644
--- a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp
+++ b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp
@@ -81,6 +81,7 @@ class mip_solver_settings_t {
   f_t time_limit       = std::numeric_limits<f_t>::infinity();
   bool heuristics_only = false;
   i_t num_cpu_threads  = -1;  // -1 means use default number of threads in branch and bound
+  i_t num_gpus         = 1;
   bool log_to_console  = true;
   std::string log_file;
   std::string sol_file;
diff --git a/cpp/include/cuopt/linear_programming/pdlp/solver_settings.hpp b/cpp/include/cuopt/linear_programming/pdlp/solver_settings.hpp
index 9a2f8aa4c9..ea697e1e10 100644
--- a/cpp/include/cuopt/linear_programming/pdlp/solver_settings.hpp
+++ b/cpp/include/cuopt/linear_programming/pdlp/solver_settings.hpp
@@ -210,8 +210,9 @@ class pdlp_solver_settings_t {
   bool dual_postsolve{true};
   int num_gpus{1};
   method_t method{method_t::Concurrent};
+  bool inside_mip{false};
   // For concurrent termination
-  volatile int* concurrent_halt;
+  volatile int* concurrent_halt{nullptr};
   static constexpr f_t minimal_absolute_tolerance = 1.0e-12;
 
  private:
diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp
index 71cf4ba465..77acca8f7d 100644
--- a/cpp/src/dual_simplex/branch_and_bound.cpp
+++ b/cpp/src/dual_simplex/branch_and_bound.cpp
@@ -9,6 +9,7 @@
 #include <algorithm>
 #include <dual_simplex/bounds_strengthening.hpp>
 #include <dual_simplex/branch_and_bound.hpp>
+#include <dual_simplex/crossover.hpp>
 #include <dual_simplex/initial_basis.hpp>
 #include <dual_simplex/logger.hpp>
 #include <dual_simplex/mip_node.hpp>
@@ -16,7 +17,6 @@
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/pseudo_costs.hpp>
 #include <dual_simplex/random.hpp>
-#include <dual_simplex/solve.hpp>
 #include <dual_simplex/tic_toc.hpp>
 #include <dual_simplex/user_problem.hpp>
 
@@ -24,6 +24,7 @@
 #include <cstdio>
 #include <cstdlib>
 #include <deque>
+#include <future>
 #include <limits>
 #include <optional>
 #include <string>
@@ -217,6 +218,7 @@ branch_and_bound_t<i_t, f_t>::branch_and_bound_t(
     original_lp_(user_problem.handle_ptr, 1, 1, 1),
     incumbent_(1),
     root_relax_soln_(1, 1),
+    root_crossover_soln_(1, 1),
     pc_(1),
     solver_status_(mip_exploration_status_t::UNSET)
 {
@@ -1203,6 +1205,81 @@ void branch_and_bound_t<i_t, f_t>::diving_thread(const csr_matrix_t<i_t, f_t>& A
   }
 }
 
+template <typename i_t, typename f_t>
+lp_status_t branch_and_bound_t<i_t, f_t>::solve_root_relaxation(
+  simplex_solver_settings_t<i_t, f_t> const& lp_settings)
+{
+  // Root node path
+  lp_status_t root_status;
+  std::future<lp_status_t> root_status_future;
+  root_status_future = std::async(std::launch::async,
+                                  &solve_linear_program_advanced<i_t, f_t>,
+                                  std::ref(original_lp_),
+                                  exploration_stats_.start_time,
+                                  std::ref(lp_settings),
+                                  std::ref(root_relax_soln_),
+                                  std::ref(root_vstatus_),
+                                  std::ref(edge_norms_));
+  // Wait for the root relaxation solution to be sent by the diversity manager or dual simplex
+  // to finish
+  while (!root_crossover_solution_set_.load(std::memory_order_acquire) &&
+         *get_root_concurrent_halt() == 0) {
+    std::this_thread::sleep_for(std::chrono::milliseconds(1));
+    continue;
+  }
+
+  if (root_crossover_solution_set_.load(std::memory_order_acquire)) {
+    // Crush the root relaxation solution on converted user problem
+    std::vector<f_t> crushed_root_x;
+    crush_primal_solution(
+      original_problem_, original_lp_, root_crossover_soln_.x, new_slacks_, crushed_root_x);
+    std::vector<f_t> crushed_root_y;
+    std::vector<f_t> crushed_root_z;
+
+    f_t dual_res_inf = crush_dual_solution(original_problem_,
+                                           original_lp_,
+                                           new_slacks_,
+                                           root_crossover_soln_.y,
+                                           root_crossover_soln_.z,
+                                           crushed_root_y,
+                                           crushed_root_z);
+
+    root_crossover_soln_.x = crushed_root_x;
+    root_crossover_soln_.y = crushed_root_y;
+    root_crossover_soln_.z = crushed_root_z;
+
+    // Call crossover on the crushed solution
+    auto root_crossover_settings            = settings_;
+    root_crossover_settings.log.log         = false;
+    root_crossover_settings.concurrent_halt = get_root_concurrent_halt();
+    crossover_status_t crossover_status     = crossover(original_lp_,
+                                                    root_crossover_settings,
+                                                    root_crossover_soln_,
+                                                    exploration_stats_.start_time,
+                                                    root_crossover_soln_,
+                                                    crossover_vstatus_);
+
+    if (crossover_status == crossover_status_t::OPTIMAL) {
+      settings_.log.printf("Crossover status: %d\n", crossover_status);
+    }
+
+    // Check if crossover was stopped by dual simplex
+    if (crossover_status == crossover_status_t::OPTIMAL) {
+      set_root_concurrent_halt(1);  // Stop dual simplex
+      root_status = root_status_future.get();
+      // Override the root relaxation solution with the crossover solution
+      root_relax_soln_ = root_crossover_soln_;
+      root_vstatus_    = crossover_vstatus_;
+      root_status      = lp_status_t::OPTIMAL;
+    } else {
+      root_status = root_status_future.get();
+    }
+  } else {
+    root_status = root_status_future.get();
+  }
+  return root_status;
+}
+
 template <typename i_t, typename f_t>
 mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solution)
 {
@@ -1233,14 +1310,24 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
   root_relax_soln_.resize(original_lp_.num_rows, original_lp_.num_cols);
 
   settings_.log.printf("Solving LP root relaxation\n");
-  simplex_solver_settings_t lp_settings  = settings_;
-  lp_settings.inside_mip                 = 1;
-  lp_status_t root_status                = solve_linear_program_advanced(original_lp_,
-                                                          exploration_stats_.start_time,
-                                                          lp_settings,
-                                                          root_relax_soln_,
-                                                          root_vstatus_,
-                                                          edge_norms_);
+
+  lp_status_t root_status;
+  simplex_solver_settings_t lp_settings = settings_;
+  lp_settings.inside_mip                = 1;
+  lp_settings.concurrent_halt           = get_root_concurrent_halt();
+  // RINS/SUBMIP path
+  if (!enable_concurrent_lp_root_solve()) {
+    root_status = solve_linear_program_advanced(original_lp_,
+                                                exploration_stats_.start_time,
+                                                lp_settings,
+                                                root_relax_soln_,
+                                                root_vstatus_,
+                                                edge_norms_);
+
+  } else {
+    root_status = solve_root_relaxation(lp_settings);
+  }
+
   exploration_stats_.total_lp_iters      = root_relax_soln_.iterations;
   exploration_stats_.total_lp_solve_time = toc(exploration_stats_.start_time);
 
diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp
index 4d0d9c284a..7891711f75 100644
--- a/cpp/src/dual_simplex/branch_and_bound.hpp
+++ b/cpp/src/dual_simplex/branch_and_bound.hpp
@@ -14,7 +14,9 @@
 #include <dual_simplex/pseudo_costs.hpp>
 #include <dual_simplex/simplex_solver_settings.hpp>
 #include <dual_simplex/solution.hpp>
+#include <dual_simplex/solve.hpp>
 #include <dual_simplex/types.hpp>
+#include <utilities/macros.cuh>
 #include <utilities/omp_helpers.hpp>
 
 #include <omp.h>
@@ -78,9 +80,29 @@ class branch_and_bound_t {
   // Set an initial guess based on the user_problem. This should be called before solve.
   void set_initial_guess(const std::vector<f_t>& user_guess) { guess_ = user_guess; }
 
+  // Set the root solution found by PDLP
+  void set_root_relaxation_solution(const std::vector<f_t>& primal,
+                                    const std::vector<f_t>& dual,
+                                    const std::vector<f_t>& reduced_costs,
+                                    f_t objective,
+                                    f_t user_objective,
+                                    i_t iterations)
+  {
+    root_crossover_soln_.x              = primal;
+    root_crossover_soln_.y              = dual;
+    root_crossover_soln_.z              = reduced_costs;
+    root_objective_                     = objective;
+    root_crossover_soln_.objective      = objective;
+    root_crossover_soln_.user_objective = user_objective;
+    root_crossover_soln_.iterations     = iterations;
+    root_crossover_solution_set_.store(true, std::memory_order_release);
+  }
+
   // Set a solution based on the user problem during the course of the solve
   void set_new_solution(const std::vector<f_t>& solution);
 
+  void set_concurrent_lp_root_solve(bool enable) { enable_concurrent_lp_root_solve_ = enable; }
+
   // Repair a low-quality solution from the heuristics.
   bool repair_solution(const std::vector<f_t>& leaf_edge_norms,
                        const std::vector<f_t>& potential_solution,
@@ -90,6 +112,10 @@ class branch_and_bound_t {
   f_t get_upper_bound();
   f_t get_lower_bound();
   i_t get_heap_size();
+  bool enable_concurrent_lp_root_solve() const { return enable_concurrent_lp_root_solve_; }
+  volatile int* get_root_concurrent_halt() { return &root_concurrent_halt_; }
+  void set_root_concurrent_halt(int value) { root_concurrent_halt_ = value; }
+  lp_status_t solve_root_relaxation(simplex_solver_settings_t<i_t, f_t> const& lp_settings);
 
   // The main entry routine. Returns the solver status and populates solution with the incumbent.
   mip_status_t solve(mip_solution_t<i_t, f_t>& solution);
@@ -137,9 +163,14 @@ class branch_and_bound_t {
 
   // Variables for the root node in the search tree.
   std::vector<variable_status_t> root_vstatus_;
+  std::vector<variable_status_t> crossover_vstatus_;
   f_t root_objective_;
   lp_solution_t<i_t, f_t> root_relax_soln_;
+  lp_solution_t<i_t, f_t> root_crossover_soln_;
   std::vector<f_t> edge_norms_;
+  std::atomic<bool> root_crossover_solution_set_{false};
+  bool enable_concurrent_lp_root_solve_{false};
+  volatile int root_concurrent_halt_{0};
 
   // Pseudocosts
   pseudo_costs_t<i_t, f_t> pc_;
diff --git a/cpp/src/dual_simplex/crossover.cpp b/cpp/src/dual_simplex/crossover.cpp
index 2e7bea111a..23d9a0e8e0 100644
--- a/cpp/src/dual_simplex/crossover.cpp
+++ b/cpp/src/dual_simplex/crossover.cpp
@@ -1204,6 +1204,7 @@ crossover_status_t crossover(const lp_problem_t<i_t, f_t>& lp,
       lp, settings, start_time, solution, ft, basic_list, nonbasic_list, superbasic_list, vstatus);
     if (primal_push_status < 0) { return return_to_status(primal_push_status); }
     print_crossover_info(lp, settings, vstatus, solution, "Primal push complete");
+    compute_dual_solution_from_basis(lp, ft, basic_list, nonbasic_list, solution.y, solution.z);
   } else {
     settings.log.printf("No primal push needed. No superbasic variables\n");
   }
@@ -1386,7 +1387,10 @@ crossover_status_t crossover(const lp_problem_t<i_t, f_t>& lp,
   crossover_status_t status = crossover_status_t::NUMERICAL_ISSUES;
   if (dual_feasible) { status = crossover_status_t::DUAL_FEASIBLE; }
   if (primal_feasible) { status = crossover_status_t::PRIMAL_FEASIBLE; }
-  if (primal_feasible && dual_feasible) { status = crossover_status_t::OPTIMAL; }
+  if (primal_feasible && dual_feasible) {
+    status = crossover_status_t::OPTIMAL;
+    if (settings.concurrent_halt != nullptr) { *settings.concurrent_halt = 1; }
+  }
   return status;
 }
 
diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp
index 4932ddae95..56298ef4dd 100644
--- a/cpp/src/dual_simplex/phase2.cpp
+++ b/cpp/src/dual_simplex/phase2.cpp
@@ -2981,6 +2981,10 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
                           100.0 * dense_delta_z / (sparse_delta_z + dense_delta_z));
       ft.print_stats();
     }
+    if (settings.inside_mip && settings.concurrent_halt != nullptr) {
+      settings.log.debug("Setting concurrent halt in Dual Simplex Phase 2\n");
+      *settings.concurrent_halt = 1;
+    }
   }
   return status;
 }
diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp
index 6ae3504c15..d247fbf67c 100644
--- a/cpp/src/dual_simplex/presolve.cpp
+++ b/cpp/src/dual_simplex/presolve.cpp
@@ -508,13 +508,15 @@ i_t find_dependent_rows(lp_problem_t<i_t, f_t>& problem,
 
 template <typename i_t, typename f_t>
 i_t add_artifical_variables(lp_problem_t<i_t, f_t>& problem,
-                            std::vector<i_t>& equality_rows,
+                            const std::vector<i_t>& range_rows,
+                            const std::vector<i_t>& equality_rows,
                             std::vector<i_t>& new_slacks)
 {
-  const i_t n        = problem.num_cols;
-  const i_t m        = problem.num_rows;
-  const i_t num_cols = n + equality_rows.size();
-  const i_t nnz      = problem.A.col_start[n] + equality_rows.size();
+  const i_t n                   = problem.num_cols;
+  const i_t m                   = problem.num_rows;
+  const i_t num_artificial_vars = equality_rows.size() - range_rows.size();
+  const i_t num_cols            = n + num_artificial_vars;
+  i_t nnz                       = problem.A.col_start[n] + num_artificial_vars;
   problem.A.col_start.resize(num_cols + 1);
   problem.A.i.resize(nnz);
   problem.A.x.resize(nnz);
@@ -522,9 +524,15 @@ i_t add_artifical_variables(lp_problem_t<i_t, f_t>& problem,
   problem.upper.resize(num_cols);
   problem.objective.resize(num_cols);
 
+  std::vector<bool> is_range_row(problem.num_rows, false);
+  for (i_t i : range_rows) {
+    is_range_row[i] = true;
+  }
+
   i_t p = problem.A.col_start[n];
   i_t j = n;
   for (i_t i : equality_rows) {
+    if (is_range_row[i]) { continue; }
     // Add an artifical variable z to the equation a_i^T x == b
     // This now becomes a_i^T x + z == b,   0 <= z =< 0
     problem.A.col_start[j] = p;
@@ -541,7 +549,7 @@ i_t add_artifical_variables(lp_problem_t<i_t, f_t>& problem,
   assert(j == num_cols);
   assert(p == nnz);
   constexpr bool verbose = false;
-  if (verbose) { printf("Added %d artificial variables\n", num_cols - n); }
+  if (verbose) { printf("Added %d artificial variables\n", num_artificial_vars); }
   problem.A.n      = num_cols;
   problem.num_cols = num_cols;
   return 0;
@@ -800,7 +808,9 @@ void convert_user_problem(const user_problem_t<i_t, f_t>& user_problem,
   }
 
   // Add artifical variables
-  if (!settings.barrier_presolve) { add_artifical_variables(problem, equality_rows, new_slacks); }
+  if (!settings.barrier_presolve) {
+    add_artifical_variables(problem, user_problem.range_rows, equality_rows, new_slacks);
+  }
 }
 
 template <typename i_t, typename f_t>
@@ -1227,13 +1237,13 @@ void crush_primal_solution_with_slack(const user_problem_t<i_t, f_t>& user_probl
 }
 
 template <typename i_t, typename f_t>
-void crush_dual_solution(const user_problem_t<i_t, f_t>& user_problem,
-                         const lp_problem_t<i_t, f_t>& problem,
-                         const std::vector<i_t>& new_slacks,
-                         const std::vector<f_t>& user_y,
-                         const std::vector<f_t>& user_z,
-                         std::vector<f_t>& y,
-                         std::vector<f_t>& z)
+f_t crush_dual_solution(const user_problem_t<i_t, f_t>& user_problem,
+                        const lp_problem_t<i_t, f_t>& problem,
+                        const std::vector<i_t>& new_slacks,
+                        const std::vector<f_t>& user_y,
+                        const std::vector<f_t>& user_z,
+                        std::vector<f_t>& y,
+                        std::vector<f_t>& z)
 {
   y.resize(problem.num_rows);
   for (i_t i = 0; i < user_problem.num_rows; i++) {
@@ -1244,6 +1254,12 @@ void crush_dual_solution(const user_problem_t<i_t, f_t>& user_problem,
     z[j] = user_z[j];
   }
 
+  std::vector<bool> is_range_row(problem.num_rows, false);
+  for (i_t i = 0; i < user_problem.range_rows.size(); i++) {
+    is_range_row[user_problem.range_rows[i]] = true;
+  }
+  assert(user_problem.num_rows == problem.num_rows);
+
   for (i_t j : new_slacks) {
     const i_t col_start = problem.A.col_start[j];
     const i_t col_end   = problem.A.col_start[j + 1];
@@ -1255,7 +1271,11 @@ void crush_dual_solution(const user_problem_t<i_t, f_t>& user_problem,
     // e_i^T y + z_j = c_j = 0
     // y_i + z_j = 0
     // z_j = - y_i;
-    z[j] = -y[i];
+    if (is_range_row[i]) {
+      z[j] = y[i];
+    } else {
+      z[j] = -y[i];
+    }
   }
 
   // A^T y + z = c or A^T y + z - c = 0
@@ -1292,6 +1312,7 @@ void crush_dual_solution(const user_problem_t<i_t, f_t>& user_problem,
   }
   const f_t dual_res_inf = vector_norm_inf<i_t, f_t>(dual_residual);
   assert(dual_res_inf < 1e-6);
+  return dual_res_inf;
 }
 
 template <typename i_t, typename f_t>
@@ -1511,6 +1532,14 @@ template void crush_primal_solution<int, double>(const user_problem_t<int, doubl
                                                  const std::vector<int>& new_slacks,
                                                  std::vector<double>& solution);
 
+template double crush_dual_solution<int, double>(const user_problem_t<int, double>& user_problem,
+                                                 const lp_problem_t<int, double>& problem,
+                                                 const std::vector<int>& new_slacks,
+                                                 const std::vector<double>& user_y,
+                                                 const std::vector<double>& user_z,
+                                                 std::vector<double>& y,
+                                                 std::vector<double>& z);
+
 template void uncrush_primal_solution<int, double>(const user_problem_t<int, double>& user_problem,
                                                    const lp_problem_t<int, double>& problem,
                                                    const std::vector<double>& solution,
diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp
index 0019c71171..557ebe648b 100644
--- a/cpp/src/dual_simplex/presolve.hpp
+++ b/cpp/src/dual_simplex/presolve.hpp
@@ -150,13 +150,13 @@ void crush_primal_solution_with_slack(const user_problem_t<i_t, f_t>& user_probl
                                       std::vector<f_t>& solution);
 
 template <typename i_t, typename f_t>
-void crush_dual_solution(const user_problem_t<i_t, f_t>& user_problem,
-                         const lp_problem_t<i_t, f_t>& problem,
-                         const std::vector<i_t>& new_slacks,
-                         const std::vector<f_t>& user_y,
-                         const std::vector<f_t>& user_z,
-                         std::vector<f_t>& y,
-                         std::vector<f_t>& z);
+f_t crush_dual_solution(const user_problem_t<i_t, f_t>& user_problem,
+                        const lp_problem_t<i_t, f_t>& problem,
+                        const std::vector<i_t>& new_slacks,
+                        const std::vector<f_t>& user_y,
+                        const std::vector<f_t>& user_z,
+                        std::vector<f_t>& y,
+                        std::vector<f_t>& z);
 
 template <typename i_t, typename f_t>
 void uncrush_primal_solution(const user_problem_t<i_t, f_t>& user_problem,
diff --git a/cpp/src/dual_simplex/right_looking_lu.cpp b/cpp/src/dual_simplex/right_looking_lu.cpp
index 59dac3ac90..a63c1181fa 100644
--- a/cpp/src/dual_simplex/right_looking_lu.cpp
+++ b/cpp/src/dual_simplex/right_looking_lu.cpp
@@ -52,12 +52,12 @@ i_t initialize_degree_data(const csc_matrix_t<i_t, f_t>& A,
   }
 
   for (i_t k = 0; k < n; ++k) {
-    assert(Cdegree[k] <= n && Cdegree[k] >= 0);
+    assert(Cdegree[k] <= m && Cdegree[k] >= 0);
     col_count[Cdegree[k]].push_back(k);
   }
 
   for (i_t k = 0; k < m; ++k) {
-    assert(Rdegree[k] <= m && Rdegree[k] >= 0);
+    assert(Rdegree[k] <= n && Rdegree[k] >= 0);
     row_count[Rdegree[k]].push_back(k);
     if (Rdegree[k] == 0) {
       constexpr bool verbose = false;
diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp
index 919b74e26e..98be9d4cbe 100644
--- a/cpp/src/dual_simplex/simplex_solver_settings.hpp
+++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp
@@ -15,6 +15,7 @@
 #include <atomic>
 #include <functional>
 #include <limits>
+#include <vector>
 
 namespace cuopt::linear_programming::dual_simplex {
 
diff --git a/cpp/src/linear_programming/solve.cu b/cpp/src/linear_programming/solve.cu
index 6fb8003a84..ed141a0c43 100644
--- a/cpp/src/linear_programming/solve.cu
+++ b/cpp/src/linear_programming/solve.cu
@@ -535,6 +535,7 @@ static optimization_problem_solution_t<i_t, f_t> run_pdlp_solver(
                                                      problem.handle_ptr->get_stream()};
   }
   detail::pdlp_solver_t<i_t, f_t> solver(problem, settings, is_batch_mode);
+  if (settings.inside_mip) { solver.set_inside_mip(true); }
   return solver.run_solver(timer);
 }
 
@@ -647,7 +648,6 @@ void run_dual_simplex_thread(
 
 template <typename i_t, typename f_t>
 optimization_problem_solution_t<i_t, f_t> run_concurrent(
-  const optimization_problem_t<i_t, f_t>& op_problem,
   detail::problem_t<i_t, f_t>& problem,
   pdlp_solver_settings_t<i_t, f_t> const& settings,
   const timer_t& timer,
@@ -666,8 +666,8 @@ optimization_problem_solution_t<i_t, f_t> run_concurrent(
   // Make sure allocations are done on the original stream
   problem.handle_ptr->sync_stream();
 
-  int device_count = raft::device_setter::get_device_count();
   if (settings.num_gpus > 1) {
+    int device_count = raft::device_setter::get_device_count();
     CUOPT_LOG_INFO("Running PDLP and Barrier on %d GPUs", device_count);
     cuopt_expects(
       device_count > 1, error_type_t::RuntimeError, "Multi-GPU mode requires at least 2 GPUs");
@@ -682,12 +682,15 @@ optimization_problem_solution_t<i_t, f_t> run_concurrent(
   std::unique_ptr<
     std::tuple<dual_simplex::lp_solution_t<i_t, f_t>, dual_simplex::lp_status_t, f_t, f_t, f_t>>
     sol_dual_simplex_ptr;
-  std::thread dual_simplex_thread(run_dual_simplex_thread<i_t, f_t>,
-                                  std::ref(dual_simplex_problem),
-                                  std::ref(settings_pdlp),
-                                  std::ref(sol_dual_simplex_ptr),
-                                  std::ref(timer));
-
+  std::thread dual_simplex_thread;
+  if (!settings.inside_mip) {
+    dual_simplex_thread = std::thread(run_dual_simplex_thread<i_t, f_t>,
+                                      std::ref(dual_simplex_problem),
+                                      std::ref(settings_pdlp),
+                                      std::ref(sol_dual_simplex_ptr),
+                                      std::ref(timer));
+  }
+  dual_simplex::user_problem_t<i_t, f_t> barrier_problem = dual_simplex_problem;
   // Create a thread for barrier
   std::unique_ptr<
     std::tuple<dual_simplex::lp_solution_t<i_t, f_t>, dual_simplex::lp_status_t, f_t, f_t, f_t>>
@@ -722,18 +725,22 @@ optimization_problem_solution_t<i_t, f_t> run_concurrent(
   auto sol_pdlp = run_pdlp(problem, settings_pdlp, timer, is_batch_mode);
 
   // Wait for dual simplex thread to finish
-  dual_simplex_thread.join();
+  if (!settings.inside_mip) { dual_simplex_thread.join(); }
 
   barrier_thread.join();
 
   // copy the dual simplex solution to the device
-  auto sol_dual_simplex = convert_dual_simplex_sol(problem,
-                                                   std::get<0>(*sol_dual_simplex_ptr),
-                                                   std::get<1>(*sol_dual_simplex_ptr),
-                                                   std::get<2>(*sol_dual_simplex_ptr),
-                                                   std::get<3>(*sol_dual_simplex_ptr),
-                                                   std::get<4>(*sol_dual_simplex_ptr),
-                                                   0);
+  auto sol_dual_simplex =
+    !settings.inside_mip
+      ? convert_dual_simplex_sol(problem,
+                                 std::get<0>(*sol_dual_simplex_ptr),
+                                 std::get<1>(*sol_dual_simplex_ptr),
+                                 std::get<2>(*sol_dual_simplex_ptr),
+                                 std::get<3>(*sol_dual_simplex_ptr),
+                                 std::get<4>(*sol_dual_simplex_ptr),
+                                 0)
+      : optimization_problem_solution_t<i_t, f_t>{pdlp_termination_status_t::ConcurrentLimit,
+                                                  problem.handle_ptr->get_stream()};
 
   // copy the barrier solution to the device
   auto sol_barrier = convert_dual_simplex_sol(problem,
@@ -748,11 +755,12 @@ optimization_problem_solution_t<i_t, f_t> run_concurrent(
   CUOPT_LOG_INFO(
     "Concurrent time:  %.3fs, total time %.3fs", timer_concurrent.elapsed_time(), end_time);
   // Check status to see if we should return the pdlp solution or the dual simplex solution
-  if (sol_dual_simplex.get_termination_status() == pdlp_termination_status_t::Optimal ||
-      sol_dual_simplex.get_termination_status() == pdlp_termination_status_t::PrimalInfeasible ||
-      sol_dual_simplex.get_termination_status() == pdlp_termination_status_t::DualInfeasible) {
+  if (!settings.inside_mip &&
+      (sol_dual_simplex.get_termination_status() == pdlp_termination_status_t::Optimal ||
+       sol_dual_simplex.get_termination_status() == pdlp_termination_status_t::PrimalInfeasible ||
+       sol_dual_simplex.get_termination_status() == pdlp_termination_status_t::DualInfeasible)) {
     CUOPT_LOG_INFO("Solved with dual simplex");
-    sol_pdlp.copy_from(op_problem.get_handle_ptr(), sol_dual_simplex);
+    sol_pdlp.copy_from(problem.handle_ptr, sol_dual_simplex);
     sol_pdlp.set_solve_time(end_time);
     CUOPT_LOG_INFO("Status: %s   Objective: %.8e  Iterations: %d  Time: %.3fs",
                    sol_pdlp.get_termination_status_string().c_str(),
@@ -762,7 +770,7 @@ optimization_problem_solution_t<i_t, f_t> run_concurrent(
     return sol_pdlp;
   } else if (sol_barrier.get_termination_status() == pdlp_termination_status_t::Optimal) {
     CUOPT_LOG_INFO("Solved with barrier");
-    sol_pdlp.copy_from(op_problem.get_handle_ptr(), sol_barrier);
+    sol_pdlp.copy_from(problem.handle_ptr, sol_barrier);
     sol_pdlp.set_solve_time(end_time);
     CUOPT_LOG_INFO("Status: %s   Objective: %.8e  Iterations: %d  Time: %.3fs",
                    sol_pdlp.get_termination_status_string().c_str(),
@@ -773,7 +781,8 @@ optimization_problem_solution_t<i_t, f_t> run_concurrent(
   } else if (sol_pdlp.get_termination_status() == pdlp_termination_status_t::Optimal) {
     CUOPT_LOG_INFO("Solved with PDLP");
     return sol_pdlp;
-  } else if (sol_pdlp.get_termination_status() == pdlp_termination_status_t::ConcurrentLimit) {
+  } else if (!settings.inside_mip &&
+             sol_pdlp.get_termination_status() == pdlp_termination_status_t::ConcurrentLimit) {
     CUOPT_LOG_INFO("Using dual simplex solve info");
     return sol_dual_simplex;
   } else {
@@ -784,7 +793,6 @@ optimization_problem_solution_t<i_t, f_t> run_concurrent(
 
 template <typename i_t, typename f_t>
 optimization_problem_solution_t<i_t, f_t> solve_lp_with_method(
-  const optimization_problem_t<i_t, f_t>& op_problem,
   detail::problem_t<i_t, f_t>& problem,
   pdlp_solver_settings_t<i_t, f_t> const& settings,
   const timer_t& timer,
@@ -795,7 +803,7 @@ optimization_problem_solution_t<i_t, f_t> solve_lp_with_method(
   } else if (settings.method == method_t::Barrier) {
     return run_barrier(problem, settings, timer);
   } else if (settings.method == method_t::Concurrent) {
-    return run_concurrent(op_problem, problem, settings, timer, is_batch_mode);
+    return run_concurrent(problem, settings, timer, is_batch_mode);
   } else {
     return run_pdlp(problem, settings, timer, is_batch_mode);
   }
@@ -900,7 +908,7 @@ optimization_problem_solution_t<i_t, f_t> solve_lp(
 
     setup_device_symbols(op_problem.get_handle_ptr()->get_stream());
 
-    auto solution = solve_lp_with_method(op_problem, problem, settings, lp_timer, is_batch_mode);
+    auto solution = solve_lp_with_method(problem, settings, lp_timer, is_batch_mode);
 
     if (run_presolve) {
       auto primal_solution = cuopt::device_copy(solution.get_primal_solution(),
@@ -1068,7 +1076,6 @@ optimization_problem_solution_t<i_t, f_t> solve_lp(
     bool use_pdlp_solver_mode);                                                        \
                                                                                        \
   template optimization_problem_solution_t<int, F_TYPE> solve_lp_with_method(          \
-    const optimization_problem_t<int, F_TYPE>& op_problem,                             \
     detail::problem_t<int, F_TYPE>& problem,                                           \
     pdlp_solver_settings_t<int, F_TYPE> const& settings,                               \
     const timer_t& timer,                                                              \
diff --git a/cpp/src/linear_programming/solve.cuh b/cpp/src/linear_programming/solve.cuh
index 83d7529c9b..a3c3240f40 100644
--- a/cpp/src/linear_programming/solve.cuh
+++ b/cpp/src/linear_programming/solve.cuh
@@ -22,7 +22,6 @@ cuopt::linear_programming::optimization_problem_t<i_t, f_t> mps_data_model_to_op
 
 template <typename i_t, typename f_t>
 cuopt::linear_programming::optimization_problem_solution_t<i_t, f_t> solve_lp_with_method(
-  const optimization_problem_t<i_t, f_t>& op_problem,
   detail::problem_t<i_t, f_t>& problem,
   pdlp_solver_settings_t<i_t, f_t> const& settings,
   const timer_t& timer,
diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu
index a8cdf9c7a2..4e3dc64650 100644
--- a/cpp/src/math_optimization/solver_settings.cu
+++ b/cpp/src/math_optimization/solver_settings.cu
@@ -87,7 +87,8 @@ solver_settings_t<i_t, f_t>::solver_settings_t() : pdlp_settings(), mip_settings
     {CUOPT_DUALIZE, &pdlp_settings.dualize, -1, 1, -1},
     {CUOPT_ORDERING, &pdlp_settings.ordering, -1, 1, -1},
     {CUOPT_BARRIER_DUAL_INITIAL_POINT, &pdlp_settings.barrier_dual_initial_point, -1, 1, -1},
-    {CUOPT_NUM_GPUS, &pdlp_settings.num_gpus, 1, 2, 1}
+    {CUOPT_NUM_GPUS, &pdlp_settings.num_gpus, 1, 2, 1},
+    {CUOPT_NUM_GPUS, &mip_settings.num_gpus, 1, 2, 1}
   };
 
     // Bool parameters
diff --git a/cpp/src/mip/diversity/diversity_manager.cu b/cpp/src/mip/diversity/diversity_manager.cu
index 1d317a6ea8..686951ebba 100644
--- a/cpp/src/mip/diversity/diversity_manager.cu
+++ b/cpp/src/mip/diversity/diversity_manager.cu
@@ -5,15 +5,17 @@
  */
 /* clang-format on */
 
+#include "cuda_profiler_api.h"
+#include "diversity_manager.cuh"
+
 #include <mip/mip_constants.hpp>
 #include <mip/presolve/probing_cache.cuh>
 #include <mip/presolve/trivial_presolve.cuh>
 #include <mip/problem/problem_helpers.cuh>
-#include "diversity_manager.cuh"
 
-#include <utilities/scope_guard.hpp>
+#include <linear_programming/solve.cuh>
 
-#include "cuda_profiler_api.h"
+#include <utilities/scope_guard.hpp>
 
 constexpr bool fj_only_run = false;
 
@@ -34,6 +36,7 @@ std::vector<recombiner_enum_t> recombiner_t<i_t, f_t>::enabled_recombiners;
 template <typename i_t, typename f_t>
 diversity_manager_t<i_t, f_t>::diversity_manager_t(mip_solver_context_t<i_t, f_t>& context_)
   : context(context_),
+    branch_and_bound_ptr(nullptr),
     problem_ptr(context.problem_ptr),
     diversity_config(),
     population("population",
@@ -336,17 +339,29 @@ solution_t<i_t, f_t> diversity_manager_t<i_t, f_t>::run_solver()
   if (bb_thread_solution_exists) {
     ls.lp_optimal_exists = true;
   } else if (!fj_only_run) {
-    relaxed_lp_settings_t lp_settings;
-    lp_settings.time_limit            = lp_time_limit;
-    lp_settings.tolerance             = context.settings.tolerances.absolute_tolerance;
-    lp_settings.return_first_feasible = false;
-    lp_settings.save_state            = true;
-    lp_settings.concurrent_halt       = &global_concurrent_halt;
-    lp_settings.has_initial_primal    = false;
+    convert_greater_to_less(*problem_ptr);
+
+    f_t tolerance_divisor =
+      problem_ptr->tolerances.absolute_tolerance / problem_ptr->tolerances.relative_tolerance;
+    if (tolerance_divisor == 0) { tolerance_divisor = 1; }
+    f_t absolute_tolerance = context.settings.tolerances.absolute_tolerance;
+
+    pdlp_solver_settings_t<i_t, f_t> pdlp_settings{};
+    pdlp_settings.tolerances.relative_primal_tolerance = absolute_tolerance / tolerance_divisor;
+    pdlp_settings.tolerances.relative_dual_tolerance   = absolute_tolerance / tolerance_divisor;
+    pdlp_settings.time_limit                           = lp_time_limit;
+    pdlp_settings.first_primal_feasible                = false;
+    pdlp_settings.concurrent_halt                      = &global_concurrent_halt;
+    pdlp_settings.method                               = method_t::Concurrent;
+    pdlp_settings.inside_mip                           = true;
+    pdlp_settings.pdlp_solver_mode                     = pdlp_solver_mode_t::Stable2;
+    pdlp_settings.num_gpus                             = context.settings.num_gpus;
+
     rmm::device_uvector<f_t> lp_optimal_solution_copy(lp_optimal_solution.size(),
                                                       problem_ptr->handle_ptr->get_stream());
-    auto lp_result =
-      get_relaxed_lp_solution(*problem_ptr, lp_optimal_solution_copy, lp_state, lp_settings);
+    timer_t lp_timer(lp_time_limit);
+    auto lp_result = solve_lp_with_method<i_t, f_t>(*problem_ptr, pdlp_settings, lp_timer);
+
     {
       std::lock_guard<std::mutex> guard(relaxed_solution_mutex);
       if (!simplex_solution_exists.load()) {
@@ -382,6 +397,41 @@ solution_t<i_t, f_t> diversity_manager_t<i_t, f_t>::run_solver()
       // note to developer, in debug mode the LP run might be too slow and it might cause PDLP not
       // to bring variables within the bounds
     }
+
+    // Send PDLP relaxed solution to branch and bound before it solves the root node
+    if (problem_ptr->set_root_relaxation_solution_callback != nullptr) {
+      auto& d_primal_solution = lp_result.get_primal_solution();
+      auto& d_dual_solution   = lp_result.get_dual_solution();
+      auto& d_reduced_costs   = lp_result.get_reduced_cost();
+
+      std::vector<f_t> host_primal(d_primal_solution.size());
+      std::vector<f_t> host_dual(d_dual_solution.size());
+      std::vector<f_t> host_reduced_costs(d_reduced_costs.size());
+      raft::copy(host_primal.data(),
+                 d_primal_solution.data(),
+                 d_primal_solution.size(),
+                 problem_ptr->handle_ptr->get_stream());
+      raft::copy(host_dual.data(),
+                 d_dual_solution.data(),
+                 d_dual_solution.size(),
+                 problem_ptr->handle_ptr->get_stream());
+      raft::copy(host_reduced_costs.data(),
+                 d_reduced_costs.data(),
+                 d_reduced_costs.size(),
+                 problem_ptr->handle_ptr->get_stream());
+      problem_ptr->handle_ptr->sync_stream();
+
+      auto user_obj   = problem_ptr->get_user_obj_from_solver_obj(lp_result.get_objective_value());
+      auto iterations = lp_result.get_additional_termination_information().number_of_steps_taken;
+      // Set for the B&B
+      problem_ptr->set_root_relaxation_solution_callback(host_primal,
+                                                         host_dual,
+                                                         host_reduced_costs,
+                                                         lp_result.get_objective_value(),
+                                                         user_obj,
+                                                         iterations);
+    }
+
     // in case the pdlp returned var boudns that are out of bounds
     clamp_within_var_bounds(lp_optimal_solution, problem_ptr, problem_ptr->handle_ptr);
   }
@@ -738,6 +788,7 @@ void diversity_manager_t<i_t, f_t>::set_simplex_solution(const std::vector<f_t>&
   std::lock_guard<std::mutex> lock(relaxed_solution_mutex);
   simplex_solution_exists.store(true, std::memory_order_release);
   global_concurrent_halt = 1;
+  CUOPT_LOG_DEBUG("Setting concurrent halt for PDLP inside diversity manager");
   // global_concurrent_halt.store(1, std::memory_order_release);
   // it is safe to use lp_optimal_solution while executing the copy operation
   // the operations are ordered as long as they are on the same stream
diff --git a/cpp/src/mip/diversity/diversity_manager.cuh b/cpp/src/mip/diversity/diversity_manager.cuh
index ab11e4cc22..4a78f6cff7 100644
--- a/cpp/src/mip/diversity/diversity_manager.cuh
+++ b/cpp/src/mip/diversity/diversity_manager.cuh
@@ -71,6 +71,7 @@ class diversity_manager_t {
                             f_t objective);
 
   mip_solver_context_t<i_t, f_t>& context;
+  dual_simplex::branch_and_bound_t<i_t, f_t>* branch_and_bound_ptr;
   problem_t<i_t, f_t>* problem_ptr;
   diversity_config_t diversity_config;
   population_t<i_t, f_t> population;
diff --git a/cpp/src/mip/diversity/lns/rins.cu b/cpp/src/mip/diversity/lns/rins.cu
index c2e46da4d5..f4381507a5 100644
--- a/cpp/src/mip/diversity/lns/rins.cu
+++ b/cpp/src/mip/diversity/lns/rins.cu
@@ -23,7 +23,6 @@
 #include <mip/presolve/trivial_presolve.cuh>
 
 namespace cuopt::linear_programming::detail {
-
 template <typename i_t, typename f_t>
 rins_t<i_t, f_t>::rins_t(mip_solver_context_t<i_t, f_t>& context_,
                          diversity_manager_t<i_t, f_t>& dm_,
diff --git a/cpp/src/mip/diversity/lns/rins.cuh b/cpp/src/mip/diversity/lns/rins.cuh
index 2194e9ee92..3d16875ce5 100644
--- a/cpp/src/mip/diversity/lns/rins.cuh
+++ b/cpp/src/mip/diversity/lns/rins.cuh
@@ -21,6 +21,7 @@
 #include <mip/solution/solution.cuh>
 #include <mip/solver.cuh>
 #include <mip/utilities/cpu_worker_thread.cuh>
+
 #include <utilities/timer.hpp>
 
 #include <atomic>
diff --git a/cpp/src/mip/problem/problem.cu b/cpp/src/mip/problem/problem.cu
index 7ca87b5802..8d251768c9 100644
--- a/cpp/src/mip/problem/problem.cu
+++ b/cpp/src/mip/problem/problem.cu
@@ -143,7 +143,8 @@ problem_t<i_t, f_t>::problem_t(
     Q_values(problem_.get_quadratic_objective_values())
 {
   op_problem_cstr_body(problem_);
-  branch_and_bound_callback = nullptr;
+  branch_and_bound_callback             = nullptr;
+  set_root_relaxation_solution_callback = nullptr;
 }
 
 template <typename i_t, typename f_t>
@@ -154,6 +155,7 @@ problem_t<i_t, f_t>::problem_t(const problem_t<i_t, f_t>& problem_)
     integer_fixed_problem(problem_.integer_fixed_problem),
     integer_fixed_variable_map(problem_.integer_fixed_variable_map, handle_ptr->get_stream()),
     branch_and_bound_callback(nullptr),
+    set_root_relaxation_solution_callback(nullptr),
     n_variables(problem_.n_variables),
     n_constraints(problem_.n_constraints),
     n_binary_vars(problem_.n_binary_vars),
diff --git a/cpp/src/mip/problem/problem.cuh b/cpp/src/mip/problem/problem.cuh
index 47d491dacb..ed0adb971f 100644
--- a/cpp/src/mip/problem/problem.cuh
+++ b/cpp/src/mip/problem/problem.cuh
@@ -207,6 +207,9 @@ class problem_t {
   rmm::device_uvector<i_t> integer_fixed_variable_map;
 
   std::function<void(const std::vector<f_t>&)> branch_and_bound_callback;
+  std::function<void(
+    const std::vector<f_t>&, const std::vector<f_t>&, const std::vector<f_t>&, f_t, f_t, i_t)>
+    set_root_relaxation_solution_callback;
 
   typename mip_solver_settings_t<i_t, f_t>::tolerances_t tolerances{};
   i_t n_variables{0};
diff --git a/cpp/src/mip/problem/problem_helpers.cuh b/cpp/src/mip/problem/problem_helpers.cuh
index 13f5d24ea4..517b52a1d9 100644
--- a/cpp/src/mip/problem/problem_helpers.cuh
+++ b/cpp/src/mip/problem/problem_helpers.cuh
@@ -16,6 +16,7 @@
 #include <utilities/copy_helpers.hpp>
 
 #include <cuda_runtime_api.h>
+#include <thrust/count.h>
 #include <thrust/functional.h>
 #include <thrust/logical.h>
 #include <thrust/sort.h>
@@ -311,4 +312,54 @@ static bool check_bounds_sanity(const detail::problem_t<i_t, f_t>& problem)
          check_constraint_bounds_sanity<i_t, f_t>(problem);
 }
 
+template <typename i_t, typename f_t>
+__global__ void kernel_convert_greater_to_less(raft::device_span<f_t> coefficients,
+                                               raft::device_span<const i_t> offsets,
+                                               raft::device_span<f_t> constraint_lower_bounds,
+                                               raft::device_span<f_t> constraint_upper_bounds)
+{
+  const i_t constraint_id = blockIdx.x;
+
+  const f_t lb = constraint_lower_bounds[constraint_id];
+  const f_t ub = constraint_upper_bounds[constraint_id];
+
+  if (!isfinite(lb) || isfinite(ub)) return;
+
+  auto row_start = offsets[constraint_id];
+  auto row_end   = offsets[constraint_id + 1];
+  auto row_size  = row_end - row_start;
+
+  for (i_t tid = threadIdx.x; tid < row_size; tid += blockDim.x) {
+    coefficients[row_start + tid] = -coefficients[row_start + tid];
+  }
+
+  if (threadIdx.x == 0) {
+    constraint_lower_bounds[constraint_id] = -ub;
+    constraint_upper_bounds[constraint_id] = -lb;
+  }
+}
+
+template <typename i_t, typename f_t>
+static void convert_greater_to_less(detail::problem_t<i_t, f_t>& problem)
+{
+  raft::common::nvtx::range scope("convert_greater_to_less");
+
+  auto* handle_ptr = problem.handle_ptr;
+
+  constexpr i_t TPB = 256;
+  kernel_convert_greater_to_less<i_t, f_t>
+    <<<problem.n_constraints, TPB, 0, handle_ptr->get_stream()>>>(
+      raft::device_span<f_t>(problem.coefficients.data(), problem.coefficients.size()),
+      raft::device_span<const i_t>(problem.offsets.data(), problem.offsets.size()),
+      raft::device_span<f_t>(problem.constraint_lower_bounds.data(),
+                             problem.constraint_lower_bounds.size()),
+      raft::device_span<f_t>(problem.constraint_upper_bounds.data(),
+                             problem.constraint_upper_bounds.size()));
+  RAFT_CHECK_CUDA(handle_ptr->get_stream());
+
+  problem.compute_transpose_of_problem();
+
+  handle_ptr->sync_stream();
+}
+
 }  // namespace cuopt::linear_programming::detail
diff --git a/cpp/src/mip/solver.cu b/cpp/src/mip/solver.cu
index f5d8c9f301..0da4c6398f 100644
--- a/cpp/src/mip/solver.cu
+++ b/cpp/src/mip/solver.cu
@@ -133,8 +133,7 @@ solution_t<i_t, f_t> mip_solver_t<i_t, f_t>::run_solver()
     auto lp_timer       = timer_t(settings.time_limit);
     settings.method     = method_t::Concurrent;
 
-    auto opt_sol = solve_lp_with_method<i_t, f_t>(
-      *context.problem_ptr->original_problem_ptr, *context.problem_ptr, settings, lp_timer);
+    auto opt_sol = solve_lp_with_method<i_t, f_t>(*context.problem_ptr, settings, lp_timer);
 
     solution_t<i_t, f_t> sol(*context.problem_ptr);
     sol.copy_new_assignment(host_copy(opt_sol.get_primal_solution()));
@@ -207,12 +206,22 @@ solution_t<i_t, f_t> mip_solver_t<i_t, f_t>::run_solver()
     branch_and_bound = std::make_unique<dual_simplex::branch_and_bound_t<i_t, f_t>>(
       branch_and_bound_problem, branch_and_bound_settings);
     context.branch_and_bound_ptr = branch_and_bound.get();
+    branch_and_bound->set_concurrent_lp_root_solve(true);
 
     // Set the primal heuristics -> branch and bound callback
     context.problem_ptr->branch_and_bound_callback =
       std::bind(&dual_simplex::branch_and_bound_t<i_t, f_t>::set_new_solution,
                 branch_and_bound.get(),
                 std::placeholders::_1);
+    context.problem_ptr->set_root_relaxation_solution_callback =
+      std::bind(&dual_simplex::branch_and_bound_t<i_t, f_t>::set_root_relaxation_solution,
+                branch_and_bound.get(),
+                std::placeholders::_1,
+                std::placeholders::_2,
+                std::placeholders::_3,
+                std::placeholders::_4,
+                std::placeholders::_5,
+                std::placeholders::_6);
 
     // Fork a thread for branch and bound
     // std::async and std::future allow us to get the return value of bb::solve()