   assert(blas_status == CUBLAS_STATUS_SUCCESS);
 
   // Allocate device memory.
+  // const auto dev_input = cuda::copy_to_device(matrix);
+  // auto dev_matrix = cuda::copy_on_device(dev_input, matrix.size());
   auto dev_matrix = cuda::copy_to_device(matrix);
   auto dev_tau = cuda::allocate<float>(n);
   auto dev_info = cuda::allocate<int>(1);

  assert(solver_status == CUSOLVER_STATUS_SUCCESS);

  auto dev_workspace = cuda::allocate<float>(workspace_size);

  // Compute QR factorization.
  solver_status = cusolverDnSgeqrf(solver_handle, n, n, dev_matrix, n, dev_tau, dev_workspace, workspace_size, dev_info);
  assert(solver_status == CUSOLVER_STATUS_SUCCESS);

  int info;
  cuda::copy_to_host(&info, dev_info, 1);
  assert(info == 0);

  const dim3 blocks(2, 2);
  const dim3 threads(2, 2);

  // Allocate device memory.
  auto dev_q = cuda::allocate<float>(n * n);

  // Execute kernel.
  construct_q_matrix<<<blocks, threads>>>(dev_q, dev_matrix, dev_tau, n);
  cuda::device_sync();

  // Allocate device memory.
  auto dev_qq = cuda::allocate<float>(n * n);

  // Copy data to device memory.
  cuda::copy_to_device(dev_matrix, matrix.data(), matrix.size());

  // --
  // Simulate a single iteration of traditional QR algorithm.
  const auto alpha = 1.0f;
  const auto beta = 0.0f;
  blas_status = cublasSgemm(blas_handle, CUBLAS_OP_T, CUBLAS_OP_N, n, n, n, &alpha, dev_q, n, dev_matrix, n, &beta, dev_qq, n);
  assert(blas_status == CUBLAS_STATUS_SUCCESS);

  blas_status = cublasSgemm(blas_handle, CUBLAS_OP_N, CUBLAS_OP_N, n, n, n, &alpha, dev_qq, n, dev_q, n, &beta, dev_qq, n);
  assert(blas_status == CUBLAS_STATUS_SUCCESS);

  std::vector<float> result(n * n);
  cuda::copy_to_host(result, dev_qq);
  std::cout << "\nResult:\n";
  print_matrix(result, n);
  // --

  std::vector<float> tau(n);
  cuda_status = cudaMemcpy(tau.data(), dev_tau, sizeof(float) * n, cudaMemcpyDeviceToHost);
  assert(cuda_status == cudaSuccess);
  std::cout << "\nTAU: ";
  for (const auto& t : tau) std::cout << t << ", ";
  std::cout << "\n";

  cuda::free(dev_qq);
  cuda::free(dev_q);
  cuda::free(dev_matrix);
  cuda::free(dev_tau);
  cuda::free(dev_workspace);;
  cuda::free(dev_info);

  if (solver_handle) cusolverDnDestroy(solver_handle);
  if (blas_handle)   cublasDestroy(blas_handle);

  cuda::device_reset();


    }
    std::cout << "\n";
  }
}

template<typename T>
void print_device_matrix(T* dev_ptr, int n)
{
  std::vector<float> v(n * n);
  cuda::copy_to_host(v, dev_ptr);
  print_matrix(v, n);
}

__device__
