// HOW TO BUILD: nvcc -O3 -std=c++17 -Xcompiler -fopenmp cub_sort_test.cu -o cub_sort_test
#include <cub/cub.cuh>
#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <thrust/system/cuda/experimental/pinned_allocator.h>
#include <algorithm>
#include <chrono>
#include <iostream>
#include <parallel/algorithm>
#include <random>
#include <vector>
#include <iostream>
#define DEBUG
#ifdef DEBUG
#define CheckCudaError(instruction) \
  { AssertNoCudaError((instruction), __FILE__, __LINE__); }
#else
#define CheckCudaError(instruction) instruction
#endif
inline void AssertNoCudaError(cudaError_t error_code, const char* file, int line) {
  if (error_code != cudaSuccess) {
    std::cout << "Error: " << cudaGetErrorString(error_code) << " " << file << " " << line << "\n";
template <typename T>
using PinnedHostVector = thrust::host_vector<T, thrust::system::cuda::experimental::pinned_allocator<T>>;
std::mt19937 SeedRandomGenerator(uint32_t distribution_seed) {
    const size_t seeds_bytes = sizeof(std::mt19937::result_type) * std::mt19937::state_size;
    const size_t seeds_length = seeds_bytes / sizeof(std::seed_seq::result_type);
    std::vector<std::seed_seq::result_type> seeds(seeds_length);
    std::generate(seeds.begin(), seeds.end(), [&]() {
        distribution_seed = (distribution_seed << 1) | (distribution_seed >> (-1 & 31));
        return distribution_seed;
    std::seed_seq seed_sequence(seeds.begin(), seeds.end());
    return std::mt19937{seed_sequence};
int main(int argc, char* argv[]) {
    if (argc != 4) {
        std::cerr << "Usage: ./cub-sort-test <num_keys> <gpu_id> <bit_entropy>" << std::endl;
        return -1;
    size_t num_keys = std::stoull(argv[1]);
    int gpu = std::stoi(argv[2]);
    size_t bit_entropy = std::stoi(argv[3]);
    cudaStream_t stream;
    CheckCudaError(cudaSetDevice(gpu));
    CheckCudaError(cudaStreamCreateWithFlags(&stream, cudaStreamNonBlocking));
    PinnedHostVector<uint32_t> keys(num_keys);
#pragma omp parallel num_threads(64)
        uint32_t max = (1 << bit_entropy) - 1;
      if (bit_entropy == sizeof(uint32_t) * 8) {
        max = std::numeric_limits<uint32_t>::max();
      } else if (bit_entropy == 1) {
        max = 2;
      std::mt19937 random_generator = SeedRandomGenerator(2147483647 + static_cast<size_t>(omp_get_thread_num()));
      std::uniform_real_distribution<double> uniform_dist(0, max);
#pragma omp for schedule(static)
      for (size_t i = 0; i < num_keys; ++i) {
        keys[i] = static_cast<uint32_t>(uniform_dist(random_generator));
    thrust::device_vector<uint32_t> device_vector(num_keys);
    thrust::copy(keys.begin(), keys.end(), device_vector.begin());
    CheckCudaError(cudaDeviceSynchronize());
    size_t num_temporary_bytes = 0;
    cub::DeviceRadixSort::SortKeys(
        NULL, num_temporary_bytes, thrust::raw_pointer_cast(device_vector.data()),
        thrust::raw_pointer_cast(device_vector.data()), num_keys, 0, bit_entropy + 1, stream); // bit subrange is [begin_bit, end_bit), thus bit_entropy + 1
    uint8_t* temporary_storage = nullptr;
    CheckCudaError(cudaMalloc(reinterpret_cast<void**>(&temporary_storage), num_temporary_bytes));
    cub::DeviceRadixSort::SortKeys(
    (void*)temporary_storage, num_temporary_bytes, thrust::raw_pointer_cast(device_vector.data()),
    thrust::raw_pointer_cast(device_vector.data()), num_keys, 0, bit_entropy + 1, stream);
    CheckCudaError(cudaStreamSynchronize(stream));
    thrust::copy(device_vector.begin(), device_vector.end(), keys.begin());
    CheckCudaError(cudaFree(temporary_storage));
    if (std::is_sorted(keys.begin(), keys.end()) == false) {
        std::cout << "Error: Sorting failed." << std::endl;
    return 0;
}