refactor: promote shared primitives to mathematical_optimization (#1481)

Moves linear algebra primitives and other shared utilities into the optimization namespace to reduce the need for `using` statements or namespace qualification introduced by #1446. Linear algebra goes into a new directory named `cpp/src/linear_algebra` and other utilities go into the pre-existing `cpp/src/math_optimization` directory.

No behavioral or API change; pure relocation.

Authors:
  - Miles Lubin (https://github.com/mlubin)

Approvers:
  - Ramakrishna Prabhu (https://github.com/ramakrishnap-nv)
  - Akif ÇÖRDÜK (https://github.com/akifcorduk)

URL: #1481

Author: mlubin

cpp/src/CMakeLists.txt

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

cpp/src/barrier/barrier.cu

@@ -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>
@@ -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,
@@ -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
@@ -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;
@@ -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
@@ -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)
@@ -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) {
@@ -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) {
@@ -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);
       }
@@ -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",
@@ -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
@@ -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
@@ -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));
     }
   }
 
@@ -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;
@@ -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)
@@ -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) {
@@ -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
   }
 
@@ -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);
@@ -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 {
@@ -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);
@@ -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) {

cpp/src/barrier/barrier.hpp

@@ -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 {
@@ -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,

cpp/src/barrier/conjugate_gradient.hpp

@@ -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>
@@ -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);
@@ -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(
@@ -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);
@@ -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",