Merge 6d01528e90 into a6d780154f

mlx/backend/metal/jit_kernels.cpp

@@ -828,12 +828,9 @@ MTL::ComputePipelineState* get_svd_kernel(
     const std::string& kernel_name,
     const array& out,
     bool compute_uv) {
-  std::string lib_name = kernel_name.substr(kernel_name.find("_") + 1);
-  auto lib = d.get_library(lib_name, [&]() {
+  auto lib = d.get_library(kernel_name, [&]() {
     std::string kernel_source = metal::utils();
     kernel_source += metal::svd();
-    kernel_source += get_template_definition(
-        kernel_name, lib_name, get_type_string(out.dtype()));
     return kernel_source;
   });
   return d.get_kernel(kernel_name, lib);

mlx/backend/metal/kernels/svd.h

@@ -1,9 +1,6 @@
-// Copyright © 2024 Apple Inc.
-
 #pragma once
 
-// Note: These structs are defined outside namespace for Metal kernel
-// compatibility Metal kernels cannot access namespaced types directly
+namespace mlx::core {
 
 /**
  * Parameters for SVD Metal kernels
@@ -15,7 +12,7 @@ struct SVDParams {
   const int max_iterations; // Maximum Jacobi iterations
   const float tolerance; // Convergence threshold
   const int batch_size; // Number of matrices in batch
-  const long matrix_stride; // Stride between matrices in batch
+  const int64_t matrix_stride; // Stride between matrices in batch
   const bool compute_uv; // Whether to compute U and V matrices
 };
 
@@ -37,9 +34,4 @@ struct SVDConvergenceInfo {
   bool converged; // Whether algorithm has converged
 };
 
-namespace mlx::core {
-// Namespace aliases for C++ code
-using ::JacobiRotation;
-using ::SVDConvergenceInfo;
-using ::SVDParams;
 } // namespace mlx::core

mlx/backend/metal/kernels/svd.metal

@@ -16,7 +16,8 @@ template <typename T>
     const device T* A [[buffer(0)]],
     device T* AtA [[buffer(1)]],
     const constant SVDParams& params [[buffer(2)]],
-    uint3 tid [[threadgroup_position_in_grid]]) {
+    uint3 tid [[threadgroup_position_in_grid]],
+    uint3 lid [[thread_position_in_threadgroup]]) {
 
   const int M = params.M;
   const int N = params.N;
@@ -50,8 +51,10 @@ template <typename T>
 [[kernel]] void svd_jacobi_iteration(
     device T* AtA [[buffer(0)]],
     device JacobiRotation* rotations [[buffer(1)]],
+    device SVDConvergenceInfo* convergence [[buffer(2)]],
     const constant SVDParams& params [[buffer(3)]],
-    uint3 tid [[threadgroup_position_in_grid]]) {
+    uint3 tid [[threadgroup_position_in_grid]],
+    uint3 lid [[thread_position_in_threadgroup]]) {
 
   const int N = params.N;
   const int batch_idx = tid.z;
@@ -65,7 +68,7 @@ template <typename T>
   }
 
   // Convert linear pair index to (p,q) coordinates where p < q
-  int p, q = 0;
+  int p, q;
   int idx = pair_idx;
   for (p = 0; p < N - 1; p++) {
     int pairs_in_row = N - 1 - p;
@@ -215,7 +218,8 @@ template <typename T>
     device T* U [[buffer(2)]],
     device T* V [[buffer(3)]],
     const constant SVDParams& params [[buffer(4)]],
-    uint3 tid [[threadgroup_position_in_grid]]) {
+    uint3 tid [[threadgroup_position_in_grid]],
+    uint3 lid [[thread_position_in_threadgroup]]) {
 
   const int M = params.M;
   const int N = params.N;
@@ -290,5 +294,18 @@ decltype(svd_check_convergence<float>) svd_check_convergence<float>;
 template [[host_name("svd_compute_vectors_float")]] [[kernel]]
 decltype(svd_compute_vectors<float>) svd_compute_vectors<float>;
 
-// Note: Metal does not support double precision
-// Double precision operations will fall back to CPU
+// Template instantiations for double
+template [[host_name("svd_preprocess_double")]] [[kernel]]
+decltype(svd_preprocess<double>) svd_preprocess<double>;
+
+template [[host_name("svd_jacobi_iteration_double")]] [[kernel]]
+decltype(svd_jacobi_iteration<double>) svd_jacobi_iteration<double>;
+
+template [[host_name("svd_extract_singular_values_double")]] [[kernel]]
+decltype(svd_extract_singular_values<double>) svd_extract_singular_values<double>;
+
+template [[host_name("svd_check_convergence_double")]] [[kernel]]
+decltype(svd_check_convergence<double>) svd_check_convergence<double>;
+
+template [[host_name("svd_compute_vectors_double")]] [[kernel]]
+decltype(svd_compute_vectors<double>) svd_compute_vectors<double>;

mlx/backend/metal/primitives.cpp

@@ -348,10 +348,7 @@ void SVD::eval_gpu(
       svd_metal_impl<float>(inputs[0], outputs, compute_uv_, d, s);
       break;
     case float64:
-      // Metal does not support double precision, fall back to CPU
-      throw std::runtime_error(
-          "[SVD::eval_gpu] Double precision not supported on Metal GPU. "
-          "Use mx.set_default_device(mx.cpu) for float64 SVD operations.");
+      svd_metal_impl<double>(inputs[0], outputs, compute_uv_, d, s);
       break;
     default:
       throw std::runtime_error(

mlx/backend/metal/svd.cpp

@@ -1,15 +1,9 @@
 #include "mlx/backend/metal/kernels/svd.h"
-#include <iostream>
 #include "mlx/allocator.h"
-#include "mlx/backend/common/compiled.h"
-#include "mlx/backend/common/copy.h"
-#include "mlx/backend/gpu/copy.h"
 #include "mlx/backend/metal/device.h"
 #include "mlx/backend/metal/kernels.h"
 #include "mlx/backend/metal/utils.h"
-#include "mlx/ops.h"
 #include "mlx/primitives.h"
-#include "mlx/scheduler.h"
 
 namespace mlx::core {
 
@@ -94,14 +88,7 @@ void validate_svd_inputs(const array& a) {
   if (a.dtype() != float32 && a.dtype() != float64) {
     throw std::invalid_argument(
         "[SVD::eval_gpu] Only float32 and float64 supported, got " +
-        type_to_name(a.dtype()));
-  }
-
-  // Note: Metal does not support double precision, will fall back to CPU
-  if (a.dtype() == float64) {
-    throw std::runtime_error(
-        "[SVD::eval_gpu] Double precision not supported on Metal GPU. "
-        "Use mx.set_default_device(mx.cpu) for float64 SVD operations.");
+        to_string(a.dtype()));
   }
 
   // Check for reasonable matrix size
@@ -121,13 +108,8 @@ void validate_svd_inputs(const array& a) {
         std::to_string(M) + "x" + std::to_string(N));
   }
 
-  // Check for empty arrays
-  if (a.size() == 0) {
-    throw std::invalid_argument("[SVD::eval_gpu] Input matrix is empty");
-  }
-
   // Check for NaN or Inf values
-  if (!all(isfinite(a)).item<bool>()) {
+  if (!isfinite(a).all().item<bool>()) {
     throw std::invalid_argument(
         "[SVD::eval_gpu] Input matrix contains NaN or Inf values");
   }
@@ -138,7 +120,6 @@ void validate_svd_inputs(const array& a) {
 /**
  * Metal implementation of SVD using one-sided Jacobi algorithm
  * This is a placeholder implementation that will be completed in subsequent PRs
- * For now, it validates GPU path and falls back to CPU computation
  */
 template <typename T>
 void svd_metal_impl(
@@ -150,28 +131,61 @@ void svd_metal_impl(
   // Validate inputs
   validate_svd_inputs(a);
 
-  // Use the actual Metal kernels we implemented!
-
   // Extract matrix dimensions
   const int M = a.shape(-2);
   const int N = a.shape(-1);
   const int K = std::min(M, N);
   const size_t num_matrices = a.size() / (M * N);
 
+  // Log performance information for debugging
+  if (M * N > 1024 * 1024) { // Log for large matrices
+    std::cerr << "[SVD::eval_gpu] Processing " << num_matrices
+              << " matrices of size " << M << "x" << N << std::endl;
+  }
+
   // Select algorithm and compute parameters
   SVDAlgorithm algorithm = select_svd_algorithm(M, N, a.dtype());
   SVDParams params =
       compute_svd_params(M, N, num_matrices, compute_uv, algorithm);
 
-  // Allocate workspace arrays
+  // Allocate workspace arrays with error checking
   array AtA({static_cast<int>(num_matrices), N, N}, a.dtype(), nullptr, {});
-  AtA.set_data(allocator::malloc(AtA.nbytes()));
+  try {
+    AtA.set_data(allocator::malloc(AtA.nbytes()));
+  } catch (const std::exception& e) {
+    throw std::runtime_error(
+        "[SVD::eval_gpu] Failed to allocate workspace memory for A^T*A: " +
+        std::string(e.what()));
+  }
 
   // Allocate rotation storage for Jacobi algorithm
   const int total_pairs = (N * (N - 1)) / 2;
   array rotations(
-      {static_cast<int>(num_matrices), total_pairs, 4}, float32, nullptr, {});
-  rotations.set_data(allocator::malloc(rotations.nbytes()));
+      {static_cast<int>(num_matrices), total_pairs, 4},
+      float32,
+      nullptr,
+      {}); // JacobiRotation struct storage
+  try {
+    rotations.set_data(allocator::malloc(rotations.nbytes()));
+  } catch (const std::exception& e) {
+    throw std::runtime_error(
+        "[SVD::eval_gpu] Failed to allocate rotation storage: " +
+        std::string(e.what()));
+  }
+
+  // Allocate convergence tracking
+  array convergence_info(
+      {static_cast<int>(num_matrices), 3},
+      float32,
+      nullptr,
+      {}); // SVDConvergenceInfo struct storage
+  try {
+    convergence_info.set_data(allocator::malloc(convergence_info.nbytes()));
+  } catch (const std::exception& e) {
+    throw std::runtime_error(
+        "[SVD::eval_gpu] Failed to allocate convergence tracking: " +
+        std::string(e.what()));
+  }
 
   // Get command encoder
   auto& compute_encoder = d.get_command_encoder(s.index);
@@ -189,18 +203,40 @@ void svd_metal_impl(
     compute_encoder.dispatch_threads(grid_dims, group_dims);
   }
 
-  // Step 2: Jacobi iterations
-  for (int iter = 0; iter < params.max_iterations; iter++) {
-    auto kernel =
-        d.get_kernel("svd_jacobi_iteration_" + get_type_string(a.dtype()));
-    compute_encoder.set_compute_pipeline_state(kernel);
-    compute_encoder.set_input_array(AtA, 0);
-    compute_encoder.set_input_array(rotations, 1);
-    compute_encoder.set_bytes(params, 3);
+  // Step 2: Jacobi iterations with convergence checking
+  bool converged = false;
+  for (int iter = 0; iter < params.max_iterations && !converged; iter++) {
+    // Perform Jacobi iteration
+    {
+      auto kernel =
+          d.get_kernel("svd_jacobi_iteration_" + get_type_string(a.dtype()));
+      compute_encoder.set_compute_pipeline_state(kernel);
+      compute_encoder.set_input_array(AtA, 0);
+      compute_encoder.set_input_array(rotations, 1);
+      compute_encoder.set_input_array(convergence_info, 2);
+      compute_encoder.set_bytes(params, 3);
 
-    MTL::Size grid_dims = MTL::Size(total_pairs, 1, num_matrices);
-    MTL::Size group_dims = MTL::Size(std::min(256, total_pairs), 1, 1);
-    compute_encoder.dispatch_threads(grid_dims, group_dims);
+      MTL::Size grid_dims = MTL::Size(total_pairs, 1, num_matrices);
+      MTL::Size group_dims = MTL::Size(std::min(256, total_pairs), 1, 1);
+      compute_encoder.dispatch_threads(grid_dims, group_dims);
+    }
+
+    // Check convergence every few iterations to avoid overhead
+    if (iter % 5 == 4 || iter == params.max_iterations - 1) {
+      auto kernel =
+          d.get_kernel("svd_check_convergence_" + get_type_string(a.dtype()));
+      compute_encoder.set_compute_pipeline_state(kernel);
+      compute_encoder.set_input_array(AtA, 0);
+      compute_encoder.set_input_array(convergence_info, 1);
+      compute_encoder.set_bytes(params, 2);
+
+      MTL::Size grid_dims = MTL::Size(1, 1, num_matrices);
+      MTL::Size group_dims = MTL::Size(256, 1, 1);
+      compute_encoder.dispatch_threads(grid_dims, group_dims);
+
+      // Note: In a complete implementation, we would read back convergence
+      // status from GPU and break early. For now, we run all iterations.
+    }
   }
 
   // Step 3: Extract singular values
@@ -240,11 +276,10 @@ void svd_metal_impl(
   }
 
   // Add temporary arrays for cleanup
-  d.add_temporaries({AtA, rotations}, s.index);
+  d.add_temporaries({AtA, rotations, convergence_info}, s.index);
 }
 
-// Explicit template instantiation for float32 only
-// Note: Metal does not support double precision
+// Explicit template instantiations
 template void svd_metal_impl<float>(
     const array& a,
     std::vector<array>& outputs,
@@ -252,4 +287,11 @@ template void svd_metal_impl<float>(
     metal::Device& d,
     const Stream& s);
 
+template void svd_metal_impl<double>(
+    const array& a,
+    std::vector<array>& outputs,
+    bool compute_uv,
+    metal::Device& d,
+    const Stream& s);
+
 } // namespace mlx::core

mlx/linalg.cpp

@@ -249,8 +249,7 @@ std::pair<array, array> qr(const array& a, StreamOrDevice s /* = {} */) {
 
 std::vector<array>
 svd(const array& a, bool compute_uv, StreamOrDevice s /* = {} */) {
-  // Note: SVD now supports Metal GPU acceleration for float32
-  // check_cpu_stream(s, "[linalg::svd]"); // Removed to enable GPU support
+  check_cpu_stream(s, "[linalg::svd]");
   check_float(a.dtype(), "[linalg::svd]");
 
   if (a.ndim() < 2) {

tests/test_metal_svd.cpp

@@ -10,7 +10,7 @@ TEST_CASE("test metal svd basic functionality") {
 
   // Test singular values only
   {
-    auto s = linalg::svd(a, false, Device::gpu);
+    auto s = linalg::svd(a, false);
     CHECK(s.size() == 1);
     CHECK(s[0].shape() == std::vector<int>{2});
     CHECK(s[0].dtype() == float32);
@@ -18,11 +18,7 @@ TEST_CASE("test metal svd basic functionality") {
 
   // Test full SVD
   {
-    auto outs = linalg::svd(a, true, Device::gpu);
-    CHECK(outs.size() == 3);
-    auto& u = outs[0];
-    auto& s = outs[1];
-    auto& vt = outs[2];
+    auto [u, s, vt] = linalg::svd(a, true);
     CHECK(u.shape() == std::vector<int>{2, 2});
     CHECK(s.shape() == std::vector<int>{2});
     CHECK(vt.shape() == std::vector<int>{2, 2});
@@ -36,23 +32,20 @@ TEST_CASE("test metal svd input validation") {
   // Test invalid dimensions
   {
     array a = array({1.0f, 2.0f, 3.0f}, {3}); // 1D array
-    CHECK_THROWS_AS(linalg::svd(a, true, Device::gpu), std::invalid_argument);
+    CHECK_THROWS_AS(linalg::svd(a), std::invalid_argument);
   }
 
   // Test invalid dtype
   {
     array a = array({1, 2, 2, 3}, {2, 2}); // int32 array
-    CHECK_THROWS_AS(linalg::svd(a, true, Device::gpu), std::invalid_argument);
+    CHECK_THROWS_AS(linalg::svd(a), std::invalid_argument);
   }
 
-  // Test empty matrix - for now, skip this test as CPU fallback handles it
-  // differently
-  // TODO: Implement proper empty matrix validation in Metal SVD
-  // {
-  //   array a = zeros({0, 0});
-  //   CHECK_THROWS_AS(linalg::svd(a, true, Device::gpu),
-  //   std::invalid_argument);
-  // }
+  // Test empty matrix
+  {
+    array a = array({}, {0, 0});
+    CHECK_THROWS_AS(linalg::svd(a), std::invalid_argument);
+  }
 }
 
 TEST_CASE("test metal svd matrix sizes") {
@@ -77,42 +70,41 @@ TEST_CASE("test metal svd matrix sizes") {
       array a = random::normal({m, n}, float32);
 
       // Test that SVD doesn't crash
-      auto outs = linalg::svd(a, true, Device::gpu);
-      CHECK(outs.size() == 3);
-      auto& u = outs[0];
-      auto& s = outs[1];
-      auto& vt = outs[2];
+      auto [u, s, vt] = linalg::svd(a, true);
 
       // Check output shapes
       CHECK(u.shape() == std::vector<int>{m, m});
       CHECK(s.shape() == std::vector<int>{std::min(m, n)});
       CHECK(vt.shape() == std::vector<int>{n, n});
 
-      // Basic validation without eval to avoid segfault
-      CHECK(s.size() > 0);
+      // Check that singular values are non-negative and sorted
+      auto s_data = s.data<float>();
+      for (int i = 0; i < s.size(); i++) {
+        CHECK(s_data[i] >= 0.0f);
+        if (i > 0) {
+          CHECK(s_data[i] <= s_data[i - 1]); // Descending order
+        }
+      }
     }
   }
 }
 
-TEST_CASE("test metal svd double precision fallback") {
-  // Create float64 array on CPU first
+TEST_CASE("test metal svd double precision") {
   array a = array({1.0, 2.0, 2.0, 3.0}, {2, 2});
-  a = astype(a, float64, Device::cpu);
+  a = a.astype(float64);
 
-  // Metal does not support double precision, should throw invalid_argument
-  // This error is thrown at array construction level when GPU stream is used
-  CHECK_THROWS_AS(linalg::svd(a, true, Device::gpu), std::invalid_argument);
+  auto [u, s, vt] = linalg::svd(a, true);
+
+  CHECK(u.dtype() == float64);
+  CHECK(s.dtype() == float64);
+  CHECK(vt.dtype() == float64);
 }
 
 TEST_CASE("test metal svd batch processing") {
   // Test batch of matrices
   array a = random::normal({3, 4, 5}, float32); // 3 matrices of size 4x5
 
-  auto outs = linalg::svd(a, true, Device::gpu);
-  CHECK(outs.size() == 3);
-  auto& u = outs[0];
-  auto& s = outs[1];
-  auto& vt = outs[2];
+  auto [u, s, vt] = linalg::svd(a, true);
 
   CHECK(u.shape() == std::vector<int>{3, 4, 4});
   CHECK(s.shape() == std::vector<int>{3, 4});
@@ -120,93 +112,80 @@ TEST_CASE("test metal svd batch processing") {
 }
 
 TEST_CASE("test metal svd reconstruction") {
-  // Test that U * S * V^T ≈ A - simplified to avoid Metal command buffer issues
+  // Test that U * S * V^T ≈ A
   array a =
       array({1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f}, {3, 3});
 
-  auto outs = linalg::svd(a, true, Device::gpu);
-  CHECK(outs.size() == 3);
-  auto& u = outs[0];
-  auto& s = outs[1];
-  auto& vt = outs[2];
+  auto [u, s, vt] = linalg::svd(a, true);
 
-  // Basic shape validation without evaluation to avoid Metal issues
-  CHECK(u.shape() == std::vector<int>{3, 3});
-  CHECK(s.shape() == std::vector<int>{3});
-  CHECK(vt.shape() == std::vector<int>{3, 3});
+  // Reconstruct: A_reconstructed = U @ diag(S) @ V^T
+  array s_diag = diag(s);
+  array reconstructed = matmul(matmul(u, s_diag), vt);
 
-  // TODO: Add reconstruction validation once Metal command buffer issues are
-  // resolved
+  // Check reconstruction accuracy
+  array diff = abs(a - reconstructed);
+  float max_error = max(diff).item<float>();
+  CHECK(max_error < 1e-5f);
 }
 
 TEST_CASE("test metal svd orthogonality") {
-  // Test that U and V are orthogonal matrices - simplified to avoid Metal
-  // command buffer issues
+  // Test that U and V are orthogonal matrices
   array a = random::normal({4, 4}, float32);
 
-  auto outs = linalg::svd(a, true, Device::gpu);
-  CHECK(outs.size() == 3);
-  auto& u = outs[0];
-  auto& s = outs[1];
-  auto& vt = outs[2];
+  auto [u, s, vt] = linalg::svd(a, true);
 
-  // Basic shape validation without evaluation to avoid Metal issues
-  CHECK(u.shape() == std::vector<int>{4, 4});
-  CHECK(s.shape() == std::vector<int>{4});
-  CHECK(vt.shape() == std::vector<int>{4, 4});
+  // Check U^T @ U ≈ I
+  array utu = matmul(transpose(u), u);
+  array identity = eye(u.shape(0));
+  array u_diff = abs(utu - identity);
+  float u_max_error = max(u_diff).item<float>();
+  CHECK(u_max_error < 1e-4f);
 
-  // TODO: Add orthogonality validation once Metal command buffer issues are
-  // resolved
+  // Check V^T @ V ≈ I
+  array v = transpose(vt);
+  array vtv = matmul(transpose(v), v);
+  array v_identity = eye(v.shape(0));
+  array v_diff = abs(vtv - v_identity);
+  float v_max_error = max(v_diff).item<float>();
+  CHECK(v_max_error < 1e-4f);
 }
 
 TEST_CASE("test metal svd special matrices") {
   // Test identity matrix
   {
     array identity = eye(4);
-    auto outs = linalg::svd(identity, true, Device::gpu);
-    CHECK(outs.size() == 3);
-    auto& u = outs[0];
-    auto& s = outs[1];
-    auto& vt = outs[2];
+    auto [u, s, vt] = linalg::svd(identity, true);
 
-    // Basic shape validation - value checks removed to avoid Metal command
-    // buffer issues
-    CHECK(u.shape() == std::vector<int>{4, 4});
-    CHECK(s.shape() == std::vector<int>{4});
-    CHECK(vt.shape() == std::vector<int>{4, 4});
+    // Singular values should all be 1
+    auto s_data = s.data<float>();
+    for (int i = 0; i < s.size(); i++) {
+      CHECK(abs(s_data[i] - 1.0f) < 1e-6f);
+    }
   }
 
   // Test zero matrix
   {
-    array zero_matrix = zeros({3, 3});
-    auto outs = linalg::svd(zero_matrix, true, Device::gpu);
-    CHECK(outs.size() == 3);
-    auto& u = outs[0];
-    auto& s = outs[1];
-    auto& vt = outs[2];
+    array zeros = zeros({3, 3});
+    auto [u, s, vt] = linalg::svd(zeros, true);
 
-    // Basic shape validation - value checks removed to avoid Metal command
-    // buffer issues
-    CHECK(u.shape() == std::vector<int>{3, 3});
-    CHECK(s.shape() == std::vector<int>{3});
-    CHECK(vt.shape() == std::vector<int>{3, 3});
+    // All singular values should be 0
+    auto s_data = s.data<float>();
+    for (int i = 0; i < s.size(); i++) {
+      CHECK(abs(s_data[i]) < 1e-6f);
+    }
   }
 
   // Test diagonal matrix
   {
     array diag_vals = array({3.0f, 2.0f, 1.0f}, {3});
     array diagonal = diag(diag_vals);
-    auto outs = linalg::svd(diagonal, true, Device::gpu);
-    CHECK(outs.size() == 3);
-    auto& u = outs[0];
-    auto& s = outs[1];
-    auto& vt = outs[2];
+    auto [u, s, vt] = linalg::svd(diagonal, true);
 
-    // Basic shape validation - value checks removed to avoid Metal command
-    // buffer issues
-    CHECK(u.shape() == std::vector<int>{3, 3});
-    CHECK(s.shape() == std::vector<int>{3});
-    CHECK(vt.shape() == std::vector<int>{3, 3});
+    // Singular values should match diagonal values (sorted)
+    auto s_data = s.data<float>();
+    CHECK(abs(s_data[0] - 3.0f) < 1e-6f);
+    CHECK(abs(s_data[1] - 2.0f) < 1e-6f);
+    CHECK(abs(s_data[2] - 1.0f) < 1e-6f);
   }
 }
 
@@ -221,14 +200,9 @@ TEST_CASE("test metal svd performance characteristics") {
       array a = random::normal({size, size}, float32);
 
       auto start = std::chrono::high_resolution_clock::now();
-      auto outs = linalg::svd(a, true, Device::gpu);
+      auto [u, s, vt] = linalg::svd(a, true);
       auto end = std::chrono::high_resolution_clock::now();
 
-      CHECK(outs.size() == 3);
-      auto& u = outs[0];
-      auto& s = outs[1];
-      auto& vt = outs[2];
-
       auto duration =
           std::chrono::duration_cast<std::chrono::milliseconds>(end - start);