diff --git a/cpp/libmps_parser/CMakeLists.txt b/cpp/libmps_parser/CMakeLists.txt index 4fe4971576..427d4ac17b 100644 --- a/cpp/libmps_parser/CMakeLists.txt +++ b/cpp/libmps_parser/CMakeLists.txt @@ -111,6 +111,7 @@ add_library(cuopt::mps_parser ALIAS mps_parser) target_include_directories(mps_parser PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/../thirdparty" + "${CMAKE_CURRENT_SOURCE_DIR}/../src" "${CMAKE_CURRENT_SOURCE_DIR}/src" PUBLIC "$" diff --git a/cpp/libmps_parser/include/mps_parser/data_model_view.hpp b/cpp/libmps_parser/include/mps_parser/data_model_view.hpp index eb34682ce2..c2a8f84980 100644 --- a/cpp/libmps_parser/include/mps_parser/data_model_view.hpp +++ b/cpp/libmps_parser/include/mps_parser/data_model_view.hpp @@ -1,6 +1,6 @@ /* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2023-2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2023-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ /* clang-format on */ @@ -249,13 +249,15 @@ class data_model_view_t { * @param[in] Q_offsets Device memory pointer to offsets of the CSR representation of the * quadratic objective matrix * @param size_offsets Size of the Q_offsets array + * @param is_symmetrized Whether the quadratic objective matrix is a symmetrized matrix */ void set_quadratic_objective_matrix(const f_t* Q_values, i_t size_values, const i_t* Q_indices, i_t size_indices, const i_t* Q_offsets, - i_t size_offsets); + i_t size_offsets, + const bool is_symmetrized = false); /** * @brief Get the sense value (false:minimize, true:maximize) @@ -406,6 +408,13 @@ class data_model_view_t { */ bool has_quadratic_objective() const noexcept; + /** + * @brief Check if the quadratic objective matrix is a symmetrized matrix + * + * @return bool + */ + bool is_Q_symmetrized() const noexcept; + private: bool maximize_{false}; span A_; @@ -434,7 +443,7 @@ class data_model_view_t { span Q_objective_; span Q_objective_indices_; span Q_objective_offsets_; - + bool is_Q_symmetrized_{false}; }; // class data_model_view_t } // namespace cuopt::mps_parser diff --git a/cpp/libmps_parser/include/mps_parser/mps_writer.hpp b/cpp/libmps_parser/include/mps_parser/mps_writer.hpp index 0d16c6add0..30f2fdf942 100644 --- a/cpp/libmps_parser/include/mps_parser/mps_writer.hpp +++ b/cpp/libmps_parser/include/mps_parser/mps_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 */ @@ -8,9 +8,11 @@ #pragma once #include +#include #include #include +#include #include #include #include @@ -31,10 +33,16 @@ class mps_writer_t { * @brief Ctor. Takes a data model view as input and writes it out as a MPS formatted file * * @param[in] problem Data model view to write - * @param[in] file Path to the MPS file to write */ mps_writer_t(const data_model_view_t& problem); + /** + * @brief Ctor. Takes a data model as input and writes it out as a MPS formatted file + * + * @param[in] problem Data model to write + */ + mps_writer_t(const mps_data_model_t& problem); + /** * @brief Writes the problem to an MPS formatted file * @@ -43,7 +51,13 @@ class mps_writer_t { void write(const std::string& mps_file_path); private: + // Owned view (created when constructing from mps_data_model_t) + std::unique_ptr> owned_view_; + // Reference to the view (either external or owned) const data_model_view_t& problem_; + + // Helper to create view from data model + static data_model_view_t create_view(const mps_data_model_t& model); }; // class mps_writer_t } // namespace cuopt::mps_parser diff --git a/cpp/libmps_parser/src/data_model_view.cpp b/cpp/libmps_parser/src/data_model_view.cpp index 7db2b390c9..62b441aa60 100644 --- a/cpp/libmps_parser/src/data_model_view.cpp +++ b/cpp/libmps_parser/src/data_model_view.cpp @@ -1,6 +1,6 @@ /* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2023-2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2023-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ /* clang-format on */ @@ -148,7 +148,8 @@ void data_model_view_t::set_quadratic_objective_matrix(const f_t* Q_va const i_t* Q_indices, i_t size_indices, const i_t* Q_offsets, - i_t size_offsets) + i_t size_offsets, + const bool is_symmetrized) { if (size_values != 0) { mps_parser_expects( @@ -167,6 +168,8 @@ void data_model_view_t::set_quadratic_objective_matrix(const f_t* Q_va mps_parser_expects( size_offsets > 0, error_type_t::ValidationError, "size_offsets cannot be empty"); Q_objective_offsets_ = span(Q_offsets, size_offsets); + + is_Q_symmetrized_ = is_symmetrized; } template @@ -346,6 +349,12 @@ bool data_model_view_t::has_quadratic_objective() const noexcept return Q_objective_.size() > 0; } +template +bool data_model_view_t::is_Q_symmetrized() const noexcept +{ + return is_Q_symmetrized_; +} + // NOTE: Explicitly instantiate all types here in order to avoid linker error template class data_model_view_t; diff --git a/cpp/libmps_parser/src/mps_writer.cpp b/cpp/libmps_parser/src/mps_writer.cpp index 4c562ec24e..3a0997774b 100644 --- a/cpp/libmps_parser/src/mps_writer.cpp +++ b/cpp/libmps_parser/src/mps_writer.cpp @@ -8,7 +8,9 @@ #include #include +#include #include +#include #include #include @@ -16,6 +18,7 @@ #include #include #include +#include namespace cuopt::mps_parser { @@ -24,6 +27,92 @@ mps_writer_t::mps_writer_t(const data_model_view_t& problem) { } +template +data_model_view_t mps_writer_t::create_view( + const mps_data_model_t& model) +{ + data_model_view_t view; + + // Set basic data + view.set_maximize(model.get_sense()); + + // Constraint matrix + const auto& A_values = model.get_constraint_matrix_values(); + const auto& A_indices = model.get_constraint_matrix_indices(); + const auto& A_offsets = model.get_constraint_matrix_offsets(); + if (!A_values.empty()) { + view.set_csr_constraint_matrix(A_values.data(), + static_cast(A_values.size()), + A_indices.data(), + static_cast(A_indices.size()), + A_offsets.data(), + static_cast(A_offsets.size())); + } + + // Constraint bounds + const auto& b = model.get_constraint_bounds(); + if (!b.empty()) { view.set_constraint_bounds(b.data(), static_cast(b.size())); } + + // Objective coefficients + const auto& c = model.get_objective_coefficients(); + if (!c.empty()) { view.set_objective_coefficients(c.data(), static_cast(c.size())); } + + view.set_objective_scaling_factor(model.get_objective_scaling_factor()); + view.set_objective_offset(model.get_objective_offset()); + + // Variable bounds + const auto& lb = model.get_variable_lower_bounds(); + const auto& ub = model.get_variable_upper_bounds(); + if (!lb.empty()) { view.set_variable_lower_bounds(lb.data(), static_cast(lb.size())); } + if (!ub.empty()) { view.set_variable_upper_bounds(ub.data(), static_cast(ub.size())); } + + // Variable types + const auto& var_types = model.get_variable_types(); + if (!var_types.empty()) { + view.set_variable_types(var_types.data(), static_cast(var_types.size())); + } + + // Row types + const auto& row_types = model.get_row_types(); + if (!row_types.empty()) { + view.set_row_types(row_types.data(), static_cast(row_types.size())); + } + + // Constraint bounds (lower/upper) + const auto& cl = model.get_constraint_lower_bounds(); + const auto& cu = model.get_constraint_upper_bounds(); + if (!cl.empty()) { view.set_constraint_lower_bounds(cl.data(), static_cast(cl.size())); } + if (!cu.empty()) { view.set_constraint_upper_bounds(cu.data(), static_cast(cu.size())); } + + // Names + view.set_problem_name(model.get_problem_name()); + view.set_objective_name(model.get_objective_name()); + view.set_variable_names(model.get_variable_names()); + view.set_row_names(model.get_row_names()); + + // Quadratic objective + const auto& Q_values = model.get_quadratic_objective_values(); + const auto& Q_indices = model.get_quadratic_objective_indices(); + const auto& Q_offsets = model.get_quadratic_objective_offsets(); + if (!Q_values.empty()) { + view.set_quadratic_objective_matrix(Q_values.data(), + static_cast(Q_values.size()), + Q_indices.data(), + static_cast(Q_indices.size()), + Q_offsets.data(), + static_cast(Q_offsets.size())); + } + + return view; +} + +template +mps_writer_t::mps_writer_t(const mps_data_model_t& problem) + : owned_view_(std::make_unique>(create_view(problem))), + problem_(*owned_view_) +{ +} + template void mps_writer_t::write(const std::string& mps_file_path) { @@ -280,6 +369,64 @@ void mps_writer_t::write(const std::string& mps_file_path) } } + // QUADOBJ section for quadratic objective terms (if present) + // MPS format: QUADOBJ stores upper triangular elements (row <= col) + // MPS uses (1/2) x^T H x, cuOpt uses x^T Q x + // For equivalence: H[i,j] = Q[i,j] + Q[j,i] (works for both diagonal and off-diagonal) + // We symmetrize Q first (H = Q + Q^T), then extract upper triangular + if (problem_.has_quadratic_objective()) { + auto Q_values_span = problem_.get_quadratic_objective_values(); + auto Q_indices_span = problem_.get_quadratic_objective_indices(); + auto Q_offsets_span = problem_.get_quadratic_objective_offsets(); + + // Copy span data to local vectors for indexed access + std::vector Q_values(Q_values_span.data(), Q_values_span.data() + Q_values_span.size()); + std::vector Q_indices(Q_indices_span.data(), + Q_indices_span.data() + Q_indices_span.size()); + std::vector Q_offsets(Q_offsets_span.data(), + Q_offsets_span.data() + Q_offsets_span.size()); + + if (Q_values.size() > 0) { + // Symmetrize Q: compute H = Q + Q^T + std::vector H_values; + std::vector H_indices; + std::vector H_offsets; + + if (problem_.is_Q_symmetrized()) { + H_values = std::move(Q_values); + H_indices = std::move(Q_indices); + H_offsets = std::move(Q_offsets); + } else { + cuopt::symmetrize_csr( + Q_values, Q_indices, Q_offsets, H_values, H_indices, H_offsets); + } + + i_t n_rows = static_cast(H_offsets.size()) - 1; + + mps_file << "QUADOBJ\n"; + + // Write upper triangular entries from symmetric H + for (i_t i = 0; i < n_rows; ++i) { + std::string row_name = static_cast(i) < problem_.get_variable_names().size() + ? problem_.get_variable_names()[i] + : "C" + std::to_string(i); + + for (i_t p = H_offsets[i]; p < H_offsets[i + 1]; ++p) { + i_t j = H_indices[p]; + f_t v = H_values[p]; + + // Only write upper triangular (i <= j) + if (i <= j && v != f_t(0)) { + std::string col_name = static_cast(j) < problem_.get_variable_names().size() + ? problem_.get_variable_names()[j] + : "C" + std::to_string(j); + mps_file << " " << row_name << " " << col_name << " " << v << "\n"; + } + } + } + } + } + mps_file << "ENDATA\n"; mps_file.close(); } diff --git a/cpp/libmps_parser/tests/mps_parser_test.cpp b/cpp/libmps_parser/tests/mps_parser_test.cpp index 510719d02a..f915fb2df5 100644 --- a/cpp/libmps_parser/tests/mps_parser_test.cpp +++ b/cpp/libmps_parser/tests/mps_parser_test.cpp @@ -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 */ @@ -8,10 +8,12 @@ #include #include +#include #include #include +#include #include #include #include @@ -892,4 +894,237 @@ TEST(qps_parser, test_qps_files) } } +// ================================================================================================ +// MPS Round-Trip Tests (Read -> Write -> Read -> Compare) +// ================================================================================================ + +// Helper function to compare two data models +template +void compare_data_models(const mps_data_model_t& original, + const mps_data_model_t& reloaded, + f_t tol = 1e-9) +{ + // Compare basic dimensions + EXPECT_EQ(original.get_n_variables(), reloaded.get_n_variables()); + EXPECT_EQ(original.get_n_constraints(), reloaded.get_n_constraints()); + + // Compare objective coefficients + auto orig_c = original.get_objective_coefficients(); + auto reload_c = reloaded.get_objective_coefficients(); + ASSERT_EQ(orig_c.size(), reload_c.size()); + for (size_t i = 0; i < orig_c.size(); ++i) { + EXPECT_NEAR(orig_c[i], reload_c[i], tol) << "Objective coefficient mismatch at index " << i; + } + + // Compare constraint matrix values + auto orig_A = original.get_constraint_matrix_values(); + auto reload_A = reloaded.get_constraint_matrix_values(); + ASSERT_EQ(orig_A.size(), reload_A.size()); + for (size_t i = 0; i < orig_A.size(); ++i) { + EXPECT_NEAR(orig_A[i], reload_A[i], tol) << "Constraint matrix value mismatch at index " << i; + } + + // Compare constraint matrix indices + auto orig_A_idx = original.get_constraint_matrix_indices(); + auto reload_A_idx = reloaded.get_constraint_matrix_indices(); + ASSERT_EQ(orig_A_idx.size(), reload_A_idx.size()); + for (size_t i = 0; i < orig_A_idx.size(); ++i) { + EXPECT_EQ(orig_A_idx[i], reload_A_idx[i]) << "Constraint matrix index mismatch at index " << i; + } + + // Compare constraint matrix offsets + auto orig_A_off = original.get_constraint_matrix_offsets(); + auto reload_A_off = reloaded.get_constraint_matrix_offsets(); + ASSERT_EQ(orig_A_off.size(), reload_A_off.size()); + for (size_t i = 0; i < orig_A_off.size(); ++i) { + EXPECT_EQ(orig_A_off[i], reload_A_off[i]) << "Constraint matrix offset mismatch at index " << i; + } + + // Compare variable bounds + auto orig_lb = original.get_variable_lower_bounds(); + auto reload_lb = reloaded.get_variable_lower_bounds(); + ASSERT_EQ(orig_lb.size(), reload_lb.size()); + for (size_t i = 0; i < orig_lb.size(); ++i) { + if (std::isinf(orig_lb[i]) && std::isinf(reload_lb[i])) { + EXPECT_EQ(std::signbit(orig_lb[i]), std::signbit(reload_lb[i])) + << "Variable lower bound infinity sign mismatch at index " << i; + } else { + EXPECT_NEAR(orig_lb[i], reload_lb[i], tol) << "Variable lower bound mismatch at index " << i; + } + } + + auto orig_ub = original.get_variable_upper_bounds(); + auto reload_ub = reloaded.get_variable_upper_bounds(); + ASSERT_EQ(orig_ub.size(), reload_ub.size()); + for (size_t i = 0; i < orig_ub.size(); ++i) { + if (std::isinf(orig_ub[i]) && std::isinf(reload_ub[i])) { + EXPECT_EQ(std::signbit(orig_ub[i]), std::signbit(reload_ub[i])) + << "Variable upper bound infinity sign mismatch at index " << i; + } else { + EXPECT_NEAR(orig_ub[i], reload_ub[i], tol) << "Variable upper bound mismatch at index " << i; + } + } + + // Compare constraint bounds + auto orig_cl = original.get_constraint_lower_bounds(); + auto reload_cl = reloaded.get_constraint_lower_bounds(); + ASSERT_EQ(orig_cl.size(), reload_cl.size()); + for (size_t i = 0; i < orig_cl.size(); ++i) { + if (std::isinf(orig_cl[i]) && std::isinf(reload_cl[i])) { + EXPECT_EQ(std::signbit(orig_cl[i]), std::signbit(reload_cl[i])) + << "Constraint lower bound infinity sign mismatch at index " << i; + } else { + EXPECT_NEAR(orig_cl[i], reload_cl[i], tol) + << "Constraint lower bound mismatch at index " << i; + } + } + + auto orig_cu = original.get_constraint_upper_bounds(); + auto reload_cu = reloaded.get_constraint_upper_bounds(); + ASSERT_EQ(orig_cu.size(), reload_cu.size()); + for (size_t i = 0; i < orig_cu.size(); ++i) { + if (std::isinf(orig_cu[i]) && std::isinf(reload_cu[i])) { + EXPECT_EQ(std::signbit(orig_cu[i]), std::signbit(reload_cu[i])) + << "Constraint upper bound infinity sign mismatch at index " << i; + } else { + EXPECT_NEAR(orig_cu[i], reload_cu[i], tol) + << "Constraint upper bound mismatch at index " << i; + } + } + + // Compare quadratic objective if present + EXPECT_EQ(original.has_quadratic_objective(), reloaded.has_quadratic_objective()); + if (original.has_quadratic_objective() && reloaded.has_quadratic_objective()) { + auto orig_Q = original.get_quadratic_objective_values(); + auto orig_Q_idx = original.get_quadratic_objective_indices(); + auto orig_Q_off = original.get_quadratic_objective_offsets(); + auto reload_Q = reloaded.get_quadratic_objective_values(); + auto reload_Q_idx = reloaded.get_quadratic_objective_indices(); + auto reload_Q_off = reloaded.get_quadratic_objective_offsets(); + + // Compare Q matrix structure and values + ASSERT_EQ(orig_Q.size(), reload_Q.size()) << "Q values size mismatch"; + ASSERT_EQ(orig_Q_idx.size(), reload_Q_idx.size()) << "Q indices size mismatch"; + ASSERT_EQ(orig_Q_off.size(), reload_Q_off.size()) << "Q offsets size mismatch"; + + for (size_t i = 0; i < orig_Q.size(); ++i) { + EXPECT_NEAR(orig_Q[i], reload_Q[i], tol) << "Q value mismatch at index " << i; + } + for (size_t i = 0; i < orig_Q_idx.size(); ++i) { + EXPECT_EQ(orig_Q_idx[i], reload_Q_idx[i]) << "Q index mismatch at index " << i; + } + for (size_t i = 0; i < orig_Q_off.size(); ++i) { + EXPECT_EQ(orig_Q_off[i], reload_Q_off[i]) << "Q offset mismatch at index " << i; + } + } +} + +TEST(mps_roundtrip, linear_programming_basic) +{ + std::string input_file = + cuopt::test::get_rapids_dataset_root_dir() + "/linear_programming/good-mps-1.mps"; + std::string temp_file = "/tmp/mps_roundtrip_lp_test.mps"; + + // Read original + auto original = parse_mps(input_file, true); + + // Write to temp file + mps_writer_t writer(original); + writer.write(temp_file); + + // Read back + auto reloaded = parse_mps(temp_file, false); + + // Compare + compare_data_models(original, reloaded); + + // Cleanup + std::filesystem::remove(temp_file); +} + +TEST(mps_roundtrip, linear_programming_with_bounds) +{ + if (!file_exists("linear_programming/lp_model_with_var_bounds.mps")) { + GTEST_SKIP() << "Test file not found"; + } + + std::string input_file = + cuopt::test::get_rapids_dataset_root_dir() + "/linear_programming/lp_model_with_var_bounds.mps"; + std::string temp_file = "/tmp/mps_roundtrip_lp_bounds_test.mps"; + + // Read original + auto original = parse_mps(input_file, false); + + // Write to temp file + mps_writer_t writer(original); + writer.write(temp_file); + + // Read back + auto reloaded = parse_mps(temp_file, false); + + // Compare + compare_data_models(original, reloaded); + + // Cleanup + std::filesystem::remove(temp_file); +} + +TEST(mps_roundtrip, quadratic_programming_qp_test_1) +{ + if (!file_exists("quadratic_programming/QP_Test_1.qps")) { + GTEST_SKIP() << "Test file not found"; + } + + std::string input_file = + cuopt::test::get_rapids_dataset_root_dir() + "/quadratic_programming/QP_Test_1.qps"; + std::string temp_file = "/tmp/mps_roundtrip_qp_test_1.mps"; + + // Read original + auto original = parse_mps(input_file, false); + ASSERT_TRUE(original.has_quadratic_objective()) << "Original should have quadratic objective"; + + // Write to temp file + mps_writer_t writer(original); + writer.write(temp_file); + + // Read back + auto reloaded = parse_mps(temp_file, false); + ASSERT_TRUE(reloaded.has_quadratic_objective()) << "Reloaded should have quadratic objective"; + + // Compare + compare_data_models(original, reloaded); + + // Cleanup + std::filesystem::remove(temp_file); +} + +TEST(mps_roundtrip, quadratic_programming_qp_test_2) +{ + if (!file_exists("quadratic_programming/QP_Test_2.qps")) { + GTEST_SKIP() << "Test file not found"; + } + + std::string input_file = + cuopt::test::get_rapids_dataset_root_dir() + "/quadratic_programming/QP_Test_2.qps"; + std::string temp_file = "/tmp/mps_roundtrip_qp_test_2.mps"; + + // Read original + auto original = parse_mps(input_file, false); + ASSERT_TRUE(original.has_quadratic_objective()) << "Original should have quadratic objective"; + + // Write to temp file + mps_writer_t writer(original); + writer.write(temp_file); + + // Read back + auto reloaded = parse_mps(temp_file, false); + ASSERT_TRUE(reloaded.has_quadratic_objective()) << "Reloaded should have quadratic objective"; + + // Compare + compare_data_models(original, reloaded); + + // Cleanup + std::filesystem::remove(temp_file); +} + } // namespace cuopt::mps_parser diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index b9ee419517..d2a68d96de 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -891,14 +891,19 @@ i_t presolve(const lp_problem_t& original, } } + std::vector kahan_compensation(problem.num_rows, 0.0); for (i_t j = 0; j < problem.num_cols; j++) { if (lower_bounds_removed[j]) { i_t col_start = problem.A.col_start[j]; i_t col_end = problem.A.col_start[j + 1]; for (i_t p = col_start; p < col_end; p++) { - i_t i = problem.A.i[p]; - f_t aij = problem.A.x[p]; - problem.rhs[i] -= aij * problem.lower[j]; + i_t i = problem.A.i[p]; + f_t aij = problem.A.x[p]; + f_t val = -aij * problem.lower[j]; + f_t y = val - kahan_compensation[i]; + f_t t = problem.rhs[i] + y; + kahan_compensation[i] = (t - problem.rhs[i]) - y; + problem.rhs[i] = t; } problem.obj_constant += old_objective[j] * problem.lower[j]; problem.upper[j] -= problem.lower[j]; @@ -1496,13 +1501,15 @@ void uncrush_solution(const presolve_info_t& presolve_info, } if (presolve_info.removed_lower_bounds.size() > 0) { - settings.log.printf("Post-solve: Handling removed lower bounds %d\n", - presolve_info.removed_lower_bounds.size()); + i_t num_lower_bounds = 0; + // We removed some lower bounds so we need to map the crushed solution back to the original // variables for (i_t j = 0; j < input_x.size(); j++) { + if (presolve_info.removed_lower_bounds[j] != 0.0) { num_lower_bounds++; } input_x[j] += presolve_info.removed_lower_bounds[j]; } + settings.log.printf("Post-solve: Handling removed lower bounds %d\n", num_lower_bounds); } assert(uncrushed_x.size() == input_x.size()); assert(uncrushed_y.size() == input_y.size()); diff --git a/cpp/src/pdlp/cpu_optimization_problem.cpp b/cpp/src/pdlp/cpu_optimization_problem.cpp index 39b23ad7c2..406b0b6541 100644 --- a/cpp/src/pdlp/cpu_optimization_problem.cpp +++ b/cpp/src/pdlp/cpu_optimization_problem.cpp @@ -122,6 +122,8 @@ void cpu_optimization_problem_t::set_quadratic_objective_matrix( cuopt_expects(Q_offsets != nullptr, error_type_t::ValidationError, "Q_offsets cannot be null"); } + // Store raw CSR Q; symmetrization (H = Q + Q^T) is applied in optimization_problem_t setter on + // GPU. Q_values_.resize(size_values); Q_indices_.resize(size_indices); Q_offsets_.resize(size_offsets); @@ -609,7 +611,7 @@ cpu_optimization_problem_t::to_optimization_problem(raft::handle_t con // Set objective coefficients if (!c_.empty()) { gpu_problem->set_objective_coefficients(c_.data(), c_.size()); } - // Set quadratic objective if present + // Set quadratic objective if present (GPU setter symmetrizes once: H = Q + Q^T) if (!Q_values_.empty()) { gpu_problem->set_quadratic_objective_matrix(Q_values_.data(), Q_values_.size(), @@ -726,6 +728,18 @@ void cpu_optimization_problem_t::write_to_mps(const std::string& mps_f data_model_view.set_variable_types(var_types_char.data(), var_types_char.size()); } + if (!Q_values_.empty()) { + // cpu optimization problem stores the raw CSR matrix, so we need to symmetrize it + const bool is_symmetrized = false; + data_model_view.set_quadratic_objective_matrix(Q_values_.data(), + static_cast(Q_values_.size()), + Q_indices_.data(), + static_cast(Q_indices_.size()), + Q_offsets_.data(), + static_cast(Q_offsets_.size()), + false); + } + cuopt::mps_parser::write_mps(data_model_view, mps_file_path); } diff --git a/cpp/src/pdlp/optimization_problem.cu b/cpp/src/pdlp/optimization_problem.cu index f1b7712445..87ff9dab08 100644 --- a/cpp/src/pdlp/optimization_problem.cu +++ b/cpp/src/pdlp/optimization_problem.cu @@ -15,6 +15,7 @@ #include #include #include +#include #include #include @@ -188,87 +189,11 @@ void optimization_problem_t::set_quadratic_objective_matrix( size_offsets >= 1, error_type_t::ValidationError, "Q_offsets must have at least 1 element"); cuopt_expects(Q_offsets != nullptr, error_type_t::ValidationError, "Q_offsets cannot be null"); - // Replace Q with Q + Q^T - i_t qn = size_offsets - 1; // Number of variables - i_t q_nnz = size_indices; - // Construct H = Q + Q^T in triplet form first - std::vector H_i; - std::vector H_j; - std::vector H_x; - - H_i.reserve(2 * q_nnz); - H_j.reserve(2 * q_nnz); - H_x.reserve(2 * q_nnz); - - for (i_t i = 0; i < qn; ++i) { - i_t row_start = Q_offsets[i]; - i_t row_end = Q_offsets[i + 1]; - for (i_t p = row_start; p < row_end; ++p) { - i_t j = Q_indices[p]; - f_t x = Q_values[p]; - // Add H(i,j) - H_i.push_back(i); - H_j.push_back(j); - if (i == j) { H_x.push_back(2 * x); } - if (i != j) { - H_x.push_back(x); - // Add H(j,i) - H_i.push_back(j); - H_j.push_back(i); - H_x.push_back(x); - } - } - } - // Convert H to CSR format - // Get row counts - i_t H_nz = H_x.size(); - std::vector H_row_counts(qn, 0); - for (i_t k = 0; k < H_nz; ++k) { - H_row_counts[H_i[k]]++; - } - std::vector H_cumulative_counts(qn + 1, 0); - for (i_t k = 0; k < qn; ++k) { - H_cumulative_counts[k + 1] = H_cumulative_counts[k] + H_row_counts[k]; - } - std::vector H_row_starts = H_cumulative_counts; - std::vector H_indices(H_nz); - std::vector H_values(H_nz); - for (i_t k = 0; k < H_nz; ++k) { - i_t p = H_cumulative_counts[H_i[k]]++; - H_indices[p] = H_j[k]; - H_values[p] = H_x[k]; - } - - // H_row_starts, H_indices, H_values are the CSR representation of H - // But this contains duplicate entries - - std::vector workspace(qn, -1); - Q_offsets_.resize(qn + 1); - std::fill(Q_offsets_.begin(), Q_offsets_.end(), 0); - Q_indices_.resize(H_nz); - Q_values_.resize(H_nz); - i_t nz = 0; - for (i_t i = 0; i < qn; ++i) { - i_t q = nz; // row i will start at q - const i_t row_start = H_row_starts[i]; - const i_t row_end = H_row_starts[i + 1]; - for (i_t p = row_start; p < row_end; ++p) { - i_t j = H_indices[p]; - if (workspace[j] >= q) { - Q_values_[workspace[j]] += H_values[p]; // H(i,j) is duplicate - } else { - workspace[j] = nz; // record where column j occurs - Q_indices_[nz] = j; // keep H(i,j) - Q_values_[nz] = H_values[p]; - nz++; - } - } - Q_offsets_[i] = q; // record start of row i - } - - Q_offsets_[qn] = nz; // finalize Q - Q_indices_.resize(nz); - Q_values_.resize(nz); + // Symmetrize Q to H = Q + Q^T (same algorithm as cuopt::symmetrize_csr in + // sparse_matrix_helpers.hpp) + const i_t qn = size_offsets - 1; + cuopt::symmetrize_csr( + Q_values, Q_indices, Q_offsets, qn, Q_values_, Q_indices_, Q_offsets_); // FIX ME:: check for positive semi definite matrix } @@ -883,6 +808,18 @@ void optimization_problem_t::write_to_mps(const std::string& mps_file_ data_model_view.set_variable_types(variable_types.data(), variable_types.size()); } + if (!Q_values_.empty()) { + // Only the symmetrized matrix is stored in the optimization_problem_t + const bool is_symmetrized = true; + data_model_view.set_quadratic_objective_matrix(Q_values_.data(), + static_cast(Q_values_.size()), + Q_indices_.data(), + static_cast(Q_indices_.size()), + Q_offsets_.data(), + static_cast(Q_offsets_.size()), + is_symmetrized); + } + cuopt::mps_parser::write_mps(data_model_view, mps_file_path); } diff --git a/cpp/src/utilities/sparse_matrix_helpers.hpp b/cpp/src/utilities/sparse_matrix_helpers.hpp new file mode 100644 index 0000000000..2bf2d96d62 --- /dev/null +++ b/cpp/src/utilities/sparse_matrix_helpers.hpp @@ -0,0 +1,160 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2022-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#pragma once + +#include +#include + +namespace cuopt { + +/** + * @brief Symmetrize a CSR matrix by computing A + A^T + * + * Given a CSR matrix A, computes the symmetric matrix H = A + A^T. + * Diagonal entries are doubled (A[i,i] + A[i,i] = 2*A[i,i]). + * Off-diagonal entries are summed (H[i,j] = A[i,j] + A[j,i]). + * + * This is useful for quadratic programming where the objective is + * (1/2) x^T Q x, and Q needs to be symmetric. + * + * @tparam i_t Integer type for indices + * @tparam f_t Floating point type for values + * + * @param[in] in_values CSR values array + * @param[in] in_indices CSR column indices array + * @param[in] in_offsets CSR row offsets array (size = n_rows + 1) + * @param[in] n_rows Number of rows (and columns, assuming square matrix) + * @param[out] out_values Output CSR values + * @param[out] out_indices Output CSR column indices + * @param[out] out_offsets Output CSR row offsets + */ +template +void symmetrize_csr(const f_t* in_values, + const i_t* in_indices, + const i_t* in_offsets, + i_t n_rows, + std::vector& out_values, + std::vector& out_indices, + std::vector& out_offsets) +{ + static_assert(std::is_integral_v && std::is_signed_v, + "symmetrize_csr: i_t must be a signed integral type (workspace uses -1 sentinel)."); + + // Optimized 3-pass algorithm (no COO intermediate) + // Memory: ~3× nnz temporary storage before deduplication + + // Pass 1: Count entries per row in A + A^T + std::vector row_counts(n_rows, 0); + for (i_t i = 0; i < n_rows; ++i) { + for (i_t p = in_offsets[i]; p < in_offsets[i + 1]; ++p) { + i_t j = in_indices[p]; + row_counts[i]++; + if (i != j) { row_counts[j]++; } + } + } + + // Build temporary offsets via prefix sum + std::vector temp_offsets(n_rows + 1); + temp_offsets[0] = 0; + for (i_t i = 0; i < n_rows; ++i) { + temp_offsets[i + 1] = temp_offsets[i] + row_counts[i]; + } + + i_t total_entries = temp_offsets[n_rows]; + std::vector temp_indices(total_entries); + std::vector temp_values(total_entries); + + // Pass 2: Fill entries directly + std::vector row_pos = temp_offsets; // Copy for tracking insertion positions + + for (i_t i = 0; i < n_rows; ++i) { + for (i_t p = in_offsets[i]; p < in_offsets[i + 1]; ++p) { + i_t j = in_indices[p]; + f_t x = in_values[p]; + + // Add entry (i, j) with value 2x for diagonal, x for off-diagonal + temp_indices[row_pos[i]] = j; + temp_values[row_pos[i]] = (i == j) ? (2 * x) : x; + row_pos[i]++; + + // Add transpose entry (j, i) if off-diagonal + if (i != j) { + temp_indices[row_pos[j]] = i; + temp_values[row_pos[j]] = x; + row_pos[j]++; + } + } + } + + // Pass 3: Deduplicate and build final CSR + std::vector workspace(n_rows, -1); + out_offsets.resize(n_rows + 1); + out_indices.resize(total_entries); + out_values.resize(total_entries); + + i_t nz = 0; + for (i_t i = 0; i < n_rows; ++i) { + i_t row_start_out = nz; + out_offsets[i] = row_start_out; + + for (i_t p = temp_offsets[i]; p < temp_offsets[i + 1]; ++p) { + i_t j = temp_indices[p]; + f_t x = temp_values[p]; + + if (workspace[j] >= row_start_out) { + out_values[workspace[j]] += x; + } else { + workspace[j] = nz; + out_indices[nz] = j; + out_values[nz] = x; + nz++; + } + } + } + + out_offsets[n_rows] = nz; + out_indices.resize(nz); + out_values.resize(nz); +} + +/** + * @brief Symmetrize a CSR matrix in-place using std::vector + * + * Convenience overload that takes and returns std::vectors. + * + * @tparam i_t Integer type for indices + * @tparam f_t Floating point type for values + * + * @param[in] in_values Input CSR values + * @param[in] in_indices Input CSR column indices + * @param[in] in_offsets Input CSR row offsets + * @param[out] out_values Output CSR values (can be same as in_values) + * @param[out] out_indices Output CSR column indices (can be same as in_indices) + * @param[out] out_offsets Output CSR row offsets (can be same as in_offsets) + */ +template +void symmetrize_csr(const std::vector& in_values, + const std::vector& in_indices, + const std::vector& in_offsets, + std::vector& out_values, + std::vector& out_indices, + std::vector& out_offsets) +{ + if (in_offsets.size() <= 1) { return; } + + i_t n_rows = static_cast(in_offsets.size()) - 1; + symmetrize_csr(in_values.data(), + in_indices.data(), + in_offsets.data(), + n_rows, + out_values, + out_indices, + out_offsets); +} + +} // namespace cuopt diff --git a/cpp/tests/qp/CMakeLists.txt b/cpp/tests/qp/CMakeLists.txt index eefa843c7c..9c3ee8923c 100644 --- a/cpp/tests/qp/CMakeLists.txt +++ b/cpp/tests/qp/CMakeLists.txt @@ -6,4 +6,5 @@ ConfigureTest(QP_UNIT_TEST ${CMAKE_CURRENT_SOURCE_DIR}/unit_tests/no_constraints.cu ${CMAKE_CURRENT_SOURCE_DIR}/unit_tests/two_variable_test.cu + ${CMAKE_CURRENT_SOURCE_DIR}/unit_tests/mps_writer_test.cpp ) diff --git a/cpp/tests/qp/unit_tests/mps_writer_test.cpp b/cpp/tests/qp/unit_tests/mps_writer_test.cpp new file mode 100644 index 0000000000..62951ccb9b --- /dev/null +++ b/cpp/tests/qp/unit_tests/mps_writer_test.cpp @@ -0,0 +1,260 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2022-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +/** + * Builds small unconstrained QPs as optimization_problem_t, writes MPS via write_to_mps, and checks + * QUADOBJ coefficients against the symmetric Hessian H = Q + Q^T stored internally (MPS uses + * (1/2) x^T H x for the quadratic part). + */ + +#include + +#include +#include + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +struct temp_mps_file_guard_t { + explicit temp_mps_file_guard_t(std::string p) : path(std::move(p)) {} + std::string path; + ~temp_mps_file_guard_t() + { + if (!path.empty()) { + std::error_code ec; + std::filesystem::remove(path, ec); + } + } +}; + +std::string read_entire_file(std::string const& file_path) +{ + std::ifstream in(file_path); + if (!in) { return {}; } + std::ostringstream ss; + ss << in.rdbuf(); + return ss.str(); +} + +/** Text between the line after "QUADOBJ" and "ENDATA". */ +std::string extract_quadobj_body(std::string const& content) +{ + auto pos = content.find("QUADOBJ"); + if (pos == std::string::npos) { return {}; } + pos = content.find('\n', pos); + if (pos == std::string::npos) { return {}; } + ++pos; + auto const end = content.find("ENDATA", pos); + if (end == std::string::npos) { return content.substr(pos); } + return content.substr(pos, end - pos); +} + +/** One QUADOBJ line: default variable names are C0, C1, ... */ +struct quadobj_entry_t { + int row{}; + int col{}; + double value{}; +}; + +bool parse_default_c_index(std::string const& name, int& index_out) +{ + if (name.size() < 2 || name[0] != 'C') { return false; } + try { + index_out = std::stoi(name.substr(1)); + return true; + } catch (...) { + return false; + } +} + +std::vector parse_quadobj_entries(std::string const& quad_body) +{ + std::vector entries; + std::istringstream stream(quad_body); + std::string line; + while (std::getline(stream, line)) { + if (line.empty()) { continue; } + std::istringstream ls(line); + std::string row_name; + std::string col_name; + if (!(ls >> row_name >> col_name)) { continue; } + std::string tok; + std::string last; + while (ls >> tok) { + last = tok; + } + if (last.empty()) { continue; } + char* endptr = nullptr; + double const v = std::strtod(last.c_str(), &endptr); + if (endptr == last.c_str()) { continue; } + int row = 0; + int col = 0; + if (!parse_default_c_index(row_name, row) || !parse_default_c_index(col_name, col)) { + continue; + } + entries.push_back({row, col, v}); + } + return entries; +} + +void sort_quadobj_entries_lex(std::vector& entries) +{ + std::sort(entries.begin(), entries.end(), [](quadobj_entry_t const& a, quadobj_entry_t const& b) { + return std::tie(a.row, a.col) < std::tie(b.row, b.col); + }); +} + +/** Two-variable QP with no constraints (same pattern as no_constraints.cu). */ +template +void setup_two_var_unconstrained_qp(Op& op) +{ + double A_values_host[] = {}; + int A_indices_host[] = {}; + int A_offsets_host[] = {0}; + op.set_csr_constraint_matrix(A_values_host, 0, A_indices_host, 0, A_offsets_host, 1); + + double lb_host[] = {0.0, 0.0}; + double ub_host[] = {std::numeric_limits::infinity(), + std::numeric_limits::infinity()}; + op.set_variable_lower_bounds(lb_host, 2); + op.set_variable_upper_bounds(ub_host, 2); + + cuopt::linear_programming::var_t const var_types_host[] = { + cuopt::linear_programming::var_t::CONTINUOUS, cuopt::linear_programming::var_t::CONTINUOUS}; + op.set_variable_types(var_types_host, 2); + + double c_host[] = {0.0, 0.0}; + op.set_objective_coefficients(c_host, 2); +} + +/** Three-variable unconstrained QP; same structure as the 2-variable helper. */ +template +void setup_three_var_unconstrained_qp(Op& op) +{ + double A_values_host[] = {}; + int A_indices_host[] = {}; + int A_offsets_host[] = {0}; + op.set_csr_constraint_matrix(A_values_host, 0, A_indices_host, 0, A_offsets_host, 1); + + double lb_host[] = {0.0, 0.0, 0.0}; + double ub_host[] = {std::numeric_limits::infinity(), + std::numeric_limits::infinity(), + std::numeric_limits::infinity()}; + op.set_variable_lower_bounds(lb_host, 3); + op.set_variable_upper_bounds(ub_host, 3); + + cuopt::linear_programming::var_t const var_types_host[] = { + cuopt::linear_programming::var_t::CONTINUOUS, + cuopt::linear_programming::var_t::CONTINUOUS, + cuopt::linear_programming::var_t::CONTINUOUS}; + op.set_variable_types(var_types_host, 3); + + double c_host[] = {0.0, 0.0, 0.0}; + op.set_objective_coefficients(c_host, 3); +} + +} // namespace + +namespace cuopt::linear_programming { + +TEST(mps_writer_op, write_to_mps_diagonal_qp_quadobj_matches_symmetrized_hessian) +{ + raft::handle_t handle; + + // minimize x1^2 + x2^2 => Q = diag(1,1) in CSR; symmetrized H = Q + Q^T has H_ii = 2. + auto op = optimization_problem_t(&handle); + setup_two_var_unconstrained_qp(op); + + double Q_values_host[] = {1.0, 1.0}; + int Q_indices_host[] = {0, 1}; + int Q_offsets_host[] = {0, 1, 2}; + op.set_quadratic_objective_matrix(Q_values_host, 2, Q_indices_host, 2, Q_offsets_host, 3); + + std::string const path = std::string(::testing::TempDir()) + "qp_diag_write.mps"; + temp_mps_file_guard_t guard(path); + op.write_to_mps(path); + + std::string const content = read_entire_file(path); + ASSERT_FALSE(content.empty()) << "MPS file was empty or could not be read"; + + std::string const quad_body = extract_quadobj_body(content); + ASSERT_FALSE(quad_body.empty()) << "No QUADOBJ section in:\n" << content; + + auto entries = parse_quadobj_entries(quad_body); + ASSERT_EQ(entries.size(), 2u) + << "Expected two upper-triangular QUADOBJ entries for 2x2 diagonal H"; + sort_quadobj_entries_lex(entries); + + EXPECT_EQ(entries[0].row, 0); + EXPECT_EQ(entries[0].col, 0); + EXPECT_NEAR(entries[0].value, 2.0, 1e-9); + EXPECT_EQ(entries[1].row, 1); + EXPECT_EQ(entries[1].col, 1); + EXPECT_NEAR(entries[1].value, 2.0, 1e-9); +} + +TEST(mps_writer_op, write_to_mps_nonsymmetric_Q_quadobj_matches_Q_plus_Q_transpose) +{ + raft::handle_t handle; + + // Non-symmetric 3x3 Q; zeros at (0,2) and (2,0) omitted from CSR. + // Dense Q: + // [ 1 2 0 ] + // [ 3 4 5 ] + // [ 0 7 8 ] + // H = Q + Q^T has upper-triangle nonzeros 2, 5, 8, 12, 16 (H(0,2)=0 is not written to QUADOBJ). + auto op = optimization_problem_t(&handle); + setup_three_var_unconstrained_qp(op); + + double Q_values_host[] = {1.0, 2.0, 3.0, 4.0, 5.0, 7.0, 8.0}; + int Q_indices_host[] = {0, 1, 0, 1, 2, 1, 2}; + int Q_offsets_host[] = {0, 2, 5, 7}; + op.set_quadratic_objective_matrix(Q_values_host, 7, Q_indices_host, 7, Q_offsets_host, 4); + + std::string const path = std::string(::testing::TempDir()) + "qp_nonsym_sparse_write.mps"; + temp_mps_file_guard_t guard(path); + op.write_to_mps(path); + + std::string const content = read_entire_file(path); + ASSERT_FALSE(content.empty()) << "MPS file was empty or could not be read"; + + std::string const quad_body = extract_quadobj_body(content); + ASSERT_FALSE(quad_body.empty()) << "No QUADOBJ section in:\n" << content; + + auto entries = parse_quadobj_entries(quad_body); + ASSERT_EQ(entries.size(), 5u) + << "Expected five nonzero upper-triangular QUADOBJ entries (H(0,2)=0 skipped)"; + sort_quadobj_entries_lex(entries); + + std::vector const expected = { + {0, 0, 2.0}, + {0, 1, 5.0}, + {1, 1, 8.0}, + {1, 2, 12.0}, + {2, 2, 16.0}, + }; + for (size_t i = 0; i < expected.size(); ++i) { + EXPECT_EQ(entries[i].row, expected[i].row) << "entry " << i; + EXPECT_EQ(entries[i].col, expected[i].col) << "entry " << i; + EXPECT_NEAR(entries[i].value, expected[i].value, 1e-9) << "entry " << i; + } +} + +} // namespace cuopt::linear_programming