cuPDLPx

benchmarks/linear_programming/cuopt/benchmark_helper.hpp

@@ -134,6 +134,8 @@ void fill_pdlp_hyper_params(const std::string& pdlp_hyper_params_path)
        handle_some_primal_gradients_on_finite_bounds_as_residuals},
     {"project_initial_primal",
      &cuopt::linear_programming::pdlp_hyper_params::project_initial_primal},
+    {"use_adaptive_step_size_strategy",
+     &cuopt::linear_programming::pdlp_hyper_params::use_adaptive_step_size_strategy},
   };
 
   while (std::getline(file, line)) {

benchmarks/linear_programming/cuopt/run_pdlp.cu

@@ -54,10 +54,11 @@ static void parse_arguments(argparse::ArgumentParser& program)
     .default_value(1e-4)
     .scan<'g', double>();
 
+  // TODO replace all comments with Stable2 with Stable3
   program.add_argument("--pdlp-solver-mode")
-    .help("Solver mode for PDLP. Possible values: Stable2 (default), Methodical1, Fast1")
-    .default_value("Stable2")
-    .choices("Stable2", "Methodical1", "Fast1");
+    .help("Solver mode for PDLP. Possible values: Stable3 (default), Methodical1, Fast1")
+    .default_value("Stable3")
+    .choices("Stable3", "Methodical1", "Fast1");
 
   program.add_argument("--method")
     .help(
@@ -90,7 +91,9 @@ static cuopt::linear_programming::pdlp_solver_mode_t string_to_pdlp_solver_mode(
     return cuopt::linear_programming::pdlp_solver_mode_t::Methodical1;
   else if (mode == "Fast1")
     return cuopt::linear_programming::pdlp_solver_mode_t::Fast1;
-  return cuopt::linear_programming::pdlp_solver_mode_t::Stable2;
+  else if (mode == "Stable3")
+    return cuopt::linear_programming::pdlp_solver_mode_t::Stable3;
+  return cuopt::linear_programming::pdlp_solver_mode_t::Stable3;
 }
 
 static cuopt::linear_programming::pdlp_solver_settings_t<int, double> create_solver_settings(

cpp/cuopt_cli.cpp

@@ -122,13 +122,13 @@ int run_single_file(const std::string& file_path,
     (op_problem.get_problem_category() == cuopt::linear_programming::problem_category_t::MIP ||
      op_problem.get_problem_category() == cuopt::linear_programming::problem_category_t::IP);
 
-  auto initial_solution =
-    initial_solution_file.empty()
-      ? std::vector<double>()
-      : cuopt::linear_programming::solution_reader_t::get_variable_values_from_sol_file(
-          initial_solution_file, mps_data_model.get_variable_names());
-
   try {
+    auto initial_solution =
+      initial_solution_file.empty()
+        ? std::vector<double>()
+        : cuopt::linear_programming::solution_reader_t::get_variable_values_from_sol_file(
+            initial_solution_file, mps_data_model.get_variable_names());
+
     if (is_mip && !solve_relaxation) {
       auto& mip_settings = settings.get_mip_settings();
       if (initial_solution.size() > 0) {

cpp/include/cuopt/linear_programming/constants.h

@@ -101,6 +101,7 @@
 #define CUOPT_PDLP_SOLVER_MODE_STABLE2     1
 #define CUOPT_PDLP_SOLVER_MODE_METHODICAL1 2
 #define CUOPT_PDLP_SOLVER_MODE_FAST1       3
+#define CUOPT_PDLP_SOLVER_MODE_STABLE3     4
 
 #define CUOPT_METHOD_CONCURRENT   0
 #define CUOPT_METHOD_PDLP         1

cpp/include/cuopt/linear_programming/pdlp/pdlp_hyper_params.cuh

@@ -56,5 +56,17 @@ extern bool update_primal_weight_on_initial_solution;
 extern bool update_step_size_on_initial_solution;
 extern bool handle_some_primal_gradients_on_finite_bounds_as_residuals;
 extern bool project_initial_primal;
+extern bool use_adaptive_step_size_strategy;
+extern bool initial_step_size_max_singular_value;
+extern bool initial_primal_weight_combined_bounds;
+extern bool bound_objective_rescaling;
+extern bool use_reflected_primal_dual;
+extern bool use_fixed_point_error;
+extern double reflection_coefficient;
+extern double restart_k_p;
+extern double restart_k_i;
+extern double restart_k_d;
+extern double restart_i_smooth;
+extern bool use_conditional_major;
 
 }  // namespace cuopt::linear_programming::pdlp_hyper_params

cpp/include/cuopt/linear_programming/pdlp/solver_settings.hpp

@@ -35,21 +35,22 @@ class solver_settings_t;
  * @brief Enum representing the different solver modes under which PDLP can
  * operate.
  *
- * Stable2: Best overall mode from experiments; balances speed and convergence
- * success. If you want to use the legacy version, use Stable1.
+ * Stable3: Best overall mode from experiments; balances speed and convergence
+ * success. If you want to use the legacy version, use Stable2.
  * Methodical1: Usually leads to slower individual steps but fewer are needed to
  * converge. It uses from 1.3x up to 1.7x times more memory.
  * Fast1: Less convergence success but usually yields the highest speed
  *
- * @note Default mode is Stable2.
+ * @note Default mode is Stable3.
  */
 // Forced to use an enum instead of an enum class for compatibility with the
 // Cython layer
 enum pdlp_solver_mode_t : int {
   Stable1     = CUOPT_PDLP_SOLVER_MODE_STABLE1,
   Stable2     = CUOPT_PDLP_SOLVER_MODE_STABLE2,
   Methodical1 = CUOPT_PDLP_SOLVER_MODE_METHODICAL1,
-  Fast1       = CUOPT_PDLP_SOLVER_MODE_FAST1
+  Fast1       = CUOPT_PDLP_SOLVER_MODE_FAST1,
+  Stable3     = CUOPT_PDLP_SOLVER_MODE_STABLE3
 };
 
 /**
@@ -200,7 +201,7 @@ class pdlp_solver_settings_t {
   bool strict_infeasibility{false};
   i_t iteration_limit{std::numeric_limits<i_t>::max()};
   double time_limit{std::numeric_limits<double>::infinity()};
-  pdlp_solver_mode_t pdlp_solver_mode{pdlp_solver_mode_t::Stable2};
+  pdlp_solver_mode_t pdlp_solver_mode{pdlp_solver_mode_t::Stable3};
   bool log_to_console{true};
   std::string log_file{""};
   std::string sol_file{""};

cpp/src/linear_programming/cusparse_view.cu

@@ -131,7 +131,8 @@ cusparse_view_t<i_t, f_t>::cusparse_view_t(
   saddle_point_state_t<i_t, f_t>& current_saddle_point_state,
   rmm::device_uvector<f_t>& _tmp_primal,
   rmm::device_uvector<f_t>& _tmp_dual,
-  rmm::device_uvector<f_t>& _potential_next_dual_solution)
+  rmm::device_uvector<f_t>& _potential_next_dual_solution,
+  rmm::device_uvector<f_t>& _reflected_primal_solution)
   : handle_ptr_(handle_ptr),
     A{},
     A_T{},
@@ -218,6 +219,12 @@ cusparse_view_t<i_t, f_t>::cusparse_view_t(
     &tmp_primal, op_problem_scaled.n_variables, _tmp_primal.data()));
   RAFT_CUSPARSE_TRY(raft::sparse::detail::cusparsecreatednvec(
     &tmp_dual, op_problem_scaled.n_constraints, _tmp_dual.data()));
+  if (pdlp_hyper_params::use_reflected_primal_dual) {
+    cuopt_assert(_reflected_primal_solution.size() > 0, "Reflected primal solution empty");
+    RAFT_CUSPARSE_TRY(raft::sparse::detail::cusparsecreatednvec(&reflected_primal_solution,
+                                                                op_problem_scaled.n_variables,
+                                                                _reflected_primal_solution.data()));
+  }
 
   const rmm::device_scalar<f_t> alpha{1, handle_ptr->get_stream()};
   const rmm::device_scalar<f_t> beta{1, handle_ptr->get_stream()};
@@ -284,6 +291,8 @@ cusparse_view_t<i_t, f_t>::cusparse_view_t(raft::handle_t const* handle_ptr,
                                            rmm::device_uvector<f_t>& _dual_solution,
                                            rmm::device_uvector<f_t>& _tmp_primal,
                                            rmm::device_uvector<f_t>& _tmp_dual,
+                                           rmm::device_uvector<f_t>& _potential_next_primal,
+                                           rmm::device_uvector<f_t>& _potential_next_dual,
                                            const rmm::device_uvector<f_t>& _A_T,
                                            const rmm::device_uvector<i_t>& _A_T_offsets,
                                            const rmm::device_uvector<i_t>& _A_T_indices)
@@ -335,10 +344,17 @@ cusparse_view_t<i_t, f_t>::cusparse_view_t(raft::handle_t const* handle_ptr,
   RAFT_CUSPARSE_TRY(raft::sparse::detail::cusparsecreatednvec(
     &c, op_problem.n_variables, const_cast<f_t*>(op_problem.objective_coefficients.data())));
 
-  RAFT_CUSPARSE_TRY(raft::sparse::detail::cusparsecreatednvec(
-    &primal_solution, op_problem.n_variables, _primal_solution.data()));
-  RAFT_CUSPARSE_TRY(raft::sparse::detail::cusparsecreatednvec(
-    &dual_solution, op_problem.n_constraints, _dual_solution.data()));
+  if (!pdlp_hyper_params::use_adaptive_step_size_strategy) {
+    RAFT_CUSPARSE_TRY(raft::sparse::detail::cusparsecreatednvec(
+      &primal_solution, op_problem.n_variables, _potential_next_primal.data()));
+    RAFT_CUSPARSE_TRY(raft::sparse::detail::cusparsecreatednvec(
+      &dual_solution, op_problem.n_constraints, _potential_next_dual.data()));
+  } else {
+    RAFT_CUSPARSE_TRY(raft::sparse::detail::cusparsecreatednvec(
+      &primal_solution, op_problem.n_variables, _primal_solution.data()));
+    RAFT_CUSPARSE_TRY(raft::sparse::detail::cusparsecreatednvec(
+      &dual_solution, op_problem.n_constraints, _dual_solution.data()));
+  }
 
   RAFT_CUSPARSE_TRY(raft::sparse::detail::cusparsecreatednvec(
     &tmp_primal, op_problem.n_variables, _tmp_primal.data()));

cpp/src/linear_programming/cusparse_view.hpp

@@ -34,14 +34,17 @@ class cusparse_view_t {
                   saddle_point_state_t<i_t, f_t>& current_saddle_point_state,
                   rmm::device_uvector<f_t>& _tmp_primal,
                   rmm::device_uvector<f_t>& _tmp_dual,
-                  rmm::device_uvector<f_t>& _potential_next_dual_solution);
+                  rmm::device_uvector<f_t>& _potential_next_dual_solution,
+                  rmm::device_uvector<f_t>& _reflected_primal_solution);
 
   cusparse_view_t(raft::handle_t const* handle_ptr,
                   const problem_t<i_t, f_t>& op_problem,
                   rmm::device_uvector<f_t>& _primal_solution,
                   rmm::device_uvector<f_t>& _dual_solution,
                   rmm::device_uvector<f_t>& _tmp_primal,
                   rmm::device_uvector<f_t>& _tmp_dual,
+                  rmm::device_uvector<f_t>& _potential_next_primal,
+                  rmm::device_uvector<f_t>& _potential_next_dual,
                   const rmm::device_uvector<f_t>& _A_T,
                   const rmm::device_uvector<i_t>& _A_T_offsets,
                   const rmm::device_uvector<i_t>& _A_T_indices);
@@ -89,6 +92,9 @@ class cusparse_view_t {
   rmm::device_uvector<uint8_t> buffer_non_transpose;
   rmm::device_uvector<uint8_t> buffer_transpose;
 
+  // Only when using reflection
+  cusparseDnVecDescr_t reflected_primal_solution;
+
   // Ref to the A_T found in either
   // Initial problem, we use it to have an unscaled A_T
   // PDLP copy of the problem which holds the scaled version