Add support for FP32 and mixed precision in PDLP

benchmarks/linear_programming/cuopt/run_pdlp.cu

@@ -76,6 +76,20 @@ static void parse_arguments(argparse::ArgumentParser& program)
     .choices("None", "Papilo", "PSLP", "Default");
 
   program.add_argument("--solution-path").help("Path where solution file will be generated");
+
+  program.add_argument("--pdlp-fp32")
+    .help(
+      "Use FP32 (float) precision instead of FP64 (double). Only supported for PDLP method without "
+      "crossover.")
+    .default_value(false)
+    .implicit_value(true);
+
+  program.add_argument("--mixed-precision-spmv")
+    .help(
+      "Enable mixed precision SpMV (FP32 matrix, FP64 vectors) during PDHG iterations. Only "
+      "supported for PDLP method in FP64.")
+    .default_value(false)
+    .implicit_value(true);
 }
 
 static cuopt::linear_programming::presolver_t string_to_presolver(const std::string& presolver)
@@ -102,41 +116,33 @@ static cuopt::linear_programming::pdlp_solver_mode_t string_to_pdlp_solver_mode(
   return cuopt::linear_programming::pdlp_solver_mode_t::Stable3;
 }
 
-static cuopt::linear_programming::pdlp_solver_settings_t<int, double> create_solver_settings(
+template <typename f_t>
+static cuopt::linear_programming::pdlp_solver_settings_t<int, f_t> create_solver_settings(
   const argparse::ArgumentParser& program)
 {
-  cuopt::linear_programming::pdlp_solver_settings_t<int, double> settings =
-    cuopt::linear_programming::pdlp_solver_settings_t<int, double>{};
+  cuopt::linear_programming::pdlp_solver_settings_t<int, f_t> settings =
+    cuopt::linear_programming::pdlp_solver_settings_t<int, f_t>{};
 
-  settings.time_limit      = program.get<double>("--time-limit");
+  settings.time_limit      = static_cast<f_t>(program.get<double>("--time-limit"));
   settings.iteration_limit = program.get<int>("--iteration-limit");
-  settings.set_optimality_tolerance(program.get<double>("--optimality-tolerance"));
+  settings.set_optimality_tolerance(
+    static_cast<f_t>(program.get<double>("--optimality-tolerance")));
   settings.pdlp_solver_mode =
     string_to_pdlp_solver_mode(program.get<std::string>("--pdlp-solver-mode"));
   settings.method = static_cast<cuopt::linear_programming::method_t>(program.get<int>("--method"));
-  settings.crossover = program.get<int>("--crossover");
-  settings.presolver = string_to_presolver(program.get<std::string>("--presolver"));
+  settings.crossover            = program.get<int>("--crossover");
+  settings.presolver            = string_to_presolver(program.get<std::string>("--presolver"));
+  settings.mixed_precision_spmv = program.get<bool>("--mixed-precision-spmv");
 
   return settings;
 }
 
-int main(int argc, char* argv[])
+template <typename f_t>
+static int run_solver(const argparse::ArgumentParser& program, const raft::handle_t& handle_)
 {
-  // Parse binary arguments
-  argparse::ArgumentParser program("solve_LP");
-  parse_arguments(program);
-
-  try {
-    program.parse_args(argc, argv);
-  } catch (const std::runtime_error& err) {
-    std::cerr << err.what() << std::endl;
-    std::cerr << program;
-    return 1;
-  }
-
   // Initialize solver settings from binary arguments
-  cuopt::linear_programming::pdlp_solver_settings_t<int, double> settings =
-    create_solver_settings(program);
+  cuopt::linear_programming::pdlp_solver_settings_t<int, f_t> settings =
+    create_solver_settings<f_t>(program);
 
   bool use_pdlp_solver_mode = true;
   if (program.is_used("--pdlp-hyper-params-path")) {
@@ -145,20 +151,13 @@ int main(int argc, char* argv[])
     use_pdlp_solver_mode = false;
   }
 
-  // Setup up RMM memory pool
-  auto memory_resource = make_pool();
-  rmm::mr::set_current_device_resource(memory_resource.get());
-
-  // Initialize raft handle and running stream
-  const raft::handle_t handle_{};
-
   // Parse MPS file
-  cuopt::mps_parser::mps_data_model_t<int, double> op_problem =
-    cuopt::mps_parser::parse_mps<int, double>(program.get<std::string>("--path"));
+  cuopt::mps_parser::mps_data_model_t<int, f_t> op_problem =
+    cuopt::mps_parser::parse_mps<int, f_t>(program.get<std::string>("--path"));
 
   // Solve LP problem
   bool problem_checking = true;
-  cuopt::linear_programming::optimization_problem_solution_t<int, double> solution =
+  cuopt::linear_programming::optimization_problem_solution_t<int, f_t> solution =
     cuopt::linear_programming::solve_lp(
       &handle_, op_problem, settings, problem_checking, use_pdlp_solver_mode);
 
@@ -168,3 +167,33 @@ int main(int argc, char* argv[])
 
   return 0;
 }
+
+int main(int argc, char* argv[])
+{
+  // Parse binary arguments
+  argparse::ArgumentParser program("solve_LP");
+  parse_arguments(program);
+
+  try {
+    program.parse_args(argc, argv);
+  } catch (const std::runtime_error& err) {
+    std::cerr << err.what() << std::endl;
+    std::cerr << program;
+    return 1;
+  }
+
+  // Setup up RMM memory pool
+  auto memory_resource = make_pool();
+  rmm::mr::set_current_device_resource(memory_resource.get());
+
+  // Initialize raft handle and running stream
+  const raft::handle_t handle_{};
+
+  // Run solver with appropriate precision
+  bool use_fp32 = program.get<bool>("--pdlp-fp32");
+  if (use_fp32) {
+    return run_solver<float>(program, handle_);
+  } else {
+    return run_solver<double>(program, handle_);
+  }
+}

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

@@ -224,7 +224,7 @@ class pdlp_solver_settings_t {
   bool detect_infeasibility{false};
   bool strict_infeasibility{false};
   i_t iteration_limit{std::numeric_limits<i_t>::max()};
-  double time_limit{std::numeric_limits<double>::infinity()};
+  f_t time_limit{std::numeric_limits<f_t>::infinity()};
   pdlp_solver_mode_t pdlp_solver_mode{pdlp_solver_mode_t::Stable3};
   bool log_to_console{true};
   std::string log_file{""};
@@ -239,6 +239,15 @@ class pdlp_solver_settings_t {
   i_t ordering{-1};
   i_t barrier_dual_initial_point{-1};
   bool eliminate_dense_columns{true};
+  /**
+   * @brief Enable mixed precision SpMV during PDHG iterations (FP64 mode only).
+   *
+   * When true, the constraint matrix A and its transpose are stored in FP32 while
+   * vectors and compute type remain in FP64, reducing memory bandwidth during SpMV.
+   * Convergence checking and restarts always use the full FP64 matrix, so this does
+   * not reduce overall memory usage.  Has no effect in FP32 mode.
+   */
+  bool mixed_precision_spmv{false};
   bool save_best_primal_so_far{false};
   bool first_primal_feasible{false};
   presolver_t presolver{presolver_t::Default};

cpp/src/dual_simplex/sparse_matrix.cpp

@@ -10,6 +10,7 @@
 #include <dual_simplex/sparse_vector.hpp>
 
 #include <dual_simplex/types.hpp>
+#include <mip_heuristics/mip_constants.hpp>
 
 // #include <thrust/for_each.h>
 // #include <thrust/iterator/counting_iterator.h>
@@ -932,6 +933,12 @@ f_t sparse_dot(const std::vector<i_t>& xind,
   return dot;
 }
 
+#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT
+// Minimal float instantiation for LP usage
+template class csc_matrix_t<int, float>;
+template class csr_matrix_t<int, float>;
+#endif
+
 #ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE
 template class csc_matrix_t<int, double>;
 

cpp/src/math_optimization/solution_writer.cu

@@ -9,15 +9,18 @@
 #include <utilities/logger.hpp>
 #include "solution_writer.hpp"
 
+#include <mip_heuristics/mip_constants.hpp>
+
 #include <fstream>
 
 namespace cuopt::linear_programming {
 
+template <typename f_t>
 void solution_writer_t::write_solution_to_sol_file(const std::string& filename,
                                                    const std::string& status,
-                                                   const double objective_value,
+                                                   const f_t objective_value,
                                                    const std::vector<std::string>& variable_names,
-                                                   const std::vector<double>& variable_values)
+                                                   const std::vector<f_t>& variable_values)
 {
   raft::common::nvtx::range fun_scope("write final solution to .sol file");
   std::ofstream file(filename.data());
@@ -27,7 +30,7 @@ void solution_writer_t::write_solution_to_sol_file(const std::string& filename,
     return;
   }
 
-  file.precision(std::numeric_limits<double>::max_digits10 + 1);
+  file.precision(std::numeric_limits<f_t>::max_digits10 + 1);
 
   file << "# Status: " << status << std::endl;
 
@@ -39,4 +42,22 @@ void solution_writer_t::write_solution_to_sol_file(const std::string& filename,
   }
 }
 
+#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT
+template void solution_writer_t::write_solution_to_sol_file<float>(
+  const std::string& filename,
+  const std::string& status,
+  const float objective_value,
+  const std::vector<std::string>& variable_names,
+  const std::vector<float>& variable_values);
+#endif
+
+#if MIP_INSTANTIATE_DOUBLE
+template void solution_writer_t::write_solution_to_sol_file<double>(
+  const std::string& filename,
+  const std::string& status,
+  const double objective_value,
+  const std::vector<std::string>& variable_names,
+  const std::vector<double>& variable_values);
+#endif
+
 }  // namespace cuopt::linear_programming

cpp/src/math_optimization/solution_writer.hpp

@@ -1,6 +1,6 @@
 /* clang-format off */
 /*
- * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+ * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
  * SPDX-License-Identifier: Apache-2.0
  */
 /* clang-format on */
@@ -23,10 +23,11 @@ namespace cuopt::linear_programming {
  */
 class solution_writer_t {
  public:
+  template <typename f_t>
   static void write_solution_to_sol_file(const std::string& sol_file_path,
                                          const std::string& status,
-                                         const double objective_value,
+                                         const f_t objective_value,
                                          const std::vector<std::string>& variable_names,
-                                         const std::vector<double>& variable_values);
+                                         const std::vector<f_t>& variable_values);
 };
 }  // namespace cuopt::linear_programming

cpp/src/math_optimization/solver_settings.cu

@@ -58,24 +58,24 @@ solver_settings_t<i_t, f_t>::solver_settings_t() : pdlp_settings(), mip_settings
   // clang-format off
   // Float parameters
   float_parameters = {
-    {CUOPT_TIME_LIMIT, &mip_settings.time_limit, 0.0, std::numeric_limits<f_t>::infinity(), std::numeric_limits<f_t>::infinity()},
-    {CUOPT_TIME_LIMIT, &pdlp_settings.time_limit, 0.0, std::numeric_limits<f_t>::infinity(), std::numeric_limits<f_t>::infinity()},
-    {CUOPT_WORK_LIMIT, &mip_settings.work_limit, 0.0, std::numeric_limits<f_t>::infinity(), std::numeric_limits<f_t>::infinity()},
-    {CUOPT_ABSOLUTE_DUAL_TOLERANCE, &pdlp_settings.tolerances.absolute_dual_tolerance, 0.0, 1e-1, 1e-4},
-    {CUOPT_RELATIVE_DUAL_TOLERANCE, &pdlp_settings.tolerances.relative_dual_tolerance, 0.0, 1e-1, 1e-4},
-    {CUOPT_ABSOLUTE_PRIMAL_TOLERANCE, &pdlp_settings.tolerances.absolute_primal_tolerance, 0.0, 1e-1, 1e-4},
-    {CUOPT_RELATIVE_PRIMAL_TOLERANCE, &pdlp_settings.tolerances.relative_primal_tolerance, 0.0, 1e-1, 1e-4},
-    {CUOPT_ABSOLUTE_GAP_TOLERANCE, &pdlp_settings.tolerances.absolute_gap_tolerance, 0.0, 1e-1, 1e-4},
-    {CUOPT_RELATIVE_GAP_TOLERANCE, &pdlp_settings.tolerances.relative_gap_tolerance, 0.0, 1e-1, 1e-4},
-    {CUOPT_MIP_ABSOLUTE_TOLERANCE, &mip_settings.tolerances.absolute_tolerance, 0.0, 1e-1, 1e-4},
-    {CUOPT_MIP_RELATIVE_TOLERANCE, &mip_settings.tolerances.relative_tolerance, 0.0, 1e-1, 1e-4},
-    {CUOPT_MIP_INTEGRALITY_TOLERANCE, &mip_settings.tolerances.integrality_tolerance, 0.0, 1e-1, 1e-5},
-    {CUOPT_MIP_ABSOLUTE_GAP, &mip_settings.tolerances.absolute_mip_gap, 0.0, CUOPT_INFINITY, 1e-10},
-    {CUOPT_MIP_RELATIVE_GAP, &mip_settings.tolerances.relative_mip_gap, 0.0, 1e-1, 1e-4},
-    {CUOPT_PRIMAL_INFEASIBLE_TOLERANCE, &pdlp_settings.tolerances.primal_infeasible_tolerance, 0.0, 1e-1, 1e-10},
-    {CUOPT_DUAL_INFEASIBLE_TOLERANCE, &pdlp_settings.tolerances.dual_infeasible_tolerance, 0.0, 1e-1, 1e-10},
-    {CUOPT_MIP_CUT_CHANGE_THRESHOLD, &mip_settings.cut_change_threshold, 0.0, std::numeric_limits<f_t>::infinity(), 1e-3},
-    {CUOPT_MIP_CUT_MIN_ORTHOGONALITY, &mip_settings.cut_min_orthogonality, 0.0, 1.0, 0.5}
+    {CUOPT_TIME_LIMIT, &mip_settings.time_limit, f_t(0.0), std::numeric_limits<f_t>::infinity(), std::numeric_limits<f_t>::infinity()},
+    {CUOPT_TIME_LIMIT, &pdlp_settings.time_limit, f_t(0.0), std::numeric_limits<f_t>::infinity(), std::numeric_limits<f_t>::infinity()},
+    {CUOPT_WORK_LIMIT, &mip_settings.work_limit, f_t(0.0), std::numeric_limits<f_t>::infinity(), std::numeric_limits<f_t>::infinity()},
+    {CUOPT_ABSOLUTE_DUAL_TOLERANCE, &pdlp_settings.tolerances.absolute_dual_tolerance, f_t(0.0), f_t(1e-1), f_t(1e-4)},
+    {CUOPT_RELATIVE_DUAL_TOLERANCE, &pdlp_settings.tolerances.relative_dual_tolerance, f_t(0.0), f_t(1e-1), f_t(1e-4)},
+    {CUOPT_ABSOLUTE_PRIMAL_TOLERANCE, &pdlp_settings.tolerances.absolute_primal_tolerance, f_t(0.0), f_t(1e-1), f_t(1e-4)},
+    {CUOPT_RELATIVE_PRIMAL_TOLERANCE, &pdlp_settings.tolerances.relative_primal_tolerance, f_t(0.0), f_t(1e-1), f_t(1e-4)},
+    {CUOPT_ABSOLUTE_GAP_TOLERANCE, &pdlp_settings.tolerances.absolute_gap_tolerance, f_t(0.0), f_t(1e-1), f_t(1e-4)},
+    {CUOPT_RELATIVE_GAP_TOLERANCE, &pdlp_settings.tolerances.relative_gap_tolerance, f_t(0.0), f_t(1e-1), f_t(1e-4)},
+    {CUOPT_MIP_ABSOLUTE_TOLERANCE, &mip_settings.tolerances.absolute_tolerance, f_t(0.0), f_t(1e-1), f_t(1e-4)},
+    {CUOPT_MIP_RELATIVE_TOLERANCE, &mip_settings.tolerances.relative_tolerance, f_t(0.0), f_t(1e-1), f_t(1e-4)},
+    {CUOPT_MIP_INTEGRALITY_TOLERANCE, &mip_settings.tolerances.integrality_tolerance, f_t(0.0), f_t(1e-1), f_t(1e-5)},
+    {CUOPT_MIP_ABSOLUTE_GAP, &mip_settings.tolerances.absolute_mip_gap, f_t(0.0), std::numeric_limits<f_t>::infinity(), std::max(f_t(1e-10), std::numeric_limits<f_t>::epsilon())},
+    {CUOPT_MIP_RELATIVE_GAP, &mip_settings.tolerances.relative_mip_gap, f_t(0.0), f_t(1e-1), f_t(1e-4)},
+    {CUOPT_PRIMAL_INFEASIBLE_TOLERANCE, &pdlp_settings.tolerances.primal_infeasible_tolerance, f_t(0.0), f_t(1e-1), std::max(f_t(1e-10), std::numeric_limits<f_t>::epsilon())},
+    {CUOPT_DUAL_INFEASIBLE_TOLERANCE, &pdlp_settings.tolerances.dual_infeasible_tolerance, f_t(0.0), f_t(1e-1), std::max(f_t(1e-10), std::numeric_limits<f_t>::epsilon())},
+    {CUOPT_MIP_CUT_CHANGE_THRESHOLD, &mip_settings.cut_change_threshold, f_t(0.0), std::numeric_limits<f_t>::infinity(), f_t(1e-3)},
+    {CUOPT_MIP_CUT_MIN_ORTHOGONALITY, &mip_settings.cut_min_orthogonality, f_t(0.0), f_t(1.0), f_t(0.5)}
    };
 
   // Int parameters

cpp/src/mip_heuristics/diversity/lns/rins.cu

@@ -186,7 +186,7 @@ void rins_t<i_t, f_t>::run_rins()
 
   total_calls++;
   node_count_at_last_rins = node_count.load();
-  time_limit              = std::min(time_limit, dm.timer.remaining_time());
+  time_limit              = std::min(time_limit, static_cast<f_t>(dm.timer.remaining_time()));
   CUOPT_LOG_DEBUG("Running RINS on solution with objective %g, fixing %d/%d",
                   best_sol.get_user_objective(),
                   vars_to_fix.size(),
@@ -288,22 +288,22 @@ void rins_t<i_t, f_t>::run_rins()
   if (branch_and_bound_status == dual_simplex::mip_status_t::OPTIMAL) {
     CUOPT_LOG_DEBUG("RINS submip optimal");
     // do goldilocks update
-    fixrate    = std::max(fixrate - 0.05, settings.min_fixrate);
-    time_limit = std::max(time_limit - 2, settings.min_time_limit);
+    fixrate    = std::max(fixrate - f_t(0.05), static_cast<f_t>(settings.min_fixrate));
+    time_limit = std::max(time_limit - f_t(2), static_cast<f_t>(settings.min_time_limit));
   } else if (branch_and_bound_status == dual_simplex::mip_status_t::TIME_LIMIT) {
     CUOPT_LOG_DEBUG("RINS submip time limit");
     // do goldilocks update
-    fixrate    = std::min(fixrate + 0.05, settings.max_fixrate);
-    time_limit = std::min(time_limit + 2, settings.max_time_limit);
+    fixrate    = std::min(fixrate + f_t(0.05), static_cast<f_t>(settings.max_fixrate));
+    time_limit = std::min(time_limit + f_t(2), static_cast<f_t>(settings.max_time_limit));
   } else if (branch_and_bound_status == dual_simplex::mip_status_t::INFEASIBLE) {
     CUOPT_LOG_DEBUG("RINS submip infeasible");
     // do goldilocks update, decreasing fixrate
-    fixrate = std::max(fixrate - 0.05, settings.min_fixrate);
+    fixrate = std::max(fixrate - f_t(0.05), static_cast<f_t>(settings.min_fixrate));
   } else {
     CUOPT_LOG_DEBUG("RINS solution not found");
     // do goldilocks update
-    fixrate    = std::min(fixrate + 0.05, settings.max_fixrate);
-    time_limit = std::min(time_limit + 2, settings.max_time_limit);
+    fixrate    = std::min(fixrate + f_t(0.05), static_cast<f_t>(settings.max_fixrate));
+    time_limit = std::min(time_limit + f_t(2), static_cast<f_t>(settings.max_time_limit));
   }
 
   cpu_fj_thread.stop_cpu_solver();

cpp/src/mip_heuristics/local_search/rounding/simple_rounding.cu

@@ -179,7 +179,7 @@ void invoke_correct_integers(solution_t<i_t, f_t>& solution, f_t tol)
   template void invoke_correct_integers<int, F_TYPE>(solution_t<int, F_TYPE> & solution,     \
                                                      F_TYPE tol);
 
-#if MIP_INSTANTIATE_FLOAT
+#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT
 INSTANTIATE(float)
 #endif
 

cpp/src/mip_heuristics/mip_constants.hpp

@@ -1,6 +1,6 @@
 /* clang-format off */
 /*
- * SPDX-FileCopyrightText: Copyright (c) 2022-2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+ * SPDX-FileCopyrightText: Copyright (c) 2022-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
  * SPDX-License-Identifier: Apache-2.0
  */
 /* clang-format on */
@@ -11,3 +11,5 @@
 
 #define MIP_INSTANTIATE_FLOAT  CUOPT_INSTANTIATE_FLOAT
 #define MIP_INSTANTIATE_DOUBLE CUOPT_INSTANTIATE_DOUBLE
+
+#define PDLP_INSTANTIATE_FLOAT 1

cpp/src/mip_heuristics/presolve/gf2_presolve.cpp

@@ -247,7 +247,7 @@ papilo::PresolveStatus GF2Presolve<f_t>::execute(const papilo::Problem<f_t>& pro
 
 #define INSTANTIATE(F_TYPE) template class GF2Presolve<F_TYPE>;
 
-#if MIP_INSTANTIATE_FLOAT
+#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT
 INSTANTIATE(float)
 #endif