diff --git a/cpp/src/CMakeLists.txt b/cpp/src/CMakeLists.txt
index 1ae6988466..928f1db76c 100644
--- a/cpp/src/CMakeLists.txt
+++ b/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)
diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu
index a1f493823d..784f6c0901 100644
--- a/cpp/src/barrier/barrier.cu
+++ b/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,19 +2269,19 @@ 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
@@ -2301,7 +2289,7 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
     // 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) {
diff --git a/cpp/src/barrier/barrier.hpp b/cpp/src/barrier/barrier.hpp
index bdf3ba3f83..07094a54a5 100644
--- a/cpp/src/barrier/barrier.hpp
+++ b/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,
diff --git a/cpp/src/barrier/conjugate_gradient.hpp b/cpp/src/barrier/conjugate_gradient.hpp
index 21fe33ba3d..39e6eadefc 100644
--- a/cpp/src/barrier/conjugate_gradient.hpp
+++ b/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",
diff --git a/cpp/src/barrier/cusparse_view.cu b/cpp/src/barrier/cusparse_view.cu
index 2585084097..637848e3eb 100644
--- a/cpp/src/barrier/cusparse_view.cu
+++ b/cpp/src/barrier/cusparse_view.cu
@@ -6,10 +6,10 @@
 /* clang-format on */
 
 #include <barrier/cusparse_view.hpp>
-#include <barrier/dense_vector.hpp>
 #include <barrier/pinned_host_allocator.hpp>
+#include <linear_algebra/dense_vector.hpp>
 
-#include <dual_simplex/sparse_matrix.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
 
 #include <utilities/copy_helpers.hpp>
 #include <utilities/macros.cuh>
@@ -126,7 +126,7 @@ static cusparseSpMVAlg_t get_spmv_alg(int num_rows)
 
 template <typename i_t, typename f_t>
 cusparse_view_t<i_t, f_t>::cusparse_view_t(raft::handle_t const* handle_ptr,
-                                           const simplex::csc_matrix_t<i_t, f_t>& A)
+                                           const csc_matrix_t<i_t, f_t>& A)
   : handle_ptr_(handle_ptr),
     A_offsets_(0, handle_ptr->get_stream()),
     A_indices_(0, handle_ptr->get_stream()),
@@ -146,7 +146,7 @@ cusparse_view_t<i_t, f_t>::cusparse_view_t(raft::handle_t const* handle_ptr,
   // TMP matrix data should already be on the GPU
   constexpr bool debug = false;
   if (debug) { printf("A hash: %zu\n", A.hash()); }
-  simplex::csr_matrix_t<i_t, f_t> A_csr(A.m, A.n, 1);
+  csr_matrix_t<i_t, f_t> A_csr(A.m, A.n, 1);
   A.to_compressed_row(A_csr);
   i_t rows                        = A_csr.m;
   i_t cols                        = A_csr.n;
diff --git a/cpp/src/barrier/cusparse_view.hpp b/cpp/src/barrier/cusparse_view.hpp
index 6af4813bc2..94aaa1ad42 100644
--- a/cpp/src/barrier/cusparse_view.hpp
+++ b/cpp/src/barrier/cusparse_view.hpp
@@ -6,7 +6,7 @@
 /* clang-format on */
 #pragma once
 
-#include <dual_simplex/sparse_matrix.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
 
 #include <pdlp/cusparse_view.hpp>
 
@@ -27,7 +27,7 @@ template <typename i_t, typename f_t>
 class cusparse_view_t {
  public:
   // TMP matrix data should already be on the GPU and in CSR not CSC
-  cusparse_view_t(raft::handle_t const* handle_ptr, const simplex::csc_matrix_t<i_t, f_t>& A);
+  cusparse_view_t(raft::handle_t const* handle_ptr, const csc_matrix_t<i_t, f_t>& A);
   ~cusparse_view_t();
 
   pdlp::cusparse_dn_vec_descr_wrapper_t<f_t> create_vector(rmm::device_uvector<f_t> const& vec);
diff --git a/cpp/src/barrier/device_sparse_matrix.cu b/cpp/src/barrier/device_sparse_matrix.cu
index 8016b3466f..29da786994 100644
--- a/cpp/src/barrier/device_sparse_matrix.cu
+++ b/cpp/src/barrier/device_sparse_matrix.cu
@@ -8,14 +8,14 @@
 #include <barrier/device_sparse_matrix.cuh>
 #include <barrier/pinned_host_allocator.hpp>
 
-#include <dual_simplex/sparse_matrix.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
 
 // This translation unit provides out-of-line definitions and explicit
-// instantiations of shared simplex sparse-matrix templates (csc_matrix_t,
+// instantiations of shared sparse-matrix templates (csc_matrix_t,
 // matrix_transpose_vector_multiply) specialized with barrier's
-// PinnedHostAllocator. They must live in the simplex namespace (where the
-// templates are declared), even though the file resides under barrier/.
-namespace cuopt::mathematical_optimization::simplex {
+// PinnedHostAllocator. They must live in the mathematical_optimization namespace
+// (where the templates are declared), even though the file resides under barrier/.
+namespace cuopt::mathematical_optimization {
 
 using cuopt::mathematical_optimization::barrier::PinnedHostAllocator;
 
@@ -72,4 +72,4 @@ template void csc_matrix_t<int, double>::scale_columns<PinnedHostAllocator<doubl
 
 #endif
 
-}  // namespace cuopt::mathematical_optimization::simplex
+}  // namespace cuopt::mathematical_optimization
diff --git a/cpp/src/barrier/device_sparse_matrix.cuh b/cpp/src/barrier/device_sparse_matrix.cuh
index 6f97db274b..90df0c86cf 100644
--- a/cpp/src/barrier/device_sparse_matrix.cuh
+++ b/cpp/src/barrier/device_sparse_matrix.cuh
@@ -7,8 +7,8 @@
 
 #pragma once
 
-#include <dual_simplex/sparse_matrix.hpp>
-#include <dual_simplex/types.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
+#include <math_optimization/types.hpp>
 
 #include <cub/cub.cuh>
 #include <rmm/device_scalar.hpp>
@@ -125,7 +125,7 @@ class device_csc_matrix_t {
   {
   }
 
-  device_csc_matrix_t(const simplex::csc_matrix_t<i_t, f_t>& A, rmm::cuda_stream_view stream)
+  device_csc_matrix_t(const csc_matrix_t<i_t, f_t>& A, rmm::cuda_stream_view stream)
     : m(A.m),
       n(A.n),
       nz_max(A.col_start[A.n]),
@@ -146,16 +146,16 @@ class device_csc_matrix_t {
     nz_max = nnz;
   }
 
-  simplex::csc_matrix_t<i_t, f_t> to_host(rmm::cuda_stream_view stream)
+  csc_matrix_t<i_t, f_t> to_host(rmm::cuda_stream_view stream)
   {
-    simplex::csc_matrix_t<i_t, f_t> A(m, n, nz_max);
+    csc_matrix_t<i_t, f_t> A(m, n, nz_max);
     A.col_start = cuopt::host_copy(col_start, stream);
     A.i         = cuopt::host_copy(i, stream);
     A.x         = cuopt::host_copy(x, stream);
     return A;
   }
 
-  void copy(simplex::csc_matrix_t<i_t, f_t>& A, rmm::cuda_stream_view stream)
+  void copy(csc_matrix_t<i_t, f_t>& A, rmm::cuda_stream_view stream)
   {
     m      = A.m;
     n      = A.n;
@@ -168,7 +168,8 @@ class device_csc_matrix_t {
     raft::copy(x.data(), A.x.data(), A.x.size(), stream);
   }
 
-  /** Same semantics as simplex::csc_matrix_t::to_compressed_row, entirely on device. */
+  /** Same semantics as csc_matrix_t::to_compressed_row, entirely on
+   * device. */
   void to_compressed_row(device_csr_matrix_t<i_t, f_t>& Arow, rmm::cuda_stream_view stream) const;
 
   void form_col_index(rmm::cuda_stream_view stream)
@@ -252,7 +253,7 @@ class device_csr_matrix_t {
   {
   }
 
-  device_csr_matrix_t(const simplex::csr_matrix_t<i_t, f_t>& A, rmm::cuda_stream_view stream)
+  device_csr_matrix_t(const csr_matrix_t<i_t, f_t>& A, rmm::cuda_stream_view stream)
     : m(A.m),
       n(A.n),
       nz_max(A.row_start[A.m]),
@@ -273,16 +274,16 @@ class device_csr_matrix_t {
     nz_max = nnz;
   }
 
-  simplex::csr_matrix_t<i_t, f_t> to_host(rmm::cuda_stream_view stream)
+  csr_matrix_t<i_t, f_t> to_host(rmm::cuda_stream_view stream)
   {
-    simplex::csr_matrix_t<i_t, f_t> A(m, n, nz_max);
+    csr_matrix_t<i_t, f_t> A(m, n, nz_max);
     A.row_start = cuopt::host_copy(row_start, stream);
     A.j         = cuopt::host_copy(j, stream);
     A.x         = cuopt::host_copy(x, stream);
     return A;
   }
 
-  void copy(simplex::csr_matrix_t<i_t, f_t>& A, rmm::cuda_stream_view stream)
+  void copy(csr_matrix_t<i_t, f_t>& A, rmm::cuda_stream_view stream)
   {
     m      = A.m;
     n      = A.n;
diff --git a/cpp/src/barrier/iterative_refinement.hpp b/cpp/src/barrier/iterative_refinement.hpp
index 93da00c933..7efe8b114d 100644
--- a/cpp/src/barrier/iterative_refinement.hpp
+++ b/cpp/src/barrier/iterative_refinement.hpp
@@ -6,11 +6,12 @@
 /* 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.cuh>
+#include <linear_algebra/vector_math.hpp>
+#include <math_optimization/types.hpp>
 
 #include <thrust/execution_policy.h>
 #include <thrust/extrema.h>
@@ -54,38 +55,6 @@ struct subtract_scaled_op {
   __host__ __device__ T operator()(T a, T b) const { return a - scale * b; }
 };
 
-template <typename f_t>
-f_t vector_norm_inf(const rmm::device_uvector<f_t>& x)
-{
-  auto begin   = x.data();
-  auto end     = x.data() + x.size();
-  auto max_abs = thrust::transform_reduce(
-    rmm::exec_policy(x.stream()),
-    begin,
-    end,
-    [] __host__ __device__(f_t val) { return abs(val); },
-    static_cast<f_t>(0),
-    thrust::maximum<f_t>{});
-  RAFT_CHECK_CUDA(x.stream());
-  return max_abs;
-}
-
-template <typename f_t>
-f_t vector_norm2(const rmm::device_uvector<f_t>& x)
-{
-  auto begin          = x.data();
-  auto end            = x.data() + x.size();
-  auto sum_of_squares = thrust::transform_reduce(
-    rmm::exec_policy(x.stream()),
-    begin,
-    end,
-    [] __host__ __device__(f_t val) { return val * val; },
-    f_t(0),
-    thrust::plus<f_t>{});
-  RAFT_CHECK_CUDA(x.stream());
-  return std::sqrt(sum_of_squares);
-}
-
 template <typename i_t, typename f_t, typename T>
 f_t iterative_refinement_simple(T& op,
                                 const rmm::device_uvector<f_t>& b,
diff --git a/cpp/src/barrier/pinned_host_allocator.cu b/cpp/src/barrier/pinned_host_allocator.cu
index 797d679e2b..619042d2d0 100644
--- a/cpp/src/barrier/pinned_host_allocator.cu
+++ b/cpp/src/barrier/pinned_host_allocator.cu
@@ -7,8 +7,8 @@
 
 #include <raft/util/cuda_utils.cuh>
 
-#include <dual_simplex/types.hpp>
-#include <dual_simplex/vector_math.hpp>
+#include <linear_algebra/vector_math.hpp>
+#include <math_optimization/types.hpp>
 
 namespace cuopt::mathematical_optimization::barrier {
 
@@ -55,10 +55,10 @@ template class PinnedHostAllocator<int>;
 
 }  // namespace cuopt::mathematical_optimization::barrier
 
-namespace cuopt::mathematical_optimization::simplex {
+namespace cuopt::mathematical_optimization {
 
-// Explicit instantiation of the shared simplex vector_math template with
-// barrier's PinnedHostAllocator must live in simplex (the template's namespace).
+// Explicit instantiation of the shared vector_math template with barrier's
+// PinnedHostAllocator must live in mathematical_optimization (the template's namespace).
 #ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE
 template double
 vector_norm_inf<int,
@@ -68,4 +68,4 @@ vector_norm_inf<int,
     x);
 #endif
 
-}  // namespace cuopt::mathematical_optimization::simplex
+}  // namespace cuopt::mathematical_optimization
diff --git a/cpp/src/barrier/pinned_host_allocator.hpp b/cpp/src/barrier/pinned_host_allocator.hpp
index 7ab452cc10..f5701db881 100644
--- a/cpp/src/barrier/pinned_host_allocator.hpp
+++ b/cpp/src/barrier/pinned_host_allocator.hpp
@@ -7,7 +7,7 @@
 
 #pragma once
 
-#include <barrier/dense_vector.hpp>
+#include <linear_algebra/dense_vector.hpp>
 
 namespace cuopt::mathematical_optimization::barrier {
 
diff --git a/cpp/src/barrier/sparse_cholesky.cuh b/cpp/src/barrier/sparse_cholesky.cuh
index 29c1c616bb..0bac84af83 100644
--- a/cpp/src/barrier/sparse_cholesky.cuh
+++ b/cpp/src/barrier/sparse_cholesky.cuh
@@ -6,12 +6,12 @@
 /* clang-format on */
 #pragma once
 
-#include <barrier/dense_vector.hpp>
 #include <barrier/device_sparse_matrix.cuh>
+#include <linear_algebra/dense_vector.hpp>
 
 #include <dual_simplex/simplex_solver_settings.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 <cuda_runtime.h>
 #include <utilities/driver_helpers.cuh>
@@ -26,8 +26,8 @@ template <typename i_t, typename f_t>
 class sparse_cholesky_base_t {
  public:
   virtual ~sparse_cholesky_base_t()                                                 = default;
-  virtual i_t analyze(const simplex::csc_matrix_t<i_t, f_t>& A_in)                  = 0;
-  virtual i_t factorize(const simplex::csc_matrix_t<i_t, f_t>& A_in)                = 0;
+  virtual i_t analyze(const csc_matrix_t<i_t, f_t>& A_in)                           = 0;
+  virtual i_t factorize(const csc_matrix_t<i_t, f_t>& A_in)                         = 0;
   virtual i_t analyze(device_csr_matrix_t<i_t, f_t>& A_in)                          = 0;
   virtual i_t factorize(device_csr_matrix_t<i_t, f_t>& A_in)                        = 0;
   virtual i_t solve(const dense_vector_t<i_t, f_t>& b, dense_vector_t<i_t, f_t>& x) = 0;
@@ -390,8 +390,8 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
 
 #ifdef WRITE_MATRIX_MARKET
     {
-      simplex::csr_matrix_t<i_t, f_t> Arow_host = Arow.to_host(Arow.row_start.stream());
-      simplex::csc_matrix_t<i_t, f_t> A_col(Arow_host.m, Arow_host.n, 1);
+      csr_matrix_t<i_t, f_t> Arow_host = Arow.to_host(Arow.row_start.stream());
+      csc_matrix_t<i_t, f_t> A_col(Arow_host.m, Arow_host.n, 1);
       Arow_host.to_compressed_col(A_col);
       FILE* fid = fopen("A_to_factorize.mtx", "w");
       settings_.log.printf("writing matrix matrix\n");
@@ -473,7 +473,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
     }
 
     // Perform symbolic analysis
-    f_t start_symbolic = simplex::tic();
+    f_t start_symbolic = tic();
     f_t start_symbolic_factor;
 
     {
@@ -490,9 +490,9 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
           status);
         return -1;
       }
-      f_t reordering_time = simplex::toc(start_symbolic);
+      f_t reordering_time = toc(start_symbolic);
       settings_.log.printf("Reordering time             : %.2fs\n", reordering_time);
-      start_symbolic_factor = simplex::tic();
+      start_symbolic_factor = tic();
 
       status = cudssExecute(
         handle, CUDSS_PHASE_SYMBOLIC_FACTORIZATION, solverConfig, solverData, A, cudss_x, cudss_b);
@@ -508,9 +508,9 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
       }
     }
     RAFT_CUDA_TRY(cudaStreamSynchronize(stream));
-    f_t symbolic_factorization_time = simplex::toc(start_symbolic_factor);
+    f_t symbolic_factorization_time = toc(start_symbolic_factor);
     settings_.log.printf("Symbolic factorization time : %.2fs\n", symbolic_factorization_time);
-    settings_.log.printf("Total symbolic time         : %.2fs\n", simplex::toc(start_symbolic));
+    settings_.log.printf("Total symbolic time         : %.2fs\n", toc(start_symbolic));
     int64_t lu_nz       = 0;
     size_t size_written = 0;
     CUDSS_CALL_AND_CHECK(
@@ -532,8 +532,8 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
 // #define PRINT_MATRIX_NORM
 #ifdef PRINT_MATRIX_NORM
     cudaStreamSynchronize(stream);
-    simplex::csr_matrix_t<i_t, f_t> Arow_host = Arow.to_host(Arow.row_start.stream());
-    simplex::csc_matrix_t<i_t, f_t> A_col(Arow_host.m, Arow_host.n, 1);
+    csr_matrix_t<i_t, f_t> Arow_host = Arow.to_host(Arow.row_start.stream());
+    csc_matrix_t<i_t, f_t> A_col(Arow_host.m, Arow_host.n, 1);
     Arow_host.to_compressed_col(A_col);
     settings_.log.printf(
       "before factorize || A to factor|| = %.16e hash: %zu\n", A_col.norm1(), A_col.hash());
@@ -551,7 +551,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
     CUDSS_CALL_AND_CHECK(
       cudssMatrixSetValues(A, Arow.x.data()), status, "cudssMatrixSetValues for A");
 
-    f_t start_numeric = simplex::tic();
+    f_t start_numeric = tic();
     status            = cudssExecute(
       handle, CUDSS_PHASE_FACTORIZATION, solverConfig, solverData, A, cudss_x, cudss_b);
     if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) {
@@ -569,7 +569,7 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
     RAFT_CUDA_TRY(cudaStreamSynchronize(stream));
 #endif
 
-    f_t numeric_time = simplex::toc(start_numeric);
+    f_t numeric_time = toc(start_numeric);
     if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) {
       return CONCURRENT_HALT_RETURN;
     }
@@ -599,9 +599,9 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
     return 0;
   }
 
-  i_t analyze(const simplex::csc_matrix_t<i_t, f_t>& A_in) override
+  i_t analyze(const csc_matrix_t<i_t, f_t>& A_in) override
   {
-    simplex::csr_matrix_t<i_t, f_t> Arow(A_in.n, A_in.m, A_in.col_start[A_in.n]);
+    csr_matrix_t<i_t, f_t> Arow(A_in.n, A_in.m, A_in.col_start[A_in.n]);
 
 #ifdef WRITE_MATRIX_MARKET
     FILE* fid = fopen("A.mtx", "w");
@@ -692,18 +692,18 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
     if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) {
       return CONCURRENT_HALT_RETURN;
     }
-    f_t start_analysis = simplex::tic();
+    f_t start_analysis = tic();
     CUDSS_CALL_AND_CHECK(
       cudssExecute(handle, CUDSS_PHASE_REORDERING, solverConfig, solverData, A, cudss_x, cudss_b),
       status,
       "cudssExecute for reordering");
 
-    f_t reorder_time = simplex::toc(start_analysis);
+    f_t reorder_time = toc(start_analysis);
     if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) {
       return CONCURRENT_HALT_RETURN;
     }
 
-    f_t start_symbolic = simplex::tic();
+    f_t start_symbolic = tic();
 
     CUDSS_CALL_AND_CHECK(
       cudssExecute(
@@ -711,8 +711,8 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
       status,
       "cudssExecute for symbolic factorization");
 
-    f_t symbolic_time = simplex::toc(start_symbolic);
-    f_t analysis_time = simplex::toc(start_analysis);
+    f_t symbolic_time = toc(start_symbolic);
+    f_t analysis_time = toc(start_analysis);
     settings_.log.printf("Symbolic factorization time : %.2fs\n", symbolic_time);
     if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) {
       RAFT_CUDA_TRY(cudaStreamSynchronize(stream));
@@ -733,9 +733,9 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
 
     return 0;
   }
-  i_t factorize(const simplex::csc_matrix_t<i_t, f_t>& A_in) override
+  i_t factorize(const csc_matrix_t<i_t, f_t>& A_in) override
   {
-    simplex::csr_matrix_t<i_t, f_t> Arow(A_in.n, A_in.m, A_in.col_start[A_in.n]);
+    csr_matrix_t<i_t, f_t> Arow(A_in.n, A_in.m, A_in.col_start[A_in.n]);
     A_in.to_compressed_row(Arow);
 
     if (A_in.n != n) { settings_.log.printf("Error A in n %d != size %d\n", A_in.n, n); }
@@ -757,14 +757,14 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
     CUDSS_CALL_AND_CHECK(
       cudssMatrixSetValues(A, csr_values_d), status, "cudssMatrixSetValues for A");
 
-    f_t start_numeric = simplex::tic();
+    f_t start_numeric = tic();
     CUDSS_CALL_AND_CHECK(
       cudssExecute(
         handle, CUDSS_PHASE_FACTORIZATION, solverConfig, solverData, A, cudss_x, cudss_b),
       status,
       "cudssExecute for factorization");
 
-    f_t numeric_time = simplex::toc(start_numeric);
+    f_t numeric_time = toc(start_numeric);
     if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) {
       return CONCURRENT_HALT_RETURN;
     }
@@ -851,9 +851,9 @@ class sparse_cholesky_cudss_t : public sparse_cholesky_base_t<i_t, f_t> {
     raft::copy(x_host.data(), x.data(), n, stream);
     cudaStreamSynchronize(stream);
     settings_.log.printf("RHS norm %.16e, hash: %zu, Solution norm %.16e, hash: %zu\n",
-                         simplex::vector_norm2<i_t, f_t>(b_host),
+                         vector_norm2<i_t, f_t>(b_host),
                          compute_hash(b_host),
-                         simplex::vector_norm2<i_t, f_t>(x_host),
+                         vector_norm2<i_t, f_t>(x_host),
                          compute_hash(x_host));
 #endif
 
diff --git a/cpp/src/barrier/translate_soc.hpp b/cpp/src/barrier/translate_soc.hpp
index aac7da2025..3b288973bc 100644
--- a/cpp/src/barrier/translate_soc.hpp
+++ b/cpp/src/barrier/translate_soc.hpp
@@ -12,9 +12,9 @@
 
 #include <dual_simplex/right_looking_lu.hpp>
 #include <dual_simplex/solution.hpp>
-#include <dual_simplex/sparse_matrix.hpp>
-#include <dual_simplex/tic_toc.hpp>
 #include <dual_simplex/user_problem.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
+#include <math_optimization/tic_toc.hpp>
 
 #include <algorithm>
 #include <cmath>
@@ -56,7 +56,7 @@ void convert_quadratic_constraints_to_second_order_cones(
   i_t n,
   const std::vector<typename optimization_problem_interface_t<i_t, f_t>::quadratic_constraint_t>&
     qcs,
-  simplex::csr_matrix_t<i_t, f_t>& csr_A,
+  csr_matrix_t<i_t, f_t>& csr_A,
   simplex::user_problem_t<i_t, f_t>& user_problem)
 {
   cuopt_expects(!qcs.empty(),
@@ -624,7 +624,7 @@ void convert_quadratic_constraints_to_second_order_cones(
         }
       }
 
-      simplex::csc_matrix_t<i_t, f_t> H_csc(n_local, n_local, h_nnz);
+      csc_matrix_t<i_t, f_t> H_csc(n_local, n_local, h_nnz);
       {
         i_t p = 0;
         for (i_t j = 0; j < n_local; j++) {
@@ -644,10 +644,10 @@ void convert_quadratic_constraints_to_second_order_cones(
       // Step 2: Factorize H = P * L * D * L^T * P^T
       simplex::simplex_solver_settings_t<i_t, f_t> ldlt_settings;
       std::vector<i_t> ldlt_perm;
-      simplex::csc_matrix_t<i_t, f_t> L_factor(n, n, 1);
+      csc_matrix_t<i_t, f_t> L_factor(n, n, 1);
       std::vector<f_t> D_factor;
       f_t ldlt_work  = 0;
-      f_t ldlt_start = simplex::tic();
+      f_t ldlt_start = tic();
 
       i_t rank = simplex::right_looking_ldlt(
         H_csc, ldlt_settings, f_t(1e-12), ldlt_start, ldlt_perm, L_factor, D_factor, ldlt_work);
@@ -706,7 +706,7 @@ void convert_quadratic_constraints_to_second_order_cones(
       user_problem.row_sense.resize(m_before + n_new_rows);
       if (!user_problem.row_names.empty()) { user_problem.row_names.resize(m_before + n_new_rows); }
 
-      simplex::sparse_vector_t<i_t, f_t> eq_row;
+      sparse_vector_t<i_t, f_t> eq_row;
       eq_row.n = csr_A.n;
 
       // y-linking rows: y_k - sqrt(D[k]) * [row k of L^T P] * x = 0
@@ -837,7 +837,7 @@ void convert_quadratic_constraints_to_second_order_cones(
     if (!user_problem.row_names.empty()) { user_problem.row_names.resize(m_aug); }
 
     csr_A.n = std::max(csr_A.n, n_aug);
-    simplex::sparse_vector_t<i_t, f_t> eq_row;
+    sparse_vector_t<i_t, f_t> eq_row;
     eq_row.n = csr_A.n;
 
     for (size_t qc_i = 0; qc_i < qcs.size(); ++qc_i) {
@@ -945,7 +945,7 @@ void convert_quadratic_constraints_to_second_order_cones(
 
     csr_A.n = n_prob;
 
-    simplex::sparse_vector_t<i_t, f_t> eq_row;
+    sparse_vector_t<i_t, f_t> eq_row;
     size_t ri      = 0;
     i_t slack_base = n_old;
     i_t row_idx    = m_old;
@@ -1084,7 +1084,7 @@ void convert_quadratic_constraints_to_second_order_cones(
       if (!user_problem.row_names.empty()) { user_problem.row_names.resize(m_new); }
 
       csr_A.n = n_new;
-      simplex::sparse_vector_t<i_t, f_t> eq_row;
+      sparse_vector_t<i_t, f_t> eq_row;
       eq_row.n    = n_new;
       i_t row_idx = m_old;
       for (const auto& [alias, original] : cone_alias_pairs) {
@@ -1167,7 +1167,7 @@ void convert_quadratic_constraints_to_second_order_cones(
       if (!user_problem.row_names.empty()) { user_problem.row_names.resize(m_new); }
 
       csr_A.n = n_new;
-      simplex::sparse_vector_t<i_t, f_t> eq_row;
+      sparse_vector_t<i_t, f_t> eq_row;
       eq_row.n    = n_new;
       i_t row_idx = m_old;
       for (const auto& [alias, original] : bound_split_pairs) {
diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp
index a17bc6716c..d8859aea2f 100644
--- a/cpp/src/branch_and_bound/branch_and_bound.cpp
+++ b/cpp/src/branch_and_bound/branch_and_bound.cpp
@@ -26,8 +26,8 @@
 #include <dual_simplex/phase2.hpp>
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/random.hpp>
-#include <dual_simplex/tic_toc.hpp>
 #include <dual_simplex/user_problem.hpp>
+#include <math_optimization/tic_toc.hpp>
 
 #include <raft/core/nvtx.hpp>
 #include <utilities/circular_deque.hpp>
@@ -52,26 +52,20 @@ using simplex::compute_objective;
 using simplex::compute_user_objective;
 using simplex::crossover_status_t;
 using simplex::crush_primal_solution;
-using simplex::csr_matrix_t;
 using simplex::decompress_vstatus;
 using simplex::dual_phase2_with_advanced_basis;
 using simplex::dual_status_t;
-using simplex::inf;
 using simplex::logger_t;
 using simplex::lp_problem_t;
 using simplex::lp_solution_t;
 using simplex::lp_status_t;
-using simplex::matrix_vector_multiply;
 using simplex::mip_solution_t;
 using simplex::simplex_solver_settings_t;
 using simplex::solve_linear_program_with_advanced_basis;
-using simplex::tic;
-using simplex::toc;
 using simplex::uncrush_primal_solution;
 using simplex::user_problem_t;
 using simplex::variable_status_t;
 using simplex::variable_type_t;
-using simplex::vector_norm_inf;
 
 namespace {
 
diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp
index 18eb38de66..12c93fcd91 100644
--- a/cpp/src/branch_and_bound/branch_and_bound.hpp
+++ b/cpp/src/branch_and_bound/branch_and_bound.hpp
@@ -22,7 +22,7 @@
 #include <dual_simplex/simplex_solver_settings.hpp>
 #include <dual_simplex/solution.hpp>
 #include <dual_simplex/solve.hpp>
-#include <dual_simplex/types.hpp>
+#include <math_optimization/types.hpp>
 
 #include <utilities/macros.cuh>
 #include <utilities/omp_helpers.hpp>
@@ -173,7 +173,7 @@ class branch_and_bound_t {
   std::vector<f_t> guess_;
 
   // LP relaxation
-  simplex::csr_matrix_t<i_t, f_t> Arow_;
+  csr_matrix_t<i_t, f_t> Arow_;
   simplex::lp_problem_t<i_t, f_t> original_lp_;
   std::vector<i_t> new_slacks_;
   std::vector<simplex::variable_type_t> var_types_;
@@ -375,7 +375,7 @@ class branch_and_bound_t {
   // ============================================================================
 
   // Main deterministic coordinator loop
-  void run_deterministic_coordinator(const simplex::csr_matrix_t<i_t, f_t>& Arow);
+  void run_deterministic_coordinator(const csr_matrix_t<i_t, f_t>& Arow);
 
   // Gather all events generated, sort by WU timestamp, apply
   void deterministic_sort_replay_events(const bb_event_batch_t<i_t, f_t>& events);
diff --git a/cpp/src/branch_and_bound/deterministic_workers.hpp b/cpp/src/branch_and_bound/deterministic_workers.hpp
index f7c960cbdd..60d3436401 100644
--- a/cpp/src/branch_and_bound/deterministic_workers.hpp
+++ b/cpp/src/branch_and_bound/deterministic_workers.hpp
@@ -86,7 +86,7 @@ class deterministic_worker_base_t : public branch_and_bound_worker_t<i_t, f_t> {
 
   deterministic_worker_base_t(int id,
                               const simplex::lp_problem_t<i_t, f_t>& original_lp,
-                              const simplex::csr_matrix_t<i_t, f_t>& Arow,
+                              const csr_matrix_t<i_t, f_t>& Arow,
                               const std::vector<simplex::variable_type_t>& var_types,
                               const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
                               const std::string& context_name)
@@ -138,7 +138,7 @@ class deterministic_bfs_worker_t
 
   explicit deterministic_bfs_worker_t(int id,
                                       const simplex::lp_problem_t<i_t, f_t>& original_lp,
-                                      const simplex::csr_matrix_t<i_t, f_t>& Arow,
+                                      const csr_matrix_t<i_t, f_t>& Arow,
                                       const std::vector<simplex::variable_type_t>& var_types,
                                       const simplex::simplex_solver_settings_t<i_t, f_t>& settings)
     : base_t(id, original_lp, Arow, var_types, settings, "BB_Worker_" + std::to_string(id))
@@ -297,7 +297,7 @@ class deterministic_diving_worker_t
     int id,
     search_strategy_t type,
     const simplex::lp_problem_t<i_t, f_t>& original_lp,
-    const simplex::csr_matrix_t<i_t, f_t>& Arow,
+    const csr_matrix_t<i_t, f_t>& Arow,
     const std::vector<simplex::variable_type_t>& var_types,
     const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
     const std::vector<f_t>* root_sol)
@@ -405,7 +405,7 @@ class deterministic_bfs_worker_pool_t
  public:
   deterministic_bfs_worker_pool_t(int num_workers,
                                   const simplex::lp_problem_t<i_t, f_t>& original_lp,
-                                  const simplex::csr_matrix_t<i_t, f_t>& Arow,
+                                  const csr_matrix_t<i_t, f_t>& Arow,
                                   const std::vector<simplex::variable_type_t>& var_types,
                                   const simplex::simplex_solver_settings_t<i_t, f_t>& settings)
   {
@@ -440,7 +440,7 @@ class deterministic_diving_worker_pool_t
   deterministic_diving_worker_pool_t(int num_workers,
                                      const std::vector<search_strategy_t>& diving_types,
                                      const simplex::lp_problem_t<i_t, f_t>& original_lp,
-                                     const simplex::csr_matrix_t<i_t, f_t>& Arow,
+                                     const csr_matrix_t<i_t, f_t>& Arow,
                                      const std::vector<simplex::variable_type_t>& var_types,
                                      const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
                                      const std::vector<f_t>* root_solution)
diff --git a/cpp/src/branch_and_bound/diving_heuristics.cpp b/cpp/src/branch_and_bound/diving_heuristics.cpp
index aebb7998cc..82843a82f4 100644
--- a/cpp/src/branch_and_bound/diving_heuristics.cpp
+++ b/cpp/src/branch_and_bound/diving_heuristics.cpp
@@ -24,7 +24,7 @@ branch_variable_t<i_t> line_search_diving(const std::vector<i_t>& fractional,
   branch_direction_t round_dir = branch_direction_t::NONE;
 
   for (i_t j : fractional) {
-    f_t score              = simplex::inf;
+    f_t score              = inf;
     branch_direction_t dir = branch_direction_t::NONE;
 
     if (solution[j] < root_solution[j] - eps) {
diff --git a/cpp/src/branch_and_bound/mip_node.hpp b/cpp/src/branch_and_bound/mip_node.hpp
index fc7b3a90e9..c763e5ada5 100644
--- a/cpp/src/branch_and_bound/mip_node.hpp
+++ b/cpp/src/branch_and_bound/mip_node.hpp
@@ -10,7 +10,7 @@
 #include <branch_and_bound/constants.hpp>
 
 #include <dual_simplex/initial_basis.hpp>
-#include <dual_simplex/types.hpp>
+#include <math_optimization/types.hpp>
 
 #include <utilities/hashing.hpp>
 #include <utilities/omp_helpers.hpp>
diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp
index baea370fa4..82aca6c458 100644
--- a/cpp/src/branch_and_bound/pseudo_costs.cpp
+++ b/cpp/src/branch_and_bound/pseudo_costs.cpp
@@ -13,7 +13,7 @@
 #include <dual_simplex/phase2.hpp>
 #include <dual_simplex/simplex_solver_settings.hpp>
 #include <dual_simplex/solve.hpp>
-#include <dual_simplex/tic_toc.hpp>
+#include <math_optimization/tic_toc.hpp>
 
 #include <mip_heuristics/mip_constants.hpp>
 
@@ -30,19 +30,13 @@ namespace cuopt::mathematical_optimization::mip {
 using simplex::basis_update_mpf_t;
 using simplex::compute_initial_nonbasic_end;
 using simplex::compute_objective;
-using simplex::csc_matrix_t;
-using simplex::csr_matrix_t;
 using simplex::dual_status_t;
 using simplex::logger_t;
 using simplex::lp_problem_t;
 using simplex::lp_solution_t;
 using simplex::simplex_solver_settings_t;
-using simplex::sparse_vector_t;
-using simplex::tic;
-using simplex::toc;
 using simplex::variable_status_t;
 using simplex::variable_type_t;
-using simplex::vector_norm_inf;
 
 namespace {
 
@@ -65,7 +59,7 @@ f_t compute_step_length(const simplex_solver_settings_t<i_t, f_t>& settings,
                         const std::vector<f_t>& delta_z,
                         const std::vector<i_t>& delta_z_indices)
 {
-  f_t step_length = simplex::inf;
+  f_t step_length = inf;
   f_t pivot_tol   = settings.pivot_tol;
   const i_t nz    = delta_z_indices.size();
   for (i_t h = 0; h < nz; h++) {
@@ -153,7 +147,7 @@ objective_change_estimate_t<f_t> single_pivot_objective_change_estimate(
     for (i_t j = 0; j < lp.num_cols; j++) {
       dual_residual[j] -= lp.objective[j];
     }
-    simplex::matrix_transpose_vector_multiply(lp.A, 1.0, lp_solution.y, 1.0, dual_residual);
+    matrix_transpose_vector_multiply(lp.A, 1.0, lp_solution.y, 1.0, dual_residual);
     f_t dual_residual_norm = vector_norm_inf<i_t, f_t>(dual_residual);
     settings.log.printf("Dual residual norm: %e\n", dual_residual_norm);
   }
diff --git a/cpp/src/branch_and_bound/pseudo_costs.hpp b/cpp/src/branch_and_bound/pseudo_costs.hpp
index 75821b79b7..eddf2cf28f 100644
--- a/cpp/src/branch_and_bound/pseudo_costs.hpp
+++ b/cpp/src/branch_and_bound/pseudo_costs.hpp
@@ -13,7 +13,7 @@
 #include <dual_simplex/basis_updates.hpp>
 #include <dual_simplex/logger.hpp>
 #include <dual_simplex/simplex_solver_settings.hpp>
-#include <dual_simplex/types.hpp>
+#include <math_optimization/types.hpp>
 
 #include <utilities/omp_helpers.hpp>
 #include <utilities/pcgenerator.hpp>
@@ -213,12 +213,12 @@ class pseudo_costs_t {
                                  f_t avg_down,
                                  f_t avg_up) const;
 
-  std::shared_ptr<simplex::csc_matrix_t<i_t, f_t>> AT;  // Transpose of the constraint matrix A
+  std::shared_ptr<csc_matrix_t<i_t, f_t>> AT;  // Transpose of the constraint matrix A
   std::shared_ptr<batch_pdlp_warm_cache_t<i_t, f_t>> pdlp_warm_cache;
 
   reliability_branching_settings_t<i_t, f_t> reliability_branching_settings;
   simplex::simplex_solver_settings_t<i_t, f_t> settings;
-  simplex::csr_matrix_t<i_t, f_t> Arow;
+  csr_matrix_t<i_t, f_t> Arow;
 
  protected:
   std::vector<omp_atomic_t<f_t>> pseudo_cost_sum_up;
diff --git a/cpp/src/branch_and_bound/symmetry.hpp b/cpp/src/branch_and_bound/symmetry.hpp
index 1848b3dad8..f8e9ce469b 100644
--- a/cpp/src/branch_and_bound/symmetry.hpp
+++ b/cpp/src/branch_and_bound/symmetry.hpp
@@ -10,9 +10,9 @@
 #include <branch_and_bound/mip_node.hpp>
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/simplex_solver_settings.hpp>
-#include <dual_simplex/tic_toc.hpp>
-#include <dual_simplex/types.hpp>
 #include <dual_simplex/user_problem.hpp>
+#include <math_optimization/tic_toc.hpp>
+#include <math_optimization/types.hpp>
 
 #include "dejavu.h"
 
@@ -688,7 +688,7 @@ std::unique_ptr<mip_symmetry_t<i_t, f_t>> detect_symmetry(
 
   has_symmetry = false;
 
-  f_t start_time = simplex::tic();
+  f_t start_time = tic();
   simplex::lp_problem_t<i_t, f_t> problem(user_problem.handle_ptr, 1, 1, 1);
   std::vector<i_t> new_slacks;
   simplex::dualize_info_t<i_t, f_t> dualize_info;
@@ -837,7 +837,7 @@ std::unique_ptr<mip_symmetry_t<i_t, f_t>> detect_symmetry(
   // Let r_i be the vertex associated with the row i
   // Let V_i,c be the set of variables in row i with the same color c
   // We create a new vertex w_i,c and edges (v_j, w_i,c) and (w_i,c, r_i) for all v_j in V_i,c
-  simplex::csr_matrix_t<i_t, f_t> A_row(problem.num_rows, problem.num_cols, 0);
+  csr_matrix_t<i_t, f_t> A_row(problem.num_rows, problem.num_cols, 0);
   problem.A.to_compressed_row(A_row);
 
   std::vector<f_t> nonzeros = A_row.x;
@@ -939,8 +939,8 @@ std::unique_ptr<mip_symmetry_t<i_t, f_t>> detect_symmetry(
 
 #endif
 
-  settings.log.printf("Graph construction time %f\n", simplex::toc(start_time));
-  f_t dejavu_start_time = simplex::tic();
+  settings.log.printf("Graph construction time %f\n", toc(start_time));
+  f_t dejavu_start_time = tic();
 
   // The graph should now be described by:
   // vertices, edge_in, edge_out, vertex_colors
@@ -1062,7 +1062,7 @@ std::unique_ptr<mip_symmetry_t<i_t, f_t>> detect_symmetry(
                       grp_size_str.str().c_str(),
                       num_dejavu_generators,
                       projected_count);
-  settings.log.printf("Dejavu time %f\n", simplex::toc(dejavu_start_time));
+  settings.log.printf("Dejavu time %f\n", toc(dejavu_start_time));
 
   result->num_generators = result->generators.num_generators();
   if (projected_count > static_cast<int>(result->num_generators)) {
@@ -1116,7 +1116,7 @@ std::unique_ptr<mip_symmetry_t<i_t, f_t>> detect_symmetry(
                    (total_vars_in_orbits >= 10);
   }
 
-  settings.log.printf("Total symmetry detection time %f\n", simplex::toc(start_time));
+  settings.log.printf("Total symmetry detection time %f\n", toc(start_time));
 
   if (!has_symmetry) {
     settings.log.printf(
diff --git a/cpp/src/branch_and_bound/worker.hpp b/cpp/src/branch_and_bound/worker.hpp
index 758bbe16d8..4fd278417e 100644
--- a/cpp/src/branch_and_bound/worker.hpp
+++ b/cpp/src/branch_and_bound/worker.hpp
@@ -91,7 +91,7 @@ class branch_and_bound_worker_t {
 
   branch_and_bound_worker_t(i_t worker_id,
                             const simplex::lp_problem_t<i_t, f_t>& original_lp,
-                            const simplex::csr_matrix_t<i_t, f_t>& Arow,
+                            const csr_matrix_t<i_t, f_t>& Arow,
                             const std::vector<simplex::variable_type_t>& var_type,
                             const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
                             uint64_t rng_offset = 0)
@@ -143,7 +143,7 @@ class bfs_worker_t : public branch_and_bound_worker_t<i_t, f_t> {
   using Base = branch_and_bound_worker_t<i_t, f_t>;
   bfs_worker_t(i_t worker_id,
                const simplex::lp_problem_t<i_t, f_t>& original_lp,
-               const simplex::csr_matrix_t<i_t, f_t>& Arow,
+               const csr_matrix_t<i_t, f_t>& Arow,
                const std::vector<simplex::variable_type_t>& var_type,
                const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
                uint64_t rng_offset = 0)
diff --git a/cpp/src/branch_and_bound/worker_pool.hpp b/cpp/src/branch_and_bound/worker_pool.hpp
index 2189fc5b2a..e9d55bafe3 100644
--- a/cpp/src/branch_and_bound/worker_pool.hpp
+++ b/cpp/src/branch_and_bound/worker_pool.hpp
@@ -20,7 +20,7 @@ class worker_pool_t {
 
   void init(i_t num_workers,
             const simplex::lp_problem_t<i_t, f_t>& original_lp,
-            const simplex::csr_matrix_t<i_t, f_t>& Arow,
+            const csr_matrix_t<i_t, f_t>& Arow,
             const std::vector<simplex::variable_type_t>& var_type,
             mip_symmetry_t<i_t, f_t>* symmetry,
             const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp
index 40efbe6185..39d75fcc97 100644
--- a/cpp/src/cuts/cuts.cpp
+++ b/cpp/src/cuts/cuts.cpp
@@ -9,7 +9,7 @@
 #include <cuts/rational.hpp>
 
 #include <dual_simplex/basis_solves.hpp>
-#include <dual_simplex/tic_toc.hpp>
+#include <math_optimization/tic_toc.hpp>
 #include <mip_heuristics/presolve/conflict_graph/clique_table.cuh>
 #include <utilities/macros.cuh>
 
@@ -22,7 +22,7 @@
 #include <tuple>
 #include <unordered_set>
 
-#include <barrier/dense_matrix.hpp>
+#include <linear_algebra/dense_matrix.hpp>
 
 #include <numeric>
 #include <queue>
@@ -30,22 +30,12 @@
 namespace cuopt::mathematical_optimization::mip {
 
 using simplex::basis_update_mpf_t;
-using simplex::csc_matrix_t;
-using simplex::csr_matrix_t;
 using simplex::form_b;
-using simplex::inf;
 using simplex::lp_problem_t;
 using simplex::lp_solution_t;
-using simplex::matrix_vector_multiply;
 using simplex::simplex_solver_settings_t;
-using simplex::sparse_dot;
-using simplex::sparse_vector_t;
-using simplex::tic;
-using simplex::toc;
 using simplex::variable_status_t;
 using simplex::variable_type_t;
-using simplex::vector_norm2;
-using simplex::vector_norm_inf;
 
 namespace {
 
@@ -2349,8 +2339,8 @@ f_t knapsack_generation_t<i_t, f_t>::solve_knapsack_problem(const std::vector<f_
   solution.assign(n, 0.0);
 
   // dp(j, v) = minimum weight using first j items to get value v
-  barrier::dense_matrix_t<i_t, i_t> dp(n + 1, sum_value + 1, INT_INF);
-  barrier::dense_matrix_t<i_t, uint8_t> take(n + 1, sum_value + 1, 0);
+  dense_matrix_t<i_t, i_t> dp(n + 1, sum_value + 1, INT_INF);
+  dense_matrix_t<i_t, uint8_t> take(n + 1, sum_value + 1, 0);
   dp(0, 0) = 0;
 
   // 4. Dynamic programming
@@ -2420,8 +2410,8 @@ f_t knapsack_generation_t<i_t, f_t>::exact_knapsack_problem_integer_values_fract
   solution.assign(n, 0.0);
 
   // dp(j, v) = minimum weight using first j items to get value v
-  barrier::dense_matrix_t<i_t, f_t> dp(n + 1, sum_value + 1, inf);
-  barrier::dense_matrix_t<i_t, uint8_t> take(n + 1, sum_value + 1, 0);
+  dense_matrix_t<i_t, f_t> dp(n + 1, sum_value + 1, inf);
+  dense_matrix_t<i_t, uint8_t> take(n + 1, sum_value + 1, 0);
   dp(0, 0) = 0;
 
   // 4. Dynamic programming
diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp
index e8dfba4ccd..c8559d4987 100644
--- a/cpp/src/cuts/cuts.hpp
+++ b/cpp/src/cuts/cuts.hpp
@@ -9,9 +9,9 @@
 #include <dual_simplex/basis_updates.hpp>
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/simplex_solver_settings.hpp>
-#include <dual_simplex/sparse_vector.hpp>
-#include <dual_simplex/types.hpp>
 #include <dual_simplex/user_problem.hpp>
+#include <linear_algebra/sparse_vector.hpp>
+#include <math_optimization/types.hpp>
 
 #include <algorithm>
 #include <array>
@@ -111,11 +111,10 @@ template <typename i_t, typename f_t>
 struct inequality_t {
   inequality_t() : vector(), rhs(0.0) {}
   inequality_t(i_t num_cols) : vector(num_cols, 0), rhs(0.0) {}
-  inequality_t(simplex::csr_matrix_t<i_t, f_t>& A, i_t row, f_t rhs_value)
-    : vector(A, row), rhs(rhs_value)
+  inequality_t(csr_matrix_t<i_t, f_t>& A, i_t row, f_t rhs_value) : vector(A, row), rhs(rhs_value)
   {
   }
-  simplex::sparse_vector_t<i_t, f_t> vector;
+  sparse_vector_t<i_t, f_t> vector;
   f_t rhs;
 
   void push_back(i_t j, f_t x)
@@ -262,7 +261,7 @@ void write_solution_for_cut_verification(const simplex::lp_problem_t<i_t, f_t>&
                                          const std::vector<f_t>& solution);
 
 template <typename i_t, typename f_t>
-void verify_cuts_against_saved_solution(const simplex::csr_matrix_t<i_t, f_t>& cuts,
+void verify_cuts_against_saved_solution(const csr_matrix_t<i_t, f_t>& cuts,
                                         const std::vector<f_t>& cut_rhs,
                                         const std::vector<f_t>& saved_solution);
 
@@ -296,7 +295,7 @@ class cut_pool_t {
   void score_cuts(std::vector<f_t>& x_relax);
 
   // We return the cuts in the form best_cuts*x <= best_rhs
-  i_t get_best_cuts(simplex::csr_matrix_t<i_t, f_t>& best_cuts,
+  i_t get_best_cuts(csr_matrix_t<i_t, f_t>& best_cuts,
                     std::vector<f_t>& best_rhs,
                     std::vector<cut_type_t>& best_cut_types);
 
@@ -318,7 +317,7 @@ class cut_pool_t {
   i_t original_vars_;
   const simplex::simplex_solver_settings_t<i_t, f_t>& settings_;
 
-  simplex::csr_matrix_t<i_t, f_t> cut_storage_;
+  csr_matrix_t<i_t, f_t> cut_storage_;
   std::vector<f_t> rhs_storage_;
   std::vector<i_t> cut_age_;
   std::vector<cut_type_t> cut_type_;
@@ -381,7 +380,7 @@ template <typename i_t, typename f_t>
 struct flow_cover_context_t {
   const simplex::lp_problem_t<i_t, f_t>& lp;
   const simplex::simplex_solver_settings_t<i_t, f_t>& settings;
-  simplex::csr_matrix_t<i_t, f_t>& Arow;
+  csr_matrix_t<i_t, f_t>& Arow;
   const variable_bounds_t<i_t, f_t>& variable_bounds;
   const std::vector<simplex::variable_type_t>& var_types;
   const std::vector<f_t>& xstar;
@@ -398,12 +397,12 @@ class flow_cover_generation_t {
  public:
   flow_cover_generation_t(const simplex::lp_problem_t<i_t, f_t>& lp,
                           const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
-                          simplex::csr_matrix_t<i_t, f_t>& Arow,
+                          csr_matrix_t<i_t, f_t>& Arow,
                           const std::vector<i_t>& new_slacks);
 
   i_t generate_cut(const simplex::lp_problem_t<i_t, f_t>& lp,
                    const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
-                   simplex::csr_matrix_t<i_t, f_t>& Arow,
+                   csr_matrix_t<i_t, f_t>& Arow,
                    const variable_bounds_t<i_t, f_t>& variable_bounds,
                    const std::vector<simplex::variable_type_t>& var_types,
                    const std::vector<f_t>& xstar,
@@ -521,13 +520,13 @@ class knapsack_generation_t {
  public:
   knapsack_generation_t(const simplex::lp_problem_t<i_t, f_t>& lp,
                         const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
-                        simplex::csr_matrix_t<i_t, f_t>& Arow,
+                        csr_matrix_t<i_t, f_t>& Arow,
                         const std::vector<i_t>& new_slacks,
                         const std::vector<simplex::variable_type_t>& var_types);
 
   i_t generate_knapsack_cut(const simplex::lp_problem_t<i_t, f_t>& lp,
                             const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
-                            simplex::csr_matrix_t<i_t, f_t>& Arow,
+                            csr_matrix_t<i_t, f_t>& Arow,
                             const std::vector<i_t>& new_slacks,
                             const std::vector<simplex::variable_type_t>& var_types,
                             const std::vector<f_t>& xstar,
@@ -589,7 +588,7 @@ class cut_generation_t {
   cut_generation_t(cut_pool_t<i_t, f_t>& cut_pool,
                    const simplex::lp_problem_t<i_t, f_t>& lp,
                    const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
-                   simplex::csr_matrix_t<i_t, f_t>& Arow,
+                   csr_matrix_t<i_t, f_t>& Arow,
                    const std::vector<i_t>& new_slacks,
                    const std::vector<simplex::variable_type_t>& var_types,
                    const simplex::user_problem_t<i_t, f_t>& user_problem,
@@ -608,7 +607,7 @@ class cut_generation_t {
 
   bool generate_cuts(const simplex::lp_problem_t<i_t, f_t>& lp,
                      const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
-                     simplex::csr_matrix_t<i_t, f_t>& Arow,
+                     csr_matrix_t<i_t, f_t>& Arow,
                      const std::vector<i_t>& new_slacks,
                      const std::vector<simplex::variable_type_t>& var_types,
                      simplex::basis_update_mpf_t<i_t, f_t>& basis_update,
@@ -624,7 +623,7 @@ class cut_generation_t {
   // Generate all mixed integer gomory cuts
   void generate_gomory_cuts(const simplex::lp_problem_t<i_t, f_t>& lp,
                             const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
-                            simplex::csr_matrix_t<i_t, f_t>& Arow,
+                            csr_matrix_t<i_t, f_t>& Arow,
                             const std::vector<i_t>& new_slacks,
                             const std::vector<simplex::variable_type_t>& var_types,
                             simplex::basis_update_mpf_t<i_t, f_t>& basis_update,
@@ -635,7 +634,7 @@ class cut_generation_t {
   // Generate all mixed integer rounding cuts
   void generate_mir_cuts(const simplex::lp_problem_t<i_t, f_t>& lp,
                          const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
-                         simplex::csr_matrix_t<i_t, f_t>& Arow,
+                         csr_matrix_t<i_t, f_t>& Arow,
                          const std::vector<i_t>& new_slacks,
                          const std::vector<simplex::variable_type_t>& var_types,
                          const std::vector<f_t>& xstar,
@@ -645,7 +644,7 @@ class cut_generation_t {
   // Generate all knapsack cuts
   void generate_knapsack_cuts(const simplex::lp_problem_t<i_t, f_t>& lp,
                               const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
-                              simplex::csr_matrix_t<i_t, f_t>& Arow,
+                              csr_matrix_t<i_t, f_t>& Arow,
                               const std::vector<i_t>& new_slacks,
                               const std::vector<simplex::variable_type_t>& var_types,
                               const std::vector<f_t>& xstar,
@@ -654,7 +653,7 @@ class cut_generation_t {
   // Generate all flow cover cuts
   void generate_flow_cover_cuts(const simplex::lp_problem_t<i_t, f_t>& lp,
                                 const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
-                                simplex::csr_matrix_t<i_t, f_t>& Arow,
+                                csr_matrix_t<i_t, f_t>& Arow,
                                 const std::vector<simplex::variable_type_t>& var_types,
                                 const std::vector<f_t>& xstar,
                                 variable_bounds_t<i_t, f_t>& variable_bounds,
@@ -763,7 +762,7 @@ class tableau_equality_t {
   // Generates the base inequalities: C*x == d that will be turned into cuts
   i_t generate_base_equality(const simplex::lp_problem_t<i_t, f_t>& lp,
                              const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
-                             simplex::csr_matrix_t<i_t, f_t>& Arow,
+                             csr_matrix_t<i_t, f_t>& Arow,
                              const std::vector<simplex::variable_type_t>& var_types,
                              simplex::basis_update_mpf_t<i_t, f_t>& basis_update,
                              const std::vector<f_t>& xstar,
@@ -786,7 +785,7 @@ class variable_bounds_t {
   variable_bounds_t(const simplex::lp_problem_t<i_t, f_t>& lp,
                     const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
                     const std::vector<simplex::variable_type_t>& var_types,
-                    const simplex::csr_matrix_t<i_t, f_t>& Arow,
+                    const csr_matrix_t<i_t, f_t>& Arow,
                     const std::vector<i_t>& new_slacks);
 
   std::vector<i_t> upper_offsets;
@@ -833,15 +832,14 @@ class variable_bounds_t {
   {
     if (num_lower_inf == 0) {
       return activity - lower_activity_i - lower_activity_j;
-    } else if (num_lower_inf == 1 && lower_activity_j == -simplex::inf) {
+    } else if (num_lower_inf == 1 && lower_activity_j == -inf) {
       return activity - lower_activity_i;
-    } else if (num_lower_inf == 1 && lower_activity_i == -simplex::inf) {
+    } else if (num_lower_inf == 1 && lower_activity_i == -inf) {
       return activity - lower_activity_j;
-    } else if (num_lower_inf == 2 && lower_activity_i == -simplex::inf &&
-               lower_activity_j == -simplex::inf) {
+    } else if (num_lower_inf == 2 && lower_activity_i == -inf && lower_activity_j == -inf) {
       return activity;
     } else {
-      return -simplex::inf;
+      return -inf;
     }
   }
 
@@ -860,15 +858,14 @@ class variable_bounds_t {
   {
     if (num_upper_inf == 0) {
       return activity - upper_activity_i - upper_activity_j;
-    } else if (num_upper_inf == 1 && upper_activity_j == simplex::inf) {
+    } else if (num_upper_inf == 1 && upper_activity_j == inf) {
       return activity - upper_activity_i;
-    } else if (num_upper_inf == 1 && upper_activity_i == simplex::inf) {
+    } else if (num_upper_inf == 1 && upper_activity_i == inf) {
       return activity - upper_activity_j;
-    } else if (num_upper_inf == 2 && upper_activity_i == simplex::inf &&
-               upper_activity_j == simplex::inf) {
+    } else if (num_upper_inf == 2 && upper_activity_i == inf && upper_activity_j == inf) {
       return activity;
     } else {
-      return simplex::inf;
+      return inf;
     }
   }
 
@@ -890,7 +887,7 @@ class complemented_mixed_integer_rounding_cut_t {
 
   void compute_initial_scores_for_rows(const simplex::lp_problem_t<i_t, f_t>& lp,
                                        const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
-                                       const simplex::csr_matrix_t<i_t, f_t>& Arow,
+                                       const csr_matrix_t<i_t, f_t>& Arow,
                                        const std::vector<f_t>& xstar,
                                        const std::vector<f_t>& ystar,
                                        std::vector<f_t>& score);
@@ -957,14 +954,14 @@ class complemented_mixed_integer_rounding_cut_t {
                                  inequality_t<i_t, f_t>& cut);
 
   void substitute_slacks(const simplex::lp_problem_t<i_t, f_t>& lp,
-                         simplex::csr_matrix_t<i_t, f_t>& Arow,
+                         csr_matrix_t<i_t, f_t>& Arow,
                          inequality_t<i_t, f_t>& cut);
 
   // Combine the pivot row with the inequality to eliminate the variable j
   // The new inequality is returned in inequality and inequality_rhs
   // The multiplier for the pivot row is returned
   f_t combine_rows(const simplex::lp_problem_t<i_t, f_t>& lp,
-                   simplex::csr_matrix_t<i_t, f_t>& Arow,
+                   csr_matrix_t<i_t, f_t>& Arow,
                    i_t j,
                    const inequality_t<i_t, f_t>& pivot_row,
                    inequality_t<i_t, f_t>& inequality);
@@ -1049,7 +1046,7 @@ class strong_cg_cut_t {
 
 template <typename i_t, typename f_t>
 i_t add_cuts(const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
-             const simplex::csr_matrix_t<i_t, f_t>& cuts,
+             const csr_matrix_t<i_t, f_t>& cuts,
              const std::vector<f_t>& cut_rhs,
              simplex::lp_problem_t<i_t, f_t>& lp,
              std::vector<i_t>& new_slacks,
@@ -1064,7 +1061,7 @@ template <typename i_t, typename f_t>
 i_t remove_cuts(simplex::lp_problem_t<i_t, f_t>& lp,
                 const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
                 f_t start_time,
-                simplex::csr_matrix_t<i_t, f_t>& Arow,
+                csr_matrix_t<i_t, f_t>& Arow,
                 std::vector<i_t>& new_slacks,
                 i_t original_rows,
                 std::vector<simplex::variable_type_t>& var_types,
diff --git a/cpp/src/dual_simplex/CMakeLists.txt b/cpp/src/dual_simplex/CMakeLists.txt
index 7341c0de30..228f2aedd7 100644
--- a/cpp/src/dual_simplex/CMakeLists.txt
+++ b/cpp/src/dual_simplex/CMakeLists.txt
@@ -19,11 +19,7 @@ set(DUAL_SIMPLEX_SRC_FILES
   ${CMAKE_CURRENT_SOURCE_DIR}/scaling.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/singletons.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/solve.cpp
-  ${CMAKE_CURRENT_SOURCE_DIR}/sparse_matrix.cpp
-  ${CMAKE_CURRENT_SOURCE_DIR}/sparse_vector.cpp
-  ${CMAKE_CURRENT_SOURCE_DIR}/tic_toc.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/triangle_solve.cpp
-  ${CMAKE_CURRENT_SOURCE_DIR}/vector_math.cpp
   )
 
 # Uncomment to enable debug info
diff --git a/cpp/src/dual_simplex/basis_solves.cpp b/cpp/src/dual_simplex/basis_solves.cpp
index a961e085aa..66541c60ff 100644
--- a/cpp/src/dual_simplex/basis_solves.cpp
+++ b/cpp/src/dual_simplex/basis_solves.cpp
@@ -10,8 +10,8 @@
 #include <dual_simplex/initial_basis.hpp>
 #include <dual_simplex/right_looking_lu.hpp>
 #include <dual_simplex/singletons.hpp>
-#include <dual_simplex/tic_toc.hpp>
 #include <dual_simplex/triangle_solve.hpp>
+#include <math_optimization/tic_toc.hpp>
 
 #include <raft/core/nvtx.hpp>
 
diff --git a/cpp/src/dual_simplex/basis_solves.hpp b/cpp/src/dual_simplex/basis_solves.hpp
index 74b213f362..dfa8f7f788 100644
--- a/cpp/src/dual_simplex/basis_solves.hpp
+++ b/cpp/src/dual_simplex/basis_solves.hpp
@@ -9,8 +9,8 @@
 
 #include <dual_simplex/initial_basis.hpp>
 #include <dual_simplex/presolve.hpp>
-#include <dual_simplex/sparse_matrix.hpp>
-#include <dual_simplex/types.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
+#include <math_optimization/types.hpp>
 
 namespace cuopt::mathematical_optimization::simplex {
 
diff --git a/cpp/src/dual_simplex/basis_updates.hpp b/cpp/src/dual_simplex/basis_updates.hpp
index 9d98a1fb43..d1c623db55 100644
--- a/cpp/src/dual_simplex/basis_updates.hpp
+++ b/cpp/src/dual_simplex/basis_updates.hpp
@@ -9,9 +9,9 @@
 
 #include <dual_simplex/initial_basis.hpp>
 #include <dual_simplex/simplex_solver_settings.hpp>
-#include <dual_simplex/sparse_matrix.hpp>
-#include <dual_simplex/sparse_vector.hpp>
-#include <dual_simplex/types.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
+#include <linear_algebra/sparse_vector.hpp>
+#include <math_optimization/types.hpp>
 
 #include <numeric>
 
diff --git a/cpp/src/dual_simplex/bound_flipping_ratio_test.cpp b/cpp/src/dual_simplex/bound_flipping_ratio_test.cpp
index 2434ea1acf..cb0964dc05 100644
--- a/cpp/src/dual_simplex/bound_flipping_ratio_test.cpp
+++ b/cpp/src/dual_simplex/bound_flipping_ratio_test.cpp
@@ -7,7 +7,7 @@
 
 #include <dual_simplex/bound_flipping_ratio_test.hpp>
 
-#include <dual_simplex/tic_toc.hpp>
+#include <math_optimization/tic_toc.hpp>
 
 #include <algorithm>
 #include <cmath>
diff --git a/cpp/src/dual_simplex/bound_flipping_ratio_test.hpp b/cpp/src/dual_simplex/bound_flipping_ratio_test.hpp
index 25ca75618f..2e73d05eff 100644
--- a/cpp/src/dual_simplex/bound_flipping_ratio_test.hpp
+++ b/cpp/src/dual_simplex/bound_flipping_ratio_test.hpp
@@ -8,7 +8,7 @@
 
 #include <dual_simplex/initial_basis.hpp>
 #include <dual_simplex/simplex_solver_settings.hpp>
-#include <dual_simplex/types.hpp>
+#include <math_optimization/types.hpp>
 
 #include <vector>
 
diff --git a/cpp/src/dual_simplex/crossover.cpp b/cpp/src/dual_simplex/crossover.cpp
index 068166cc5e..e1ba272adf 100644
--- a/cpp/src/dual_simplex/crossover.cpp
+++ b/cpp/src/dual_simplex/crossover.cpp
@@ -13,7 +13,7 @@
 #include <dual_simplex/phase2.hpp>
 #include <dual_simplex/primal.hpp>
 #include <dual_simplex/solve.hpp>
-#include <dual_simplex/tic_toc.hpp>
+#include <math_optimization/tic_toc.hpp>
 
 #include <raft/core/nvtx.hpp>
 
diff --git a/cpp/src/dual_simplex/crossover.hpp b/cpp/src/dual_simplex/crossover.hpp
index 1b0b0ba823..13317107a9 100644
--- a/cpp/src/dual_simplex/crossover.hpp
+++ b/cpp/src/dual_simplex/crossover.hpp
@@ -10,8 +10,8 @@
 #include <dual_simplex/initial_basis.hpp>
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/solution.hpp>
-#include <dual_simplex/types.hpp>
 #include <dual_simplex/user_problem.hpp>
+#include <math_optimization/types.hpp>
 
 namespace cuopt::mathematical_optimization::simplex {
 
diff --git a/cpp/src/dual_simplex/folding.cpp b/cpp/src/dual_simplex/folding.cpp
index 8e6724a07f..a56659b3ff 100644
--- a/cpp/src/dual_simplex/folding.cpp
+++ b/cpp/src/dual_simplex/folding.cpp
@@ -7,7 +7,7 @@
 
 #include <dual_simplex/folding.hpp>
 
-#include <dual_simplex/tic_toc.hpp>
+#include <math_optimization/tic_toc.hpp>
 
 #include <map>
 #include <numeric>
diff --git a/cpp/src/dual_simplex/folding.hpp b/cpp/src/dual_simplex/folding.hpp
index 74dd2e8d61..adb88a0239 100644
--- a/cpp/src/dual_simplex/folding.hpp
+++ b/cpp/src/dual_simplex/folding.hpp
@@ -9,7 +9,7 @@
 
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/simplex_solver_settings.hpp>
-#include <dual_simplex/types.hpp>
+#include <math_optimization/types.hpp>
 
 namespace cuopt::mathematical_optimization::simplex {
 
diff --git a/cpp/src/dual_simplex/initial_basis.cpp b/cpp/src/dual_simplex/initial_basis.cpp
index 19426853b5..e75bad0e67 100644
--- a/cpp/src/dual_simplex/initial_basis.cpp
+++ b/cpp/src/dual_simplex/initial_basis.cpp
@@ -9,7 +9,7 @@
 
 #include <dual_simplex/right_looking_lu.hpp>
 #include <dual_simplex/singletons.hpp>
-#include <dual_simplex/tic_toc.hpp>
+#include <math_optimization/tic_toc.hpp>
 
 #include <raft/core/nvtx.hpp>
 
diff --git a/cpp/src/dual_simplex/initial_basis.hpp b/cpp/src/dual_simplex/initial_basis.hpp
index ef7230946b..d31da0f9e4 100644
--- a/cpp/src/dual_simplex/initial_basis.hpp
+++ b/cpp/src/dual_simplex/initial_basis.hpp
@@ -9,8 +9,8 @@
 
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/simplex_solver_settings.hpp>
-#include <dual_simplex/sparse_matrix.hpp>
-#include <dual_simplex/types.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
+#include <math_optimization/types.hpp>
 
 namespace cuopt::mathematical_optimization::simplex {
 
diff --git a/cpp/src/dual_simplex/phase1.cpp b/cpp/src/dual_simplex/phase1.cpp
index 96b941d992..1b9bcdf0f2 100644
--- a/cpp/src/dual_simplex/phase1.cpp
+++ b/cpp/src/dual_simplex/phase1.cpp
@@ -9,7 +9,7 @@
 
 #include <dual_simplex/basis_solves.hpp>
 #include <dual_simplex/initial_basis.hpp>
-#include <dual_simplex/sparse_matrix.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
 
 #include <cassert>
 #include <cstdio>
diff --git a/cpp/src/dual_simplex/phase1.hpp b/cpp/src/dual_simplex/phase1.hpp
index e180746bef..119d39f48b 100644
--- a/cpp/src/dual_simplex/phase1.hpp
+++ b/cpp/src/dual_simplex/phase1.hpp
@@ -9,7 +9,7 @@
 
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/solution.hpp>
-#include <dual_simplex/types.hpp>
+#include <math_optimization/types.hpp>
 
 #include <limits>
 #include <vector>
diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp
index 36daed4c8d..2e3c1e05c5 100644
--- a/cpp/src/dual_simplex/phase2.cpp
+++ b/cpp/src/dual_simplex/phase2.cpp
@@ -13,8 +13,8 @@
 #include <dual_simplex/phase2.hpp>
 #include <dual_simplex/random.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 <utilities/scope_guard.hpp>
 #include <utilities/timing_utils.hpp>
diff --git a/cpp/src/dual_simplex/phase2.hpp b/cpp/src/dual_simplex/phase2.hpp
index 294c2b63bc..daa946e019 100644
--- a/cpp/src/dual_simplex/phase2.hpp
+++ b/cpp/src/dual_simplex/phase2.hpp
@@ -12,7 +12,7 @@
 #include <dual_simplex/logger.hpp>
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/simplex_solver_settings.hpp>
-#include <dual_simplex/types.hpp>
+#include <math_optimization/types.hpp>
 #include <utilities/memory_instrumentation.hpp>
 
 #include <vector>
diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp
index 9c55185c6d..30a7623bbe 100644
--- a/cpp/src/dual_simplex/presolve.cpp
+++ b/cpp/src/dual_simplex/presolve.cpp
@@ -11,7 +11,7 @@
 #include <dual_simplex/folding.hpp>
 #include <dual_simplex/right_looking_lu.hpp>
 #include <dual_simplex/solve.hpp>
-#include <dual_simplex/tic_toc.hpp>
+#include <math_optimization/tic_toc.hpp>
 
 #include <algorithm>
 #include <cmath>
diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp
index 855a4e6145..5dd1df0304 100644
--- a/cpp/src/dual_simplex/presolve.hpp
+++ b/cpp/src/dual_simplex/presolve.hpp
@@ -9,9 +9,9 @@
 
 #include <dual_simplex/simplex_solver_settings.hpp>
 #include <dual_simplex/solution.hpp>
-#include <dual_simplex/sparse_matrix.hpp>
-#include <dual_simplex/types.hpp>
 #include <dual_simplex/user_problem.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
+#include <math_optimization/types.hpp>
 
 #include <fstream>
 #include <iomanip>
diff --git a/cpp/src/dual_simplex/primal.cpp b/cpp/src/dual_simplex/primal.cpp
index 5451a57f08..78c7107ca3 100644
--- a/cpp/src/dual_simplex/primal.cpp
+++ b/cpp/src/dual_simplex/primal.cpp
@@ -12,7 +12,7 @@
 #include <dual_simplex/initial_basis.hpp>
 #include <dual_simplex/phase1.hpp>
 #include <dual_simplex/solve.hpp>
-#include <dual_simplex/tic_toc.hpp>
+#include <math_optimization/tic_toc.hpp>
 
 namespace cuopt::mathematical_optimization::simplex {
 
diff --git a/cpp/src/dual_simplex/primal.hpp b/cpp/src/dual_simplex/primal.hpp
index 33bcfa722d..930958a802 100644
--- a/cpp/src/dual_simplex/primal.hpp
+++ b/cpp/src/dual_simplex/primal.hpp
@@ -11,7 +11,7 @@
 #include <dual_simplex/phase1.hpp>
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/simplex_solver_settings.hpp>
-#include <dual_simplex/types.hpp>
+#include <math_optimization/types.hpp>
 
 #include <vector>
 
diff --git a/cpp/src/dual_simplex/right_looking_lu.cpp b/cpp/src/dual_simplex/right_looking_lu.cpp
index 2882662f17..b982ee35a0 100644
--- a/cpp/src/dual_simplex/right_looking_lu.cpp
+++ b/cpp/src/dual_simplex/right_looking_lu.cpp
@@ -6,7 +6,7 @@
 /* clang-format on */
 
 #include <dual_simplex/right_looking_lu.hpp>
-#include <dual_simplex/tic_toc.hpp>
+#include <math_optimization/tic_toc.hpp>
 #include <utilities/memory_instrumentation.hpp>
 
 #include <raft/core/nvtx.hpp>
diff --git a/cpp/src/dual_simplex/right_looking_lu.hpp b/cpp/src/dual_simplex/right_looking_lu.hpp
index 58f0b4ee7d..abe148cf20 100644
--- a/cpp/src/dual_simplex/right_looking_lu.hpp
+++ b/cpp/src/dual_simplex/right_looking_lu.hpp
@@ -8,7 +8,7 @@
 #pragma once
 
 #include <dual_simplex/simplex_solver_settings.hpp>
-#include <dual_simplex/sparse_matrix.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
 
 #include <optional>
 
diff --git a/cpp/src/dual_simplex/scaling.cpp b/cpp/src/dual_simplex/scaling.cpp
index ac96496ec6..21f42a29dd 100644
--- a/cpp/src/dual_simplex/scaling.cpp
+++ b/cpp/src/dual_simplex/scaling.cpp
@@ -6,7 +6,7 @@
 /* clang-format on */
 
 #include <dual_simplex/scaling.hpp>
-#include <dual_simplex/sparse_matrix.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
 
 #include <cmath>
 
diff --git a/cpp/src/dual_simplex/scaling.hpp b/cpp/src/dual_simplex/scaling.hpp
index eab60c486b..9ddc670315 100644
--- a/cpp/src/dual_simplex/scaling.hpp
+++ b/cpp/src/dual_simplex/scaling.hpp
@@ -9,8 +9,8 @@
 
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/simplex_solver_settings.hpp>
-#include <dual_simplex/sparse_matrix.hpp>
-#include <dual_simplex/types.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
+#include <math_optimization/types.hpp>
 
 #include <vector>
 
diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp
index 45168ced8a..7cf185cefc 100644
--- a/cpp/src/dual_simplex/simplex_solver_settings.hpp
+++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp
@@ -10,7 +10,7 @@
 #include <cuopt/mathematical_optimization/mip/diving_hyper_params.hpp>
 
 #include <dual_simplex/logger.hpp>
-#include <dual_simplex/types.hpp>
+#include <math_optimization/types.hpp>
 
 #include <omp.h>
 #include <algorithm>
diff --git a/cpp/src/dual_simplex/singletons.cpp b/cpp/src/dual_simplex/singletons.cpp
index 4d686b2ac1..8b151337b7 100644
--- a/cpp/src/dual_simplex/singletons.cpp
+++ b/cpp/src/dual_simplex/singletons.cpp
@@ -7,7 +7,7 @@
 
 #include <dual_simplex/singletons.hpp>
 #include <dual_simplex/triangle_solve.hpp>
-#include <dual_simplex/types.hpp>
+#include <math_optimization/types.hpp>
 
 #include <cstdio>
 
diff --git a/cpp/src/dual_simplex/singletons.hpp b/cpp/src/dual_simplex/singletons.hpp
index 59964fc428..cc6ac95f9d 100644
--- a/cpp/src/dual_simplex/singletons.hpp
+++ b/cpp/src/dual_simplex/singletons.hpp
@@ -7,7 +7,7 @@
 
 #pragma once
 
-#include <dual_simplex/sparse_matrix.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
 
 #include <numeric>
 #include <queue>
diff --git a/cpp/src/dual_simplex/solution.hpp b/cpp/src/dual_simplex/solution.hpp
index 435e9a9fce..59541dcce8 100644
--- a/cpp/src/dual_simplex/solution.hpp
+++ b/cpp/src/dual_simplex/solution.hpp
@@ -7,7 +7,7 @@
 
 #pragma once
 
-#include <dual_simplex/types.hpp>
+#include <math_optimization/types.hpp>
 
 #include <utilities/omp_helpers.hpp>
 #include <vector>
diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp
index 7c4f13d8c5..697af9e869 100644
--- a/cpp/src/dual_simplex/solve.cpp
+++ b/cpp/src/dual_simplex/solve.cpp
@@ -20,11 +20,11 @@
 #include <dual_simplex/primal.hpp>
 #include <dual_simplex/scaling.hpp>
 #include <dual_simplex/singletons.hpp>
-#include <dual_simplex/sparse_matrix.hpp>
-#include <dual_simplex/tic_toc.hpp>
 #include <dual_simplex/triangle_solve.hpp>
-#include <dual_simplex/types.hpp>
 #include <dual_simplex/user_problem.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
+#include <math_optimization/tic_toc.hpp>
+#include <math_optimization/types.hpp>
 
 #include <raft/core/nvtx.hpp>
 
diff --git a/cpp/src/dual_simplex/solve.hpp b/cpp/src/dual_simplex/solve.hpp
index 60a6ead66b..7cc9a9f5cf 100644
--- a/cpp/src/dual_simplex/solve.hpp
+++ b/cpp/src/dual_simplex/solve.hpp
@@ -11,7 +11,7 @@
 #include <dual_simplex/initial_basis.hpp>
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/simplex_solver_settings.hpp>
-#include <dual_simplex/types.hpp>
+#include <math_optimization/types.hpp>
 
 namespace cuopt {
 struct work_limit_context_t;
diff --git a/cpp/src/dual_simplex/triangle_solve.hpp b/cpp/src/dual_simplex/triangle_solve.hpp
index 7e907512c2..e22e206634 100644
--- a/cpp/src/dual_simplex/triangle_solve.hpp
+++ b/cpp/src/dual_simplex/triangle_solve.hpp
@@ -7,9 +7,9 @@
 
 #pragma once
 
-#include <dual_simplex/sparse_matrix.hpp>
-#include <dual_simplex/sparse_vector.hpp>
-#include <dual_simplex/types.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
+#include <linear_algebra/sparse_vector.hpp>
+#include <math_optimization/types.hpp>
 
 #include <optional>
 
diff --git a/cpp/src/dual_simplex/user_problem.hpp b/cpp/src/dual_simplex/user_problem.hpp
index fe19b4d456..419b33e1e1 100644
--- a/cpp/src/dual_simplex/user_problem.hpp
+++ b/cpp/src/dual_simplex/user_problem.hpp
@@ -8,8 +8,8 @@
 #pragma once
 
 #include <dual_simplex/solution.hpp>
-#include <dual_simplex/sparse_matrix.hpp>
-#include <dual_simplex/types.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
+#include <math_optimization/types.hpp>
 
 #include <raft/core/handle.hpp>
 
diff --git a/cpp/src/linear_algebra/CMakeLists.txt b/cpp/src/linear_algebra/CMakeLists.txt
new file mode 100644
index 0000000000..875a016544
--- /dev/null
+++ b/cpp/src/linear_algebra/CMakeLists.txt
@@ -0,0 +1,13 @@
+# cmake-format: off
+# SPDX-FileCopyrightText: Copyright (c) 2024-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+# SPDX-License-Identifier: Apache-2.0
+# cmake-format: on
+
+set(LINEAR_ALGEBRA_SRC_FILES
+  ${CMAKE_CURRENT_SOURCE_DIR}/sparse_matrix.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/sparse_vector.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/vector_math.cpp
+  )
+
+set(CUOPT_SRC_FILES ${CUOPT_SRC_FILES}
+  ${LINEAR_ALGEBRA_SRC_FILES} PARENT_SCOPE)
diff --git a/cpp/src/barrier/dense_matrix.hpp b/cpp/src/linear_algebra/dense_matrix.hpp
similarity index 95%
rename from cpp/src/barrier/dense_matrix.hpp
rename to cpp/src/linear_algebra/dense_matrix.hpp
index 260c4a2957..f991fc3bdc 100644
--- a/cpp/src/barrier/dense_matrix.hpp
+++ b/cpp/src/linear_algebra/dense_matrix.hpp
@@ -5,14 +5,14 @@
  */
 /* clang-format on */
 
-#include <barrier/dense_vector.hpp>
+#include <linear_algebra/dense_vector.hpp>
 
-#include <dual_simplex/sparse_matrix.hpp>
-#include <dual_simplex/types.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
+#include <math_optimization/types.hpp>
 
 #pragma once
 
-namespace cuopt::mathematical_optimization::barrier {
+namespace cuopt::mathematical_optimization {
 
 template <typename i_t, typename f_t>
 class dense_matrix_t {
@@ -32,7 +32,7 @@ class dense_matrix_t {
 
   f_t operator()(i_t row, i_t col) const { return values[col * m + row]; }
 
-  void from_sparse(const simplex::csc_matrix_t<i_t, f_t>& A, i_t sparse_column, i_t dense_column)
+  void from_sparse(const csc_matrix_t<i_t, f_t>& A, i_t sparse_column, i_t dense_column)
   {
     for (i_t i = 0; i < m; i++) {
       this->operator()(i, dense_column) = 0.0;
@@ -242,4 +242,4 @@ class dense_matrix_t {
   std::vector<f_t> values;
 };
 
-}  // namespace cuopt::mathematical_optimization::barrier
+}  // namespace cuopt::mathematical_optimization
diff --git a/cpp/src/barrier/dense_vector.hpp b/cpp/src/linear_algebra/dense_vector.hpp
similarity index 96%
rename from cpp/src/barrier/dense_vector.hpp
rename to cpp/src/linear_algebra/dense_vector.hpp
index c444bbc8f0..61284b65bd 100644
--- a/cpp/src/barrier/dense_vector.hpp
+++ b/cpp/src/linear_algebra/dense_vector.hpp
@@ -6,13 +6,13 @@
 /* clang-format on */
 #pragma once
 
-#include <dual_simplex/types.hpp>
+#include <math_optimization/types.hpp>
 
 #include <cmath>
 #include <cstdio>
 #include <vector>
 
-namespace cuopt::mathematical_optimization::barrier {
+namespace cuopt::mathematical_optimization {
 
 template <typename i_t, typename f_t, typename Allocator = std::allocator<f_t>>
 class dense_vector_t : public std::vector<f_t, Allocator> {
@@ -56,7 +56,7 @@ class dense_vector_t : public std::vector<f_t, Allocator> {
   f_t minimum() const
   {
     const i_t n = this->size();
-    f_t min_x   = simplex::inf;
+    f_t min_x   = inf;
     for (i_t i = 0; i < n; i++) {
       min_x = std::min(min_x, (*this)[i]);
     }
@@ -66,7 +66,7 @@ class dense_vector_t : public std::vector<f_t, Allocator> {
   f_t maximum() const
   {
     const i_t n = this->size();
-    f_t max_x   = -simplex::inf;
+    f_t max_x   = -inf;
     for (i_t i = 0; i < n; i++) {
       max_x = std::max(max_x, (*this)[i]);
     }
@@ -186,7 +186,7 @@ class dense_vector_t : public std::vector<f_t, Allocator> {
 
   void ensure_positive(f_t epsilon_adjust, const std::vector<i_t>& mask)
   {
-    f_t min_x   = simplex::inf;
+    f_t min_x   = inf;
     const i_t n = this->size();
     for (i_t i = 0; i < n; i++) {
       if (mask[i]) { min_x = std::min(min_x, (*this)[i]); }
@@ -246,4 +246,4 @@ std::vector<T> copy(const std::vector<T, Alloc>& src)
   return std::vector<T>(src.begin(), src.end());
 }
 
-}  // namespace cuopt::mathematical_optimization::barrier
+}  // namespace cuopt::mathematical_optimization
diff --git a/cpp/src/mip_heuristics/utilities/sort_csr.cuh b/cpp/src/linear_algebra/sort_csr.cuh
similarity index 93%
rename from cpp/src/mip_heuristics/utilities/sort_csr.cuh
rename to cpp/src/linear_algebra/sort_csr.cuh
index dffc2e004f..23b9fd2d57 100644
--- a/cpp/src/mip_heuristics/utilities/sort_csr.cuh
+++ b/cpp/src/linear_algebra/sort_csr.cuh
@@ -7,15 +7,16 @@
 
 #pragma once
 
-#include <mip_heuristics/problem/problem.cuh>
+#include <cuopt/mathematical_optimization/optimization_problem.hpp>
 
 #include <cub/cub.cuh>
 
+#include <raft/core/nvtx.hpp>
 #include <rmm/device_uvector.hpp>
 
 namespace cuopt {
 
-namespace mathematical_optimization::mip {
+namespace mathematical_optimization {
 
 template <typename i_t, typename f_t>
 void sort_csr(optimization_problem_t<i_t, f_t>& op_problem)
@@ -53,5 +54,5 @@ void sort_csr(optimization_problem_t<i_t, f_t>& op_problem)
   RAFT_CUDA_TRY(cudaStreamSynchronize(stream_view));
 }
 
-}  // namespace mathematical_optimization::mip
+}  // namespace mathematical_optimization
 }  // namespace cuopt
diff --git a/cpp/src/dual_simplex/sparse_matrix.cpp b/cpp/src/linear_algebra/sparse_matrix.cpp
similarity index 99%
rename from cpp/src/dual_simplex/sparse_matrix.cpp
rename to cpp/src/linear_algebra/sparse_matrix.cpp
index 4509b97fb1..158fecefcf 100644
--- a/cpp/src/dual_simplex/sparse_matrix.cpp
+++ b/cpp/src/linear_algebra/sparse_matrix.cpp
@@ -6,10 +6,10 @@
 /* clang-format on */
 
 // #include <dual_simplex/dense_vector.hpp>
-#include <dual_simplex/sparse_matrix.hpp>
-#include <dual_simplex/sparse_vector.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
+#include <linear_algebra/sparse_vector.hpp>
 
-#include <dual_simplex/types.hpp>
+#include <math_optimization/types.hpp>
 #include <mip_heuristics/mip_constants.hpp>
 
 // #include <thrust/for_each.h>
@@ -25,7 +25,7 @@
 #include <cmath>
 #include <cstdio>
 
-namespace cuopt::mathematical_optimization::simplex {
+namespace cuopt::mathematical_optimization {
 
 template <typename i_t, typename f_t>
 void csc_matrix_t<i_t, f_t>::reallocate(i_t new_nz)
@@ -1006,4 +1006,4 @@ matrix_transpose_vector_multiply<int, double, std::allocator<double>, std::alloc
 
 #endif
 
-}  // namespace cuopt::mathematical_optimization::simplex
+}  // namespace cuopt::mathematical_optimization
diff --git a/cpp/src/dual_simplex/sparse_matrix.hpp b/cpp/src/linear_algebra/sparse_matrix.hpp
similarity index 96%
rename from cpp/src/dual_simplex/sparse_matrix.hpp
rename to cpp/src/linear_algebra/sparse_matrix.hpp
index 2aad545ec8..95ee0f0f32 100644
--- a/cpp/src/dual_simplex/sparse_matrix.hpp
+++ b/cpp/src/linear_algebra/sparse_matrix.hpp
@@ -7,8 +7,8 @@
 
 #pragma once
 
-#include <dual_simplex/types.hpp>
-#include <dual_simplex/vector_math.hpp>
+#include <linear_algebra/vector_math.hpp>
+#include <math_optimization/types.hpp>
 
 #include <algorithm>
 #include <cassert>
@@ -16,7 +16,7 @@
 #include <string>
 #include <vector>
 
-namespace cuopt::mathematical_optimization::simplex {
+namespace cuopt::mathematical_optimization {
 
 template <typename i_t, typename f_t>
 class csr_matrix_t;  // Forward declaration of CSR matrix needed to define CSC
@@ -55,8 +55,7 @@ class csc_matrix_t {
   i_t col_length(i_t j) const { return col_start[j + 1] - col_start[j]; }
 
   // Convert the CSC matrix to a CSR matrix
-  i_t to_compressed_row(
-    cuopt::mathematical_optimization::simplex::csr_matrix_t<i_t, f_t>& Arow) const;
+  i_t to_compressed_row(cuopt::mathematical_optimization::csr_matrix_t<i_t, f_t>& Arow) const;
 
   // Permutes rows of a sparse matrix A. Computes C = A(p, :)
   i_t permute_rows(const std::vector<i_t>& pinv, csc_matrix_t<i_t, f_t>& C) const;
@@ -312,4 +311,4 @@ i_t matrix_vector_multiply(
   return 0;
 }
 
-}  // namespace cuopt::mathematical_optimization::simplex
+}  // namespace cuopt::mathematical_optimization
diff --git a/cpp/src/dual_simplex/sparse_vector.cpp b/cpp/src/linear_algebra/sparse_vector.cpp
similarity index 97%
rename from cpp/src/dual_simplex/sparse_vector.cpp
rename to cpp/src/linear_algebra/sparse_vector.cpp
index 08821ea719..17839b342f 100644
--- a/cpp/src/dual_simplex/sparse_vector.cpp
+++ b/cpp/src/linear_algebra/sparse_vector.cpp
@@ -5,13 +5,13 @@
  */
 /* clang-format on */
 
-#include <dual_simplex/sparse_vector.hpp>
+#include <linear_algebra/sparse_vector.hpp>
 
 #include <algorithm>
 #include <cassert>
 #include <cstdio>
 
-namespace cuopt::mathematical_optimization::simplex {
+namespace cuopt::mathematical_optimization {
 
 template <typename i_t, typename f_t>
 sparse_vector_t<i_t, f_t>::sparse_vector_t(const csc_matrix_t<i_t, f_t>& A, i_t col)
@@ -285,4 +285,4 @@ void sparse_vector_t<i_t, f_t>::squeeze(sparse_vector_t<i_t, f_t>& y) const
 template class sparse_vector_t<int, double>;
 #endif
 
-}  // namespace cuopt::mathematical_optimization::simplex
+}  // namespace cuopt::mathematical_optimization
diff --git a/cpp/src/dual_simplex/sparse_vector.hpp b/cpp/src/linear_algebra/sparse_vector.hpp
similarity index 93%
rename from cpp/src/dual_simplex/sparse_vector.hpp
rename to cpp/src/linear_algebra/sparse_vector.hpp
index 03d9d9f601..4d95b0bc82 100644
--- a/cpp/src/dual_simplex/sparse_vector.hpp
+++ b/cpp/src/linear_algebra/sparse_vector.hpp
@@ -7,12 +7,12 @@
 
 #pragma once
 
-#include <dual_simplex/sparse_matrix.hpp>
-#include <dual_simplex/types.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
+#include <math_optimization/types.hpp>
 
 #include <vector>
 
-namespace cuopt::mathematical_optimization::simplex {
+namespace cuopt::mathematical_optimization {
 
 // A sparse vector stored as a list of nonzero coefficients and their indices
 template <typename i_t, typename f_t>
@@ -71,4 +71,4 @@ class sparse_vector_t {
   std::vector<f_t> x;
 };
 
-}  // namespace cuopt::mathematical_optimization::simplex
+}  // namespace cuopt::mathematical_optimization
diff --git a/cpp/src/dual_simplex/vector_math.cpp b/cpp/src/linear_algebra/vector_math.cpp
similarity index 96%
rename from cpp/src/dual_simplex/vector_math.cpp
rename to cpp/src/linear_algebra/vector_math.cpp
index fceeec3d97..97f8e644a9 100644
--- a/cpp/src/dual_simplex/vector_math.cpp
+++ b/cpp/src/linear_algebra/vector_math.cpp
@@ -7,15 +7,15 @@
 
 #include <barrier/pinned_host_allocator.hpp>
 
-#include <dual_simplex/types.hpp>
-#include <dual_simplex/vector_math.hpp>
+#include <linear_algebra/vector_math.hpp>
+#include <math_optimization/types.hpp>
 
 #include <cassert>
 #include <cmath>
 #include <cstdio>
 #include <vector>
 
-namespace cuopt::mathematical_optimization::simplex {
+namespace cuopt::mathematical_optimization {
 
 template <typename i_t, typename f_t, typename Allocator>
 f_t vector_norm2_squared(const std::vector<f_t, Allocator>& x)
@@ -211,4 +211,4 @@ template int inverse_permutation<int>(const std::vector<int>& p, std::vector<int
 
 #endif
 
-}  // namespace cuopt::mathematical_optimization::simplex
+}  // namespace cuopt::mathematical_optimization
diff --git a/cpp/src/dual_simplex/vector_math.cuh b/cpp/src/linear_algebra/vector_math.cuh
similarity index 72%
rename from cpp/src/dual_simplex/vector_math.cuh
rename to cpp/src/linear_algebra/vector_math.cuh
index 6bb80fc981..ac9d24001b 100644
--- a/cpp/src/dual_simplex/vector_math.cuh
+++ b/cpp/src/linear_algebra/vector_math.cuh
@@ -12,12 +12,19 @@
 #include <raft/core/copy.hpp>
 #include <raft/core/device_span.hpp>
 #include <raft/core/host_span.hpp>
+#include <raft/util/cuda_rt_essentials.hpp>
 
 #include <rmm/device_scalar.hpp>
 #include <rmm/device_uvector.hpp>
+#include <rmm/exec_policy.hpp>
 
+#include <thrust/execution_policy.h>
+#include <thrust/functional.h>
+#include <thrust/transform_reduce.h>
+
+#include <cmath>
 #include <vector>
-namespace cuopt::mathematical_optimization::simplex {
+namespace cuopt::mathematical_optimization {
 
 struct norm_inf_max {
   template <typename f_t>
@@ -90,4 +97,36 @@ f_t vector_norm_inf(raft::host_span<const f_t> x, rmm::cuda_stream_view stream_v
   return device_vector_norm_inf<i_t, f_t>(d_x, stream_view);
 }
 
-}  // namespace cuopt::mathematical_optimization::simplex
+template <typename f_t>
+f_t vector_norm_inf(const rmm::device_uvector<f_t>& x)
+{
+  auto begin   = x.data();
+  auto end     = x.data() + x.size();
+  auto max_abs = thrust::transform_reduce(
+    rmm::exec_policy(x.stream()),
+    begin,
+    end,
+    [] __host__ __device__(f_t val) { return abs(val); },
+    static_cast<f_t>(0),
+    thrust::maximum<f_t>{});
+  RAFT_CHECK_CUDA(x.stream());
+  return max_abs;
+}
+
+template <typename f_t>
+f_t vector_norm2(const rmm::device_uvector<f_t>& x)
+{
+  auto begin          = x.data();
+  auto end            = x.data() + x.size();
+  auto sum_of_squares = thrust::transform_reduce(
+    rmm::exec_policy(x.stream()),
+    begin,
+    end,
+    [] __host__ __device__(f_t val) { return val * val; },
+    f_t(0),
+    thrust::plus<f_t>{});
+  RAFT_CHECK_CUDA(x.stream());
+  return std::sqrt(sum_of_squares);
+}
+
+}  // namespace cuopt::mathematical_optimization
diff --git a/cpp/src/dual_simplex/vector_math.hpp b/cpp/src/linear_algebra/vector_math.hpp
similarity index 95%
rename from cpp/src/dual_simplex/vector_math.hpp
rename to cpp/src/linear_algebra/vector_math.hpp
index a7e481c4b2..8063bc735b 100644
--- a/cpp/src/dual_simplex/vector_math.hpp
+++ b/cpp/src/linear_algebra/vector_math.hpp
@@ -11,7 +11,7 @@
 #include <cmath>
 #include <vector>
 
-namespace cuopt::mathematical_optimization::simplex {
+namespace cuopt::mathematical_optimization {
 
 // Computes || x ||_inf = max_j | x |_j
 template <typename i_t, typename f_t, typename Allocator>
@@ -69,4 +69,4 @@ i_t inverse_permute_vector(const std::vector<i_t>& p,
 template <typename i_t>
 i_t inverse_permutation(const std::vector<i_t>& p, std::vector<i_t>& pinv);
 
-}  // namespace cuopt::mathematical_optimization::simplex
+}  // namespace cuopt::mathematical_optimization
diff --git a/cpp/src/math_optimization/CMakeLists.txt b/cpp/src/math_optimization/CMakeLists.txt
index bdd644e6a6..efa1600c54 100644
--- a/cpp/src/math_optimization/CMakeLists.txt
+++ b/cpp/src/math_optimization/CMakeLists.txt
@@ -1,5 +1,5 @@
 # cmake-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
 # cmake-format: on
 
@@ -8,6 +8,7 @@ list(PREPEND
   ${CMAKE_CURRENT_SOURCE_DIR}/solver_settings.cu
   ${CMAKE_CURRENT_SOURCE_DIR}/solution_reader.cu
   ${CMAKE_CURRENT_SOURCE_DIR}/solution_writer.cu
+  ${CMAKE_CURRENT_SOURCE_DIR}/tic_toc.cpp
   )
 
 set(CUOPT_SRC_FILES ${CUOPT_SRC_FILES}
diff --git a/cpp/src/dual_simplex/tic_toc.cpp b/cpp/src/math_optimization/tic_toc.cpp
similarity index 73%
rename from cpp/src/dual_simplex/tic_toc.cpp
rename to cpp/src/math_optimization/tic_toc.cpp
index 5c27c4bc65..658536bb00 100644
--- a/cpp/src/dual_simplex/tic_toc.cpp
+++ b/cpp/src/math_optimization/tic_toc.cpp
@@ -5,11 +5,11 @@
  */
 /* clang-format on */
 
-#include <dual_simplex/tic_toc.hpp>
+#include <math_optimization/tic_toc.hpp>
 
 #include <sys/time.h>
 
-namespace cuopt::mathematical_optimization::simplex {
+namespace cuopt::mathematical_optimization {
 
 double tic()
 {
@@ -24,4 +24,4 @@ double toc(double start)
   return (now - start);
 }
 
-}  // namespace cuopt::mathematical_optimization::simplex
+}  // namespace cuopt::mathematical_optimization
diff --git a/cpp/src/dual_simplex/tic_toc.hpp b/cpp/src/math_optimization/tic_toc.hpp
similarity index 63%
rename from cpp/src/dual_simplex/tic_toc.hpp
rename to cpp/src/math_optimization/tic_toc.hpp
index 23463d63ae..e5ce719c74 100644
--- a/cpp/src/dual_simplex/tic_toc.hpp
+++ b/cpp/src/math_optimization/tic_toc.hpp
@@ -7,10 +7,10 @@
 
 #pragma once
 
-#include <dual_simplex/types.hpp>
+#include <math_optimization/types.hpp>
 
-namespace cuopt::mathematical_optimization::simplex {
+namespace cuopt::mathematical_optimization {
 double tic();
 double toc(double start);
 
-}  // namespace cuopt::mathematical_optimization::simplex
+}  // namespace cuopt::mathematical_optimization
diff --git a/cpp/src/dual_simplex/types.hpp b/cpp/src/math_optimization/types.hpp
similarity index 86%
rename from cpp/src/dual_simplex/types.hpp
rename to cpp/src/math_optimization/types.hpp
index df45c733c6..b3c4f408a6 100644
--- a/cpp/src/dual_simplex/types.hpp
+++ b/cpp/src/math_optimization/types.hpp
@@ -10,7 +10,7 @@
 #include <cstdint>
 #include <limits>
 
-namespace cuopt::mathematical_optimization::simplex {
+namespace cuopt::mathematical_optimization {
 
 #define DUAL_SIMPLEX_INSTANTIATE_DOUBLE
 
@@ -26,4 +26,4 @@ constexpr float64_t inf = std::numeric_limits<float64_t>::infinity();
 // We return this constant to signal that a matrix is indefinite (has a negative pivot)
 #define INDEFINITE_MATRIX_RETURN -4
 
-}  // namespace cuopt::mathematical_optimization::simplex
+}  // namespace cuopt::mathematical_optimization
diff --git a/cpp/src/mip_heuristics/diversity/lns/rins.cu b/cpp/src/mip_heuristics/diversity/lns/rins.cu
index a670322cb0..5ff96e00f0 100644
--- a/cpp/src/mip_heuristics/diversity/lns/rins.cu
+++ b/cpp/src/mip_heuristics/diversity/lns/rins.cu
@@ -23,7 +23,7 @@
 #include <mip_heuristics/presolve/trivial_presolve.cuh>
 
 #include <branch_and_bound/branch_and_bound.hpp>
-#include <dual_simplex/tic_toc.hpp>
+#include <math_optimization/tic_toc.hpp>
 #include <utilities/scope_guard.hpp>
 
 namespace cuopt::mathematical_optimization::mip {
@@ -264,7 +264,7 @@ void rins_t<i_t, f_t>::run_rins()
   };
   mip::probing_implied_bound_t<i_t, f_t> empty_probing(branch_and_bound_problem.num_cols);
   mip::branch_and_bound_t<i_t, f_t> branch_and_bound(
-    branch_and_bound_problem, branch_and_bound_settings, simplex::tic(), empty_probing);
+    branch_and_bound_problem, branch_and_bound_settings, tic(), empty_probing);
   branch_and_bound.set_initial_guess(cuopt::host_copy(fixed_assignment, rins_handle.get_stream()));
   branch_and_bound_status = branch_and_bound.solve(branch_and_bound_solution);
 
diff --git a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh
index 0f3625d9aa..7b8408f518 100644
--- a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh
+++ b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh
@@ -13,7 +13,7 @@
 #include <branch_and_bound/branch_and_bound.hpp>
 #include <dual_simplex/simplex_solver_settings.hpp>
 #include <dual_simplex/solve.hpp>
-#include <dual_simplex/tic_toc.hpp>
+#include <math_optimization/tic_toc.hpp>
 #include <pdlp/initial_scaling_strategy/initial_scaling.cuh>
 
 namespace cuopt::mathematical_optimization::mip {
@@ -122,7 +122,7 @@ class sub_mip_recombiner_t : public recombiner_t<i_t, f_t> {
       branch_and_bound_settings.log.log = false;
       mip::probing_implied_bound_t<i_t, f_t> empty_probing(branch_and_bound_problem.num_cols);
       mip::branch_and_bound_t<i_t, f_t> branch_and_bound(
-        branch_and_bound_problem, branch_and_bound_settings, simplex::tic(), empty_probing);
+        branch_and_bound_problem, branch_and_bound_settings, tic(), empty_probing);
       branch_and_bound_status = branch_and_bound.solve(branch_and_bound_solution);
       if (solution_vector.size() > 0) {
         cuopt_assert(fixed_assignment.size() == branch_and_bound_solution.x.size(),
diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu
index 04a312766c..ad4f584e60 100644
--- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu
+++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu
@@ -1442,7 +1442,7 @@ static std::unique_ptr<fj_cpu_climber_t<i_t, f_t>> init_fj_cpu_from_host_lp(
   const i_t n_variables   = problem.num_cols;
   const i_t n_constraints = problem.num_rows;
 
-  simplex::csr_matrix_t<i_t, f_t> csr_A(problem.num_rows, problem.num_cols, problem.A.nnz());
+  csr_matrix_t<i_t, f_t> csr_A(problem.num_rows, problem.num_cols, problem.A.nnz());
   problem.A.to_compressed_row(csr_A);
   std::vector<f_t> coefficients            = csr_A.x;
   std::vector<i_t> variables               = csr_A.j;
@@ -1469,7 +1469,7 @@ static std::unique_ptr<fj_cpu_climber_t<i_t, f_t>> init_fj_cpu_from_host_lp(
   }
 
   const i_t nnz = static_cast<i_t>(variables.size());
-  simplex::csc_matrix_t<i_t, f_t> reverse_csc(n_constraints, n_variables, nnz);
+  csc_matrix_t<i_t, f_t> reverse_csc(n_constraints, n_variables, nnz);
   csr_A.to_compressed_col(reverse_csc);
   std::vector<f_t> reverse_coefficients = std::move(reverse_csc.x);
   std::vector<i_t> reverse_constraints  = std::move(reverse_csc.i);
diff --git a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu
index bcf2d1a9dd..7ed6879296 100644
--- a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu
+++ b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu
@@ -21,9 +21,9 @@
 
 #include <algorithm>
 #include <cmath>
-#include <dual_simplex/sparse_matrix.hpp>
-#include <dual_simplex/sparse_vector.hpp>
 #include <limits>
+#include <linear_algebra/sparse_matrix.hpp>
+#include <linear_algebra/sparse_vector.hpp>
 #include <mip_heuristics/mip_constants.hpp>
 #include <mip_heuristics/utils.cuh>
 #include <utilities/logger.hpp>
@@ -32,7 +32,6 @@
 
 namespace cuopt::mathematical_optimization::mip {
 
-using simplex::csr_matrix_t;
 using simplex::user_problem_t;
 
 // do constraints with only binary variables.
diff --git a/cpp/src/mip_heuristics/presolve/semi_continuous.cu b/cpp/src/mip_heuristics/presolve/semi_continuous.cu
index ebd8d24a7a..33b7efff0e 100644
--- a/cpp/src/mip_heuristics/presolve/semi_continuous.cu
+++ b/cpp/src/mip_heuristics/presolve/semi_continuous.cu
@@ -87,7 +87,7 @@ std::vector<f_t> call_host_bounds_strengthening(const optimization_problem_t<i_t
   auto var_types = user_problem.var_types;
   var_types.resize(lp_problem.num_cols, simplex::variable_type_t::CONTINUOUS);
 
-  simplex::csr_matrix_t<i_t, f_t> Arow(1, 1, 1);
+  csr_matrix_t<i_t, f_t> Arow(1, 1, 1);
   lp_problem.A.to_compressed_row(Arow);
 
   // convert_user_problem returns an equality-form LP. Empty row_sense makes
diff --git a/cpp/src/mip_heuristics/problem/problem.cu b/cpp/src/mip_heuristics/problem/problem.cu
index 770c8c4623..e38e889495 100644
--- a/cpp/src/mip_heuristics/problem/problem.cu
+++ b/cpp/src/mip_heuristics/problem/problem.cu
@@ -19,7 +19,7 @@
 
 #include <cuts/objective_step.hpp>
 #include <cuts/rational.hpp>
-#include <dual_simplex/tic_toc.hpp>
+#include <math_optimization/tic_toc.hpp>
 #include <mip_heuristics/presolve/third_party_presolve.hpp>
 #include <mip_heuristics/presolve/trivial_presolve.cuh>
 #include <mip_heuristics/utils.cuh>
@@ -50,7 +50,6 @@
 
 namespace cuopt::mathematical_optimization::mip {
 
-using simplex::csr_matrix_t;
 using simplex::user_problem_t;
 using simplex::variable_type_t;
 
@@ -1419,7 +1418,7 @@ void problem_t<i_t, f_t>::recompute_objective_integrality()
 template <typename i_t, typename f_t>
 void problem_t<i_t, f_t>::compute_objective_step()
 {
-  f_t start_time = simplex::tic();
+  f_t start_time = tic();
   // Copy info from device to host
   auto h_obj_coefs = cuopt::host_copy(objective_coefficients, handle_ptr->get_stream());
   auto h_var_types = cuopt::host_copy(variable_types, handle_ptr->get_stream());
diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu
index a7cb4884e4..80f68ee500 100644
--- a/cpp/src/mip_heuristics/solve.cu
+++ b/cpp/src/mip_heuristics/solve.cu
@@ -8,6 +8,7 @@
 #include <cuopt/error.hpp>
 #include <cuopt/mathematical_optimization/solve_remote.hpp>
 
+#include <linear_algebra/sort_csr.cuh>
 #include <mip_heuristics/feasibility_jump/early_cpufj.cuh>
 #include <mip_heuristics/feasibility_jump/early_gpufj.cuh>
 #include <mip_heuristics/mip_constants.hpp>
@@ -16,7 +17,6 @@
 #include <mip_heuristics/presolve/third_party_presolve.hpp>
 #include <mip_heuristics/presolve/trivial_presolve.cuh>
 #include <mip_heuristics/solver.cuh>
-#include <mip_heuristics/utilities/sort_csr.cuh>
 #include <mip_heuristics/utils.cuh>
 
 #include <pdlp/pdlp.cuh>
@@ -545,7 +545,7 @@ mip_solution_t<i_t, f_t> solve_mip_helper(optimization_problem_t<i_t, f_t>& op_p
 
     auto constexpr const dual_postsolve = false;
     if (run_presolve) {
-      mip::sort_csr(op_problem);
+      sort_csr(op_problem);
       // allocate not more than 10% of the time limit to presolve.
       // Note that this is not the presolve time, but the time limit for presolve.
       const auto& hp = settings.heuristic_params;
diff --git a/cpp/src/pdlp/solve.cu b/cpp/src/pdlp/solve.cu
index 08b580b6a6..057fd22045 100644
--- a/cpp/src/pdlp/solve.cu
+++ b/cpp/src/pdlp/solve.cu
@@ -19,11 +19,11 @@
 #include <pdlp/utils.cuh>
 #include <utilities/logger.hpp>
 
+#include <linear_algebra/sort_csr.cuh>
 #include <mip_heuristics/mip_constants.hpp>
 #include <mip_heuristics/presolve/third_party_presolve.hpp>
 #include <mip_heuristics/presolve/trivial_presolve.cuh>
 #include <mip_heuristics/solver.cuh>
-#include <mip_heuristics/utilities/sort_csr.cuh>
 
 #include <cuopt/mathematical_optimization/backend_selection.hpp>
 #include <cuopt/mathematical_optimization/cpu_optimization_problem.hpp>
@@ -44,7 +44,7 @@
 
 #include <dual_simplex/crossover.hpp>
 #include <dual_simplex/solve.hpp>
-#include <dual_simplex/tic_toc.hpp>
+#include <math_optimization/tic_toc.hpp>
 #include <pdlp/utilities/problem_checking.cuh>
 
 #include <raft/sparse/detail/cusparse_wrappers.h>
@@ -491,8 +491,8 @@ std::tuple<simplex::lp_solution_t<i_t, f_t>, simplex::lp_status_t, f_t, f_t, f_t
   pdlp_solver_settings_t<i_t, f_t> const& settings,
   const timer_t& timer)
 {
-  f_t norm_user_objective = simplex::vector_norm2<i_t, f_t>(user_problem.objective);
-  f_t norm_rhs            = simplex::vector_norm2<i_t, f_t>(user_problem.rhs);
+  f_t norm_user_objective = vector_norm2<i_t, f_t>(user_problem.objective);
+  f_t norm_rhs            = vector_norm2<i_t, f_t>(user_problem.rhs);
 
   simplex::simplex_solver_settings_t<i_t, f_t> barrier_settings;
   barrier_settings.num_gpus                        = settings.num_gpus;
@@ -581,8 +581,8 @@ std::tuple<simplex::lp_solution_t<i_t, f_t>, simplex::lp_status_t, f_t, f_t, f_t
   pdlp_solver_settings_t<i_t, f_t> const& settings,
   const timer_t& timer)
 {
-  f_t norm_user_objective = simplex::vector_norm2<i_t, f_t>(user_problem.objective);
-  f_t norm_rhs            = simplex::vector_norm2<i_t, f_t>(user_problem.rhs);
+  f_t norm_user_objective = vector_norm2<i_t, f_t>(user_problem.objective);
+  f_t norm_rhs            = vector_norm2<i_t, f_t>(user_problem.rhs);
 
   simplex::simplex_solver_settings_t<i_t, f_t> dual_simplex_settings;
   dual_simplex_settings.time_limit      = settings.time_limit;
@@ -1947,7 +1947,7 @@ optimization_problem_solution_t<i_t, f_t> solve_lp(
     std::optional<mip::third_party_presolve_result_t<i_t, f_t>> result;
 
     if (run_presolve) {
-      mip::sort_csr(op_problem);
+      sort_csr(op_problem);
       // allocate no more than 10% of the time limit to presolve.
       // Note that this is not the presolve time, but the time limit for presolve.
       // But no less than 1 second, to avoid early timeout triggering known crashes
diff --git a/cpp/src/pdlp/translate.hpp b/cpp/src/pdlp/translate.hpp
index 3c6119e164..ab0d0e0c24 100644
--- a/cpp/src/pdlp/translate.hpp
+++ b/cpp/src/pdlp/translate.hpp
@@ -11,8 +11,8 @@
 #include <cuopt/mathematical_optimization/optimization_problem_interface.hpp>
 
 #include <dual_simplex/presolve.hpp>
-#include <dual_simplex/sparse_matrix.hpp>
-#include <dual_simplex/sparse_vector.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
+#include <linear_algebra/sparse_vector.hpp>
 
 #include <barrier/translate_soc.hpp>
 
@@ -40,7 +40,7 @@ static simplex::user_problem_t<i_t, f_t> cuopt_problem_to_user_problem(
   user_problem.num_cols  = n;
   user_problem.objective = problem.get_objective_coefficients_host();
 
-  simplex::csr_matrix_t<i_t, f_t> csr_A(m, n, static_cast<i_t>(A_values.size()));
+  csr_matrix_t<i_t, f_t> csr_A(m, n, static_cast<i_t>(A_values.size()));
   csr_A.x         = std::move(A_values);
   csr_A.j         = std::move(A_indices);
   csr_A.row_start = std::move(A_offsets);
@@ -112,7 +112,7 @@ static simplex::user_problem_t<i_t, f_t> cuopt_problem_to_user_problem(
   user_problem.num_cols  = n;
   user_problem.objective = cuopt::host_copy(model.objective_coefficients, handle_ptr->get_stream());
 
-  simplex::csr_matrix_t<i_t, f_t> csr_A(m, n, nz);
+  csr_matrix_t<i_t, f_t> csr_A(m, n, nz);
   csr_A.x = std::vector<f_t>(cuopt::host_copy(model.coefficients, handle_ptr->get_stream()));
   csr_A.j = std::vector<i_t>(cuopt::host_copy(model.variables, handle_ptr->get_stream()));
   csr_A.row_start = std::vector<i_t>(cuopt::host_copy(model.offsets, handle_ptr->get_stream()));
@@ -224,7 +224,7 @@ static simplex::user_problem_t<i_t, f_t> cuopt_optimization_problem_to_user_prob
     }
   }
 
-  simplex::csr_matrix_t<i_t, f_t> csr_A(m, n, nz);
+  csr_matrix_t<i_t, f_t> csr_A(m, n, nz);
   csr_A.x         = model.get_constraint_matrix_values_host();
   csr_A.j         = model.get_constraint_matrix_indices_host();
   csr_A.row_start = model.get_constraint_matrix_offsets_host();
@@ -333,7 +333,7 @@ void translate_to_crossover_problem(const mip::problem_t<i_t, f_t>& problem,
   auto stream                     = problem.handle_ptr->get_stream();
   std::vector<f_t> pdlp_objective = cuopt::host_copy(problem.objective_coefficients, stream);
 
-  simplex::csr_matrix_t<i_t, f_t> csr_A(problem.n_constraints, problem.n_variables, problem.nnz);
+  csr_matrix_t<i_t, f_t> csr_A(problem.n_constraints, problem.n_variables, problem.nnz);
   csr_A.x         = std::vector<f_t>(cuopt::host_copy(problem.coefficients, stream));
   csr_A.j         = std::vector<i_t>(cuopt::host_copy(problem.variables, stream));
   csr_A.row_start = std::vector<i_t>(cuopt::host_copy(problem.offsets, stream));
@@ -346,7 +346,7 @@ void translate_to_crossover_problem(const mip::problem_t<i_t, f_t>& problem,
   std::vector<f_t> slack(problem.n_constraints);
   std::vector<f_t> tmp_x = cuopt::host_copy(sol.get_primal_solution(), stream);
   stream.synchronize();
-  simplex::matrix_vector_multiply(lp.A, f_t(1.0), tmp_x, f_t(0.0), slack);
+  matrix_vector_multiply(lp.A, f_t(1.0), tmp_x, f_t(0.0), slack);
   CUOPT_LOG_DEBUG("Multiplied A and x");
 
   lp.A.col_start.resize(problem.n_variables + problem.n_constraints + 1);
diff --git a/cpp/tests/dual_simplex/unit_tests/right_looking_ldlt.cpp b/cpp/tests/dual_simplex/unit_tests/right_looking_ldlt.cpp
index 0c80dbe472..b03b773f73 100644
--- a/cpp/tests/dual_simplex/unit_tests/right_looking_ldlt.cpp
+++ b/cpp/tests/dual_simplex/unit_tests/right_looking_ldlt.cpp
@@ -9,8 +9,8 @@
 
 #include <dual_simplex/right_looking_lu.hpp>
 #include <dual_simplex/simplex_solver_settings.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 <cmath>
 #include <numeric>
diff --git a/cpp/tests/dual_simplex/unit_tests/solve.cpp b/cpp/tests/dual_simplex/unit_tests/solve.cpp
index 33d5392a9b..2e44442599 100644
--- a/cpp/tests/dual_simplex/unit_tests/solve.cpp
+++ b/cpp/tests/dual_simplex/unit_tests/solve.cpp
@@ -13,8 +13,8 @@
 
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/solve.hpp>
-#include <dual_simplex/tic_toc.hpp>
 #include <dual_simplex/user_problem.hpp>
+#include <math_optimization/tic_toc.hpp>
 
 #include <cuopt/mathematical_optimization/io/parser.hpp>
 #include <utilities/logger.hpp>
@@ -65,8 +65,8 @@ TEST(dual_simplex, chess_set)
   user_problem.lower[0] = 0;
   user_problem.lower[1] = 0.0;
   user_problem.upper.resize(n);
-  user_problem.upper[0]       = simplex::inf;
-  user_problem.upper[1]       = simplex::inf;
+  user_problem.upper[0]       = inf;
+  user_problem.upper[1]       = inf;
   user_problem.num_range_rows = 0;
   user_problem.problem_name   = "chess set";
   user_problem.row_names.resize(m);
@@ -313,8 +313,8 @@ TEST(dual_simplex, dual_variable_greater_than)
   user_problem.lower[1] = 0.0;
 
   user_problem.upper.resize(n);
-  user_problem.upper[0] = simplex::inf;
-  user_problem.upper[1] = simplex::inf;
+  user_problem.upper[0] = inf;
+  user_problem.upper[1] = inf;
 
   user_problem.num_range_rows = 0;
   user_problem.problem_name   = "dual_variable_greater_than";
diff --git a/cpp/tests/dual_simplex/unit_tests/solve_barrier.cu b/cpp/tests/dual_simplex/unit_tests/solve_barrier.cu
index 87ce592b82..76ad2a3e58 100644
--- a/cpp/tests/dual_simplex/unit_tests/solve_barrier.cu
+++ b/cpp/tests/dual_simplex/unit_tests/solve_barrier.cu
@@ -18,8 +18,8 @@
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/scaling.hpp>
 #include <dual_simplex/solve.hpp>
-#include <dual_simplex/tic_toc.hpp>
 #include <dual_simplex/user_problem.hpp>
+#include <math_optimization/tic_toc.hpp>
 
 #include <raft/sparse/detail/cusparse_wrappers.h>
 #include <raft/core/cusparse_macros.hpp>
@@ -84,8 +84,8 @@ TEST(barrier, chess_set)
   user_problem.lower[0] = 0;
   user_problem.lower[1] = 0.0;
   user_problem.upper.resize(n);
-  user_problem.upper[0]       = simplex::inf;
-  user_problem.upper[1]       = simplex::inf;
+  user_problem.upper[0]       = inf;
+  user_problem.upper[1]       = inf;
   user_problem.num_range_rows = 0;
   user_problem.problem_name   = "chess set";
   user_problem.row_names.resize(m);
@@ -163,8 +163,8 @@ TEST(barrier, dual_variable_greater_than)
   user_problem.lower[1] = 0.0;
 
   user_problem.upper.resize(n);
-  user_problem.upper[0] = simplex::inf;
-  user_problem.upper[1] = simplex::inf;
+  user_problem.upper[0] = inf;
+  user_problem.upper[1] = inf;
 
   user_problem.num_range_rows = 0;
   user_problem.problem_name   = "dual_variable_greater_than";
diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu
index 642f69a04e..4df3b788df 100644
--- a/cpp/tests/mip/cuts_test.cu
+++ b/cpp/tests/mip/cuts_test.cu
@@ -1349,7 +1349,7 @@ struct flow_cover_test_problem_t {
   raft::handle_t handle;
   simplex::simplex_solver_settings_t<int, double> settings;
   simplex::lp_problem_t<int, double> lp;
-  simplex::csr_matrix_t<int, double> Arow;
+  csr_matrix_t<int, double> Arow;
   std::vector<int> new_slacks;
   std::vector<simplex::variable_type_t> var_types;
 
diff --git a/cpp/tests/socp/general_quadratic_test.cu b/cpp/tests/socp/general_quadratic_test.cu
index 870c1b62ae..21c4e151c6 100644
--- a/cpp/tests/socp/general_quadratic_test.cu
+++ b/cpp/tests/socp/general_quadratic_test.cu
@@ -10,8 +10,8 @@
 #include <barrier/translate_soc.hpp>
 #include <cuopt/mathematical_optimization/optimization_problem_interface.hpp>
 #include <dual_simplex/solve.hpp>
-#include <dual_simplex/sparse_matrix.hpp>
 #include <dual_simplex/user_problem.hpp>
+#include <linear_algebra/sparse_matrix.hpp>
 
 #include <raft/sparse/detail/cusparse_wrappers.h>
 #include <raft/core/cusparse_macros.hpp>
@@ -21,8 +21,6 @@
 
 namespace cuopt::mathematical_optimization::barrier::test {
 
-using simplex::csr_matrix_t;
-using simplex::inf;
 using simplex::lp_solution_t;
 using simplex::lp_status_t;
 using simplex::simplex_solver_settings_t;