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

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

__device__

  assert(solver_status == CUSOLVER_STATUS_SUCCESS);

  auto dev_workspace = cuda::allocate<float>(workspace_size);
  auto dev_qr = cuda::allocate<float>(matrix.size());
  auto dev_q = cuda::allocate<float>(matrix.size());
  auto dev_temp = cuda::allocate<float>(matrix.size());

  cuda::copy_on_device(dev_qr, dev_matrix, matrix.size());

  for (auto iter = 0; iter < 50; ++iter) {
    // Compute QR factorization.
    solver_status = cusolverDnSgeqrf(solver_handle, n, n, dev_qr, 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);

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

    // --
    // Simulate a single iteration of traditional QR algorithm.
    constexpr auto alpha = 1.0f;
    constexpr 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_temp, n);
    assert(blas_status == CUBLAS_STATUS_SUCCESS);

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

    cuda::copy_on_device(dev_qr, dev_matrix, matrix.size());
  }

  std::cout << "\nEigenvalue matrix:\n";
  print_device_matrix(dev_matrix, n);
  // --

  cuda::free(dev_temp);
  cuda::free(dev_q);
  cuda::free(dev_qr);
  cuda::free(dev_matrix);
  cuda::free(dev_tau);
  cuda::free(dev_workspace);;
  cuda::free(dev_info);

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

  cuda::device_reset();

