diff --git a/cpp/include/cuopt/mathematical_optimization/cpu_optimization_problem.hpp b/cpp/include/cuopt/mathematical_optimization/cpu_optimization_problem.hpp index fd4c50fa75..7277487f65 100644 --- a/cpp/include/cuopt/mathematical_optimization/cpu_optimization_problem.hpp +++ b/cpp/include/cuopt/mathematical_optimization/cpu_optimization_problem.hpp @@ -7,6 +7,7 @@ #pragma once +#include #include #include @@ -45,6 +46,11 @@ class cpu_optimization_problem_t : public optimization_problem_interface_t::quadratic_constraint_t; cpu_optimization_problem_t(); + cpu_optimization_problem_t(const cpu_optimization_problem_t& other); + cpu_optimization_problem_t(cpu_optimization_problem_t&& other) noexcept; + cpu_optimization_problem_t& operator=(const cpu_optimization_problem_t& other); + cpu_optimization_problem_t& operator=(cpu_optimization_problem_t&& other) noexcept; + ~cpu_optimization_problem_t(); // Setters void set_maximize(bool maximize) override; @@ -212,6 +218,8 @@ class cpu_optimization_problem_t : public optimization_problem_interface_t Q_values_; std::vector quadratic_constraints_{}; + // Scratch storage reused across add_quadratic_constraint calls. Not thread-safe. + std::unique_ptr> qc_coo_workspace_; std::vector variable_lower_bounds_; std::vector variable_upper_bounds_; diff --git a/cpp/include/cuopt/mathematical_optimization/io/mps_data_model.hpp b/cpp/include/cuopt/mathematical_optimization/io/mps_data_model.hpp index 4ed1d7244f..238a72ff33 100644 --- a/cpp/include/cuopt/mathematical_optimization/io/mps_data_model.hpp +++ b/cpp/include/cuopt/mathematical_optimization/io/mps_data_model.hpp @@ -15,6 +15,19 @@ namespace cuopt::mathematical_optimization::io { +template +class coo_canonicalization_workspace_t; + +/** + * @brief Canonicalize a symmetric matrix in COO form to upper-triangular storage + * using reusable workspace storage. + */ +template +void canonicalize_coo_matrix(std::vector& rows, + std::vector& cols, + std::vector& vals, + coo_canonicalization_workspace_t& workspace); + /** * @brief A representation of a linear programming (LP) optimization problem * @@ -278,7 +291,8 @@ class mps_data_model_t { f_t rhs_value, std::span vals, std::span rows, - std::span cols); + std::span cols, + coo_canonicalization_workspace_t& workspace); const std::vector& get_quadratic_constraints() const; @@ -397,8 +411,7 @@ class mps_data_model_t { * * @param constraints Quadratic constraints whose rows, cols, and vals are updated. */ -template -void canonicalize_quadratic_constraints( - std::vector::quadratic_constraint_t>& constraints); +template +void canonicalize_quadratic_constraints(std::vector& constraints); } // namespace cuopt::mathematical_optimization::io diff --git a/cpp/include/cuopt/mathematical_optimization/optimization_problem.hpp b/cpp/include/cuopt/mathematical_optimization/optimization_problem.hpp index 5c755281ca..cab0d7522c 100644 --- a/cpp/include/cuopt/mathematical_optimization/optimization_problem.hpp +++ b/cpp/include/cuopt/mathematical_optimization/optimization_problem.hpp @@ -7,6 +7,7 @@ #pragma once +#include #include #include @@ -117,8 +118,9 @@ class optimization_problem_t : public optimization_problem_interface_t explicit optimization_problem_t(raft::handle_t const* handle_ptr); optimization_problem_t(const optimization_problem_t& other); - optimization_problem_t(optimization_problem_t&&) noexcept = default; - optimization_problem_t& operator=(optimization_problem_t&&) noexcept = default; + optimization_problem_t(optimization_problem_t&& other) noexcept; + optimization_problem_t& operator=(optimization_problem_t&& other) noexcept; + ~optimization_problem_t(); std::vector mip_callbacks_; @@ -411,6 +413,8 @@ class optimization_problem_t : public optimization_problem_interface_t /** QCQP: quadratic constraints **/ std::vector quadratic_constraints_{}; + // Scratch storage reused across add_quadratic_constraint calls. Not thread-safe. + std::unique_ptr> qc_coo_workspace_; rmm::device_uvector variable_lower_bounds_; rmm::device_uvector variable_upper_bounds_; diff --git a/cpp/include/cuopt/mathematical_optimization/optimization_problem_utils.hpp b/cpp/include/cuopt/mathematical_optimization/optimization_problem_utils.hpp index 7d511ef6f9..097b0b0403 100644 --- a/cpp/include/cuopt/mathematical_optimization/optimization_problem_utils.hpp +++ b/cpp/include/cuopt/mathematical_optimization/optimization_problem_utils.hpp @@ -294,7 +294,7 @@ void populate_from_data_model_view( // Raw Q COO from data_model_view is canonicalized here before solver storage. if (data_model->has_quadratic_constraints()) { auto qcs = data_model->get_quadratic_constraints(); - io::canonicalize_quadratic_constraints(qcs); + io::canonicalize_quadratic_constraints(qcs); problem->set_quadratic_constraints(std::move(qcs)); } } diff --git a/cpp/src/grpc/codegen/generate_conversions.py b/cpp/src/grpc/codegen/generate_conversions.py index 1cd4458230..1286418d1a 100644 --- a/cpp/src/grpc/codegen/generate_conversions.py +++ b/cpp/src/grpc/codegen/generate_conversions.py @@ -2215,12 +2215,12 @@ def _gen_proto_to_problem(registry, indent=" "): f'{b["name"]} size mismatch");' ) lines.append(f"{ind} }}") + lines.append(f"{ind} _entries.push_back(std::move(_entry));") + lines.append(f"{ind} }}") if name == "quadratic_constraints": lines.append( - f"{ind} io::canonicalize_coo_matrix(_entry.rows, _entry.cols, _entry.vals);" + f"{ind} io::canonicalize_quadratic_constraints(_entries);" ) - lines.append(f"{ind} _entries.push_back(std::move(_entry));") - lines.append(f"{ind} }}") lines.append(f"{ind} cpu_problem.{setter}(std::move(_entries));") lines.append(f"{ind}}}") @@ -2815,12 +2815,12 @@ def _gen_chunked_arrays_to_problem(registry, indent=" "): f'{b["name"]} size mismatch");' ) lines.append(f"{ind} }}") + lines.append(f"{ind} _entries.push_back(std::move(_entry));") + lines.append(f"{ind} }}") if name == "quadratic_constraints": lines.append( - f"{ind} io::canonicalize_coo_matrix(_entry.rows, _entry.cols, _entry.vals);" + f"{ind} io::canonicalize_quadratic_constraints(_entries);" ) - lines.append(f"{ind} _entries.push_back(std::move(_entry));") - lines.append(f"{ind} }}") lines.append(f"{ind} cpu_problem.{setter}(std::move(_entries));") lines.append(f"{ind}}}") diff --git a/cpp/src/grpc/codegen/generated/generated_chunked_arrays_to_problem.inc b/cpp/src/grpc/codegen/generated/generated_chunked_arrays_to_problem.inc index e71cc81179..514fd6f0d5 100644 --- a/cpp/src/grpc/codegen/generated/generated_chunked_arrays_to_problem.inc +++ b/cpp/src/grpc/codegen/generated/generated_chunked_arrays_to_problem.inc @@ -194,8 +194,8 @@ if (_entry.cols.size() != _entry.vals.size()) { throw std::invalid_argument("set_quadratic_constraints: cols/vals size mismatch"); } - io::canonicalize_coo_matrix(_entry.rows, _entry.cols, _entry.vals); _entries.push_back(std::move(_entry)); } + io::canonicalize_quadratic_constraints(_entries); cpu_problem.set_quadratic_constraints(std::move(_entries)); } diff --git a/cpp/src/grpc/codegen/generated/generated_proto_to_problem.inc b/cpp/src/grpc/codegen/generated/generated_proto_to_problem.inc index 7f99650caf..e0becdade5 100644 --- a/cpp/src/grpc/codegen/generated/generated_proto_to_problem.inc +++ b/cpp/src/grpc/codegen/generated/generated_proto_to_problem.inc @@ -96,8 +96,8 @@ if (_entry.cols.size() != _entry.vals.size()) { throw std::invalid_argument("set_quadratic_constraints: cols/vals size mismatch"); } - io::canonicalize_coo_matrix(_entry.rows, _entry.cols, _entry.vals); _entries.push_back(std::move(_entry)); } + io::canonicalize_quadratic_constraints(_entries); cpu_problem.set_quadratic_constraints(std::move(_entries)); } diff --git a/cpp/src/io/lp_parser.cpp b/cpp/src/io/lp_parser.cpp index cf9d8d6538..9a2adec236 100644 --- a/cpp/src/io/lp_parser.cpp +++ b/cpp/src/io/lp_parser.cpp @@ -9,6 +9,7 @@ #include #include +#include #include #include @@ -1471,6 +1472,7 @@ void flush_quadratic_constraints(mps_data_model_t& problem, const lp_parser_t& parser) { const i_t linear_row_count = static_cast(parser.row_names.size()); + coo_canonicalization_workspace_t qc_workspace; for (i_t k = 0; k < static_cast(parser.quadratic_constraint_blocks.size()); k++) { const auto& block = parser.quadratic_constraint_blocks[k]; problem.append_quadratic_constraint(linear_row_count + k, @@ -1481,7 +1483,8 @@ void flush_quadratic_constraints(mps_data_model_t& problem, block.rhs_value, block.quad_triples.vals, block.quad_triples.rows, - block.quad_triples.cols); + block.quad_triples.cols, + qc_workspace); } } diff --git a/cpp/src/io/mps_data_model.cpp b/cpp/src/io/mps_data_model.cpp index 941f16a96a..822cccf366 100644 --- a/cpp/src/io/mps_data_model.cpp +++ b/cpp/src/io/mps_data_model.cpp @@ -6,16 +6,147 @@ /* clang-format on */ #include +#include #include #include #include +#include +#include #include #include namespace cuopt::mathematical_optimization::io { +template +void canonicalize_coo_matrix(std::vector& rows, + std::vector& cols, + std::vector& vals, + coo_canonicalization_workspace_t& workspace) +{ + using workspace_t = coo_canonicalization_workspace_t; + + const size_t nnz = vals.size(); + mps_parser_expects(rows.size() == nnz && cols.size() == nnz, + error_type_t::ValidationError, + "COO rows/cols/vals length mismatch"); + if (nnz == 0) { + rows.clear(); + cols.clear(); + vals.clear(); + return; + } + + i_t max_index = 0; + for (size_t k = 0; k < nnz; ++k) { + mps_parser_expects(rows[k] >= 0 && cols[k] >= 0, + error_type_t::ValidationError, + "COO indices must be non-negative"); + max_index = std::max(max_index, std::max(rows[k], cols[k])); + } + workspace.ensure_index_capacity(max_index); + + std::vector& active_rows = workspace.active_rows_; + std::vector& offsets = workspace.offsets_; + std::vector& cursor = workspace.cursor_; + std::vector& permutation = workspace.permutation_; + std::vector& entries = workspace.entries_; + + active_rows.clear(); + offsets.clear(); + cursor.clear(); + permutation.clear(); + entries.clear(); + + const auto reset_lookup_entries = [&]() noexcept { + for (const auto& active : active_rows) { + workspace.row_perm_[active.row] = workspace_t::invalid_index; + } + for (const auto& entry : entries) { + workspace.col_position_[entry.col] = workspace_t::invalid_index; + } + }; + + // Discover active canonical rows and count entries in first-appearance order. + for (size_t k = 0; k < nnz; ++k) { + const i_t row = std::min(rows[k], cols[k]); + if (workspace.row_perm_[row] == workspace_t::invalid_index) { + const i_t id = active_rows.size(); + active_rows.push_back(typename workspace_t::active_row_t{row, 0}); + workspace.row_perm_[row] = id; + } + active_rows[workspace.row_perm_[row]].count++; + } + + std::sort(active_rows.begin(), active_rows.end(), [](const auto& lhs, const auto& rhs) { + return lhs.row < rhs.row; + }); + const size_t n_active_rows = active_rows.size(); + for (size_t id = 0; id < n_active_rows; ++id) { + workspace.row_perm_[active_rows[id].row] = id; + } + + offsets.resize(n_active_rows + 1); + offsets[0] = 0; + for (size_t id = 0; id < n_active_rows; ++id) { + offsets[id + 1] = offsets[id] + active_rows[id].count; + } + + cursor.assign(offsets.begin(), offsets.end() - 1); + permutation.resize(nnz); + for (size_t k = 0; k < nnz; ++k) { + const i_t row = std::min(rows[k], cols[k]); + const i_t dense_row = workspace.row_perm_[row]; + permutation[cursor[dense_row]] = k; + cursor[dense_row]++; + } + + entries.reserve(nnz); + for (size_t dense_row = 0; dense_row < n_active_rows; ++dense_row) { + const size_t row_output_start = entries.size(); + const i_t row = active_rows[dense_row].row; + + for (i_t p = offsets[dense_row]; p < offsets[dense_row + 1]; ++p) { + const i_t k = permutation[p]; + const i_t col = std::max(rows[k], cols[k]); + + if (workspace.col_position_[col] == workspace_t::invalid_index) { + const i_t position = entries.size(); + entries.push_back(typename workspace_t::entry_t{row, col, vals[k]}); + workspace.col_position_[col] = position; + } else { + entries[workspace.col_position_[col]].val += vals[k]; + } + } + + const size_t row_output_end = entries.size(); + for (size_t p = row_output_start; p < row_output_end; ++p) { + workspace.col_position_[entries[p].col] = workspace_t::invalid_index; + } + std::sort(entries.begin() + row_output_start, + entries.begin() + row_output_end, + [](const auto& lhs, const auto& rhs) { return lhs.col < rhs.col; }); + } + + const f_t eps = std::numeric_limits::epsilon(); + rows.clear(); + cols.clear(); + vals.clear(); + rows.reserve(entries.size()); + cols.reserve(entries.size()); + vals.reserve(entries.size()); + for (const auto& entry : entries) { + if (std::abs(entry.val) > eps) { + rows.push_back(entry.row); + cols.push_back(entry.col); + vals.push_back(entry.val); + } + } + + reset_lookup_entries(); +} + template void mps_data_model_t::set_csr_constraint_matrix(std::span A_values, std::span A_indices, @@ -142,15 +273,17 @@ void mps_data_model_t::set_quadratic_objective_matrix(std::span -void mps_data_model_t::append_quadratic_constraint(i_t constraint_row_index, - const std::string& constraint_row_name, - char constraint_row_type, - std::span linear_values, - std::span linear_indices, - f_t rhs_value, - std::span vals, - std::span rows, - std::span cols) +void mps_data_model_t::append_quadratic_constraint( + i_t constraint_row_index, + const std::string& constraint_row_name, + char constraint_row_type, + std::span linear_values, + std::span linear_indices, + f_t rhs_value, + std::span vals, + std::span rows, + std::span cols, + coo_canonicalization_workspace_t& workspace) { mps_parser_expects(constraint_row_index >= 0, error_type_t::ValidationError, @@ -199,7 +332,7 @@ void mps_data_model_t::append_quadratic_constraint(i_t constraint_row_ qc.rows.assign(rows.begin(), rows.end()); qc.cols.assign(cols.begin(), cols.end()); qc.vals.assign(vals.begin(), vals.end()); - canonicalize_coo_matrix(qc.rows, qc.cols, qc.vals); + canonicalize_coo_matrix(qc.rows, qc.cols, qc.vals, workspace); } quadratic_constraints_.push_back(std::move(qc)); @@ -459,24 +592,44 @@ bool mps_data_model_t::has_quadratic_constraints() const noexcept return !quadratic_constraints_.empty(); } -template -void canonicalize_quadratic_constraints( - std::vector::quadratic_constraint_t>& constraints) +template +void canonicalize_quadratic_constraints(std::vector& constraints) { + using i_t = typename decltype(qc_t::rows)::value_type; + using f_t = typename decltype(qc_t::vals)::value_type; + coo_canonicalization_workspace_t workspace; for (auto& qc : constraints) { - canonicalize_coo_matrix(qc.rows, qc.cols, qc.vals); + canonicalize_coo_matrix(qc.rows, qc.cols, qc.vals, workspace); } } // NOTE: Explicitly instantiate all types here in order to avoid linker error +template void +canonicalize_quadratic_constraints::quadratic_constraint_t>( + std::vector::quadratic_constraint_t>&); +template void +canonicalize_quadratic_constraints::quadratic_constraint_t>( + std::vector::quadratic_constraint_t>&); +template void canonicalize_quadratic_constraints< + optimization_problem_interface_t::quadratic_constraint_t>( + std::vector::quadratic_constraint_t>&); +template void canonicalize_quadratic_constraints< + optimization_problem_interface_t::quadratic_constraint_t>( + std::vector::quadratic_constraint_t>&); + +template void canonicalize_coo_matrix(std::vector&, + std::vector&, + std::vector&, + coo_canonicalization_workspace_t&); +template void canonicalize_coo_matrix(std::vector&, + std::vector&, + std::vector&, + coo_canonicalization_workspace_t&); + template class mps_data_model_t; template class mps_data_model_t; -template void canonicalize_quadratic_constraints( - std::vector::quadratic_constraint_t>&); -template void canonicalize_quadratic_constraints( - std::vector::quadratic_constraint_t>&); // TODO current raft to cusparse wrappers only support int64_t // can be CUSPARSE_INDEX_16U, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_64I diff --git a/cpp/src/io/mps_parser.cpp b/cpp/src/io/mps_parser.cpp index 481df1ca5b..a77ca0d53f 100644 --- a/cpp/src/io/mps_parser.cpp +++ b/cpp/src/io/mps_parser.cpp @@ -242,94 +242,6 @@ void check_symmetric_offdiagonal_pairs(const std::vector& rows, } } -template -void canonicalize_coo_matrix(std::vector& rows, std::vector& cols, std::vector& vals) -{ - const size_t n = vals.size(); - mps_parser_expects(rows.size() == n && cols.size() == n, - error_type_t::ValidationError, - "COO rows/cols/vals length mismatch"); - if (n == 0) { - rows.clear(); - cols.clear(); - vals.clear(); - return; - } - - // Canonical coordinates place each entry in the upper triangle: - // row = min(r, c), col = max(r, c). - i_t min_idx = std::min(rows[0], cols[0]); - i_t max_row = 0; - i_t max_col = 0; - for (size_t k = 0; k < n; ++k) { - const i_t r = std::min(rows[k], cols[k]); - const i_t c = std::max(rows[k], cols[k]); - min_idx = std::min(min_idx, r); - max_row = std::max(max_row, r); - max_col = std::max(max_col, c); - } - - // Bucket entry indices by canonical row via a stable counting sort. - std::vector offsets(static_cast(max_row - min_idx) + 2, 0); - for (size_t k = 0; k < n; ++k) { - offsets[static_cast(std::min(rows[k], cols[k]) - min_idx) + 1]++; - } - for (i_t r = 0; r <= max_row - min_idx; ++r) { - offsets[static_cast(r) + 1] += offsets[static_cast(r)]; - } - std::vector perm(n); - std::vector cursor(offsets.begin(), offsets.end() - 1); - for (size_t k = 0; k < n; ++k) { - const i_t r = std::min(rows[k], cols[k]) - min_idx; - perm[static_cast(cursor[r])] = static_cast(k); - cursor[static_cast(r)]++; - } - - // Merge duplicate (row, col) entries within each canonical row using a marker array - // (marker[c] holds 1 + the output position of column c in the current row), then sort - // each row's entries by column. - struct entry_t { - i_t row; - i_t col; - f_t val; - }; - std::vector entries; - entries.reserve(n); - std::vector marker(static_cast(max_col - min_idx) + 1, 0); - for (i_t r = 0; r <= max_row - min_idx; ++r) { - const i_t row_start = static_cast(entries.size()); - for (i_t p = offsets[static_cast(r)]; p < offsets[static_cast(r) + 1]; ++p) { - const i_t k = perm[static_cast(p)]; - const i_t row = std::min(rows[static_cast(k)], cols[static_cast(k)]); - const i_t col = std::max(rows[static_cast(k)], cols[static_cast(k)]); - const i_t c = col - min_idx; - if (marker[static_cast(c)] <= row_start) { - entries.push_back(entry_t{row, col, vals[static_cast(k)]}); - marker[static_cast(c)] = static_cast(entries.size()); - } else { - entries[static_cast(marker[static_cast(c)] - 1)].val += - vals[static_cast(k)]; - } - } - std::sort(entries.begin() + row_start, - entries.end(), - [](const entry_t& lhs, const entry_t& rhs) { return lhs.col < rhs.col; }); - } - - // Emit canonical COO, dropping entries whose coefficients cancelled to ~0. - const f_t eps = std::numeric_limits::epsilon(); - rows.clear(); - cols.clear(); - vals.clear(); - for (const entry_t& e : entries) { - if (std::abs(e.val) > eps) { - rows.push_back(e.row); - cols.push_back(e.col); - vals.push_back(e.val); - } - } -} - BoundType convert(std::string_view str) { if (str == "LO") { @@ -619,6 +531,7 @@ void mps_parser_t::fill_problem(mps_data_model_t& problem) } // QCMATRIX: one symmetric Q per constraint row (no extra ½ factor vs file coeffs). + coo_canonicalization_workspace_t qc_workspace; for (const auto& block : qcmatrix_blocks_) { const i_t row_id = block.constraint_row_id; mps_parser_expects(row_id >= 0 && row_id < static_cast(row_types.size()), @@ -633,7 +546,8 @@ void mps_parser_t::fill_problem(mps_data_model_t& problem) b_values[row_id], block.entries.vals, block.entries.rows, - block.entries.cols); + block.entries.cols, + qc_workspace); } if (!quadratic_row_ids.empty()) { @@ -1712,11 +1626,4 @@ template class mps_parser_t; template class mps_parser_t; -template void canonicalize_coo_matrix(std::vector&, - std::vector&, - std::vector&); -template void canonicalize_coo_matrix(std::vector&, - std::vector&, - std::vector&); - } // namespace cuopt::mathematical_optimization::io diff --git a/cpp/src/io/mps_parser_internal.hpp b/cpp/src/io/mps_parser_internal.hpp index 2bc5443f11..5199a00cd2 100644 --- a/cpp/src/io/mps_parser_internal.hpp +++ b/cpp/src/io/mps_parser_internal.hpp @@ -12,12 +12,66 @@ #include #include #include +#include #include #include #include namespace cuopt::mathematical_optimization::io { +/** + * @brief Reusable workspace for COO canonicalization. + * + * Lookup vectors grow to fit the largest variable index seen so far. Each call resets only + * entries that were touched. A workspace must not be shared across concurrent calls. + */ +template +class coo_canonicalization_workspace_t { + public: + static_assert(std::is_integral_v && std::is_signed_v, + "coo_canonicalization_workspace_t: i_t must be signed integral"); + + [[nodiscard]] i_t max_index_capacity() const noexcept { return row_perm_.size(); } + + private: + static constexpr i_t invalid_index = -1; + + struct active_row_t { + i_t row; + i_t count; + }; + + struct entry_t { + i_t row; + i_t col; + f_t val; + }; + + std::vector row_perm_; + std::vector col_position_; + std::vector active_rows_; + std::vector offsets_; + std::vector cursor_; + std::vector permutation_; + std::vector entries_; + + void ensure_index_capacity(i_t max_index) + { + const size_t needed = static_cast(max_index) + 1; + if (row_perm_.size() < needed) { + row_perm_.resize(needed, invalid_index); + col_position_.resize(needed, invalid_index); + } + } + + template + friend void canonicalize_coo_matrix( + std::vector& rows, + std::vector& cols, + std::vector& vals, + coo_canonicalization_workspace_t& workspace); +}; + /** * Sparse COO (coordinate) entries for a matrix: parallel row/col/val vectors. */ @@ -54,19 +108,6 @@ void check_symmetric_offdiagonal_pairs(const std::vector& rows, const std::vector& cols, const std::vector& vals); -/** - * @brief Canonicalize a symmetric matrix in COO form to upper-triangular storage. - * - * - Merges duplicate (row, col) indices by summing coefficients. - * - Merges each off-diagonal variable pair into one (min,max) entry with coefficient - * Q[min,max]+Q[max,min] (x^T Q x semantics). - * - Sorts output by (row, col). - */ -template -void canonicalize_coo_matrix(std::vector& rows, - std::vector& cols, - std::vector& vals); - /** * @brief Different possible types of 'ROWS' */ diff --git a/cpp/src/io/mps_writer.cpp b/cpp/src/io/mps_writer.cpp index d269d6ec8a..39c3d8f3c2 100644 --- a/cpp/src/io/mps_writer.cpp +++ b/cpp/src/io/mps_writer.cpp @@ -497,10 +497,11 @@ void mps_writer_t::write(const std::string& mps_file_path) // QCMATRIX sections for quadratic constraints (QCQP) if (problem_.has_quadratic_constraints()) { + coo_canonicalization_workspace_t workspace; for (const auto& qc : problem_.get_quadratic_constraints()) { mps_file << "QCMATRIX " << qc.constraint_row_name << "\n"; typename mps_data_model_t::quadratic_constraint_t qc_canon = qc; - canonicalize_coo_matrix(qc_canon.rows, qc_canon.cols, qc_canon.vals); + canonicalize_coo_matrix(qc_canon.rows, qc_canon.cols, qc_canon.vals, workspace); const i_t nnz = static_cast(qc_canon.vals.size()); for (i_t p = 0; p < nnz; ++p) { const i_t i = qc_canon.rows[p]; diff --git a/cpp/src/pdlp/cpu_optimization_problem.cpp b/cpp/src/pdlp/cpu_optimization_problem.cpp index 84f03b5e85..f7feeade9b 100644 --- a/cpp/src/pdlp/cpu_optimization_problem.cpp +++ b/cpp/src/pdlp/cpu_optimization_problem.cpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -31,8 +32,87 @@ namespace cuopt::mathematical_optimization { template cpu_optimization_problem_t::cpu_optimization_problem_t() -{ -} + : qc_coo_workspace_(std::make_unique>()) +{ +} + +template +cpu_optimization_problem_t::cpu_optimization_problem_t( + const cpu_optimization_problem_t& other) + : problem_category_(other.problem_category_), + maximize_(other.maximize_), + n_vars_(other.n_vars_), + n_constraints_(other.n_constraints_), + A_(other.A_), + A_indices_(other.A_indices_), + A_offsets_(other.A_offsets_), + b_(other.b_), + c_(other.c_), + objective_scaling_factor_(other.objective_scaling_factor_), + objective_offset_(other.objective_offset_), + Q_offsets_(other.Q_offsets_), + Q_indices_(other.Q_indices_), + Q_values_(other.Q_values_), + quadratic_constraints_(other.quadratic_constraints_), + qc_coo_workspace_(std::make_unique>()), + variable_lower_bounds_(other.variable_lower_bounds_), + variable_upper_bounds_(other.variable_upper_bounds_), + constraint_lower_bounds_(other.constraint_lower_bounds_), + constraint_upper_bounds_(other.constraint_upper_bounds_), + row_types_(other.row_types_), + variable_types_(other.variable_types_), + objective_name_(other.objective_name_), + problem_name_(other.problem_name_), + var_names_(other.var_names_), + row_names_(other.row_names_) +{ +} + +template +cpu_optimization_problem_t::cpu_optimization_problem_t( + cpu_optimization_problem_t&& other) noexcept = default; + +template +cpu_optimization_problem_t& cpu_optimization_problem_t::operator=( + const cpu_optimization_problem_t& other) +{ + if (this == &other) { return *this; } + + problem_category_ = other.problem_category_; + maximize_ = other.maximize_; + n_vars_ = other.n_vars_; + n_constraints_ = other.n_constraints_; + A_ = other.A_; + A_indices_ = other.A_indices_; + A_offsets_ = other.A_offsets_; + b_ = other.b_; + c_ = other.c_; + objective_scaling_factor_ = other.objective_scaling_factor_; + objective_offset_ = other.objective_offset_; + Q_offsets_ = other.Q_offsets_; + Q_indices_ = other.Q_indices_; + Q_values_ = other.Q_values_; + quadratic_constraints_ = other.quadratic_constraints_; + qc_coo_workspace_ = std::make_unique>(); + variable_lower_bounds_ = other.variable_lower_bounds_; + variable_upper_bounds_ = other.variable_upper_bounds_; + constraint_lower_bounds_ = other.constraint_lower_bounds_; + constraint_upper_bounds_ = other.constraint_upper_bounds_; + row_types_ = other.row_types_; + variable_types_ = other.variable_types_; + objective_name_ = other.objective_name_; + problem_name_ = other.problem_name_; + var_names_ = other.var_names_; + row_names_ = other.row_names_; + return *this; +} + +template +cpu_optimization_problem_t& cpu_optimization_problem_t::operator=( + cpu_optimization_problem_t&& other) noexcept = default; + +template +cpu_optimization_problem_t::~cpu_optimization_problem_t() = default; // ============================================================================== // Setters @@ -172,7 +252,10 @@ void cpu_optimization_problem_t::add_quadratic_constraint( qc.vals.assign(coeff.begin(), coeff.end()); qc.linear_values.assign(linear_values.begin(), linear_values.end()); qc.linear_indices.assign(linear_indices.begin(), linear_indices.end()); - io::canonicalize_coo_matrix(qc.rows, qc.cols, qc.vals); + if (qc_coo_workspace_ == nullptr) { + qc_coo_workspace_ = std::make_unique>(); + } + io::canonicalize_coo_matrix(qc.rows, qc.cols, qc.vals, *qc_coo_workspace_); quadratic_constraints_.push_back(std::move(qc)); } diff --git a/cpp/src/pdlp/optimization_problem.cu b/cpp/src/pdlp/optimization_problem.cu index 85d3ef4df4..dd222d9ee5 100644 --- a/cpp/src/pdlp/optimization_problem.cu +++ b/cpp/src/pdlp/optimization_problem.cu @@ -48,6 +48,7 @@ #include #include +#include #include #include @@ -62,6 +63,7 @@ optimization_problem_t::optimization_problem_t(raft::handle_t const* h A_offsets_(0, stream_view_), b_(0, stream_view_), c_(0, stream_view_), + qc_coo_workspace_(std::make_unique>()), variable_lower_bounds_(0, stream_view_), variable_upper_bounds_(0, stream_view_), constraint_lower_bounds_(0, stream_view_), @@ -102,10 +104,22 @@ optimization_problem_t::optimization_problem_t( problem_category_{other.get_problem_category()}, var_names_{other.get_variable_names()}, row_names_{other.get_row_names()}, - quadratic_constraints_{other.get_quadratic_constraints()} + quadratic_constraints_{other.get_quadratic_constraints()}, + qc_coo_workspace_(std::make_unique>()) { } +template +optimization_problem_t::optimization_problem_t( + optimization_problem_t&& other) noexcept = default; + +template +optimization_problem_t& optimization_problem_t::operator=( + optimization_problem_t&& other) noexcept = default; + +template +optimization_problem_t::~optimization_problem_t() = default; + // ============================================================================== // Setters // ============================================================================== @@ -244,7 +258,10 @@ void optimization_problem_t::add_quadratic_constraint(char constraint_ qc.vals.assign(coeff.begin(), coeff.end()); qc.linear_values.assign(linear_values.begin(), linear_values.end()); qc.linear_indices.assign(linear_indices.begin(), linear_indices.end()); - io::canonicalize_coo_matrix(qc.rows, qc.cols, qc.vals); + if (qc_coo_workspace_ == nullptr) { + qc_coo_workspace_ = std::make_unique>(); + } + io::canonicalize_coo_matrix(qc.rows, qc.cols, qc.vals, *qc_coo_workspace_); quadratic_constraints_.push_back(std::move(qc)); } diff --git a/cpp/tests/linear_programming/parser_test.cpp b/cpp/tests/linear_programming/parser_test.cpp index a0bb3483ef..a11e7587da 100644 --- a/cpp/tests/linear_programming/parser_test.cpp +++ b/cpp/tests/linear_programming/parser_test.cpp @@ -2728,10 +2728,11 @@ TEST(append_quadratic_constraint, merges_duplicate_entries) { using model_t = mps_data_model_t; model_t model; + coo_canonicalization_workspace_t workspace; const std::vector vals = {2.0, 3.0}; const std::vector rows = {0, 0}; const std::vector cols = {1, 1}; - model.append_quadratic_constraint(0, "QC0", 'L', {}, {}, 0.0, vals, rows, cols); + model.append_quadratic_constraint(0, "QC0", 'L', {}, {}, 0.0, vals, rows, cols, workspace); ASSERT_TRUE(model.has_quadratic_constraints()); const auto& qc = model.get_quadratic_constraints().back(); @@ -2745,10 +2746,11 @@ TEST(append_quadratic_constraint, collapses_symmetric_mps_halves) { using model_t = mps_data_model_t; model_t model; + coo_canonicalization_workspace_t workspace; const std::vector vals = {2.0, 2.0}; const std::vector rows = {0, 1}; const std::vector cols = {1, 0}; - model.append_quadratic_constraint(0, "QC0", 'L', {}, {}, 0.0, vals, rows, cols); + model.append_quadratic_constraint(0, "QC0", 'L', {}, {}, 0.0, vals, rows, cols, workspace); ASSERT_TRUE(model.has_quadratic_constraints()); const auto& qc = model.get_quadratic_constraints().back(); @@ -2762,10 +2764,11 @@ TEST(append_quadratic_constraint, sums_both_orientations_for_off_diagonal_pair) { using model_t = mps_data_model_t; model_t model; + coo_canonicalization_workspace_t workspace; const std::vector vals = {2.0, 3.0}; const std::vector rows = {0, 1}; const std::vector cols = {1, 0}; - model.append_quadratic_constraint(0, "QC0", 'L', {}, {}, 0.0, vals, rows, cols); + model.append_quadratic_constraint(0, "QC0", 'L', {}, {}, 0.0, vals, rows, cols, workspace); ASSERT_TRUE(model.has_quadratic_constraints()); const auto& qc = model.get_quadratic_constraints().back(); @@ -2779,6 +2782,7 @@ TEST(qps_parser, qcmatrix_append_api) { using model_t = mps_data_model_t; model_t model; + coo_canonicalization_workspace_t workspace; // Validate default-constructed struct shape. model_t::quadratic_constraint_t default_qcm; @@ -2804,7 +2808,8 @@ TEST(qps_parser, qcmatrix_append_api) 5.0, mps_qc0_values, mps_qc0_row_indices, - mps_qc0_col_indices); + mps_qc0_col_indices, + workspace); // API-style canonical COO [[4, 2], [2, 6]] -> stored unchanged after merge/sort const std::vector api_qc1_values = {4.0, 2.0, 6.0}; @@ -2820,7 +2825,8 @@ TEST(qps_parser, qcmatrix_append_api) 10.0, api_qc1_values, api_qc1_row_indices, - api_qc1_col_indices); + api_qc1_col_indices, + workspace); ASSERT_TRUE(model.has_quadratic_constraints()); const auto& qcs = model.get_quadratic_constraints();