Skip to content

refactor: promote shared primitives to mathematical_optimization #1481

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

Merged
rapids-bot merged 1 commit into from
Jun 30, 2026
Merged

refactor: promote shared primitives to mathematical_optimization #1481

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
1 change: 1 addition & 0 deletions cpp/src/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ set(UTIL_SRC_FILES ${CMAKE_CURRENT_SOURCE_DIR}/utilities/seed_generator.cu
${CMAKE_CURRENT_SOURCE_DIR}/utilities/timestamp_utils.cpp
${CMAKE_CURRENT_SOURCE_DIR}/utilities/work_unit_scheduler.cpp)

add_subdirectory(linear_algebra)
add_subdirectory(pdlp)
add_subdirectory(math_optimization)
add_subdirectory(mip_heuristics)
Expand Down
90 changes: 38 additions & 52 deletions cpp/src/barrier/barrier.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,23 +10,23 @@
#include <barrier/conjugate_gradient.hpp>
#include <barrier/cusparse_info.hpp>
#include <barrier/cusparse_view.hpp>
#include <barrier/dense_matrix.hpp>
#include <barrier/dense_vector.hpp>
#include <barrier/device_sparse_matrix.cuh>
#include <barrier/iterative_refinement.hpp>
#include <barrier/pinned_host_allocator.hpp>
#include <barrier/second_order_cone_kernels.cuh>
#include <barrier/sparse_cholesky.cuh>
#include <barrier/sparse_matrix_kernels.cuh>
#include <linear_algebra/dense_matrix.hpp>
#include <linear_algebra/dense_vector.hpp>

#include <dual_simplex/presolve.hpp>
#include <dual_simplex/solve.hpp>

#include <dual_simplex/sparse_matrix.hpp>
#include <dual_simplex/tic_toc.hpp>
#include <dual_simplex/types.hpp>
#include <linear_algebra/sparse_matrix.hpp>
#include <math_optimization/tic_toc.hpp>
#include <math_optimization/types.hpp>

#include <dual_simplex/vector_math.cuh>
#include <linear_algebra/vector_math.cuh>

#include <rmm/device_scalar.hpp>
#include <rmm/device_uvector.hpp>
Expand All @@ -53,20 +53,10 @@
namespace cuopt::mathematical_optimization::barrier {

using simplex::compute_user_objective;
using simplex::csc_matrix_t;
using simplex::csr_matrix_t;
using simplex::device_vector_norm_inf;
using simplex::float64_t;
using simplex::inf;
using simplex::lp_problem_t;
using simplex::lp_solution_t;
using simplex::lp_status_t;
using simplex::matrix_vector_multiply;
using simplex::multiply;
using simplex::simplex_solver_settings_t;
using simplex::tic;
using simplex::toc;
using simplex::vector_norm1;

template <typename i_t, typename f_t>
bool validate_barrier_cone_layout(const lp_problem_t<i_t, f_t>& problem,
Expand Down Expand Up @@ -1024,9 +1014,9 @@ class iteration_data_t {
{
if (n_dense_columns == 0) {
// Solve ADAT * x = b
if (debug) { settings_.log.printf("||b|| = %.16e\n", simplex::vector_norm2<i_t, f_t>(b)); }
if (debug) { settings_.log.printf("||b|| = %.16e\n", vector_norm2<i_t, f_t>(b)); }
i_t solve_status = chol->solve(b, x);
if (debug) { settings_.log.printf("||x|| = %.16e\n", simplex::vector_norm2<i_t, f_t>(x)); }
if (debug) { settings_.log.printf("||x|| = %.16e\n", vector_norm2<i_t, f_t>(x)); }
return solve_status;
} else {
// Use Sherman Morrison followed by PCG
Expand Down Expand Up @@ -1062,9 +1052,9 @@ class iteration_data_t {
dense_vector_t<i_t, f_t> w(AD.m);
const bool debug = false;
const bool full_debug = false;
if (debug) { settings_.log.printf("||b|| = %.16e\n", simplex::vector_norm2<i_t, f_t>(b)); }
if (debug) { settings_.log.printf("||b|| = %.16e\n", vector_norm2<i_t, f_t>(b)); }
i_t solve_status = chol->solve(b, w);
if (debug) { settings_.log.printf("||w|| = %.16e\n", simplex::vector_norm2<i_t, f_t>(w)); }
if (debug) { settings_.log.printf("||w|| = %.16e\n", vector_norm2<i_t, f_t>(w)); }
if (solve_status != 0) {
settings_.log.printf("Linear solve failed in Sherman Morrison after ADAT solve\n");
return solve_status;
Expand Down Expand Up @@ -1102,7 +1092,7 @@ class iteration_data_t {
matrix_vector_multiply(ADAT, 1.0, M_col, -1.0, M_residual);
settings_.log.printf(
"|| A_sparse * D_sparse * A_sparse^T * M(:, k) - AD_dense(:, k) ||_2 = %e\n",
simplex::vector_norm2<i_t, f_t>(M_residual));
vector_norm2<i_t, f_t>(M_residual));
}
}
// A_sparse * D_sparse * A_sparse^T * M = U = AD_dense
Expand Down Expand Up @@ -1154,8 +1144,7 @@ class iteration_data_t {
if (debug) {
dense_vector_t<i_t, f_t> H_residual = g;
H.matrix_vector_multiply(1.0, y, -1.0, H_residual);
settings_.log.printf("|| H * y - g ||_2 = %e\n",
simplex::vector_norm2<i_t, f_t>(H_residual));
settings_.log.printf("|| H * y - g ||_2 = %e\n", vector_norm2<i_t, f_t>(H_residual));
}

// x = (A_sparse * D_sparse * A_sparse^T)^{-1} * (b - U * y)
Expand All @@ -1174,15 +1163,14 @@ class iteration_data_t {
dense_vector_t<i_t, f_t> solve_residual = v;
matrix_vector_multiply(ADAT, 1.0, x, -1.0, solve_residual);
settings_.log.printf("|| A_sparse * D * A_sparse^T * x - v ||_2 = %e\n",
simplex::vector_norm2<i_t, f_t>(solve_residual));
vector_norm2<i_t, f_t>(solve_residual));
}

if (debug) {
// Check U^T * x - y = 0;
dense_vector_t<i_t, f_t> residual_2 = y;
AD_dense.transpose_multiply(1.0, x, -1.0, residual_2);
settings_.log.printf("|| U^T * x - y ||_2 = %e\n",
simplex::vector_norm2<i_t, f_t>(residual_2));
settings_.log.printf("|| U^T * x - y ||_2 = %e\n", vector_norm2<i_t, f_t>(residual_2));
}

if (debug) {
Expand All @@ -1191,7 +1179,7 @@ class iteration_data_t {
AD_dense.matrix_vector_multiply(1.0, y, -1.0, residual_1);
matrix_vector_multiply(ADAT, 1.0, x, 1.0, residual_1);
settings_.log.printf("|| A_sparse * D_sparse * A_sparse^T * x + U * y - b ||_2 = %e\n",
simplex::vector_norm2<i_t, f_t>(residual_1));
vector_norm2<i_t, f_t>(residual_1));
}

if (full_debug && debug) {
Expand All @@ -1216,7 +1204,7 @@ class iteration_data_t {

adat_multiply(-1.0, ei, 1.0, u);

max_error = std::max(max_error, simplex::vector_norm2<i_t, f_t>(u));
max_error = std::max(max_error, vector_norm2<i_t, f_t>(u));
}
settings_.log.printf("|| ADAT(e_i) - ADA^T * e_i ||_2 = %e\n", max_error);
}
Expand Down Expand Up @@ -1356,7 +1344,7 @@ class iteration_data_t {

adat_multiply(-1.0, ei, 1.0, u);

max_error = std::max(max_error, simplex::vector_norm2<i_t, f_t>(u));
max_error = std::max(max_error, vector_norm2<i_t, f_t>(u));
}
settings_.log.printf(
"|| (A_sparse * D_sparse * A_sparse^T + U * V^T) * e_i - ADA^T * e_i ||_2 = %e\n",
Expand All @@ -1367,7 +1355,7 @@ class iteration_data_t {
dense_vector_t<i_t, f_t> total_residual = b;
adat_multiply(1.0, x, -1.0, total_residual);
settings_.log.printf("|| A * D * A^T * x - b ||_2 = %e\n",
simplex::vector_norm2<i_t, f_t>(total_residual));
vector_norm2<i_t, f_t>(total_residual));
}

// Now do some rounds of PCG
Expand Down Expand Up @@ -1436,7 +1424,7 @@ class iteration_data_t {
dense_vector_t<i_t, f_t> dual_res = z_tilde;
dual_res.axpy(-1.0, lp.objective, 1.0);
cusparse_view.transpose_spmv(1.0, solution.y, 1.0, dual_res);
f_t dual_residual_norm = simplex::vector_norm_inf<i_t, f_t>(dual_res, stream_view_);
f_t dual_residual_norm = vector_norm_inf<i_t, f_t>(dual_res, stream_view_);
#ifdef PRINT_INFO
settings_.log.printf("Solution Dual residual: %e\n", dual_residual_norm);
#endif
Expand Down Expand Up @@ -1794,20 +1782,20 @@ class iteration_data_t {

// u = A^T * y
dense_vector_t<i_t, f_t> u(n);
simplex::matrix_transpose_vector_multiply(A, 1.0, y, 0.0, u);
if (debug) { printf("||u|| = %.16e\n", simplex::vector_norm2<i_t, f_t>(u)); }
matrix_transpose_vector_multiply(A, 1.0, y, 0.0, u);
if (debug) { printf("||u|| = %.16e\n", vector_norm2<i_t, f_t>(u)); }

// w = Dinv * u
dense_vector_t<i_t, f_t> w(n);
inv_diag.pairwise_product(u, w);
if (debug) { printf("||inv_diag|| = %.16e\n", simplex::vector_norm2<i_t, f_t>(inv_diag)); }
if (debug) { printf("||inv_diag|| = %.16e\n", vector_norm2<i_t, f_t>(inv_diag)); }

// v = alpha * A * w + beta * v = alpha * A * Dinv * A^T * y + beta * v
matrix_vector_multiply(A, alpha, w, beta, v);
if (debug) {
printf("||A|| = %.16e\n", simplex::vector_norm2<i_t, f_t>(A.x));
printf("||w|| = %.16e\n", simplex::vector_norm2<i_t, f_t>(w));
printf("||v|| = %.16e\n", simplex::vector_norm2<i_t, f_t>(v));
printf("||A|| = %.16e\n", vector_norm2<i_t, f_t>(A.x));
printf("||w|| = %.16e\n", vector_norm2<i_t, f_t>(w));
printf("||v|| = %.16e\n", vector_norm2<i_t, f_t>(v));
}
}

Expand Down Expand Up @@ -2174,8 +2162,8 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
// LP block: e = 1, SOC block: e = (sqrt(2), 0, ..., 0)
if (data.has_cones()) {
const i_t cs = data.cone_start();
const f_t norm_b = simplex::vector_norm_inf<i_t, f_t>(lp.rhs);
const f_t norm_c = simplex::vector_norm_inf<i_t, f_t>(lp.objective);
const f_t norm_b = vector_norm_inf<i_t, f_t>(lp.rhs);
const f_t norm_c = vector_norm_inf<i_t, f_t>(lp.objective);
const f_t mu = std::sqrt((1.0 + norm_b) * (1.0 + norm_c));
const f_t sqrt2 = std::sqrt(2.0);
const f_t x_soc = mu * sqrt2;
Expand Down Expand Up @@ -2281,27 +2269,27 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
// rhs_x <- A * Dinv * F * u - b
data.cusparse_view_.spmv(1.0, DinvFu, -1.0, rhs_x);
#ifdef PRINT_INFO
settings.log.printf("||DinvFu|| = %e\n", simplex::vector_norm2<i_t, f_t>(DinvFu));
settings.log.printf("||DinvFu|| = %e\n", vector_norm2<i_t, f_t>(DinvFu));
#endif

// Solve A*Dinv*A'*q = A*Dinv*F*u - b
#ifdef PRINT_INFO
settings.log.printf("||rhs_x|| = %.16e\n", simplex::vector_norm2<i_t, f_t>(rhs_x));
settings.log.printf("||rhs_x|| = %.16e\n", vector_norm2<i_t, f_t>(rhs_x));
#endif
// i_t solve_status = data.chol->solve(rhs_x, q);
i_t solve_status = data.solve_adat(rhs_x, q);
if (solve_status != 0) { return status; }
#ifdef PRINT_INFO
settings.log.printf("Initial solve status %d\n", solve_status);
settings.log.printf("||q|| = %.16e\n", simplex::vector_norm2<i_t, f_t>(q));
settings.log.printf("||q|| = %.16e\n", vector_norm2<i_t, f_t>(q));
#endif

// rhs_x <- A*Dinv*A'*q - rhs_x
data.adat_multiply(1.0, q, -1.0, rhs_x);
// matrix_vector_multiply(data.ADAT, 1.0, q, -1.0, rhs_x);
#ifdef PRINT_INFO
settings.log.printf("|| A*Dinv*A'*q - (A*Dinv*F*u - b) || = %.16e\n",
simplex::vector_norm2<i_t, f_t>(rhs_x));
vector_norm2<i_t, f_t>(rhs_x));
#endif

// x = Dinv*(F*u - A'*q)
Expand All @@ -2327,8 +2315,7 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
data.cusparse_view_.spmv(1.0, data.x, -1.0, init_primal_residual);
data.handle_ptr->get_stream().synchronize();
#ifdef PRINT_INFO
settings.log.printf("||b - A * x||: %.16e\n",
simplex::vector_norm2<i_t, f_t>(init_primal_residual));
settings.log.printf("||b - A * x||: %.16e\n", vector_norm2<i_t, f_t>(init_primal_residual));
#endif

if (data.n_upper_bounds > 0) {
Expand All @@ -2338,8 +2325,7 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
init_bound_residual[k] = lp.upper[j] - data.w[k] - data.x[j];
}
#ifdef PRINT_INFO
settings.log.printf("|| u - w - x||: %e\n",
simplex::vector_norm2<i_t, f_t>(init_bound_residual));
settings.log.printf("|| u - w - x||: %e\n", vector_norm2<i_t, f_t>(init_bound_residual));
#endif
}

Expand Down Expand Up @@ -2453,7 +2439,7 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
}
#ifdef PRINT_INFO
settings.log.printf("||A^T y + z - E*v - Q*x - c ||: %e\n",
simplex::vector_norm2<i_t, f_t>(init_dual_residual));
vector_norm2<i_t, f_t>(init_dual_residual));
#endif
// Make sure (w, x, v, z) > 0. Skip free variables being handled directly.
data.w.ensure_positive(epsilon_adjust);
Expand Down Expand Up @@ -3017,7 +3003,7 @@ i_t barrier_solver_t<i_t, f_t>::gpu_compute_search_direction(iteration_data_t<i_
raft::common::nvtx::range fun_scope("Barrier: dx_residual_2 GPU");

// norm_inf(D^-1 * (A'*dy - r1) - dx)
const f_t dx_residual_2_norm = simplex::device_custom_vector_norm_inf<i_t, f_t>(
const f_t dx_residual_2_norm = device_custom_vector_norm_inf<i_t, f_t>(
thrust::make_transform_iterator(
thrust::make_zip_iterator(data.d_inv_diag.data(), data.d_r1_.data(), data.d_dx_.data()),
[] HD(thrust::tuple<f_t, f_t, f_t> t) -> f_t {
Expand Down Expand Up @@ -3096,7 +3082,7 @@ i_t barrier_solver_t<i_t, f_t>::gpu_compute_search_direction(iteration_data_t<i_
lp.A.transpose(Atranspose);
multiply(ADinv, Atranspose, ADinvAT);
matrix_vector_multiply(ADinvAT, 1.0, dy, -1.0, dx_residual_4);
const f_t dx_residual_4_norm = simplex::vector_norm_inf<i_t, f_t>(dx_residual_4, stream_view_);
const f_t dx_residual_4_norm = vector_norm_inf<i_t, f_t>(dx_residual_4, stream_view_);
max_residual = std::max(max_residual, dx_residual_4_norm);
if (dx_residual_4_norm > 1e-2) {
settings.log.printf("|| ADAT * dy - A * D^-1 * r1 - A * dx || = %.2e\n", dx_residual_4_norm);
Expand Down Expand Up @@ -4127,8 +4113,8 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time, lp_solution_t<i_t,
f_t mu;
compute_mu(data, mu);

f_t norm_b = simplex::vector_norm_inf<i_t, f_t>(data.b, stream_view_);
f_t norm_c = simplex::vector_norm_inf<i_t, f_t>(data.c, stream_view_);
f_t norm_b = vector_norm_inf<i_t, f_t>(data.b, stream_view_);
f_t norm_c = vector_norm_inf<i_t, f_t>(data.c, stream_view_);

f_t quad_objective = 0.0;
if (data.Q.n > 0) {
Expand Down
8 changes: 4 additions & 4 deletions cpp/src/barrier/barrier.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,14 +6,14 @@
/* clang-format on */
#pragma once

#include <barrier/dense_vector.hpp>
#include <linear_algebra/dense_vector.hpp>

#include <dual_simplex/presolve.hpp>
#include <dual_simplex/simplex_solver_settings.hpp>
#include <dual_simplex/solution.hpp>
#include <dual_simplex/solve.hpp>
#include <dual_simplex/sparse_matrix.hpp>
#include <dual_simplex/tic_toc.hpp>
#include <linear_algebra/sparse_matrix.hpp>
#include <math_optimization/tic_toc.hpp>

#include <rmm/device_uvector.hpp>
namespace cuopt::mathematical_optimization::barrier {
Expand All @@ -36,7 +36,7 @@ class barrier_solver_t {

private:
void my_pop_range(bool debug) const;
void create_Q(const simplex::lp_problem_t<i_t, f_t>& lp, simplex::csc_matrix_t<i_t, f_t>& Q);
void create_Q(const simplex::lp_problem_t<i_t, f_t>& lp, csc_matrix_t<i_t, f_t>& Q);
int initial_point(iteration_data_t<i_t, f_t>& data);
void compute_residual_norms(const dense_vector_t<i_t, f_t>& w,
const dense_vector_t<i_t, f_t>& x,
Expand Down
18 changes: 9 additions & 9 deletions cpp/src/barrier/conjugate_gradient.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,11 +6,11 @@
/* clang-format on */
#pragma once

#include <barrier/dense_vector.hpp>
#include <linear_algebra/dense_vector.hpp>

#include <dual_simplex/simplex_solver_settings.hpp>
#include <dual_simplex/types.hpp>
#include <dual_simplex/vector_math.hpp>
#include <linear_algebra/vector_math.hpp>
#include <math_optimization/types.hpp>

#include <cmath>
#include <cstdio>
Expand Down Expand Up @@ -42,12 +42,12 @@ i_t preconditioned_conjugate_gradient(const T& op,

dense_vector_t<i_t, f_t> Ap(b.size());
i_t iter = 0;
f_t norm_residual = simplex::vector_norm2<i_t, f_t>(residual);
f_t norm_residual = vector_norm2<i_t, f_t>(residual);
f_t initial_norm_residual = norm_residual;
if (show_pcg_info) {
settings.log.printf("PCG initial residual 2-norm %e inf-norm %e\n",
norm_residual,
simplex::vector_norm_inf<i_t, f_t>(residual));
vector_norm_inf<i_t, f_t>(residual));
}

f_t rTy = residual.inner_product(y);
Expand All @@ -62,7 +62,7 @@ i_t preconditioned_conjugate_gradient(const T& op,
// Update residual = residual + alpha * Ap
residual.axpy(alpha, Ap, 1.0);

f_t new_residual = simplex::vector_norm2<i_t, f_t>(residual);
f_t new_residual = vector_norm2<i_t, f_t>(residual);
if (new_residual > 1.1 * norm_residual || new_residual > 1.1 * initial_norm_residual) {
if (show_pcg_info) {
settings.log.printf(
Expand All @@ -78,7 +78,7 @@ i_t preconditioned_conjugate_gradient(const T& op,
// residual = A*x - b
residual = b;
op.a_multiply(1.0, x, -1.0, residual);
norm_residual = simplex::vector_norm2<i_t, f_t>(residual);
norm_residual = vector_norm2<i_t, f_t>(residual);

// Solve M y = r for y
op.m_solve(residual, y);
Expand All @@ -98,13 +98,13 @@ i_t preconditioned_conjugate_gradient(const T& op,
settings.log.printf("PCG iter %3d 2-norm_residual %.2e inf-norm_residual %.2e\n",
iter,
norm_residual,
simplex::vector_norm_inf<i_t, f_t>(residual));
vector_norm_inf<i_t, f_t>(residual));
}
}

residual = b;
op.a_multiply(1.0, x, -1.0, residual);
norm_residual = simplex::vector_norm2<i_t, f_t>(residual);
norm_residual = vector_norm2<i_t, f_t>(residual);
if (norm_residual < initial_norm_residual) {
if (show_pcg_info) {
settings.log.printf("PCG improved residual 2-norm %.2e/%.2e in %d iterations\n",
Expand Down
Loading