diff --git a/benchmarks/linear_programming/cuopt/run_pdlp.cu b/benchmarks/linear_programming/cuopt/run_pdlp.cu index 18a473d64e..a7838d773e 100644 --- a/benchmarks/linear_programming/cuopt/run_pdlp.cu +++ b/benchmarks/linear_programming/cuopt/run_pdlp.cu @@ -76,6 +76,13 @@ 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-precision") + .help( + "PDLP precision mode. default: native type, single: FP32 internally, " + "double: FP64 explicitly, mixed: mixed-precision SpMV (FP32 matrix, FP64 vectors).") + .default_value(std::string("default")) + .choices("default", "single", "double", "mixed"); } static cuopt::linear_programming::presolver_t string_to_presolver(const std::string& presolver) @@ -87,6 +94,15 @@ static cuopt::linear_programming::presolver_t string_to_presolver(const std::str return cuopt::linear_programming::presolver_t::Default; } +static cuopt::linear_programming::pdlp_precision_t string_to_pdlp_precision( + const std::string& precision) +{ + if (precision == "single") return cuopt::linear_programming::pdlp_precision_t::SinglePrecision; + if (precision == "double") return cuopt::linear_programming::pdlp_precision_t::DoublePrecision; + if (precision == "mixed") return cuopt::linear_programming::pdlp_precision_t::MixedPrecision; + return cuopt::linear_programming::pdlp_precision_t::DefaultPrecision; +} + static cuopt::linear_programming::pdlp_solver_mode_t string_to_pdlp_solver_mode( const std::string& mode) { @@ -105,8 +121,7 @@ static cuopt::linear_programming::pdlp_solver_mode_t string_to_pdlp_solver_mode( static cuopt::linear_programming::pdlp_solver_settings_t create_solver_settings( const argparse::ArgumentParser& program) { - cuopt::linear_programming::pdlp_solver_settings_t settings = - cuopt::linear_programming::pdlp_solver_settings_t{}; + cuopt::linear_programming::pdlp_solver_settings_t settings{}; settings.time_limit = program.get("--time-limit"); settings.iteration_limit = program.get("--iteration-limit"); @@ -114,29 +129,16 @@ static cuopt::linear_programming::pdlp_solver_settings_t create_sol settings.pdlp_solver_mode = string_to_pdlp_solver_mode(program.get("--pdlp-solver-mode")); settings.method = static_cast(program.get("--method")); - settings.crossover = program.get("--crossover"); - settings.presolver = string_to_presolver(program.get("--presolver")); + settings.crossover = program.get("--crossover"); + settings.presolver = string_to_presolver(program.get("--presolver")); + settings.pdlp_precision = string_to_pdlp_precision(program.get("--pdlp-precision")); return settings; } -int main(int argc, char* argv[]) +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 settings = - create_solver_settings(program); + auto settings = create_solver_settings(program); bool use_pdlp_solver_mode = true; if (program.is_used("--pdlp-hyper-params-path")) { @@ -145,13 +147,6 @@ 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 op_problem = cuopt::mps_parser::parse_mps(program.get("--path")); @@ -168,3 +163,27 @@ 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_{}; + + return run_solver(program, handle_); +} diff --git a/cpp/cuopt_cli.cpp b/cpp/cuopt_cli.cpp index e9b1ee3719..899a3118b3 100644 --- a/cpp/cuopt_cli.cpp +++ b/cpp/cuopt_cli.cpp @@ -77,8 +77,8 @@ inline auto make_async() { return std::make_shared& settings) { - return cuopt::init_logger_t(settings.get_parameter(CUOPT_LOG_FILE), - settings.get_parameter(CUOPT_LOG_TO_CONSOLE)); + return cuopt::init_logger_t(settings.template get_parameter(CUOPT_LOG_FILE), + settings.template get_parameter(CUOPT_LOG_TO_CONSOLE)); } /** @@ -287,6 +287,17 @@ int main(int argc, char* argv[]) .implicit_value(true); std::map arg_name_to_param_name; + + // Register --pdlp-precision with string-to-int mapping so that it flows + // through the settings_strings map like other settings. + program.add_argument("--pdlp-precision") + .help( + "PDLP precision mode. default: native type, single: FP32 internally, " + "double: FP64 explicitly, mixed: mixed-precision SpMV (FP32 matrix, FP64 vectors).") + .default_value(std::string("-1")) + .choices("default", "single", "double", "mixed", "-1", "0", "1", "2"); + arg_name_to_param_name["--pdlp-precision"] = CUOPT_PDLP_PRECISION; + { // Add all solver settings as arguments cuopt::linear_programming::solver_settings_t dummy_settings; @@ -341,11 +352,20 @@ int main(int argc, char* argv[]) return 1; } + // Map symbolic pdlp-precision names to integer values + static const std::map precision_name_to_value = { + {"default", "-1"}, {"single", "0"}, {"double", "1"}, {"mixed", "2"}}; + // Read everything as a string std::map settings_strings; for (auto& [arg_name, param_name] : arg_name_to_param_name) { if (program.is_used(arg_name.c_str())) { - settings_strings[param_name] = program.get(arg_name.c_str()); + auto val = program.get(arg_name.c_str()); + if (param_name == CUOPT_PDLP_PRECISION) { + auto it = precision_name_to_value.find(val); + if (it != precision_name_to_value.end()) { val = it->second; } + } + settings_strings[param_name] = val; } } // Get the values diff --git a/cpp/include/cuopt/linear_programming/constants.h b/cpp/include/cuopt/linear_programming/constants.h index 7eb0aa07d6..d9dfbce16d 100644 --- a/cpp/include/cuopt/linear_programming/constants.h +++ b/cpp/include/cuopt/linear_programming/constants.h @@ -74,6 +74,7 @@ #define CUOPT_NUM_GPUS "num_gpus" #define CUOPT_USER_PROBLEM_FILE "user_problem_file" #define CUOPT_RANDOM_SEED "random_seed" +#define CUOPT_PDLP_PRECISION "pdlp_precision" /* @brief MIP determinism mode constants */ #define CUOPT_MODE_OPPORTUNISTIC 0 @@ -125,6 +126,12 @@ #define CUOPT_METHOD_DUAL_SIMPLEX 2 #define CUOPT_METHOD_BARRIER 3 +/* @brief PDLP precision mode constants */ +#define CUOPT_PDLP_DEFAULT_PRECISION -1 +#define CUOPT_PDLP_SINGLE_PRECISION 0 +#define CUOPT_PDLP_DOUBLE_PRECISION 1 +#define CUOPT_PDLP_MIXED_PRECISION 2 + /* @brief File format constants for problem I/O */ #define CUOPT_FILE_FORMAT_MPS 0 diff --git a/cpp/include/cuopt/linear_programming/optimization_problem.hpp b/cpp/include/cuopt/linear_programming/optimization_problem.hpp index d0f624ebdf..df78dd17c7 100644 --- a/cpp/include/cuopt/linear_programming/optimization_problem.hpp +++ b/cpp/include/cuopt/linear_programming/optimization_problem.hpp @@ -312,6 +312,14 @@ class optimization_problem_t : public optimization_problem_interface_t // Conversion // ============================================================================ + /** + * @brief Convert this problem to a different floating-point precision. + * + * @tparam other_f_t Target floating-point type (e.g. float when this is double) + */ + template + optimization_problem_t convert_to_other_prec(rmm::cuda_stream_view stream) const; + /** * @brief Returns nullptr since this is already a GPU problem. * @return nullptr diff --git a/cpp/include/cuopt/linear_programming/pdlp/solver_settings.hpp b/cpp/include/cuopt/linear_programming/pdlp/solver_settings.hpp index f6ad4c8619..d3f59144cc 100644 --- a/cpp/include/cuopt/linear_programming/pdlp/solver_settings.hpp +++ b/cpp/include/cuopt/linear_programming/pdlp/solver_settings.hpp @@ -63,6 +63,21 @@ enum method_t : int { Barrier = CUOPT_METHOD_BARRIER }; +/** + * @brief Enum representing the PDLP precision modes. + * + * DefaultPrecision: Use the type of the problem (FP64 for double problems). + * SinglePrecision: Run PDLP internally in FP32, converting inputs and outputs. + * DoublePrecision: Explicitly run in FP64 (same as default for double problems). + * MixedPrecision: Use mixed precision SpMV (FP32 matrix with FP64 vectors/compute). + */ +enum pdlp_precision_t : int { + DefaultPrecision = CUOPT_PDLP_DEFAULT_PRECISION, + SinglePrecision = CUOPT_PDLP_SINGLE_PRECISION, + DoublePrecision = CUOPT_PDLP_DOUBLE_PRECISION, + MixedPrecision = CUOPT_PDLP_MIXED_PRECISION +}; + template class pdlp_solver_settings_t { public: @@ -224,7 +239,7 @@ class pdlp_solver_settings_t { bool detect_infeasibility{false}; bool strict_infeasibility{false}; i_t iteration_limit{std::numeric_limits::max()}; - double time_limit{std::numeric_limits::infinity()}; + f_t time_limit{std::numeric_limits::infinity()}; pdlp_solver_mode_t pdlp_solver_mode{pdlp_solver_mode_t::Stable3}; bool log_to_console{true}; std::string log_file{""}; @@ -239,6 +254,7 @@ class pdlp_solver_settings_t { i_t ordering{-1}; i_t barrier_dual_initial_point{-1}; bool eliminate_dense_columns{true}; + pdlp_precision_t pdlp_precision{pdlp_precision_t::DefaultPrecision}; bool save_best_primal_so_far{false}; bool first_primal_feasible{false}; presolver_t presolver{presolver_t::Default}; diff --git a/cpp/src/dual_simplex/sparse_matrix.cpp b/cpp/src/dual_simplex/sparse_matrix.cpp index 8ccccd57cf..63004be72b 100644 --- a/cpp/src/dual_simplex/sparse_matrix.cpp +++ b/cpp/src/dual_simplex/sparse_matrix.cpp @@ -10,6 +10,7 @@ #include #include +#include // #include // #include @@ -938,6 +939,12 @@ f_t sparse_dot(const std::vector& xind, return dot; } +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT +// Minimal float instantiation for LP usage +template class csc_matrix_t; +template class csr_matrix_t; +#endif + #ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE template class csc_matrix_t; diff --git a/cpp/src/math_optimization/solution_writer.cu b/cpp/src/math_optimization/solution_writer.cu index 273b8e989c..880127546d 100644 --- a/cpp/src/math_optimization/solution_writer.cu +++ b/cpp/src/math_optimization/solution_writer.cu @@ -9,15 +9,18 @@ #include #include "solution_writer.hpp" +#include + #include namespace cuopt::linear_programming { +template 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& variable_names, - const std::vector& variable_values) + const std::vector& 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::max_digits10 + 1); + file.precision(std::numeric_limits::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( + const std::string& filename, + const std::string& status, + const float objective_value, + const std::vector& variable_names, + const std::vector& variable_values); +#endif + +#if MIP_INSTANTIATE_DOUBLE +template void solution_writer_t::write_solution_to_sol_file( + const std::string& filename, + const std::string& status, + const double objective_value, + const std::vector& variable_names, + const std::vector& variable_values); +#endif + } // namespace cuopt::linear_programming diff --git a/cpp/src/math_optimization/solution_writer.hpp b/cpp/src/math_optimization/solution_writer.hpp index 0890bf260b..0ac1b64464 100644 --- a/cpp/src/math_optimization/solution_writer.hpp +++ b/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 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& variable_names, - const std::vector& variable_values); + const std::vector& variable_values); }; } // namespace cuopt::linear_programming diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index f1350ca432..7435bb37fa 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -58,24 +58,24 @@ solver_settings_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::infinity(), std::numeric_limits::infinity()}, - {CUOPT_TIME_LIMIT, &pdlp_settings.time_limit, 0.0, std::numeric_limits::infinity(), std::numeric_limits::infinity()}, - {CUOPT_WORK_LIMIT, &mip_settings.work_limit, 0.0, std::numeric_limits::infinity(), std::numeric_limits::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::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::infinity(), std::numeric_limits::infinity()}, + {CUOPT_TIME_LIMIT, &pdlp_settings.time_limit, f_t(0.0), std::numeric_limits::infinity(), std::numeric_limits::infinity()}, + {CUOPT_WORK_LIMIT, &mip_settings.work_limit, f_t(0.0), std::numeric_limits::infinity(), std::numeric_limits::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::infinity(), std::max(f_t(1e-10), std::numeric_limits::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::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::epsilon())}, + {CUOPT_MIP_CUT_CHANGE_THRESHOLD, &mip_settings.cut_change_threshold, f_t(0.0), std::numeric_limits::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 @@ -103,7 +103,8 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_PRESOLVE, reinterpret_cast(&mip_settings.presolver), CUOPT_PRESOLVE_DEFAULT, CUOPT_PRESOLVE_PSLP, CUOPT_PRESOLVE_DEFAULT}, {CUOPT_MIP_DETERMINISM_MODE, &mip_settings.determinism_mode, CUOPT_MODE_OPPORTUNISTIC, CUOPT_MODE_DETERMINISTIC, CUOPT_MODE_OPPORTUNISTIC}, {CUOPT_RANDOM_SEED, &mip_settings.seed, -1, std::numeric_limits::max(), -1}, - {CUOPT_MIP_RELIABILITY_BRANCHING, &mip_settings.reliability_branching, -1, std::numeric_limits::max(), -1} + {CUOPT_MIP_RELIABILITY_BRANCHING, &mip_settings.reliability_branching, -1, std::numeric_limits::max(), -1}, + {CUOPT_PDLP_PRECISION, reinterpret_cast(&pdlp_settings.pdlp_precision), CUOPT_PDLP_DEFAULT_PRECISION, CUOPT_PDLP_MIXED_PRECISION, CUOPT_PDLP_DEFAULT_PRECISION} }; // Bool parameters @@ -120,7 +121,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_CROSSOVER, &pdlp_settings.crossover, false}, {CUOPT_ELIMINATE_DENSE_COLUMNS, &pdlp_settings.eliminate_dense_columns, true}, {CUOPT_CUDSS_DETERMINISTIC, &pdlp_settings.cudss_deterministic, false}, - {CUOPT_DUAL_POSTSOLVE, &pdlp_settings.dual_postsolve, true} + {CUOPT_DUAL_POSTSOLVE, &pdlp_settings.dual_postsolve, true}, }; // String parameters string_parameters = { diff --git a/cpp/src/mip_heuristics/diversity/lns/rins.cu b/cpp/src/mip_heuristics/diversity/lns/rins.cu index 7fd8533f82..2fbe79ba34 100644 --- a/cpp/src/mip_heuristics/diversity/lns/rins.cu +++ b/cpp/src/mip_heuristics/diversity/lns/rins.cu @@ -186,7 +186,7 @@ void rins_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(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::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(settings.min_fixrate)); + time_limit = std::max(time_limit - f_t(2), static_cast(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(settings.max_fixrate)); + time_limit = std::min(time_limit + f_t(2), static_cast(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(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(settings.max_fixrate)); + time_limit = std::min(time_limit + f_t(2), static_cast(settings.max_time_limit)); } cpu_fj_thread.stop_cpu_solver(); diff --git a/cpp/src/mip_heuristics/local_search/rounding/simple_rounding.cu b/cpp/src/mip_heuristics/local_search/rounding/simple_rounding.cu index c9a6dd0eda..4f3a015a6c 100644 --- a/cpp/src/mip_heuristics/local_search/rounding/simple_rounding.cu +++ b/cpp/src/mip_heuristics/local_search/rounding/simple_rounding.cu @@ -179,7 +179,7 @@ void invoke_correct_integers(solution_t& solution, f_t tol) template void invoke_correct_integers(solution_t & solution, \ F_TYPE tol); -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT INSTANTIATE(float) #endif diff --git a/cpp/src/mip_heuristics/mip_constants.hpp b/cpp/src/mip_heuristics/mip_constants.hpp index 66f5ebd273..47d3d22de4 100644 --- a/cpp/src/mip_heuristics/mip_constants.hpp +++ b/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 diff --git a/cpp/src/mip_heuristics/presolve/gf2_presolve.cpp b/cpp/src/mip_heuristics/presolve/gf2_presolve.cpp index 45ea4e420f..8ab0176cc4 100644 --- a/cpp/src/mip_heuristics/presolve/gf2_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/gf2_presolve.cpp @@ -247,7 +247,7 @@ papilo::PresolveStatus GF2Presolve::execute(const papilo::Problem& pro #define INSTANTIATE(F_TYPE) template class GF2Presolve; -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT INSTANTIATE(float) #endif diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index 5a89393a6a..20a586f6fb 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -432,6 +432,7 @@ optimization_problem_t build_optimization_problem( const int* cols = constraint_matrix.getConstraintMatrix().getColumns(); const f_t* coeffs = constraint_matrix.getConstraintMatrix().getValues(); + op_problem.set_csr_constraint_matrix( &(coeffs[start]), nnz, &(cols[start]), nnz, offsets.data(), nrows + 1); @@ -535,7 +536,7 @@ void set_presolve_options(papilo::Presolve& presolver, problem_category_t category, f_t absolute_tolerance, f_t relative_tolerance, - double time_limit, + f_t time_limit, bool dual_postsolve, i_t num_cpu_threads) { @@ -572,24 +573,30 @@ template std::optional> third_party_presolve_t::apply_pslp( optimization_problem_t const& op_problem, const double time_limit) { - f_t original_obj_offset = op_problem.get_objective_offset(); - auto ctx = build_and_run_pslp_presolver(op_problem, maximize_, time_limit); + if constexpr (std::is_same_v) { + double original_obj_offset = op_problem.get_objective_offset(); + auto ctx = build_and_run_pslp_presolver(op_problem, maximize_, time_limit); - // Free previously allocated presolver and settings - if (pslp_presolver_ != nullptr) { free_presolver(pslp_presolver_); } - if (pslp_stgs_ != nullptr) { free_settings(pslp_stgs_); } + // Free previously allocated presolver and settings if they exist + if (pslp_presolver_ != nullptr) { free_presolver(pslp_presolver_); } + if (pslp_stgs_ != nullptr) { free_settings(pslp_stgs_); } - pslp_presolver_ = ctx.presolver; - pslp_stgs_ = ctx.settings; + pslp_presolver_ = ctx.presolver; + pslp_stgs_ = ctx.settings; - if (ctx.status == PresolveStatus_::INFEASIBLE || ctx.status == PresolveStatus_::UNBNDORINFEAS) { - return std::nullopt; - } + if (ctx.status == PresolveStatus_::INFEASIBLE || ctx.status == PresolveStatus_::UNBNDORINFEAS) { + return std::nullopt; + } - auto opt_problem = build_optimization_problem_from_pslp( - pslp_presolver_, op_problem.get_handle_ptr(), maximize_, original_obj_offset); + auto opt_problem = build_optimization_problem_from_pslp( + pslp_presolver_, op_problem.get_handle_ptr(), maximize_, original_obj_offset); - return std::make_optional(third_party_presolve_result_t{opt_problem, {}}); + return std::make_optional(third_party_presolve_result_t{opt_problem, {}}); + } else { + cuopt_expects( + false, error_type_t::ValidationError, "PSLP presolver only supports double precision"); + return std::nullopt; + } } template @@ -625,7 +632,7 @@ std::optional> third_party_presolve_t papilo_presolver; - set_presolve_methods(papilo_presolver, category, dual_postsolve); + set_presolve_methods(papilo_presolver, category, dual_postsolve); set_presolve_options(papilo_presolver, category, absolute_tolerance, @@ -633,7 +640,7 @@ std::optional> third_party_presolve_t( + set_presolve_parameters( papilo_presolver, category, op_problem.get_n_constraints(), op_problem.get_n_variables()); // Disable papilo logs @@ -697,12 +704,14 @@ void third_party_presolve_t::undo(rmm::device_uvector& primal_sol } if (status_to_skip) { return; } + std::vector primal_sol_vec_h(primal_solution.size()); raft::copy(primal_sol_vec_h.data(), primal_solution.data(), primal_solution.size(), stream_view); std::vector dual_sol_vec_h(dual_solution.size()); raft::copy(dual_sol_vec_h.data(), dual_solution.data(), dual_solution.size(), stream_view); std::vector reduced_costs_vec_h(reduced_costs.size()); raft::copy(reduced_costs_vec_h.data(), reduced_costs.data(), reduced_costs.size(), stream_view); + papilo::Solution reduced_sol(primal_sol_vec_h); if (dual_postsolve) { reduced_sol.dual = dual_sol_vec_h; @@ -734,26 +743,34 @@ void third_party_presolve_t::undo_pslp(rmm::device_uvector& prima rmm::device_uvector& reduced_costs, rmm::cuda_stream_view stream_view) { - std::vector h_primal_solution(primal_solution.size()); - std::vector h_dual_solution(dual_solution.size()); - std::vector h_reduced_costs(reduced_costs.size()); - raft::copy(h_primal_solution.data(), primal_solution.data(), primal_solution.size(), stream_view); - raft::copy(h_dual_solution.data(), dual_solution.data(), dual_solution.size(), stream_view); - raft::copy(h_reduced_costs.data(), reduced_costs.data(), reduced_costs.size(), stream_view); - - postsolve( - pslp_presolver_, h_primal_solution.data(), h_dual_solution.data(), h_reduced_costs.data()); - - auto uncrushed_sol = pslp_presolver_->sol; - int n_cols = uncrushed_sol->dim_x; - int n_rows = uncrushed_sol->dim_y; - - primal_solution.resize(n_cols, stream_view); - dual_solution.resize(n_rows, stream_view); - reduced_costs.resize(n_cols, stream_view); - raft::copy(primal_solution.data(), uncrushed_sol->x, n_cols, stream_view); - raft::copy(dual_solution.data(), uncrushed_sol->y, n_rows, stream_view); - raft::copy(reduced_costs.data(), uncrushed_sol->z, n_cols, stream_view); + if constexpr (std::is_same_v) { + // PSLP uses double internally, so we can use the data directly + std::vector h_primal_solution(primal_solution.size()); + std::vector h_dual_solution(dual_solution.size()); + std::vector h_reduced_costs(reduced_costs.size()); + raft::copy( + h_primal_solution.data(), primal_solution.data(), primal_solution.size(), stream_view); + raft::copy(h_dual_solution.data(), dual_solution.data(), dual_solution.size(), stream_view); + raft::copy(h_reduced_costs.data(), reduced_costs.data(), reduced_costs.size(), stream_view); + stream_view.synchronize(); + + postsolve( + pslp_presolver_, h_primal_solution.data(), h_dual_solution.data(), h_reduced_costs.data()); + + auto uncrushed_sol = pslp_presolver_->sol; + int n_cols = uncrushed_sol->dim_x; + int n_rows = uncrushed_sol->dim_y; + + primal_solution.resize(n_cols, stream_view); + dual_solution.resize(n_rows, stream_view); + reduced_costs.resize(n_cols, stream_view); + raft::copy(primal_solution.data(), uncrushed_sol->x, n_cols, stream_view); + raft::copy(dual_solution.data(), uncrushed_sol->y, n_rows, stream_view); + raft::copy(reduced_costs.data(), uncrushed_sol->z, n_cols, stream_view); + } else { + cuopt_expects( + false, error_type_t::ValidationError, "PSLP postsolve only supports double precision"); + } stream_view.synchronize(); } @@ -795,7 +812,7 @@ void papilo_postsolve_deleter::operator()(papilo::PostsolveStorage* pt delete ptr; } -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template struct papilo_postsolve_deleter; template class third_party_presolve_t; #endif diff --git a/cpp/src/mip_heuristics/problem/presolve_data.cu b/cpp/src/mip_heuristics/problem/presolve_data.cu index b11f7b108a..bf05efa875 100644 --- a/cpp/src/mip_heuristics/problem/presolve_data.cu +++ b/cpp/src/mip_heuristics/problem/presolve_data.cu @@ -245,7 +245,7 @@ void presolve_data_t::papilo_uncrush_assignment( problem.handle_ptr->sync_stream(); } -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template class presolve_data_t; #endif diff --git a/cpp/src/mip_heuristics/problem/problem.cu b/cpp/src/mip_heuristics/problem/problem.cu index bc93a9d988..fadc00a850 100644 --- a/cpp/src/mip_heuristics/problem/problem.cu +++ b/cpp/src/mip_heuristics/problem/problem.cu @@ -2380,7 +2380,7 @@ void problem_t::update_variable_bounds(const std::vector& var_ind RAFT_CHECK_CUDA(handle_ptr->get_stream()); } -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template class problem_t; #endif diff --git a/cpp/src/mip_heuristics/solution/solution.cu b/cpp/src/mip_heuristics/solution/solution.cu index 5f1c13199b..531d54372c 100644 --- a/cpp/src/mip_heuristics/solution/solution.cu +++ b/cpp/src/mip_heuristics/solution/solution.cu @@ -660,7 +660,7 @@ mip_solution_t solution_t::get_solution(bool output_feasible } } -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template class solution_t; #endif diff --git a/cpp/src/mip_heuristics/solver_solution.cu b/cpp/src/mip_heuristics/solver_solution.cu index 60556884c9..e497a21c8f 100644 --- a/cpp/src/mip_heuristics/solver_solution.cu +++ b/cpp/src/mip_heuristics/solver_solution.cu @@ -209,8 +209,8 @@ void mip_solution_t::write_to_sol_file(std::string_view filename, status = "Infeasible"; } - double objective_value = get_objective_value(); - auto& var_names = get_variable_names(); + f_t objective_value = get_objective_value(); + auto& var_names = get_variable_names(); std::vector solution; solution.resize(solution_.size()); raft::copy(solution.data(), solution_.data(), solution_.size(), stream_view.value()); @@ -234,7 +234,7 @@ void mip_solution_t::log_summary() const CUOPT_LOG_INFO("Total Solve Time: %f", get_total_solve_time()); } -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template class mip_solution_t; #endif diff --git a/cpp/src/pdlp/cpu_pdlp_warm_start_data.cu b/cpp/src/pdlp/cpu_pdlp_warm_start_data.cu index f9a73dff06..b078bc4779 100644 --- a/cpp/src/pdlp/cpu_pdlp_warm_start_data.cu +++ b/cpp/src/pdlp/cpu_pdlp_warm_start_data.cu @@ -108,14 +108,14 @@ pdlp_warm_start_data_t convert_to_gpu_warmstart( return gpu_data; } -// Explicit template instantiations +#if MIP_INSTANTIATE_DOUBLE template cpu_pdlp_warm_start_data_t convert_to_cpu_warmstart( const pdlp_warm_start_data_t&, rmm::cuda_stream_view); - template pdlp_warm_start_data_t convert_to_gpu_warmstart( const cpu_pdlp_warm_start_data_t&, rmm::cuda_stream_view); +#endif -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template cpu_pdlp_warm_start_data_t convert_to_cpu_warmstart( const pdlp_warm_start_data_t&, rmm::cuda_stream_view); diff --git a/cpp/src/pdlp/cusparse_view.cu b/cpp/src/pdlp/cusparse_view.cu index ca36dde421..64ec44f5ef 100644 --- a/cpp/src/pdlp/cusparse_view.cu +++ b/cpp/src/pdlp/cusparse_view.cu @@ -21,6 +21,12 @@ #include #include +#include + +struct double_to_float_functor { + __host__ __device__ float operator()(double val) const { return static_cast(val); } +}; + namespace cuopt::linear_programming::detail { // cusparse_sp_mat_descr_wrapper_t implementation @@ -277,7 +283,8 @@ cusparse_view_t::cusparse_view_t( rmm::device_uvector& _potential_next_dual_solution, rmm::device_uvector& _reflected_primal_solution, const std::vector& climber_strategies, - const pdlp_hyper_params::pdlp_hyper_params_t& hyper_params) + const pdlp_hyper_params::pdlp_hyper_params_t& hyper_params, + bool enable_mixed_precision_spmv) : batch_mode_(climber_strategies.size() > 1), handle_ptr_(handle_ptr), A{}, @@ -304,7 +311,12 @@ cusparse_view_t::cusparse_view_t( A_{op_problem_scaled.coefficients}, A_offsets_{op_problem_scaled.offsets}, A_indices_{op_problem_scaled.variables}, - climber_strategies_(climber_strategies) + climber_strategies_(climber_strategies), + A_float_{0, handle_ptr->get_stream()}, + A_T_float_{0, handle_ptr->get_stream()}, + buffer_non_transpose_mixed_{0, handle_ptr->get_stream()}, + buffer_transpose_mixed_{0, handle_ptr->get_stream()}, + mixed_precision_enabled_{false} { raft::common::nvtx::range fun_scope("Initializing cuSparse view"); @@ -583,6 +595,92 @@ cusparse_view_t::cusparse_view_t( handle_ptr->get_stream()); } #endif + + if constexpr (std::is_same_v) { + if (enable_mixed_precision_spmv && !batch_mode_) { + mixed_precision_enabled_ = true; + + A_float_.resize(op_problem_scaled.nnz, handle_ptr->get_stream()); + A_T_float_.resize(op_problem_scaled.nnz, handle_ptr->get_stream()); + + RAFT_CUDA_TRY(cub::DeviceTransform::Transform(op_problem_scaled.coefficients.data(), + A_float_.data(), + op_problem_scaled.nnz, + double_to_float_functor{}, + handle_ptr->get_stream().value())); + + RAFT_CUDA_TRY(cub::DeviceTransform::Transform(A_T_.data(), + A_T_float_.data(), + op_problem_scaled.nnz, + double_to_float_functor{}, + handle_ptr->get_stream().value())); + + A_mixed_.create(op_problem_scaled.n_constraints, + op_problem_scaled.n_variables, + op_problem_scaled.nnz, + const_cast(op_problem_scaled.offsets.data()), + const_cast(op_problem_scaled.variables.data()), + A_float_.data()); + + A_T_mixed_.create(op_problem_scaled.n_variables, + op_problem_scaled.n_constraints, + op_problem_scaled.nnz, + const_cast(A_T_offsets_.data()), + const_cast(A_T_indices_.data()), + A_T_float_.data()); + + const rmm::device_scalar alpha_d{1.0, handle_ptr->get_stream()}; + const rmm::device_scalar beta_d{0.0, handle_ptr->get_stream()}; + + size_t buffer_size_non_transpose_mixed = + mixed_precision_spmv_buffersize(handle_ptr_->get_cusparse_handle(), + CUSPARSE_OPERATION_NON_TRANSPOSE, + alpha_d.data(), + A_mixed_, + c, + beta_d.data(), + dual_solution, + CUSPARSE_SPMV_CSR_ALG2, + handle_ptr->get_stream()); + buffer_non_transpose_mixed_.resize(buffer_size_non_transpose_mixed, handle_ptr->get_stream()); + + size_t buffer_size_transpose_mixed = + mixed_precision_spmv_buffersize(handle_ptr_->get_cusparse_handle(), + CUSPARSE_OPERATION_NON_TRANSPOSE, + alpha_d.data(), + A_T_mixed_, + dual_solution, + beta_d.data(), + c, + CUSPARSE_SPMV_CSR_ALG2, + handle_ptr->get_stream()); + buffer_transpose_mixed_.resize(buffer_size_transpose_mixed, handle_ptr->get_stream()); + +#if CUDA_VER_12_4_UP + mixed_precision_spmv_preprocess(handle_ptr_->get_cusparse_handle(), + CUSPARSE_OPERATION_NON_TRANSPOSE, + alpha_d.data(), + A_mixed_, + c, + beta_d.data(), + dual_solution, + CUSPARSE_SPMV_CSR_ALG2, + buffer_non_transpose_mixed_.data(), + handle_ptr->get_stream()); + + mixed_precision_spmv_preprocess(handle_ptr_->get_cusparse_handle(), + CUSPARSE_OPERATION_NON_TRANSPOSE, + alpha_d.data(), + A_T_mixed_, + dual_solution, + beta_d.data(), + c, + CUSPARSE_SPMV_CSR_ALG2, + buffer_transpose_mixed_.data(), + handle_ptr->get_stream()); +#endif + } + } } // Used by pdlp object for current and average termination condition @@ -625,7 +723,12 @@ cusparse_view_t::cusparse_view_t( A_{op_problem.coefficients}, A_offsets_{op_problem.offsets}, A_indices_{op_problem.variables}, - climber_strategies_(climber_strategies) + climber_strategies_(climber_strategies), + A_float_{0, handle_ptr->get_stream()}, + A_T_float_{0, handle_ptr->get_stream()}, + buffer_non_transpose_mixed_{0, handle_ptr->get_stream()}, + buffer_transpose_mixed_{0, handle_ptr->get_stream()}, + mixed_precision_enabled_{false} { #ifdef PDLP_DEBUG_MODE RAFT_CUDA_TRY(cudaDeviceSynchronize()); @@ -832,7 +935,12 @@ cusparse_view_t::cusparse_view_t( A_{existing_cusparse_view.A_}, A_offsets_{existing_cusparse_view.A_offsets_}, A_indices_{existing_cusparse_view.A_indices_}, - climber_strategies_(existing_cusparse_view.climber_strategies_) + climber_strategies_(existing_cusparse_view.climber_strategies_), + A_float_{0, handle_ptr->get_stream()}, + A_T_float_{0, handle_ptr->get_stream()}, + buffer_non_transpose_mixed_{0, handle_ptr->get_stream()}, + buffer_transpose_mixed_{0, handle_ptr->get_stream()}, + mixed_precision_enabled_{false} { #ifdef PDLP_DEBUG_MODE RAFT_CUDA_TRY(cudaDeviceSynchronize()); @@ -942,11 +1050,105 @@ cusparse_view_t::cusparse_view_t( A_(dummy_float), A_offsets_(dummy_int), A_indices_(dummy_int), - climber_strategies_(climber_strategies) + climber_strategies_(climber_strategies), + A_float_{0, handle_ptr->get_stream()}, + A_T_float_{0, handle_ptr->get_stream()}, + buffer_non_transpose_mixed_{0, handle_ptr->get_stream()}, + buffer_transpose_mixed_{0, handle_ptr->get_stream()}, + mixed_precision_enabled_{false} +{ +} + +// Update FP32 matrix copies after scaling (must be called after scale_problem()) +template +void cusparse_view_t::update_mixed_precision_matrices() +{ + if constexpr (std::is_same_v) { + if (!mixed_precision_enabled_) { return; } + + RAFT_CUDA_TRY(cub::DeviceTransform::Transform(A_.data(), + A_float_.data(), + A_.size(), + double_to_float_functor{}, + handle_ptr_->get_stream().value())); + + RAFT_CUDA_TRY(cub::DeviceTransform::Transform(A_T_.data(), + A_T_float_.data(), + A_T_.size(), + double_to_float_functor{}, + handle_ptr_->get_stream().value())); + + handle_ptr_->get_stream().synchronize(); + } +} + +// Mixed precision SpMV implementation: FP32 matrix with FP64 vectors and FP64 compute type +size_t mixed_precision_spmv_buffersize(cusparseHandle_t handle, + cusparseOperation_t opA, + const double* alpha, + cusparseSpMatDescr_t matA, // FP32 matrix + cusparseDnVecDescr_t vecX, // FP64 vector + const double* beta, + cusparseDnVecDescr_t vecY, // FP64 vector + cusparseSpMVAlg_t alg, + cudaStream_t stream) +{ + size_t bufferSize = 0; + RAFT_CUSPARSE_TRY(cusparseSetStream(handle, stream)); + RAFT_CUSPARSE_TRY(cusparseSpMV_bufferSize( + handle, opA, alpha, matA, vecX, beta, vecY, CUDA_R_64F, alg, &bufferSize)); + return bufferSize; +} + +void mixed_precision_spmv(cusparseHandle_t handle, + cusparseOperation_t opA, + const double* alpha, + cusparseSpMatDescr_t matA, // FP32 matrix + cusparseDnVecDescr_t vecX, // FP64 vector + const double* beta, + cusparseDnVecDescr_t vecY, // FP64 vector + cusparseSpMVAlg_t alg, + void* externalBuffer, + cudaStream_t stream) +{ + RAFT_CUSPARSE_TRY(cusparseSetStream(handle, stream)); + RAFT_CUSPARSE_TRY( + cusparseSpMV(handle, opA, alpha, matA, vecX, beta, vecY, CUDA_R_64F, alg, externalBuffer)); +} + +#if CUDA_VER_12_4_UP +void mixed_precision_spmv_preprocess(cusparseHandle_t handle, + cusparseOperation_t opA, + const double* alpha, + cusparseSpMatDescr_t matA, // FP32 matrix + cusparseDnVecDescr_t vecX, // FP64 vector + const double* beta, + cusparseDnVecDescr_t vecY, // FP64 vector + cusparseSpMVAlg_t alg, + void* externalBuffer, + cudaStream_t stream) +{ + static const auto func = + dynamic_load_runtime::function("cusparseSpMV_preprocess"); + if (func.has_value()) { + RAFT_CUSPARSE_TRY(cusparseSetStream(handle, stream)); + RAFT_CUSPARSE_TRY( + (*func)(handle, opA, alpha, matA, vecX, beta, vecY, CUDA_R_64F, alg, externalBuffer)); + } +} +#endif + +bool is_cusparse_runtime_mixed_precision_supported() { + int major = 0, minor = 0; + auto status = cusparseGetProperty(libraryPropertyType_t::MAJOR_VERSION, &major); + if (status != CUSPARSE_STATUS_SUCCESS) return false; + status = cusparseGetProperty(libraryPropertyType_t::MINOR_VERSION, &minor); + if (status != CUSPARSE_STATUS_SUCCESS) return false; + return (major > 12) || (major == 12 && minor >= 5); } -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template class cusparse_sp_mat_descr_wrapper_t; template class cusparse_dn_vec_descr_wrapper_t; template class cusparse_dn_mat_descr_wrapper_t; @@ -960,7 +1162,7 @@ template class cusparse_view_t; #endif #if CUDA_VER_12_4_UP -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template void my_cusparsespmm_preprocess(cusparseHandle_t, cusparseOperation_t, cusparseOperation_t, diff --git a/cpp/src/pdlp/cusparse_view.hpp b/cpp/src/pdlp/cusparse_view.hpp index cbbc856924..416a0b1e5f 100644 --- a/cpp/src/pdlp/cusparse_view.hpp +++ b/cpp/src/pdlp/cusparse_view.hpp @@ -90,7 +90,8 @@ class cusparse_view_t { rmm::device_uvector& _potential_next_dual_solution, rmm::device_uvector& _reflected_primal_solution, const std::vector& climber_strategies, - const pdlp_hyper_params::pdlp_hyper_params_t& hyper_params); + const pdlp_hyper_params::pdlp_hyper_params_t& hyper_params, + bool enable_mixed_precision_spmv); cusparse_view_t(raft::handle_t const* handle_ptr, const problem_t& op_problem, @@ -194,8 +195,56 @@ class cusparse_view_t { const rmm::device_uvector& A_indices_; const std::vector& climber_strategies_; + + // Mixed precision SpMV support (FP32 matrix with FP64 vectors/compute) + // Only used when mixed_precision_enabled_ is true and f_t = double + rmm::device_uvector A_float_; // FP32 copy of A values + rmm::device_uvector A_T_float_; // FP32 copy of A_T values + cusparse_sp_mat_descr_wrapper_t A_mixed_; // FP32 matrix descriptor for A + cusparse_sp_mat_descr_wrapper_t A_T_mixed_; // FP32 matrix descriptor for A_T + rmm::device_uvector buffer_non_transpose_mixed_; // SpMV buffer for mixed precision A + rmm::device_uvector buffer_transpose_mixed_; // SpMV buffer for mixed precision A_T + bool mixed_precision_enabled_{false}; + + // Update FP32 matrix copies after scaling (must be called after scale_problem()) + void update_mixed_precision_matrices(); }; +// Mixed precision SpMV: FP32 matrix with FP64 vectors and FP64 compute type +void mixed_precision_spmv(cusparseHandle_t handle, + cusparseOperation_t opA, + const double* alpha, + cusparseSpMatDescr_t matA, // FP32 matrix + cusparseDnVecDescr_t vecX, // FP64 vector + const double* beta, + cusparseDnVecDescr_t vecY, // FP64 vector + cusparseSpMVAlg_t alg, + void* externalBuffer, + cudaStream_t stream); + +size_t mixed_precision_spmv_buffersize(cusparseHandle_t handle, + cusparseOperation_t opA, + const double* alpha, + cusparseSpMatDescr_t matA, // FP32 matrix + cusparseDnVecDescr_t vecX, // FP64 vector + const double* beta, + cusparseDnVecDescr_t vecY, // FP64 vector + cusparseSpMVAlg_t alg, + cudaStream_t stream); + +#if CUDA_VER_12_4_UP +void mixed_precision_spmv_preprocess(cusparseHandle_t handle, + cusparseOperation_t opA, + const double* alpha, + cusparseSpMatDescr_t matA, // FP32 matrix + cusparseDnVecDescr_t vecX, // FP64 vector + const double* beta, + cusparseDnVecDescr_t vecY, // FP64 vector + cusparseSpMVAlg_t alg, + void* externalBuffer, + cudaStream_t stream); +#endif + #if CUDA_VER_12_4_UP template < typename T, @@ -213,4 +262,6 @@ void my_cusparsespmm_preprocess(cusparseHandle_t handle, cudaStream_t stream); #endif +bool is_cusparse_runtime_mixed_precision_supported(); + } // namespace cuopt::linear_programming::detail diff --git a/cpp/src/pdlp/initial_scaling_strategy/initial_scaling.cu b/cpp/src/pdlp/initial_scaling_strategy/initial_scaling.cu index afa3ee5fb7..b618550f6e 100644 --- a/cpp/src/pdlp/initial_scaling_strategy/initial_scaling.cu +++ b/cpp/src/pdlp/initial_scaling_strategy/initial_scaling.cu @@ -858,7 +858,7 @@ pdlp_initial_scaling_strategy_t::view() int* A_T_offsets, \ int* A_T_indices); -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT INSTANTIATE(float) #endif diff --git a/cpp/src/pdlp/optimal_batch_size_handler/optimal_batch_size_handler.cu b/cpp/src/pdlp/optimal_batch_size_handler/optimal_batch_size_handler.cu index deb6b759aa..cbfb03618d 100644 --- a/cpp/src/pdlp/optimal_batch_size_handler/optimal_batch_size_handler.cu +++ b/cpp/src/pdlp/optimal_batch_size_handler/optimal_batch_size_handler.cu @@ -434,7 +434,7 @@ int optimal_batch_size_handler(const optimization_problem_t& op_proble return 0; } -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template int optimal_batch_size_handler( const optimization_problem_t& op_problem, int max_batch_size); #endif diff --git a/cpp/src/pdlp/optimization_problem.cu b/cpp/src/pdlp/optimization_problem.cu index d0888dd3ac..9b3016a113 100644 --- a/cpp/src/pdlp/optimization_problem.cu +++ b/cpp/src/pdlp/optimization_problem.cu @@ -40,6 +40,7 @@ #include #include +#include #include #include @@ -1505,15 +1506,105 @@ void optimization_problem_t::copy_variable_types_to_host(var_t* output cudaMemcpy(output, variable_types_.data(), size * sizeof(var_t), cudaMemcpyDeviceToHost)); } +template +struct cast_op { + HDI To operator()(From val) const { return static_cast(val); } +}; + +template +rmm::device_uvector gpu_cast(const rmm::device_uvector& src, rmm::cuda_stream_view stream) +{ + rmm::device_uvector dst(src.size(), stream); + if (src.size() > 0) { + RAFT_CUDA_TRY(cub::DeviceTransform::Transform( + src.data(), dst.data(), src.size(), cast_op{}, stream.value())); + } + return dst; +} + +template rmm::device_uvector gpu_cast(const rmm::device_uvector&, + rmm::cuda_stream_view); +template rmm::device_uvector gpu_cast(const rmm::device_uvector&, + rmm::cuda_stream_view); + +template +template +optimization_problem_t optimization_problem_t::convert_to_other_prec( + rmm::cuda_stream_view stream) const +{ + optimization_problem_t other(handle_ptr_); + + other.set_maximize(maximize_); + other.set_objective_offset(static_cast(objective_offset_)); + other.set_objective_scaling_factor(static_cast(objective_scaling_factor_)); + + if (A_.size() > 0) { + auto other_A = gpu_cast(A_, stream); + other.set_csr_constraint_matrix(other_A.data(), + static_cast(other_A.size()), + A_indices_.data(), + static_cast(A_indices_.size()), + A_offsets_.data(), + static_cast(A_offsets_.size())); + } + + if (c_.size() > 0) { + auto other_c = gpu_cast(c_, stream); + other.set_objective_coefficients(other_c.data(), static_cast(other_c.size())); + } + + if (b_.size() > 0) { + auto other_b = gpu_cast(b_, stream); + other.set_constraint_bounds(other_b.data(), static_cast(other_b.size())); + } + + if (constraint_lower_bounds_.size() > 0) { + auto other_clb = gpu_cast(constraint_lower_bounds_, stream); + other.set_constraint_lower_bounds(other_clb.data(), static_cast(other_clb.size())); + } + + if (constraint_upper_bounds_.size() > 0) { + auto other_cub = gpu_cast(constraint_upper_bounds_, stream); + other.set_constraint_upper_bounds(other_cub.data(), static_cast(other_cub.size())); + } + + if (variable_lower_bounds_.size() > 0) { + auto other_vlb = gpu_cast(variable_lower_bounds_, stream); + other.set_variable_lower_bounds(other_vlb.data(), static_cast(other_vlb.size())); + } + + if (variable_upper_bounds_.size() > 0) { + auto other_vub = gpu_cast(variable_upper_bounds_, stream); + other.set_variable_upper_bounds(other_vub.data(), static_cast(other_vub.size())); + } + + if (variable_types_.size() > 0) { + other.set_variable_types(variable_types_.data(), static_cast(variable_types_.size())); + } + + other.set_variable_names(var_names_); + other.set_row_names(row_names_); + other.set_objective_name(objective_name_); + other.set_problem_category(problem_category_); + + return other; +} + // ============================================================================== // Template instantiations // ============================================================================== // Explicit template instantiations matching MIP constants -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template class optimization_problem_t; #endif #if MIP_INSTANTIATE_DOUBLE template class optimization_problem_t; #endif +#if PDLP_INSTANTIATE_FLOAT || MIP_INSTANTIATE_FLOAT +template optimization_problem_t + optimization_problem_t::convert_to_other_prec( + rmm::cuda_stream_view) const; +#endif + } // namespace cuopt::linear_programming diff --git a/cpp/src/pdlp/pdhg.cu b/cpp/src/pdlp/pdhg.cu index 51e0b29381..74df7fee01 100644 --- a/cpp/src/pdlp/pdhg.cu +++ b/cpp/src/pdlp/pdhg.cu @@ -41,7 +41,8 @@ pdhg_solver_t::pdhg_solver_t( bool is_legacy_batch_mode, // Batch mode with streams const std::vector& climber_strategies, const pdlp_hyper_params::pdlp_hyper_params_t& hyper_params, - const std::vector>& new_bounds) + const std::vector>& new_bounds, + bool enable_mixed_precision_spmv) : batch_mode_(climber_strategies.size() > 1), handle_ptr_(handle_ptr), stream_view_(handle_ptr_->get_stream()), @@ -77,7 +78,8 @@ pdhg_solver_t::pdhg_solver_t( potential_next_dual_solution_, reflected_primal_, climber_strategies, - hyper_params}, + hyper_params, + enable_mixed_precision_spmv}, reusable_device_scalar_value_1_{1.0, stream_view_}, reusable_device_scalar_value_0_{0.0, stream_view_}, reusable_device_scalar_value_neg_1_{f_t(-1.0), stream_view_}, @@ -249,17 +251,33 @@ void pdhg_solver_t::compute_next_dual_solution(rmm::device_uvectorget_cusparse_handle(), - CUSPARSE_OPERATION_NON_TRANSPOSE, - reusable_device_scalar_value_1_.data(), // 1 - cusparse_view_.A, - cusparse_view_.tmp_primal, - reusable_device_scalar_value_0_.data(), // 1 - cusparse_view_.dual_gradient, - CUSPARSE_SPMV_CSR_ALG2, - (f_t*)cusparse_view_.buffer_non_transpose.data(), - stream_view_)); + if constexpr (std::is_same_v) { + if (cusparse_view_.mixed_precision_enabled_) { + mixed_precision_spmv(handle_ptr_->get_cusparse_handle(), + CUSPARSE_OPERATION_NON_TRANSPOSE, + reusable_device_scalar_value_1_.data(), + cusparse_view_.A_mixed_, + cusparse_view_.tmp_primal, + reusable_device_scalar_value_0_.data(), + cusparse_view_.dual_gradient, + CUSPARSE_SPMV_CSR_ALG2, + cusparse_view_.buffer_non_transpose_mixed_.data(), + stream_view_); + } + } + if (!cusparse_view_.mixed_precision_enabled_) { + RAFT_CUSPARSE_TRY( + raft::sparse::detail::cusparsespmv(handle_ptr_->get_cusparse_handle(), + CUSPARSE_OPERATION_NON_TRANSPOSE, + reusable_device_scalar_value_1_.data(), + cusparse_view_.A, + cusparse_view_.tmp_primal, + reusable_device_scalar_value_0_.data(), + cusparse_view_.dual_gradient, + CUSPARSE_SPMV_CSR_ALG2, + (f_t*)cusparse_view_.buffer_non_transpose.data(), + stream_view_)); + } // y - (sigma*dual_gradient) // max(min(0, sigma*constraint_upper+primal_product), sigma*constraint_lower+primal_product) @@ -287,17 +305,33 @@ void pdhg_solver_t::compute_At_y() // A_t @ y if (!batch_mode_) { - RAFT_CUSPARSE_TRY( - raft::sparse::detail::cusparsespmv(handle_ptr_->get_cusparse_handle(), - CUSPARSE_OPERATION_NON_TRANSPOSE, - reusable_device_scalar_value_1_.data(), - cusparse_view_.A_T, - cusparse_view_.dual_solution, - reusable_device_scalar_value_0_.data(), - cusparse_view_.current_AtY, - CUSPARSE_SPMV_CSR_ALG2, - (f_t*)cusparse_view_.buffer_transpose.data(), - stream_view_)); + if constexpr (std::is_same_v) { + if (cusparse_view_.mixed_precision_enabled_) { + mixed_precision_spmv(handle_ptr_->get_cusparse_handle(), + CUSPARSE_OPERATION_NON_TRANSPOSE, + reusable_device_scalar_value_1_.data(), + cusparse_view_.A_T_mixed_, + cusparse_view_.dual_solution, + reusable_device_scalar_value_0_.data(), + cusparse_view_.current_AtY, + CUSPARSE_SPMV_CSR_ALG2, + cusparse_view_.buffer_transpose_mixed_.data(), + stream_view_); + } + } + if (!cusparse_view_.mixed_precision_enabled_) { + RAFT_CUSPARSE_TRY( + raft::sparse::detail::cusparsespmv(handle_ptr_->get_cusparse_handle(), + CUSPARSE_OPERATION_NON_TRANSPOSE, + reusable_device_scalar_value_1_.data(), + cusparse_view_.A_T, + cusparse_view_.dual_solution, + reusable_device_scalar_value_0_.data(), + cusparse_view_.current_AtY, + CUSPARSE_SPMV_CSR_ALG2, + (f_t*)cusparse_view_.buffer_transpose.data(), + stream_view_)); + } } else { RAFT_CUSPARSE_TRY(raft::sparse::detail::cusparsespmm( handle_ptr_->get_cusparse_handle(), @@ -319,17 +353,33 @@ void pdhg_solver_t::compute_A_x() { // A @ x if (!batch_mode_) { - RAFT_CUSPARSE_TRY( - raft::sparse::detail::cusparsespmv(handle_ptr_->get_cusparse_handle(), - CUSPARSE_OPERATION_NON_TRANSPOSE, - reusable_device_scalar_value_1_.data(), - cusparse_view_.A, - cusparse_view_.reflected_primal_solution, - reusable_device_scalar_value_0_.data(), - cusparse_view_.dual_gradient, - CUSPARSE_SPMV_CSR_ALG2, - (f_t*)cusparse_view_.buffer_non_transpose.data(), - stream_view_)); + if constexpr (std::is_same_v) { + if (cusparse_view_.mixed_precision_enabled_) { + mixed_precision_spmv(handle_ptr_->get_cusparse_handle(), + CUSPARSE_OPERATION_NON_TRANSPOSE, + reusable_device_scalar_value_1_.data(), + cusparse_view_.A_mixed_, + cusparse_view_.reflected_primal_solution, + reusable_device_scalar_value_0_.data(), + cusparse_view_.dual_gradient, + CUSPARSE_SPMV_CSR_ALG2, + cusparse_view_.buffer_non_transpose_mixed_.data(), + stream_view_); + } + } + if (!cusparse_view_.mixed_precision_enabled_) { + RAFT_CUSPARSE_TRY( + raft::sparse::detail::cusparsespmv(handle_ptr_->get_cusparse_handle(), + CUSPARSE_OPERATION_NON_TRANSPOSE, + reusable_device_scalar_value_1_.data(), + cusparse_view_.A, + cusparse_view_.reflected_primal_solution, + reusable_device_scalar_value_0_.data(), + cusparse_view_.dual_gradient, + CUSPARSE_SPMV_CSR_ALG2, + (f_t*)cusparse_view_.buffer_non_transpose.data(), + stream_view_)); + } } else { RAFT_CUSPARSE_TRY(raft::sparse::detail::cusparsespmm( handle_ptr_->get_cusparse_handle(), @@ -1196,7 +1246,7 @@ rmm::device_uvector& pdhg_solver_t::get_dual_solution() return current_saddle_point_state_.get_dual_solution(); } -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template class pdhg_solver_t; #endif #if MIP_INSTANTIATE_DOUBLE diff --git a/cpp/src/pdlp/pdhg.hpp b/cpp/src/pdlp/pdhg.hpp index 8ff45ac0ce..0a64e49efb 100644 --- a/cpp/src/pdlp/pdhg.hpp +++ b/cpp/src/pdlp/pdhg.hpp @@ -29,7 +29,8 @@ class pdhg_solver_t { bool is_legacy_batch_mode, const std::vector& climber_strategies, const pdlp_hyper_params::pdlp_hyper_params_t& hyper_params, - const std::vector>& new_bounds); + const std::vector>& new_bounds, + bool enable_mixed_precision_spmv = false); saddle_point_state_t& get_saddle_point_state(); cusparse_view_t& get_cusparse_view(); diff --git a/cpp/src/pdlp/pdlp.cu b/cpp/src/pdlp/pdlp.cu index cda60cf5ff..82e79098a7 100644 --- a/cpp/src/pdlp/pdlp.cu +++ b/cpp/src/pdlp/pdlp.cu @@ -43,6 +43,59 @@ namespace cuopt::linear_programming::detail { +// Templated wrapper for cuBLAS geam function +// cublasSgeam for float, cublasDgeam for double +template +inline cublasStatus_t cublasGeam(cublasHandle_t handle, + cublasOperation_t transa, + cublasOperation_t transb, + int m, + int n, + const T* alpha, + const T* A, + int lda, + const T* beta, + const T* B, + int ldb, + T* C, + int ldc); + +template <> +inline cublasStatus_t cublasGeam(cublasHandle_t handle, + cublasOperation_t transa, + cublasOperation_t transb, + int m, + int n, + const float* alpha, + const float* A, + int lda, + const float* beta, + const float* B, + int ldb, + float* C, + int ldc) +{ + return cublasSgeam(handle, transa, transb, m, n, alpha, A, lda, beta, B, ldb, C, ldc); +} + +template <> +inline cublasStatus_t cublasGeam(cublasHandle_t handle, + cublasOperation_t transa, + cublasOperation_t transb, + int m, + int n, + const double* alpha, + const double* A, + int lda, + const double* beta, + const double* B, + int ldb, + double* C, + int ldc) +{ + return cublasDgeam(handle, transa, transb, m, n, alpha, A, lda, beta, B, ldb, C, ldc); +} + template static size_t batch_size_handler(const problem_t& op_problem, const pdlp_solver_settings_t& settings) @@ -88,7 +141,8 @@ pdlp_solver_t::pdlp_solver_t(problem_t& op_problem, is_legacy_batch_mode, climber_strategies_, settings_.hyper_params, - settings_.new_bounds}, + settings_.new_bounds, + settings_.pdlp_precision == pdlp_precision_t::MixedPrecision}, initial_scaling_strategy_{handle_ptr_, op_problem_scaled_, settings_.hyper_params.default_l_inf_ruiz_iterations, @@ -1925,48 +1979,48 @@ void pdlp_solver_t::transpose_primal_dual_to_row( is_dual_slack_empty ? 0 : primal_size_h_ * climber_strategies_.size(), stream_view_); RAFT_CUBLAS_TRY(cublasSetStream(handle_ptr_->get_cublas_handle(), stream_view_)); - CUBLAS_CHECK(cublasDgeam(handle_ptr_->get_cublas_handle(), - CUBLAS_OP_T, - CUBLAS_OP_N, - climber_strategies_.size(), - primal_size_h_, - reusable_device_scalar_value_1_.data(), - primal_to_transpose.data(), - primal_size_h_, - reusable_device_scalar_value_0_.data(), - nullptr, - climber_strategies_.size(), - primal_transposed.data(), - climber_strategies_.size())); + CUBLAS_CHECK(cublasGeam(handle_ptr_->get_cublas_handle(), + CUBLAS_OP_T, + CUBLAS_OP_N, + climber_strategies_.size(), + primal_size_h_, + reusable_device_scalar_value_1_.data(), + primal_to_transpose.data(), + primal_size_h_, + reusable_device_scalar_value_0_.data(), + nullptr, + climber_strategies_.size(), + primal_transposed.data(), + climber_strategies_.size())); if (!is_dual_slack_empty) { - CUBLAS_CHECK(cublasDgeam(handle_ptr_->get_cublas_handle(), - CUBLAS_OP_T, - CUBLAS_OP_N, - climber_strategies_.size(), - primal_size_h_, - reusable_device_scalar_value_1_.data(), - dual_slack_to_transpose.data(), - primal_size_h_, - reusable_device_scalar_value_0_.data(), - nullptr, - climber_strategies_.size(), - dual_slack_transposed.data(), - climber_strategies_.size())); + CUBLAS_CHECK(cublasGeam(handle_ptr_->get_cublas_handle(), + CUBLAS_OP_T, + CUBLAS_OP_N, + climber_strategies_.size(), + primal_size_h_, + reusable_device_scalar_value_1_.data(), + dual_slack_to_transpose.data(), + primal_size_h_, + reusable_device_scalar_value_0_.data(), + nullptr, + climber_strategies_.size(), + dual_slack_transposed.data(), + climber_strategies_.size())); } - CUBLAS_CHECK(cublasDgeam(handle_ptr_->get_cublas_handle(), - CUBLAS_OP_T, - CUBLAS_OP_N, - climber_strategies_.size(), - dual_size_h_, - reusable_device_scalar_value_1_.data(), - dual_to_transpose.data(), - dual_size_h_, - reusable_device_scalar_value_0_.data(), - nullptr, - climber_strategies_.size(), - dual_transposed.data(), - climber_strategies_.size())); + CUBLAS_CHECK(cublasGeam(handle_ptr_->get_cublas_handle(), + CUBLAS_OP_T, + CUBLAS_OP_N, + climber_strategies_.size(), + dual_size_h_, + reusable_device_scalar_value_1_.data(), + dual_to_transpose.data(), + dual_size_h_, + reusable_device_scalar_value_0_.data(), + nullptr, + climber_strategies_.size(), + dual_transposed.data(), + climber_strategies_.size())); // Copy that holds the tranpose to the original vector raft::copy(primal_to_transpose.data(), @@ -2002,49 +2056,49 @@ void pdlp_solver_t::transpose_primal_dual_back_to_col( is_dual_slack_empty ? 0 : primal_size_h_ * climber_strategies_.size(), stream_view_); RAFT_CUBLAS_TRY(cublasSetStream(handle_ptr_->get_cublas_handle(), stream_view_)); - CUBLAS_CHECK(cublasDgeam(handle_ptr_->get_cublas_handle(), - CUBLAS_OP_T, - CUBLAS_OP_N, - primal_size_h_, - climber_strategies_.size(), - reusable_device_scalar_value_1_.data(), - primal_to_transpose.data(), - climber_strategies_.size(), - reusable_device_scalar_value_0_.data(), - nullptr, - primal_size_h_, - primal_transposed.data(), - primal_size_h_)); + CUBLAS_CHECK(cublasGeam(handle_ptr_->get_cublas_handle(), + CUBLAS_OP_T, + CUBLAS_OP_N, + primal_size_h_, + climber_strategies_.size(), + reusable_device_scalar_value_1_.data(), + primal_to_transpose.data(), + climber_strategies_.size(), + reusable_device_scalar_value_0_.data(), + nullptr, + primal_size_h_, + primal_transposed.data(), + primal_size_h_)); if (!is_dual_slack_empty) { - CUBLAS_CHECK(cublasDgeam(handle_ptr_->get_cublas_handle(), - CUBLAS_OP_T, - CUBLAS_OP_N, - primal_size_h_, - climber_strategies_.size(), - reusable_device_scalar_value_1_.data(), - dual_slack_to_transpose.data(), - climber_strategies_.size(), - reusable_device_scalar_value_0_.data(), - nullptr, - primal_size_h_, - dual_slack_transposed.data(), - primal_size_h_)); + CUBLAS_CHECK(cublasGeam(handle_ptr_->get_cublas_handle(), + CUBLAS_OP_T, + CUBLAS_OP_N, + primal_size_h_, + climber_strategies_.size(), + reusable_device_scalar_value_1_.data(), + dual_slack_to_transpose.data(), + climber_strategies_.size(), + reusable_device_scalar_value_0_.data(), + nullptr, + primal_size_h_, + dual_slack_transposed.data(), + primal_size_h_)); } - CUBLAS_CHECK(cublasDgeam(handle_ptr_->get_cublas_handle(), - CUBLAS_OP_T, - CUBLAS_OP_N, - dual_size_h_, - climber_strategies_.size(), - reusable_device_scalar_value_1_.data(), - dual_to_transpose.data(), - climber_strategies_.size(), - reusable_device_scalar_value_0_.data(), - nullptr, - dual_size_h_, - dual_transposed.data(), - dual_size_h_)); + CUBLAS_CHECK(cublasGeam(handle_ptr_->get_cublas_handle(), + CUBLAS_OP_T, + CUBLAS_OP_N, + dual_size_h_, + climber_strategies_.size(), + reusable_device_scalar_value_1_.data(), + dual_to_transpose.data(), + climber_strategies_.size(), + reusable_device_scalar_value_0_.data(), + nullptr, + dual_size_h_, + dual_transposed.data(), + dual_size_h_)); // Copy that holds the tranpose to the original vector raft::copy(primal_to_transpose.data(), @@ -2090,6 +2144,9 @@ optimization_problem_solution_t pdlp_solver_t::run_solver(co initial_scaling_strategy_.scale_problem(); + // Update FP32 matrix copies for mixed precision SpMV after scaling + pdhg_solver_.get_cusparse_view().update_mixed_precision_matrices(); + if (!settings_.hyper_params.compute_initial_step_size_before_scaling && !settings_.get_initial_step_size().has_value()) compute_initial_step_size(); @@ -2914,7 +2971,7 @@ pdlp_solver_t::get_current_termination_strategy() return current_termination_strategy_; } -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template class pdlp_solver_t; template __global__ void compute_weights_initial_primal_weight_from_squared_norms( diff --git a/cpp/src/pdlp/pdlp_warm_start_data.cu b/cpp/src/pdlp/pdlp_warm_start_data.cu index 66bfe66914..80abf015d8 100644 --- a/cpp/src/pdlp/pdlp_warm_start_data.cu +++ b/cpp/src/pdlp/pdlp_warm_start_data.cu @@ -178,7 +178,7 @@ void pdlp_warm_start_data_t::check_sizes() "All dual vectors should be of same size"); } -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template class pdlp_warm_start_data_t; #endif diff --git a/cpp/src/pdlp/restart_strategy/localized_duality_gap_container.cu b/cpp/src/pdlp/restart_strategy/localized_duality_gap_container.cu index 2f2e2c7333..bb79e5b6e6 100644 --- a/cpp/src/pdlp/restart_strategy/localized_duality_gap_container.cu +++ b/cpp/src/pdlp/restart_strategy/localized_duality_gap_container.cu @@ -144,7 +144,7 @@ localized_duality_gap_container_t::view() return v; } -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template struct localized_duality_gap_container_t; #endif #if MIP_INSTANTIATE_DOUBLE diff --git a/cpp/src/pdlp/restart_strategy/pdlp_restart_strategy.cu b/cpp/src/pdlp/restart_strategy/pdlp_restart_strategy.cu index 8eacd4d246..149e99a431 100644 --- a/cpp/src/pdlp/restart_strategy/pdlp_restart_strategy.cu +++ b/cpp/src/pdlp/restart_strategy/pdlp_restart_strategy.cu @@ -2523,7 +2523,7 @@ bool pdlp_restart_strategy_t::get_last_restart_was_average() const const typename localized_duality_gap_container_t::view_t duality_gap_view, \ F_TYPE* primal_product); -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT INSTANTIATE(float) #endif diff --git a/cpp/src/pdlp/restart_strategy/weighted_average_solution.cu b/cpp/src/pdlp/restart_strategy/weighted_average_solution.cu index 03c76d79ae..70a448a9de 100644 --- a/cpp/src/pdlp/restart_strategy/weighted_average_solution.cu +++ b/cpp/src/pdlp/restart_strategy/weighted_average_solution.cu @@ -139,7 +139,7 @@ i_t weighted_average_solution_t::get_iterations_since_last_restart() c return iterations_since_last_restart_; } -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template __global__ void add_weight_sums(const float* primal_weight, const float* dual_weight, float* sum_primal_solution_weights, diff --git a/cpp/src/pdlp/saddle_point.cu b/cpp/src/pdlp/saddle_point.cu index c516ab7355..157e7fa389 100644 --- a/cpp/src/pdlp/saddle_point.cu +++ b/cpp/src/pdlp/saddle_point.cu @@ -166,7 +166,7 @@ rmm::device_uvector& saddle_point_state_t::get_next_AtY() return next_AtY_; } -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template class saddle_point_state_t; #endif diff --git a/cpp/src/pdlp/solve.cu b/cpp/src/pdlp/solve.cu index 5e1e25bbee..22fff31906 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -60,6 +60,10 @@ namespace cuopt::linear_programming { +template +extern rmm::device_uvector gpu_cast(const rmm::device_uvector& src, + rmm::cuda_stream_view stream); + // This serves as both a warm up but also a mandatory initial call to setup cuSparse and cuBLAS static void init_handler(const raft::handle_t* handle_ptr) { @@ -560,6 +564,122 @@ optimization_problem_solution_t run_dual_simplex( 0); } +#if PDLP_INSTANTIATE_FLOAT || CUOPT_INSTANTIATE_FLOAT + +template +static optimization_problem_solution_t run_pdlp_solver_in_fp32( + detail::problem_t& problem, + pdlp_solver_settings_t const& settings, + const timer_t& timer, + bool is_batch_mode) +{ + CUOPT_LOG_CONDITIONAL_INFO(!settings.inside_mip, "Running PDLP in FP32 precision"); + auto stream = problem.handle_ptr->get_stream(); + + // Convert the optimization problem stored inside problem_t to float + auto float_op = problem.original_problem_ptr->template convert_to_other_prec(stream); + float_op.set_objective_offset(static_cast(problem.presolve_data.objective_offset)); + float_op.set_objective_scaling_factor( + static_cast(problem.presolve_data.objective_scaling_factor)); + + detail::problem_t float_problem(float_op); + + auto objective_name = problem.objective_name; + auto var_names = problem.var_names; + auto row_names = problem.row_names; + // When crossover is off, free double-precision GPU memory to reduce peak usage. + // When crossover is on, run_pdlp needs the problem data after we return. + if (!settings.crossover) { + { + [[maybe_unused]] auto discard = detail::problem_t(std::move(problem)); + } + } + + // Create float settings from double settings + pdlp_solver_settings_t fs; + fs.tolerances.absolute_dual_tolerance = + static_cast(settings.tolerances.absolute_dual_tolerance); + fs.tolerances.relative_dual_tolerance = + static_cast(settings.tolerances.relative_dual_tolerance); + fs.tolerances.absolute_primal_tolerance = + static_cast(settings.tolerances.absolute_primal_tolerance); + fs.tolerances.relative_primal_tolerance = + static_cast(settings.tolerances.relative_primal_tolerance); + fs.tolerances.absolute_gap_tolerance = + static_cast(settings.tolerances.absolute_gap_tolerance); + fs.tolerances.relative_gap_tolerance = + static_cast(settings.tolerances.relative_gap_tolerance); + fs.tolerances.primal_infeasible_tolerance = + static_cast(settings.tolerances.primal_infeasible_tolerance); + fs.tolerances.dual_infeasible_tolerance = + static_cast(settings.tolerances.dual_infeasible_tolerance); + fs.detect_infeasibility = settings.detect_infeasibility; + fs.strict_infeasibility = settings.strict_infeasibility; + fs.iteration_limit = settings.iteration_limit; + fs.time_limit = static_cast(settings.time_limit); + fs.pdlp_solver_mode = settings.pdlp_solver_mode; + fs.log_to_console = settings.log_to_console; + fs.log_file = settings.log_file; + fs.per_constraint_residual = settings.per_constraint_residual; + fs.save_best_primal_so_far = settings.save_best_primal_so_far; + fs.first_primal_feasible = settings.first_primal_feasible; + fs.eliminate_dense_columns = settings.eliminate_dense_columns; + fs.pdlp_precision = pdlp_precision_t::DefaultPrecision; + fs.method = method_t::PDLP; + fs.inside_mip = settings.inside_mip; + fs.hyper_params = settings.hyper_params; + fs.presolver = settings.presolver; + fs.num_gpus = settings.num_gpus; + fs.concurrent_halt = settings.concurrent_halt; + + detail::pdlp_solver_t solver(float_problem, fs, is_batch_mode); + if (settings.inside_mip) { solver.set_inside_mip(true); } + auto float_sol = solver.run_solver(timer); + + // Convert float solution back to double on GPU (gpu_cast defined in optimization_problem.cu) + auto dev_primal = gpu_cast(float_sol.get_primal_solution(), stream); + auto dev_dual = gpu_cast(float_sol.get_dual_solution(), stream); + auto dev_reduced = gpu_cast(float_sol.get_reduced_cost(), stream); + + // Convert termination info (small host-side struct, stays on CPU) + auto float_term_infos = float_sol.get_additional_termination_informations(); + using double_term_info_t = + typename optimization_problem_solution_t::additional_termination_information_t; + std::vector term_infos; + for (auto& fi : float_term_infos) { + double_term_info_t di; + di.number_of_steps_taken = fi.number_of_steps_taken; + di.total_number_of_attempted_steps = fi.total_number_of_attempted_steps; + di.l2_primal_residual = static_cast(fi.l2_primal_residual); + di.l2_relative_primal_residual = static_cast(fi.l2_relative_primal_residual); + di.l2_dual_residual = static_cast(fi.l2_dual_residual); + di.l2_relative_dual_residual = static_cast(fi.l2_relative_dual_residual); + di.primal_objective = static_cast(fi.primal_objective); + di.dual_objective = static_cast(fi.dual_objective); + di.gap = static_cast(fi.gap); + di.relative_gap = static_cast(fi.relative_gap); + di.max_primal_ray_infeasibility = static_cast(fi.max_primal_ray_infeasibility); + di.primal_ray_linear_objective = static_cast(fi.primal_ray_linear_objective); + di.max_dual_ray_infeasibility = static_cast(fi.max_dual_ray_infeasibility); + di.dual_ray_linear_objective = static_cast(fi.dual_ray_linear_objective); + di.solve_time = fi.solve_time; + di.solved_by_pdlp = fi.solved_by_pdlp; + term_infos.push_back(di); + } + + auto status_vec = float_sol.get_terminations_status(); + + return optimization_problem_solution_t(dev_primal, + dev_dual, + dev_reduced, + objective_name, + var_names, + row_names, + std::move(term_infos), + std::move(status_vec)); +} +#endif + template static optimization_problem_solution_t run_pdlp_solver( detail::problem_t& problem, @@ -574,6 +694,13 @@ static optimization_problem_solution_t run_pdlp_solver( return optimization_problem_solution_t{pdlp_termination_status_t::NumericalError, problem.handle_ptr->get_stream()}; } +#if PDLP_INSTANTIATE_FLOAT || CUOPT_INSTANTIATE_FLOAT + if constexpr (std::is_same_v) { + if (settings.pdlp_precision == pdlp_precision_t::SinglePrecision) { + return run_pdlp_solver_in_fp32(problem, settings, timer, is_batch_mode); + } + } +#endif detail::pdlp_solver_t solver(problem, settings, is_batch_mode); if (settings.inside_mip) { solver.set_inside_mip(true); } return solver.run_solver(timer); @@ -585,6 +712,24 @@ optimization_problem_solution_t run_pdlp(detail::problem_t& const timer_t& timer, bool is_batch_mode) { + if constexpr (!std::is_same_v) { + cuopt_expects(!is_batch_mode, + error_type_t::ValidationError, + "PDLP batch mode is not supported for float precision. Use double precision."); + } + cuopt_expects(!(settings.pdlp_precision == pdlp_precision_t::MixedPrecision && + !detail::is_cusparse_runtime_mixed_precision_supported()), + error_type_t::ValidationError, + "Mixed-precision SpMV requires cuSPARSE runtime 12.5 or later."); + cuopt_expects( + !(is_batch_mode && settings.pdlp_precision == pdlp_precision_t::MixedPrecision), + error_type_t::ValidationError, + "Mixed-precision SpMV is not supported in batch mode. Set pdlp_precision=-1 (default) " + "or disable batch mode."); + cuopt_expects(!(settings.pdlp_precision == pdlp_precision_t::SinglePrecision && is_batch_mode), + error_type_t::ValidationError, + "Single-precision PDLP is not supported in batch mode."); + auto start_solver = std::chrono::high_resolution_clock::now(); timer_t timer_pdlp(timer.remaining_time()); auto sol = run_pdlp_solver(problem, settings, timer, is_batch_mode); @@ -606,82 +751,89 @@ optimization_problem_solution_t run_pdlp(detail::problem_t& sol.get_solve_time()); } - const bool do_crossover = settings.crossover; - i_t crossover_info = 0; - if (do_crossover && sol.get_termination_status() == pdlp_termination_status_t::Optimal) { - crossover_info = -1; - - dual_simplex::lp_problem_t lp(problem.handle_ptr, 1, 1, 1); - dual_simplex::lp_solution_t initial_solution(1, 1); - translate_to_crossover_problem(problem, sol, lp, initial_solution); - dual_simplex::simplex_solver_settings_t dual_simplex_settings; - dual_simplex_settings.time_limit = settings.time_limit; - dual_simplex_settings.iteration_limit = settings.iteration_limit; - dual_simplex_settings.concurrent_halt = settings.concurrent_halt; - dual_simplex::lp_solution_t vertex_solution(lp.num_rows, lp.num_cols); - std::vector vstatus(lp.num_cols); - dual_simplex::crossover_status_t crossover_status = dual_simplex::crossover( - lp, dual_simplex_settings, initial_solution, timer.get_tic_start(), vertex_solution, vstatus); - pdlp_termination_status_t termination_status = pdlp_termination_status_t::TimeLimit; - auto to_termination_status = [](dual_simplex::crossover_status_t status) { - switch (status) { - case dual_simplex::crossover_status_t::OPTIMAL: return pdlp_termination_status_t::Optimal; - case dual_simplex::crossover_status_t::PRIMAL_FEASIBLE: - return pdlp_termination_status_t::PrimalFeasible; - case dual_simplex::crossover_status_t::DUAL_FEASIBLE: - return pdlp_termination_status_t::NumericalError; - case dual_simplex::crossover_status_t::NUMERICAL_ISSUES: - return pdlp_termination_status_t::NumericalError; - case dual_simplex::crossover_status_t::CONCURRENT_LIMIT: - return pdlp_termination_status_t::ConcurrentLimit; - case dual_simplex::crossover_status_t::TIME_LIMIT: - return pdlp_termination_status_t::TimeLimit; - default: return pdlp_termination_status_t::NumericalError; - } - }; - termination_status = to_termination_status(crossover_status); - if (crossover_status == dual_simplex::crossover_status_t::OPTIMAL) { crossover_info = 0; } - rmm::device_uvector final_primal_solution = - cuopt::device_copy(vertex_solution.x, problem.handle_ptr->get_stream()); - rmm::device_uvector final_dual_solution = - cuopt::device_copy(vertex_solution.y, problem.handle_ptr->get_stream()); - rmm::device_uvector final_reduced_cost = - cuopt::device_copy(vertex_solution.z, problem.handle_ptr->get_stream()); - problem.handle_ptr->sync_stream(); - // Negate dual variables and reduced costs for maximization problems - if (problem.maximize) { - adjust_dual_solution_and_reduced_cost( - final_dual_solution, final_reduced_cost, problem.handle_ptr->get_stream()); + if constexpr (std::is_same_v) { + const bool do_crossover = settings.crossover; + i_t crossover_info = 0; + if (do_crossover && sol.get_termination_status() == pdlp_termination_status_t::Optimal) { + crossover_info = -1; + + dual_simplex::lp_problem_t lp(problem.handle_ptr, 1, 1, 1); + dual_simplex::lp_solution_t initial_solution(1, 1); + translate_to_crossover_problem(problem, sol, lp, initial_solution); + dual_simplex::simplex_solver_settings_t dual_simplex_settings; + dual_simplex_settings.time_limit = settings.time_limit; + dual_simplex_settings.iteration_limit = settings.iteration_limit; + dual_simplex_settings.concurrent_halt = settings.concurrent_halt; + dual_simplex::lp_solution_t vertex_solution(lp.num_rows, lp.num_cols); + std::vector vstatus(lp.num_cols); + dual_simplex::crossover_status_t crossover_status = + dual_simplex::crossover(lp, + dual_simplex_settings, + initial_solution, + timer.get_tic_start(), + vertex_solution, + vstatus); + pdlp_termination_status_t termination_status = pdlp_termination_status_t::TimeLimit; + auto to_termination_status = [](dual_simplex::crossover_status_t status) { + switch (status) { + case dual_simplex::crossover_status_t::OPTIMAL: return pdlp_termination_status_t::Optimal; + case dual_simplex::crossover_status_t::PRIMAL_FEASIBLE: + return pdlp_termination_status_t::PrimalFeasible; + case dual_simplex::crossover_status_t::DUAL_FEASIBLE: + return pdlp_termination_status_t::NumericalError; + case dual_simplex::crossover_status_t::NUMERICAL_ISSUES: + return pdlp_termination_status_t::NumericalError; + case dual_simplex::crossover_status_t::CONCURRENT_LIMIT: + return pdlp_termination_status_t::ConcurrentLimit; + case dual_simplex::crossover_status_t::TIME_LIMIT: + return pdlp_termination_status_t::TimeLimit; + default: return pdlp_termination_status_t::NumericalError; + } + }; + termination_status = to_termination_status(crossover_status); + if (crossover_status == dual_simplex::crossover_status_t::OPTIMAL) { crossover_info = 0; } + rmm::device_uvector final_primal_solution = + cuopt::device_copy(vertex_solution.x, problem.handle_ptr->get_stream()); + rmm::device_uvector final_dual_solution = + cuopt::device_copy(vertex_solution.y, problem.handle_ptr->get_stream()); + rmm::device_uvector final_reduced_cost = + cuopt::device_copy(vertex_solution.z, problem.handle_ptr->get_stream()); problem.handle_ptr->sync_stream(); - } + // Negate dual variables and reduced costs for maximization problems + if (problem.maximize) { + adjust_dual_solution_and_reduced_cost( + final_dual_solution, final_reduced_cost, problem.handle_ptr->get_stream()); + problem.handle_ptr->sync_stream(); + } - // Should be filled with more information from dual simplex - std::vector< - typename optimization_problem_solution_t::additional_termination_information_t> - info(1); - info[0].primal_objective = vertex_solution.user_objective; - info[0].number_of_steps_taken = vertex_solution.iterations; - auto crossover_end = std::chrono::high_resolution_clock::now(); - auto crossover_duration = - std::chrono::duration_cast(crossover_end - start_solver); - info[0].solve_time = crossover_duration.count() / 1000.0; - auto sol_crossover = optimization_problem_solution_t(final_primal_solution, - final_dual_solution, - final_reduced_cost, - problem.objective_name, - problem.var_names, - problem.row_names, - std::move(info), - {termination_status}); - sol.copy_from(problem.handle_ptr, sol_crossover); - CUOPT_LOG_CONDITIONAL_INFO( - !settings.inside_mip, "Crossover status %s", sol.get_termination_status_string().c_str()); - } - if (settings.method == method_t::Concurrent && settings.concurrent_halt != nullptr && - crossover_info == 0 && sol.get_termination_status() == pdlp_termination_status_t::Optimal) { - // We finished. Tell dual simplex to stop if it is still running. - CUOPT_LOG_CONDITIONAL_INFO(!settings.inside_mip, "PDLP finished. Telling others to stop"); - *settings.concurrent_halt = 1; + // Should be filled with more information from dual simplex + std::vector< + typename optimization_problem_solution_t::additional_termination_information_t> + info(1); + info[0].primal_objective = vertex_solution.user_objective; + info[0].number_of_steps_taken = vertex_solution.iterations; + auto crossover_end = std::chrono::high_resolution_clock::now(); + auto crossover_duration = + std::chrono::duration_cast(crossover_end - start_solver); + info[0].solve_time = crossover_duration.count() / 1000.0; + auto sol_crossover = optimization_problem_solution_t(final_primal_solution, + final_dual_solution, + final_reduced_cost, + problem.objective_name, + problem.var_names, + problem.row_names, + std::move(info), + {termination_status}); + sol.copy_from(problem.handle_ptr, sol_crossover); + CUOPT_LOG_CONDITIONAL_INFO( + !settings.inside_mip, "Crossover status %s", sol.get_termination_status_string().c_str()); + } + if (settings.method == method_t::Concurrent && settings.concurrent_halt != nullptr && + crossover_info == 0 && sol.get_termination_status() == pdlp_termination_status_t::Optimal) { + // We finished. Tell dual simplex to stop if it is still running. + CUOPT_LOG_CONDITIONAL_INFO(!settings.inside_mip, "PDLP finished. Telling others to stop"); + *settings.concurrent_halt = 1; + } } return sol; } @@ -1117,13 +1269,22 @@ optimization_problem_solution_t solve_lp_with_method( const timer_t& timer, bool is_batch_mode) { - if (settings.method == method_t::DualSimplex) { - return run_dual_simplex(problem, settings, timer); - } else if (settings.method == method_t::Barrier) { - return run_barrier(problem, settings, timer); - } else if (settings.method == method_t::Concurrent) { - return run_concurrent(problem, settings, timer, is_batch_mode); + if constexpr (std::is_same_v) { + if (settings.method == method_t::DualSimplex) { + return run_dual_simplex(problem, settings, timer); + } else if (settings.method == method_t::Barrier) { + return run_barrier(problem, settings, timer); + } else if (settings.method == method_t::Concurrent) { + return run_concurrent(problem, settings, timer, is_batch_mode); + } else { + return run_pdlp(problem, settings, timer, is_batch_mode); + } } else { + // Float precision only supports PDLP without presolve/crossover + cuopt_expects(settings.method == method_t::PDLP, + error_type_t::ValidationError, + "Float precision only supports PDLP method. DualSimplex, Barrier, and Concurrent " + "require double precision."); return run_pdlp(problem, settings, timer, is_batch_mode); } } @@ -1199,9 +1360,6 @@ optimization_problem_solution_t solve_lp( std::unique_ptr> presolver; auto run_presolve = settings.presolver != presolver_t::None; run_presolve = run_presolve && settings.get_pdlp_warm_start_data().total_pdlp_iterations_ == -1; - if (!run_presolve && !settings_const.inside_mip) { - CUOPT_LOG_INFO("Third-party presolve is disabled, skipping"); - } // Declare result at outer scope so that result->reduced_problem (which may be // referenced by problem.original_problem_ptr) remains alive through the solve. diff --git a/cpp/src/pdlp/solver_settings.cu b/cpp/src/pdlp/solver_settings.cu index 560e40f302..7acfc7481c 100644 --- a/cpp/src/pdlp/solver_settings.cu +++ b/cpp/src/pdlp/solver_settings.cu @@ -382,7 +382,7 @@ pdlp_solver_settings_t::get_pdlp_warm_start_data_view() const noexcept return pdlp_warm_start_data_view_; } -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template class pdlp_solver_settings_t; #endif diff --git a/cpp/src/pdlp/solver_solution.cu b/cpp/src/pdlp/solver_solution.cu index a8001b91c1..10e6a80593 100644 --- a/cpp/src/pdlp/solver_solution.cu +++ b/cpp/src/pdlp/solver_solution.cu @@ -448,7 +448,7 @@ void optimization_problem_solution_t::write_to_sol_file( std::string(filename), status, objective_value, var_names_, solution); } -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template class optimization_problem_solution_t; #endif diff --git a/cpp/src/pdlp/step_size_strategy/adaptive_step_size_strategy.cu b/cpp/src/pdlp/step_size_strategy/adaptive_step_size_strategy.cu index 24ef29b243..d17a88dd29 100644 --- a/cpp/src/pdlp/step_size_strategy/adaptive_step_size_strategy.cu +++ b/cpp/src/pdlp/step_size_strategy/adaptive_step_size_strategy.cu @@ -578,7 +578,7 @@ adaptive_step_size_strategy_t::view() F_TYPE * dual_step_size, \ int* pdhg_iteration); -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT INSTANTIATE(float) #endif diff --git a/cpp/src/pdlp/termination_strategy/convergence_information.cu b/cpp/src/pdlp/termination_strategy/convergence_information.cu index 9b01608b47..ab0c921cc7 100644 --- a/cpp/src/pdlp/termination_strategy/convergence_information.cu +++ b/cpp/src/pdlp/termination_strategy/convergence_information.cu @@ -996,7 +996,7 @@ convergence_information_t::to_primal_quality_adapter( primal_objective_.element(0, stream_view_)}; } -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template class convergence_information_t; template __global__ void compute_remaining_stats_kernel( diff --git a/cpp/src/pdlp/termination_strategy/infeasibility_information.cu b/cpp/src/pdlp/termination_strategy/infeasibility_information.cu index 14114d306f..dbb35b732d 100644 --- a/cpp/src/pdlp/termination_strategy/infeasibility_information.cu +++ b/cpp/src/pdlp/termination_strategy/infeasibility_information.cu @@ -745,7 +745,7 @@ typename infeasibility_information_t::view_t infeasibility_information return v; } -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template class infeasibility_information_t; template __global__ void compute_remaining_stats_kernel( diff --git a/cpp/src/pdlp/termination_strategy/termination_strategy.cu b/cpp/src/pdlp/termination_strategy/termination_strategy.cu index 033cbdbfda..7179df6a49 100644 --- a/cpp/src/pdlp/termination_strategy/termination_strategy.cu +++ b/cpp/src/pdlp/termination_strategy/termination_strategy.cu @@ -681,7 +681,7 @@ void pdlp_termination_strategy_t::print_termination_criteria(i_t itera bool per_constraint_residual, \ int batch_size); -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT INSTANTIATE(float) #endif diff --git a/cpp/src/pdlp/translate.hpp b/cpp/src/pdlp/translate.hpp index b8e0075733..b143a206d4 100644 --- a/cpp/src/pdlp/translate.hpp +++ b/cpp/src/pdlp/translate.hpp @@ -133,7 +133,7 @@ void translate_to_crossover_problem(const detail::problem_t& problem, std::vector slack(problem.n_constraints); std::vector tmp_x = cuopt::host_copy(sol.get_primal_solution(), stream); stream.synchronize(); - dual_simplex::matrix_vector_multiply(lp.A, 1.0, tmp_x, 0.0, slack); + dual_simplex::matrix_vector_multiply(lp.A, f_t(1.0), tmp_x, f_t(0.0), slack); CUOPT_LOG_DEBUG("Multiplied A and x"); lp.A.col_start.resize(problem.n_variables + problem.n_constraints + 1); diff --git a/cpp/src/pdlp/utilities/problem_checking.cu b/cpp/src/pdlp/utilities/problem_checking.cu index f970f8740d..b10850de27 100644 --- a/cpp/src/pdlp/utilities/problem_checking.cu +++ b/cpp/src/pdlp/utilities/problem_checking.cu @@ -340,7 +340,7 @@ bool problem_checking_t::has_crossing_bounds( #define INSTANTIATE(F_TYPE) template class problem_checking_t; -#if MIP_INSTANTIATE_FLOAT +#if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT INSTANTIATE(float) #endif diff --git a/cpp/tests/linear_programming/c_api_tests/c_api_test.c b/cpp/tests/linear_programming/c_api_tests/c_api_test.c index ecf610041c..996d60deae 100644 --- a/cpp/tests/linear_programming/c_api_tests/c_api_test.c +++ b/cpp/tests/linear_programming/c_api_tests/c_api_test.c @@ -2122,3 +2122,135 @@ cuopt_int_t test_cpu_only_mip_execution(const char* filename) cuOptDestroySolution(&solution); return status; } + +cuopt_int_t test_pdlp_precision_mixed(const char* filename, + cuopt_int_t* termination_status_ptr, + cuopt_float_t* objective_ptr) +{ + cuOptOptimizationProblem problem = NULL; + cuOptSolverSettings settings = NULL; + cuOptSolution solution = NULL; + cuopt_int_t status; + cuopt_int_t termination_status = -1; + cuopt_float_t objective_value; + + status = cuOptReadProblem(filename, &problem); + if (status != CUOPT_SUCCESS) { + printf("Error reading problem\n"); + goto DONE; + } + + status = cuOptCreateSolverSettings(&settings); + if (status != CUOPT_SUCCESS) { + printf("Error creating solver settings\n"); + goto DONE; + } + + status = cuOptSetIntegerParameter(settings, CUOPT_METHOD, CUOPT_METHOD_PDLP); + if (status != CUOPT_SUCCESS) { + printf("Error setting method\n"); + goto DONE; + } + + status = cuOptSetIntegerParameter(settings, CUOPT_PDLP_PRECISION, CUOPT_PDLP_MIXED_PRECISION); + if (status != CUOPT_SUCCESS) { + printf("Error setting pdlp_precision\n"); + goto DONE; + } + + status = cuOptSolve(problem, settings, &solution); + if (status != CUOPT_SUCCESS) { + printf("Error solving problem with pdlp_precision=mixed\n"); + goto DONE; + } + + status = cuOptGetTerminationStatus(solution, &termination_status); + if (status != CUOPT_SUCCESS) { + printf("Error getting termination status\n"); + goto DONE; + } + *termination_status_ptr = termination_status; + + status = cuOptGetObjectiveValue(solution, &objective_value); + if (status != CUOPT_SUCCESS) { + printf("Error getting objective value\n"); + goto DONE; + } + *objective_ptr = objective_value; + + printf("PDLP precision=mixed test passed: status=%s, objective=%f\n", + termination_status_to_string(termination_status), + objective_value); + +DONE: + cuOptDestroyProblem(&problem); + cuOptDestroySolverSettings(&settings); + cuOptDestroySolution(&solution); + return status; +} + +cuopt_int_t test_pdlp_precision_single(const char* filename, + cuopt_int_t* termination_status_ptr, + cuopt_float_t* objective_ptr) +{ + cuOptOptimizationProblem problem = NULL; + cuOptSolverSettings settings = NULL; + cuOptSolution solution = NULL; + cuopt_int_t status; + cuopt_int_t termination_status = -1; + cuopt_float_t objective_value; + + status = cuOptReadProblem(filename, &problem); + if (status != CUOPT_SUCCESS) { + printf("Error reading problem\n"); + goto DONE; + } + + status = cuOptCreateSolverSettings(&settings); + if (status != CUOPT_SUCCESS) { + printf("Error creating solver settings\n"); + goto DONE; + } + + status = cuOptSetIntegerParameter(settings, CUOPT_METHOD, CUOPT_METHOD_PDLP); + if (status != CUOPT_SUCCESS) { + printf("Error setting method\n"); + goto DONE; + } + + status = cuOptSetIntegerParameter(settings, CUOPT_PDLP_PRECISION, CUOPT_PDLP_SINGLE_PRECISION); + if (status != CUOPT_SUCCESS) { + printf("Error setting pdlp_precision\n"); + goto DONE; + } + + status = cuOptSolve(problem, settings, &solution); + if (status != CUOPT_SUCCESS) { + printf("Error solving problem with pdlp_precision=single\n"); + goto DONE; + } + + status = cuOptGetTerminationStatus(solution, &termination_status); + if (status != CUOPT_SUCCESS) { + printf("Error getting termination status\n"); + goto DONE; + } + *termination_status_ptr = termination_status; + + status = cuOptGetObjectiveValue(solution, &objective_value); + if (status != CUOPT_SUCCESS) { + printf("Error getting objective value\n"); + goto DONE; + } + *objective_ptr = objective_value; + + printf("PDLP precision=single test passed: status=%s, objective=%f\n", + termination_status_to_string(termination_status), + objective_value); + +DONE: + cuOptDestroyProblem(&problem); + cuOptDestroySolverSettings(&settings); + cuOptDestroySolution(&solution); + return status; +} diff --git a/cpp/tests/linear_programming/c_api_tests/c_api_tests.cpp b/cpp/tests/linear_programming/c_api_tests/c_api_tests.cpp index 33fb42cc9d..d39a970763 100644 --- a/cpp/tests/linear_programming/c_api_tests/c_api_tests.cpp +++ b/cpp/tests/linear_programming/c_api_tests/c_api_tests.cpp @@ -18,6 +18,10 @@ #include #include +namespace cuopt::linear_programming::detail { +bool is_cusparse_runtime_mixed_precision_supported(); +} + #include TEST(c_api, int_size) { EXPECT_EQ(test_int_size(), sizeof(int32_t)); } @@ -271,6 +275,43 @@ INSTANTIATE_TEST_SUITE_P(c_api, // Different instance std::make_tuple("/mip/bb_optimality.mps", 8, 60.0, 2))); +// ============================================================================= +// PDLP Precision Tests +// ============================================================================= + +TEST(c_api, pdlp_precision_single) +{ + const std::string& rapidsDatasetRootDir = cuopt::test::get_rapids_dataset_root_dir(); + std::string filename = rapidsDatasetRootDir + "/linear_programming/afiro_original.mps"; + cuopt_int_t termination_status; + cuopt_float_t objective; + EXPECT_EQ(test_pdlp_precision_single(filename.c_str(), &termination_status, &objective), + CUOPT_SUCCESS); + EXPECT_EQ(termination_status, CUOPT_TERIMINATION_STATUS_OPTIMAL); + EXPECT_NEAR(objective, -464.7531, 1e-1); +} + +TEST(c_api, pdlp_precision_mixed) +{ + using namespace cuopt::linear_programming::detail; + const std::string& rapidsDatasetRootDir = cuopt::test::get_rapids_dataset_root_dir(); + std::string filename = rapidsDatasetRootDir + "/linear_programming/afiro_original.mps"; + cuopt_int_t termination_status = -1; + cuopt_float_t objective; + if (!is_cusparse_runtime_mixed_precision_supported()) { + auto status = test_pdlp_precision_mixed(filename.c_str(), &termination_status, &objective); + bool solve_returned_error = (status != CUOPT_SUCCESS); + bool solve_returned_non_optimal = + (status == CUOPT_SUCCESS && termination_status != CUOPT_TERIMINATION_STATUS_OPTIMAL); + EXPECT_TRUE(solve_returned_error || solve_returned_non_optimal); + return; + } + EXPECT_EQ(test_pdlp_precision_mixed(filename.c_str(), &termination_status, &objective), + CUOPT_SUCCESS); + EXPECT_EQ(termination_status, CUOPT_TERIMINATION_STATUS_OPTIMAL); + EXPECT_NEAR(objective, -464.7531, 1e-1); +} + // ============================================================================= // Solution Interface Polymorphism Tests // ============================================================================= diff --git a/cpp/tests/linear_programming/c_api_tests/c_api_tests.h b/cpp/tests/linear_programming/c_api_tests/c_api_tests.h index e541316567..402c7d06a5 100644 --- a/cpp/tests/linear_programming/c_api_tests/c_api_tests.h +++ b/cpp/tests/linear_programming/c_api_tests/c_api_tests.h @@ -53,6 +53,14 @@ cuopt_int_t test_deterministic_bb(const char* filename, cuopt_int_t test_lp_solution_mip_methods(); cuopt_int_t test_mip_solution_lp_methods(); +cuopt_int_t test_pdlp_precision_single(const char* filename, + cuopt_int_t* termination_status_ptr, + cuopt_float_t* objective_ptr); + +cuopt_int_t test_pdlp_precision_mixed(const char* filename, + cuopt_int_t* termination_status_ptr, + cuopt_float_t* objective_ptr); + /* CPU-only execution tests (require env vars CUDA_VISIBLE_DEVICES="" and CUOPT_REMOTE_HOST) */ cuopt_int_t test_cpu_only_execution(const char* filename); cuopt_int_t test_cpu_only_mip_execution(const char* filename); diff --git a/cpp/tests/linear_programming/pdlp_test.cu b/cpp/tests/linear_programming/pdlp_test.cu index 8bf759367e..d5a8d69008 100644 --- a/cpp/tests/linear_programming/pdlp_test.cu +++ b/cpp/tests/linear_programming/pdlp_test.cu @@ -6,6 +6,7 @@ /* clang-format on */ #include +#include #include #include #include @@ -20,6 +21,7 @@ #include #include #include +#include #include #include @@ -45,10 +47,10 @@ namespace cuopt::linear_programming::test { -constexpr double afiro_primal_objective = -464; - +constexpr double afiro_primal_objective = -464.0; // Accept a 1% error -static bool is_incorrect_objective(double reference, double objective) +template +static bool is_incorrect_objective(f_t reference, f_t objective) { if (reference == 0) { return std::abs(objective) > 0.01; } if (objective == 0) { return std::abs(reference) > 0.01; } @@ -73,6 +75,58 @@ TEST(pdlp_class, run_double) afiro_primal_objective, solution.get_additional_termination_information().primal_objective)); } +TEST(pdlp_class, precision_mixed) +{ + using namespace cuopt::linear_programming::detail; + if (!is_cusparse_runtime_mixed_precision_supported()) { + const raft::handle_t handle_{}; + auto path = make_path_absolute("linear_programming/afiro_original.mps"); + cuopt::mps_parser::mps_data_model_t op_problem = + cuopt::mps_parser::parse_mps(path, true); + + auto settings = pdlp_solver_settings_t{}; + settings.method = cuopt::linear_programming::method_t::PDLP; + settings.pdlp_precision = cuopt::linear_programming::pdlp_precision_t::MixedPrecision; + + optimization_problem_solution_t solution = + solve_lp(&handle_, op_problem, settings); + EXPECT_EQ(solution.get_error_status().get_error_type(), cuopt::error_type_t::ValidationError); + return; + } + + const raft::handle_t handle_{}; + + auto path = make_path_absolute("linear_programming/afiro_original.mps"); + cuopt::mps_parser::mps_data_model_t op_problem = + cuopt::mps_parser::parse_mps(path, true); + + auto settings_mixed = pdlp_solver_settings_t{}; + settings_mixed.method = cuopt::linear_programming::method_t::PDLP; + settings_mixed.pdlp_precision = cuopt::linear_programming::pdlp_precision_t::MixedPrecision; + + optimization_problem_solution_t solution_mixed = + solve_lp(&handle_, op_problem, settings_mixed); + EXPECT_EQ((int)solution_mixed.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL); + EXPECT_FALSE(is_incorrect_objective( + afiro_primal_objective, + solution_mixed.get_additional_termination_information().primal_objective)); + + auto settings_full = pdlp_solver_settings_t{}; + settings_full.method = cuopt::linear_programming::method_t::PDLP; + settings_full.pdlp_precision = cuopt::linear_programming::pdlp_precision_t::DefaultPrecision; + + optimization_problem_solution_t solution_full = + solve_lp(&handle_, op_problem, settings_full); + EXPECT_EQ((int)solution_full.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL); + EXPECT_FALSE(is_incorrect_objective( + afiro_primal_objective, + solution_full.get_additional_termination_information().primal_objective)); + + EXPECT_NEAR(solution_mixed.get_additional_termination_information().primal_objective, + solution_full.get_additional_termination_information().primal_objective, + 1e-2); +} + TEST(pdlp_class, run_double_very_low_accuracy) { const raft::handle_t handle_{}; @@ -1888,6 +1942,107 @@ TEST(pdlp_class, some_climber_hit_iteration_limit) } } +TEST(pdlp_class, precision_single) +{ + const raft::handle_t handle_{}; + + auto path = make_path_absolute("linear_programming/afiro_original.mps"); + cuopt::mps_parser::mps_data_model_t op_problem = + cuopt::mps_parser::parse_mps(path, true); + + auto solver_settings = pdlp_solver_settings_t{}; + solver_settings.method = cuopt::linear_programming::method_t::PDLP; + solver_settings.pdlp_precision = cuopt::linear_programming::pdlp_precision_t::SinglePrecision; + + optimization_problem_solution_t solution = + solve_lp(&handle_, op_problem, solver_settings); + EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL); + + EXPECT_FALSE(is_incorrect_objective( + afiro_primal_objective, solution.get_additional_termination_information().primal_objective)); +} + +TEST(pdlp_class, precision_single_crossover) +{ + const raft::handle_t handle_{}; + + auto path = make_path_absolute("linear_programming/afiro_original.mps"); + cuopt::mps_parser::mps_data_model_t op_problem = + cuopt::mps_parser::parse_mps(path, true); + + auto solver_settings = pdlp_solver_settings_t{}; + solver_settings.method = cuopt::linear_programming::method_t::PDLP; + solver_settings.pdlp_precision = cuopt::linear_programming::pdlp_precision_t::SinglePrecision; + solver_settings.crossover = true; + + optimization_problem_solution_t solution = + solve_lp(&handle_, op_problem, solver_settings); + EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL); + + EXPECT_FALSE(is_incorrect_objective( + afiro_primal_objective, solution.get_additional_termination_information().primal_objective)); +} + +TEST(pdlp_class, precision_single_concurrent) +{ + const raft::handle_t handle_{}; + + auto path = make_path_absolute("linear_programming/afiro_original.mps"); + cuopt::mps_parser::mps_data_model_t op_problem = + cuopt::mps_parser::parse_mps(path, true); + + auto solver_settings = pdlp_solver_settings_t{}; + solver_settings.method = cuopt::linear_programming::method_t::Concurrent; + solver_settings.pdlp_precision = cuopt::linear_programming::pdlp_precision_t::SinglePrecision; + + optimization_problem_solution_t solution = + solve_lp(&handle_, op_problem, solver_settings); + EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL); + + EXPECT_FALSE(is_incorrect_objective( + afiro_primal_objective, solution.get_additional_termination_information().primal_objective)); +} + +TEST(pdlp_class, precision_single_papilo_presolve) +{ + const raft::handle_t handle_{}; + + auto path = make_path_absolute("linear_programming/afiro_original.mps"); + cuopt::mps_parser::mps_data_model_t op_problem = + cuopt::mps_parser::parse_mps(path, true); + + auto solver_settings = pdlp_solver_settings_t{}; + solver_settings.method = cuopt::linear_programming::method_t::PDLP; + solver_settings.pdlp_precision = cuopt::linear_programming::pdlp_precision_t::SinglePrecision; + solver_settings.presolver = cuopt::linear_programming::presolver_t::Papilo; + + optimization_problem_solution_t solution = + solve_lp(&handle_, op_problem, solver_settings); + EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL); + EXPECT_FALSE(is_incorrect_objective( + afiro_primal_objective, solution.get_additional_termination_information().primal_objective)); +} + +TEST(pdlp_class, precision_single_pslp_presolve) +{ + const raft::handle_t handle_{}; + + auto path = make_path_absolute("linear_programming/afiro_original.mps"); + cuopt::mps_parser::mps_data_model_t op_problem = + cuopt::mps_parser::parse_mps(path, true); + + auto solver_settings = pdlp_solver_settings_t{}; + solver_settings.method = cuopt::linear_programming::method_t::PDLP; + solver_settings.pdlp_precision = cuopt::linear_programming::pdlp_precision_t::SinglePrecision; + solver_settings.presolver = cuopt::linear_programming::presolver_t::PSLP; + + optimization_problem_solution_t solution = + solve_lp(&handle_, op_problem, solver_settings); + EXPECT_EQ((int)solution.get_termination_status(), CUOPT_TERIMINATION_STATUS_OPTIMAL); + EXPECT_FALSE(is_incorrect_objective( + afiro_primal_objective, solution.get_additional_termination_information().primal_objective)); +} + } // namespace cuopt::linear_programming::test CUOPT_TEST_PROGRAM_MAIN() diff --git a/docs/cuopt/source/cuopt-c/lp-qp-milp/lp-qp-milp-c-api.rst b/docs/cuopt/source/cuopt-c/lp-qp-milp/lp-qp-milp-c-api.rst index 43d15eca64..d9a42301cb 100644 --- a/docs/cuopt/source/cuopt-c/lp-qp-milp/lp-qp-milp-c-api.rst +++ b/docs/cuopt/source/cuopt-c/lp-qp-milp/lp-qp-milp-c-api.rst @@ -187,6 +187,7 @@ These constants are used as parameter names in the :c:func:`cuOptSetParameter`, .. doxygendefine:: CUOPT_SOLUTION_FILE .. doxygendefine:: CUOPT_NUM_CPU_THREADS .. doxygendefine:: CUOPT_USER_PROBLEM_FILE +.. doxygendefine:: CUOPT_PDLP_PRECISION .. _pdlp-solver-mode-constants: @@ -201,6 +202,18 @@ These constants are used to configure `CUOPT_PDLP_SOLVER_MODE` via :c:func:`cuOp .. doxygendefine:: CUOPT_PDLP_SOLVER_MODE_METHODICAL1 .. doxygendefine:: CUOPT_PDLP_SOLVER_MODE_FAST1 +.. _pdlp-precision-constants: + +PDLP Precision Constants +------------------------ + +These constants are used to configure `CUOPT_PDLP_PRECISION` via :c:func:`cuOptSetIntegerParameter`. + +.. doxygendefine:: CUOPT_PDLP_DEFAULT_PRECISION +.. doxygendefine:: CUOPT_PDLP_SINGLE_PRECISION +.. doxygendefine:: CUOPT_PDLP_DOUBLE_PRECISION +.. doxygendefine:: CUOPT_PDLP_MIXED_PRECISION + .. _method-constants: Method Constants diff --git a/docs/cuopt/source/lp-qp-features.rst b/docs/cuopt/source/lp-qp-features.rst index 4bd178ed53..e3cbddbb05 100644 --- a/docs/cuopt/source/lp-qp-features.rst +++ b/docs/cuopt/source/lp-qp-features.rst @@ -157,6 +157,17 @@ Batch Mode Users can submit a set of problems which will be solved in a batch. Problems will be solved at the same time in parallel to fully utilize the GPU. Checkout :ref:`self-hosted client ` example in thin client. +PDLP Precision Modes +-------------------- + +By default, PDLP operates in the native precision of the problem type (FP64 for double-precision problems). The ``pdlp_precision`` parameter provides several modes: + +- **single**: Run PDLP internally in FP32, with automatic conversion of inputs and outputs. FP32 uses half the memory and allows PDHG iterations to be on average twice as fast, but may require more iterations to converge. Compatible with crossover (solution is converted back to FP64 before crossover) and concurrent mode (PDLP runs in FP32 while other solvers run in FP64). +- **mixed**: Use mixed precision SpMV during PDHG iterations. The constraint matrix is stored in FP32 while vectors and compute type remain in FP64, improving SpMV performance with limited impact on convergence. Convergence checking and restart logic always use the full FP64 matrix. +- **double**: Explicitly run in FP64 (same as default for double-precision problems). + +.. note:: The default precision is the native type of the problem (FP64 for double). + Multi-GPU Mode -------------- diff --git a/docs/cuopt/source/lp-qp-milp-settings.rst b/docs/cuopt/source/lp-qp-milp-settings.rst index bd1372f70e..29c27a4ac2 100644 --- a/docs/cuopt/source/lp-qp-milp-settings.rst +++ b/docs/cuopt/source/lp-qp-milp-settings.rst @@ -192,6 +192,27 @@ Per Constraint Residual .. note:: The default value is false. +PDLP Precision +^^^^^^^^^^^^^^ + +``CUOPT_PDLP_PRECISION`` controls the precision mode used by the PDLP solver. The following modes are +available: + +- **default** (-1): Use the native precision of the problem type (FP64 for double-precision problems). +- **single** (0): Run PDLP internally in FP32 (float). Inputs are converted from FP64 to FP32 before + solving and outputs are converted back to FP64. FP32 uses half the memory and allows PDHG iterations + to be on average twice as fast, but may require more iterations to converge due to reduced numerical + accuracy. Compatible with crossover (solution is converted back to FP64 before crossover runs) and + concurrent mode (the PDLP leg runs in FP32 while Dual Simplex and Barrier run in FP64). +- **double** (1): Explicitly run in FP64 (same as default for double-precision problems). +- **mixed** (2): Use mixed precision sparse matrix-vector products (SpMV) during PDHG iterations. The + constraint matrix and its transpose are stored in FP32 while vectors and the compute type remain in + FP64, improving SpMV performance. Convergence checking and restart logic always use the + full FP64 matrix, so this mode does not reduce overall memory usage. This provides a middle ground + between full FP64 and FP32: faster PDHG iterations with limited impact on convergence. + +.. note:: The default value is 0 (default precision). + Barrier Solver Settings ^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/python/cuopt/cuopt/tests/linear_programming/test_lp_solver.py b/python/cuopt/cuopt/tests/linear_programming/test_lp_solver.py index 9f94916ff0..e284ffc0ab 100644 --- a/python/cuopt/cuopt/tests/linear_programming/test_lp_solver.py +++ b/python/cuopt/cuopt/tests/linear_programming/test_lp_solver.py @@ -18,6 +18,7 @@ CUOPT_ITERATION_LIMIT, CUOPT_METHOD, CUOPT_MIP_HEURISTICS_ONLY, + CUOPT_PDLP_PRECISION, CUOPT_PDLP_SOLVER_MODE, CUOPT_PRIMAL_INFEASIBLE_TOLERANCE, CUOPT_RELATIVE_DUAL_TOLERANCE, @@ -604,10 +605,10 @@ def test_barrier(): A_offsets = np.array([0, 2, 4]) data_model_obj.set_csr_constraint_matrix(A_values, A_indices, A_offsets) - b = np.array([200, 160]) + b = np.array([200.0, 160.0]) data_model_obj.set_constraint_bounds(b) - c = np.array([5, 20]) + c = np.array([5.0, 20.0]) data_model_obj.set_objective_coefficients(c) row_types = np.array(["L", "L"]) @@ -722,3 +723,43 @@ def test_write_files(): assert float(line.split()[-1]) == pytest.approx(80) os.remove("afiro.sol") + + +def test_pdlp_precision_single(): + file_path = ( + RAPIDS_DATASET_ROOT_DIR + "/linear_programming/afiro_original.mps" + ) + data_model_obj = cuopt_mps_parser.ParseMps(file_path) + + settings = solver_settings.SolverSettings() + settings.set_parameter(CUOPT_METHOD, SolverMethod.PDLP) + settings.set_parameter(CUOPT_PDLP_PRECISION, 0) # Single + settings.set_optimality_tolerance(1e-4) + + solution = solver.Solve(data_model_obj, settings) + + assert solution.get_termination_status() == LPTerminationStatus.Optimal + assert solution.get_primal_objective() == pytest.approx( + -464.7531, rel=1e-1 + ) + assert solution.get_solved_by_pdlp() + + +def test_pdlp_precision_single_crossover(): + file_path = ( + RAPIDS_DATASET_ROOT_DIR + "/linear_programming/afiro_original.mps" + ) + data_model_obj = cuopt_mps_parser.ParseMps(file_path) + + settings = solver_settings.SolverSettings() + settings.set_parameter(CUOPT_METHOD, SolverMethod.PDLP) + settings.set_parameter(CUOPT_PDLP_PRECISION, 1) # Single + settings.set_parameter("crossover", True) + settings.set_optimality_tolerance(1e-4) + + solution = solver.Solve(data_model_obj, settings) + + assert solution.get_termination_status() == LPTerminationStatus.Optimal + assert solution.get_primal_objective() == pytest.approx( + -464.7531, rel=1e-1 + )