Skip to content

Write quadratic terms to mps file #949

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
rg20 wants to merge 7 commits into
base: release/26.04
Choose a base branch
from
Open

Write quadratic terms to mps file #949

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
85902c3
speedup symmetrization
rg20 Feb 10, 2026
bbdbbdb
separate out heper function
rg20 Feb 10, 2026
e9f7a79
Support quadobj writing to mps
rg20 Feb 12, 2026
5b5a2bb
add missing header
rg20 Feb 13, 2026
eb79da4
Use kahan summation for post solve residual computation
rg20 Mar 5, 2026
1b94b55
Handle core review comment
rg20 Mar 16, 2026
0f3ca8e
Merge remote-tracking branch 'upstream/release/26.04' into qp_improve…
rg20 Mar 20, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 16 additions & 2 deletions cpp/libmps_parser/include/mps_parser/mps_writer.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,16 +1,18 @@
/* 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 */

#pragma once

#include <mps_parser/data_model_view.hpp>
#include <mps_parser/mps_data_model.hpp>

#include <stdarg.h>
#include <limits>
#include <memory>
#include <string>
#include <unordered_map>
#include <unordered_set>
Expand All @@ -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<i_t, f_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<i_t, f_t>& problem);

/**
* @brief Writes the problem to an MPS formatted file
*
Expand All @@ -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<data_model_view_t<i_t, f_t>> owned_view_;
// Reference to the view (either external or owned)
const data_model_view_t<i_t, f_t>& problem_;

// Helper to create view from data model
static data_model_view_t<i_t, f_t> create_view(const mps_data_model_t<i_t, f_t>& model);
}; // class mps_writer_t

} // namespace cuopt::mps_parser
139 changes: 139 additions & 0 deletions cpp/libmps_parser/src/mps_writer.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,14 +8,17 @@
#include <mps_parser/mps_writer.hpp>

#include <mps_parser/data_model_view.hpp>
#include <mps_parser/mps_data_model.hpp>
#include <utilities/error.hpp>
#include <utilities/sparse_matrix_utils.hpp>

#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <limits>
#include <map>
#include <memory>

namespace cuopt::mps_parser {

Expand All @@ -24,6 +27,92 @@ mps_writer_t<i_t, f_t>::mps_writer_t(const data_model_view_t<i_t, f_t>& problem)
{
}

template <typename i_t, typename f_t>
data_model_view_t<i_t, f_t> mps_writer_t<i_t, f_t>::create_view(
const mps_data_model_t<i_t, f_t>& model)
{
data_model_view_t<i_t, f_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<i_t>(A_values.size()),
A_indices.data(),
static_cast<i_t>(A_indices.size()),
A_offsets.data(),
static_cast<i_t>(A_offsets.size()));
}
Comment on lines +43 to +50
Copy link

@coderabbitai coderabbitai bot Mar 16, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🔴 Critical

Don’t skip CSR setup when A_values is empty

Line 43 should not gate on values. A valid zero-nnz CSR still needs offsets; otherwise write() can read constraint_matrix_offsets[row_id + 1] out of bounds for constrained-but-empty matrices.

✅ Proposed fix 
-  if (!A_values.empty()) {
+  if (!A_offsets.empty()) {
     view.set_csr_constraint_matrix(A_values.data(),
                                    static_cast<i_t>(A_values.size()),
                                    A_indices.data(),
                                    static_cast<i_t>(A_indices.size()),
                                    A_offsets.data(),
                                    static_cast<i_t>(A_offsets.size()));
   }
As per coding guidelines, "Validate algorithm correctness in optimization logic: simplex pivots, branch-and-bound decisions, routing heuristics, and constraint/objective handling must produce correct results."
🤖 Prompt for AI Agents 
Verify each finding against the current code and only fix it if needed.

In `@cpp/libmps_parser/src/mps_writer.cpp` around lines 43 - 50, The code
currently skips calling view.set_csr_constraint_matrix when A_values.empty(),
which breaks zero-nnz CSR matrices because write() later accesses
constraint_matrix_offsets[row_id + 1]; instead ensure set_csr_constraint_matrix
is invoked even when A_values is empty by removing the A_values.empty() guard
and calling view.set_csr_constraint_matrix with the existing A_offsets (and
A_values/A_indices pointers with size 0 as appropriate) whenever the constraint
matrix structure (A_offsets) is present or the matrix is constrained; update the
call site (the block around view.set_csr_constraint_matrix) to pass correct zero
sizes (static_cast<i_t>(0)) for values/indices when A_values/A_indices are empty
so constraint_matrix_offsets is always initialized for write().


// Constraint bounds
const auto& b = model.get_constraint_bounds();
if (!b.empty()) { view.set_constraint_bounds(b.data(), static_cast<i_t>(b.size())); }

// Objective coefficients
const auto& c = model.get_objective_coefficients();
if (!c.empty()) { view.set_objective_coefficients(c.data(), static_cast<i_t>(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<i_t>(lb.size())); }
if (!ub.empty()) { view.set_variable_upper_bounds(ub.data(), static_cast<i_t>(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<i_t>(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<i_t>(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<i_t>(cl.size())); }
if (!cu.empty()) { view.set_constraint_upper_bounds(cu.data(), static_cast<i_t>(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<i_t>(Q_values.size()),
Q_indices.data(),
static_cast<i_t>(Q_indices.size()),
Q_offsets.data(),
static_cast<i_t>(Q_offsets.size()));
}

return view;
}

template <typename i_t, typename f_t>
mps_writer_t<i_t, f_t>::mps_writer_t(const mps_data_model_t<i_t, f_t>& problem)
: owned_view_(std::make_unique<data_model_view_t<i_t, f_t>>(create_view(problem))),
problem_(*owned_view_)
{
}

template <typename i_t, typename f_t>
void mps_writer_t<i_t, f_t>::write(const std::string& mps_file_path)
{
Expand Down Expand Up @@ -280,6 +369,56 @@ void mps_writer_t<i_t, f_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<f_t> Q_values(Q_values_span.data(), Q_values_span.data() + Q_values_span.size());
std::vector<i_t> Q_indices(Q_indices_span.data(),
Q_indices_span.data() + Q_indices_span.size());
std::vector<i_t> 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<f_t> H_values;
std::vector<i_t> H_indices;
std::vector<i_t> H_offsets;
symmetrize_csr(Q_values, Q_indices, Q_offsets, H_values, H_indices, H_offsets);

i_t n_rows = static_cast<i_t>(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<size_t>(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<size_t>(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();
}
Expand Down
137 changes: 137 additions & 0 deletions cpp/libmps_parser/src/utilities/sparse_matrix_utils.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,137 @@
/* clang-format off */
/*
* SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
* SPDX-License-Identifier: Apache-2.0
*/
/* clang-format on */

#pragma once

#include <vector>

namespace cuopt::mps_parser {

/**
* @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]).
*
* @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 <typename i_t, typename f_t>
void symmetrize_csr(const f_t* in_values,
const i_t* in_indices,
const i_t* in_offsets,
i_t n_rows,
std::vector<f_t>& out_values,
std::vector<i_t>& out_indices,
std::vector<i_t>& out_offsets)
{
// Optimized 3-pass algorithm
// Pass 1: Count entries per row in A + A^T
std::vector<i_t> 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<i_t> 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<i_t> temp_indices(total_entries);
std::vector<f_t> temp_values(total_entries);

// Pass 2: Fill entries directly
std::vector<i_t> 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<i_t> 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++;
}
}
Comment on lines +86 to +108
Copy link

@coderabbitai coderabbitai bot Mar 16, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟠 Major

Workspace sentinel -1 breaks for unsigned index types.

The workspace is initialized with -1 (line 86), and the check at line 100 uses workspace[j] >= row_start_out to detect whether a column was already visited in the current row. If i_t is an unsigned type (e.g., uint32_t), -1 wraps to std::numeric_limits<i_t>::max(), which will always satisfy >= row_start_out, causing incorrect value accumulation on unvisited entries.

Either enforce signed types via static_assert or use a sentinel that works for both:

🛠️ Proposed fix using explicit sentinel 
+#include <limits>
+#include <type_traits>
+
   // Pass 3: Deduplicate and build final CSR
-  std::vector<i_t> workspace(n_rows, -1);
+  constexpr i_t SENTINEL = std::is_signed<i_t>::value 
+                           ? static_cast<i_t>(-1) 
+                           : std::numeric_limits<i_t>::max();
+  std::vector<i_t> workspace(n_rows, SENTINEL);
   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) {
+      if (workspace[j] != SENTINEL && workspace[j] >= row_start_out) {
         out_values[workspace[j]] += x;
       } else {
         workspace[j]    = nz;
         out_indices[nz] = j;
         out_values[nz]  = x;
         nz++;
       }
     }
   }

Alternatively, add a static_assert to enforce signed index types:

static_assert(std::is_signed<i_t>::value, "i_t must be a signed integer type");
🤖 Prompt for AI Agents 
Verify each finding against the current code and only fix it if needed.

In `@cpp/libmps_parser/src/utilities/sparse_matrix_utils.hpp` around lines 86 -
108, The workspace sentinel uses -1 which breaks when i_t is unsigned; add a
compile-time guard to require signed index types by inserting
static_assert(std::is_signed<i_t>::value, "i_t must be a signed integer type");
near the top of the function/template (before workspace is initialized) so
workspace, temp_offsets, temp_indices and related logic remain correct for all
builds; alternatively, replace the -1 sentinel with a sentinel value based on
std::numeric_limits<i_t>::min() and adjust the visited check accordingly if you
prefer a runtime-compatible sentinel instead of the static_assert.

}

out_offsets[n_rows] = nz;
out_indices.resize(nz);
out_values.resize(nz);
}

/**
* @brief Symmetrize a CSR matrix using std::vector (convenience overload)
*/
template <typename i_t, typename f_t>
void symmetrize_csr(const std::vector<f_t>& in_values,
const std::vector<i_t>& in_indices,
const std::vector<i_t>& in_offsets,
std::vector<f_t>& out_values,
std::vector<i_t>& out_indices,
std::vector<i_t>& out_offsets)
{
i_t n_rows = static_cast<i_t>(in_offsets.size()) - 1;
symmetrize_csr(in_values.data(),
in_indices.data(),
in_offsets.data(),
n_rows,
out_values,
out_indices,
out_offsets);
}
Comment on lines +119 to +135
Copy link

@coderabbitai coderabbitai bot Mar 10, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟡 Minor

Handle edge case: empty in_offsets vector.

When in_offsets is empty, in_offsets.size() - 1 will underflow (for unsigned) or produce an unexpected value. Add a guard for this edge case.

🛡️ Proposed fix 
 template <typename i_t, typename f_t>
 void symmetrize_csr(const std::vector<f_t>& in_values,
                     const std::vector<i_t>& in_indices,
                     const std::vector<i_t>& in_offsets,
                     std::vector<f_t>& out_values,
                     std::vector<i_t>& out_indices,
                     std::vector<i_t>& out_offsets)
 {
+  if (in_offsets.empty()) {
+    out_values.clear();
+    out_indices.clear();
+    out_offsets.clear();
+    return;
+  }
   i_t n_rows = static_cast<i_t>(in_offsets.size()) - 1;
   symmetrize_csr(in_values.data(),
                  in_indices.data(),
                  in_offsets.data(),
                  n_rows,
                  out_values,
                  out_indices,
                  out_offsets);
 }
🤖 Prompt for AI Agents 
Verify each finding against the current code and only fix it if needed.

In `@cpp/libmps_parser/src/utilities/sparse_matrix_utils.hpp` around lines 119 -
135, The wrapper symmetrize_csr must handle an empty in_offsets to avoid
underflow when computing n_rows; add a guard at the start of the function that
checks if in_offsets.empty() and if so clear out_values and out_indices and set
out_offsets to a single element {0} (or otherwise initialize outputs to
represent an empty CSR) then return early; otherwise proceed to compute n_rows =
in_offsets.size() - 1 and call the underlying symmetrize_csr overload as before.


} // namespace cuopt::mps_parser
Loading
Loading