Fix bug on consecutive runs caused by static lp_state_t object

cpp/src/mip/diversity/diversity_manager.cu

@@ -337,7 +337,7 @@ solution_t<i_t, f_t> diversity_manager_t<i_t, f_t>::run_solver()
     ls.constraint_prop.bounds_update.probing_cache.probing_cache;
 
   if (check_b_b_preemption()) { return population.best_feasible(); }
-  lp_state_t<i_t, f_t>& lp_state = lp_state_t<i_t, f_t>::get_default_lp_state(*problem_ptr);
+  lp_state_t<i_t, f_t>& lp_state = context.lp_state;
   // resize because some constructor might be called before the presolve
   lp_state.resize(*problem_ptr, problem_ptr->handle_ptr->get_stream());
   auto lp_result = get_relaxed_lp_solution(*problem_ptr,
@@ -546,6 +546,7 @@ diversity_manager_t<i_t, f_t>::recombine_and_local_search(solution_t<i_t, f_t>&
                          lp_offspring,
                          lp_offspring.problem_ptr->integer_indices,
                          context.settings.get_tolerances(),
+                         context.lp_state,
                          lp_run_time);
   cuopt_assert(population.test_invariant(), "");
   cuopt_assert(lp_offspring.test_number_all_integer(), "All must be integers after LP");

cpp/src/mip/diversity/recombiners/fp_recombiner.cuh

@@ -99,7 +99,7 @@ class fp_recombiner_t : public recombiner_t<i_t, f_t> {
     const bool return_first_feasible = false;
     const bool save_state            = false;
     // every sub problem is different,so it is very hard to find a valid initial solution
-    lp_state_t<i_t, f_t> lp_state = lp_state_t<i_t, f_t>::get_default_lp_state(fixed_problem);
+    lp_state_t<i_t, f_t> lp_state = this->context.lp_state;
     auto solver_response          = get_relaxed_lp_solution(fixed_problem,
                                                    fixed_assignment,
                                                    lp_state,

cpp/src/mip/local_search/feasibility_pump/feasibility_pump.cu

@@ -667,6 +667,7 @@ bool feasibility_pump_t<i_t, f_t>::run_single_fp_descent(solution_t<i_t, f_t>& s
                                solution,
                                solution.problem_ptr->integer_indices,
                                context.settings.get_tolerances(),
+                               context.lp_state,
                                lp_verify_time_limit,
                                return_first_feasible,
                                &constraint_prop.bounds_update);

cpp/src/mip/local_search/local_search.cu

@@ -198,6 +198,7 @@ bool local_search_t<i_t, f_t>::check_fj_on_lp_optimal(solution_t<i_t, f_t>& solu
                            solution,
                            solution.problem_ptr->integer_indices,
                            solution.problem_ptr->tolerances,
+                           context.lp_state,
                            lp_run_time);
   } else {
     return is_feasible;

cpp/src/mip/local_search/rounding/constraint_prop.cu

@@ -896,6 +896,7 @@ bool constraint_prop_t<i_t, f_t>::find_integer(
                            orig_sol,
                            orig_sol.problem_ptr->integer_indices,
                            context.settings.get_tolerances(),
+                           context.lp_state,
                            lp_run_time_after_feasible,
                            true);
   }

cpp/src/mip/local_search/rounding/lb_constraint_prop.cu

@@ -949,6 +949,7 @@ bool lb_constraint_prop_t<i_t, f_t>::find_integer(
                            orig_sol,
                            orig_sol.problem_ptr->integer_indices,
                            context.settings.get_tolerances(),
+                           context.lp_state,
                            lp_run_time_after_feasible,
                            true);
   }

cpp/src/mip/relaxed_lp/lp_state.cuh

@@ -24,15 +24,23 @@ namespace cuopt::linear_programming::detail {
 
 template <typename i_t, typename f_t>
 class lp_state_t {
-  // this constructor should be used only once by get_default_lp_state
- private:
+ public:
   lp_state_t(problem_t<i_t, f_t>& problem, rmm::cuda_stream_view stream)
     : prev_primal(problem.n_variables, stream), prev_dual(problem.n_constraints, stream)
   {
+    thrust::fill(problem.handle_ptr->get_thrust_policy(),
+                 prev_primal.data(),
+                 prev_primal.data() + problem.n_variables,
+                 0);
+    thrust::fill(problem.handle_ptr->get_thrust_policy(),
+                 prev_dual.data(),
+                 prev_dual.data() + problem.n_constraints,
+                 0);
   }
 
- public:
-  lp_state_t(problem_t<i_t, f_t>& problem) : lp_state_t(get_default_lp_state(problem)) {}
+  lp_state_t(problem_t<i_t, f_t>& problem) : lp_state_t(problem, problem.handle_ptr->get_stream())
+  {
+  }
 
   lp_state_t(const lp_state_t<i_t, f_t>& other)
     : prev_primal(other.prev_primal, other.prev_primal.stream()),
@@ -43,23 +51,6 @@ class lp_state_t {
   lp_state_t(lp_state_t<i_t, f_t>&& other) noexcept            = default;
   lp_state_t& operator=(lp_state_t<i_t, f_t>&& other) noexcept = default;
 
-  static lp_state_t<i_t, f_t>& get_default_lp_state(problem_t<i_t, f_t>& problem)
-  {
-    if (!default_lp_state) {
-      default_lp_state.reset(new lp_state_t<i_t, f_t>(problem, problem.handle_ptr->get_stream()));
-      thrust::fill(problem.handle_ptr->get_thrust_policy(),
-                   default_lp_state->prev_primal.data(),
-                   default_lp_state->prev_primal.data() + problem.n_variables,
-                   0);
-      thrust::fill(problem.handle_ptr->get_thrust_policy(),
-                   default_lp_state->prev_dual.data(),
-                   default_lp_state->prev_dual.data() + problem.n_constraints,
-                   0);
-    }
-    if (!root_is_initialized) { root_is_initialized = true; }
-    return *default_lp_state;
-  }
-
   void resize(problem_t<i_t, f_t>& problem, rmm::cuda_stream_view stream)
   {
     prev_primal.resize(problem.n_variables, stream);
@@ -79,14 +70,6 @@ class lp_state_t {
   }
   rmm::device_uvector<f_t> prev_primal;
   rmm::device_uvector<f_t> prev_dual;
-  static std::unique_ptr<lp_state_t<i_t, f_t>> default_lp_state;
-  static bool root_is_initialized;
 };
 
-template <typename i_t, typename f_t>
-std::unique_ptr<lp_state_t<i_t, f_t>> lp_state_t<i_t, f_t>::default_lp_state = nullptr;
-
-template <typename i_t, typename f_t>
-bool lp_state_t<i_t, f_t>::root_is_initialized = false;
-
 }  // namespace cuopt::linear_programming::detail

cpp/src/mip/relaxed_lp/relaxed_lp.cu

@@ -135,6 +135,7 @@ bool run_lp_with_vars_fixed(problem_t<i_t, f_t>& op_problem,
                             solution_t<i_t, f_t>& solution,
                             const rmm::device_uvector<i_t>& variables_to_fix,
                             typename mip_solver_settings_t<i_t, f_t>::tolerances_t tols,
+                            lp_state_t<i_t, f_t>& lp_state,
                             f_t time_limit,
                             bool return_first_feasible,
                             bound_presolve_t<i_t, f_t>* bound_presolve)
@@ -160,7 +161,6 @@ bool run_lp_with_vars_fixed(problem_t<i_t, f_t>& op_problem,
   // if we are on the original problem and fixing the integers, save the state
   // if we are in recombiners and on a smaller problem, don't update the state with integers fixed
   bool save_state      = false;
-  auto& lp_state       = lp_state_t<i_t, f_t>::get_default_lp_state(fixed_problem);
   auto solver_response = get_relaxed_lp_solution(fixed_problem,
                                                  fixed_assignment,
                                                  lp_state,
@@ -197,6 +197,7 @@ bool run_lp_with_vars_fixed(problem_t<i_t, f_t>& op_problem,
     solution_t<int, F_TYPE> & solution,                                                       \
     const rmm::device_uvector<int>& variables_to_fix,                                         \
     typename mip_solver_settings_t<int, F_TYPE>::tolerances_t tols,                           \
+    lp_state_t<int, F_TYPE>& lp_state,                                                        \
     F_TYPE time_limit,                                                                        \
     bool return_first_feasible,                                                               \
     bound_presolve_t<int, F_TYPE>* bound_presolve);

cpp/src/mip/relaxed_lp/relaxed_lp.cuh

@@ -51,6 +51,7 @@ bool run_lp_with_vars_fixed(problem_t<i_t, f_t>& op_problem,
                             solution_t<i_t, f_t>& solution,
                             const rmm::device_uvector<i_t>& variables_to_fix,
                             typename mip_solver_settings_t<i_t, f_t>::tolerances_t tols,
+                            lp_state_t<i_t, f_t>& lp_state,
                             f_t time_limit                             = 20.,
                             bool return_first_feasible                 = false,
                             bound_presolve_t<i_t, f_t>* bound_presolve = nullptr);

cpp/src/mip/solver_context.cuh

@@ -17,6 +17,7 @@
 
 #include <linear_programming/initial_scaling_strategy/initial_scaling.cuh>
 #include <mip/problem/problem.cuh>
+#include <mip/relaxed_lp/lp_state.cuh>
 #include <mip/solver_stats.cuh>
 
 #pragma once
@@ -31,14 +32,20 @@ struct mip_solver_context_t {
                                 problem_t<i_t, f_t>* problem_ptr_,
                                 mip_solver_settings_t<i_t, f_t> settings_,
                                 pdlp_initial_scaling_strategy_t<i_t, f_t>& scaling)
-    : handle_ptr(handle_ptr_), problem_ptr(problem_ptr_), settings(settings_), scaling(scaling)
+    : handle_ptr(handle_ptr_),
+      problem_ptr(problem_ptr_),
+      settings(settings_),
+      scaling(scaling),
+      lp_state(*problem_ptr)
   {
+    cuopt_assert(problem_ptr != nullptr, "problem_ptr is nullptr");
     stats.solution_bound = problem_ptr->maximize ? std::numeric_limits<f_t>::infinity()
                                                  : -std::numeric_limits<f_t>::infinity();
   }
 
   raft::handle_t const* const handle_ptr;
   problem_t<i_t, f_t>* problem_ptr;
+  lp_state_t<i_t, f_t> lp_state;
   const mip_solver_settings_t<i_t, f_t> settings;
   pdlp_initial_scaling_strategy_t<i_t, f_t>& scaling;
   solver_stats_t<i_t, f_t> stats;