 }
 
 void egblas_par_shuffle_seed(size_t n, void* x, size_t incx, void* y, size_t incy, size_t seed){
-    // Allocate memory for temporary for swapping
-    void* tx;
-    void* ty;
-    cuda_check(cudaMalloc((void**)&tx, incx));
-    cuda_check(cudaMalloc((void**)&ty, incy));
-
     std::default_random_engine g(seed);
 
     using distribution_t = typename std::uniform_int_distribution<size_t>;


    distribution_t dist;

    uint8_t* x_flat = reinterpret_cast<uint8_t*>(x);
    uint8_t* y_flat = reinterpret_cast<uint8_t*>(y);

    for (auto i = n - 1; i > 0; --i) {
        auto new_i = dist(g, param_t(0, i));

        // Swap x[i] and x[new_i]
        cuda_check(cudaMemcpy(tx, x_flat + i * incx, incx, cudaMemcpyDeviceToDevice));
        cuda_check(cudaMemcpy(x_flat + i * incx, x_flat + new_i * incx, incx, cudaMemcpyDeviceToDevice));
        cuda_check(cudaMemcpy(x_flat + new_i * incx, tx, incx, cudaMemcpyDeviceToDevice));

        // Swap y[i] and y[new_i]
        cuda_check(cudaMemcpy(ty, y_flat + i * incy, incy, cudaMemcpyDeviceToDevice));
        cuda_check(cudaMemcpy(y_flat + i * incy, y_flat + new_i * incy, incy, cudaMemcpyDeviceToDevice));
        cuda_check(cudaMemcpy(y_flat + new_i * incy, ty, incy, cudaMemcpyDeviceToDevice));
    }

    // Release the tmp memory
    cuda_check(cudaFree(tx));
    cuda_check(cudaFree(ty));
}

void egblas_par_shuffle(size_t n, void* x, size_t incx, void* y, size_t incy){


    distribution_t dist;

    // Optimized version for arrays of 4B
    if (incx == 4) {
        uint32_t* x_flat = reinterpret_cast<uint32_t*>(x);

        for (auto i = n - 1; i > 0; --i) {
            auto new_i = dist(g, param_t(0, i));

            shuffle_one_4_kernel<<<1, 1, 1>>>(x_flat, i, new_i);
        }

        cudaDeviceSynchronize();

        return;
    }

    if (incx % 8 == 0) {
        uint64_t* x_flat = reinterpret_cast<uint64_t*>(x);

        for (auto i = n - 1; i > 0; --i) {
            auto new_i = dist(g, param_t(0, i));

            shuffle_kernel_run(x_flat, size_t(i), new_i, incx / 8);
        }
    } else if (incx % 4 == 0) {
        uint32_t* x_flat = reinterpret_cast<uint32_t*>(x);

        for (auto i = n - 1; i > 0; --i) {
            auto new_i = dist(g, param_t(0, i));

            shuffle_kernel_run(x_flat, size_t(i), new_i, incx / 4);
        }
    } else {
        uint8_t* x_flat = reinterpret_cast<uint8_t*>(x);

        for (auto i = n - 1; i > 0; --i) {
            auto new_i = dist(g, param_t(0, i));

            shuffle_kernel_run(x_flat, size_t(i), new_i, incx);
        }
    }

    cudaDeviceSynchronize();
}

void egblas_shuffle(size_t n, void* x, size_t incx){
