Write quadratic terms to mps file

cpp/libmps_parser/include/mps_parser/mps_writer.hpp

@@ -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>
@@ -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
    *
@@ -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

cpp/libmps_parser/src/mps_writer.cpp

@@ -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 {
 
@@ -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()));
+  }
+
+  // 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)
 {
@@ -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();
 }

cpp/libmps_parser/src/utilities/sparse_matrix_utils.hpp

@@ -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++;
+      }
+    }
+  }
+
+  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);
+}
+
+}  // namespace cuopt::mps_parser