Why CuFFT throughput increases as the transform size gets larger?

Question

(Updated) I am trying to understand how CUDA parallelism works in CuFFT while learning CUDA coding.

I wrote my version of 1-D FFT in CUDA C++ and compared it with cuFFT. Below are the throughputs I got when I compared the two methods for a single block FFT tests.

1-D FFT Size: throughput-1 (time) in my version, throughput-2 (time) in CuFFT. Gf = GFlops (using cudaEventElapsedTime, GTX 1050, includes cufftPlan1d in cuFFT - see code below)

(1)2^16: 1.3 Gf(3.8 msec), 0.19 Gf(27 msec)
(2)2^17: 1.6 Gf(6.7 msec), 0.5 Gf(23 msec)
(3)2^18: 1.8 Gf(13 msec), 1.1 Gf(22 msec)
(4)2^19: 1.8 Gf(27 msec), 2.1 Gf(23 msec)
(5)2^20: 1.8 Gf(56 msec), 4.5 Gf(23 msec)
(6)2^21: 1.8 Gf(121 msec), 9.6 Gf(22 msec)
(7)2^23: 1.8 Gf(523 msec), 31 Gf (30 msec)
(8)2^25: 1.8 Gf(2298 msec), 106 Gf(39 msec)

I am wondering why the throughput from CuFFT increases almost linearly as the transform size? Since the throughput is computed by computational complexity of the transform divided by computation time, I thought it should stay almost constant as my version shows if everything remain same. What technique should I explore to improve my version behave like 1-D CuFFT (grid size strategy)? My code is included below ( I think the code works fine since FFT-IFFT shows the errors less than 1e-8 where the computation is double precision with the input vector data[i] = Complex(i.0, 0))

#include <iostream>
#include <vector>
#include <cuda_runtime.h>
#include <thrust/complex.h>
#include <cmath>
#include <iomanip>  // For std::fixed and std::setprecision
#include <cufft.h>

using Complex = thrust::complex<double>;
#define USE_SHARED_MEM 0 // disabled: using shared memory didn't improve performance

#define BLOCK_SIZE 1024  
// #define BLOCK_SIZE 512
// #define BLOCK_SIZE 32
#if USE_SHARED_MEM
    #define SHARED_MEM_SIZE_FACTOR 4 // relative to BLOCK_SIZE
#endif

unsigned int mixedRadixReverse(unsigned int x, int numBits) {
    unsigned int rev = 0;
    int numRadix4Stages = numBits / 2;  // Number of radix-4 stages
    int useRadix2 = numBits % 2;  // If numBits is odd, we have an extra radix-2 stage

    // First, reverse the radix-2 bit if present (apply to LSB first)
    if (useRadix2) {
        rev = (rev << 1) | (x & 1);  // Extract LSB for radix-2
        x >>= 1;
    }

    // Then, reverse the radix-4 groups (higher bits)
    for (int i = 0; i < numRadix4Stages; i++) {
        rev = (rev << 2) | (x & 3);  // Extract 2 bits at a time
        x >>= 2;
    }

    return rev;
}
__global__ void scale(Complex* d_data, float scale, int N) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N) {
        d_data[idx] *= scale;
    }
}

__global__ void fft_radix2(Complex* data, int sizeFFT, int inFFTSize, int direction) {
    int k = blockIdx.x * blockDim.x + threadIdx.x;
    if (k >= sizeFFT / 2) return;
    int outFFTSize = 2*inFFTSize;
    float sign = direction == 1 ? -1.0: 1.0;
    Complex tw = Complex(cos(2*M_PI*k/outFFTSize), sign*sin(2*M_PI*k/outFFTSize));
    Complex t = data[k + inFFTSize]*tw;
    Complex tmp = data[k];
    data[k] = tmp + t;
    data[k + inFFTSize] = tmp - t;
}

__global__ void fft_radix4(Complex* data, int sizeFFT, int inFFTSize, int direction) {
    float sign = direction == 1 ? -1.0: 1.0;
    int outFFTSize = 4*inFFTSize;
    int k = blockIdx.x * blockDim.x + threadIdx.x;
    if (k >= sizeFFT / 4) return; // max number of threads is sizeFFT /4 for radix-4 stage
    // k = 0, 1, 2, ..., 7 when sizeFFT = 32
    int k_offset_inFFTSize = k%inFFTSize; // k_inFFTSize= 0, 1 ..., 3 when inFFTSize = 4
    
    // printf("k=%d: inFFTSize=%d, outFFTSize=%d\n",k, inFFTSize,outFFTSize);
    
    int blockOffset = (k/inFFTSize)*outFFTSize; // replaces outer loop
    int offset = blockOffset + k_offset_inFFTSize; // replaces inner loop
    // dragonfly inputs
    int idx0 = offset;
    int idx1 = offset + inFFTSize;
    int idx2 = offset + 2 * inFFTSize;
    int idx3 = offset + 3 * inFFTSize;

    Complex tw = Complex(cos(2.0*M_PI*k_offset_inFFTSize/outFFTSize), 
                    sign*sin(2.0*M_PI*k_offset_inFFTSize/outFFTSize));
    Complex t0 = data[idx0];
    Complex t1 = data[idx1]*tw;
    Complex t2 = data[idx2]*tw*tw;
    Complex t3 = data[idx3]*tw*tw*tw;

    Complex a = t0 + t2;
    Complex b = t1 + t3;
    Complex c = t0 - t2;
    Complex d = (t1 - t3) * Complex(0, 1.0*sign);
    // in-place computation
    data[idx0] = a + b;
    data[idx1] = c + d;
    data[idx2] = a - b;
    data[idx3] = c - d;
}
__global__ void fft_radix2_shared(Complex* data, int sizeFFT, int inFFTSize, int direction) {
    extern __shared__ Complex s_data[];
    int k = blockIdx.x * blockDim.x + threadIdx.x;
    if (k >= sizeFFT / 2) return;

    int outFFTSize = 2 * inFFTSize;
    float sign = (direction == 1) ? -1.0 : 1.0;
    Complex tw = Complex(cos(2 * M_PI * k / outFFTSize), sign * sin(2 * M_PI * k / outFFTSize));

    // Load data into shared memory
    int idx1 = k + inFFTSize;
    s_data[threadIdx.x] = data[k];
    s_data[threadIdx.x + blockDim.x] = data[idx1];
    __syncthreads();

    // Perform FFT operation in shared memory
    Complex t = s_data[threadIdx.x + blockDim.x] * tw;
    Complex tmp = s_data[threadIdx.x];

    s_data[threadIdx.x] = tmp + t;
    s_data[threadIdx.x + blockDim.x] = tmp - t;

    __syncthreads();

    // Write back to global memory
    data[k] = s_data[threadIdx.x];
    data[idx1] = s_data[threadIdx.x + blockDim.x];
}

__global__ void fft_radix4_shared(Complex* data, int sizeFFT, int inFFTSize, int direction) {
    extern __shared__ Complex s_data[];
    float sign = (direction == 1) ? -1.0 : 1.0;
    int outFFTSize = 4 * inFFTSize;
    int k = blockIdx.x * blockDim.x + threadIdx.x;
    if (k >= sizeFFT / 4) return;

    int k_offset_inFFTSize = k % inFFTSize;
    int blockOffset = (k / inFFTSize) * outFFTSize;
    int offset = blockOffset + k_offset_inFFTSize;

    // Indices for dragonfly inputs
    int idx0 = offset;
    int idx1 = offset + inFFTSize;
    int idx2 = offset + 2 * inFFTSize;
    int idx3 = offset + 3 * inFFTSize;

    // Load into shared memory
    s_data[threadIdx.x] = data[idx0];
    s_data[threadIdx.x + blockDim.x] = data[idx1];
    s_data[threadIdx.x + 2 * blockDim.x] = data[idx2];
    s_data[threadIdx.x + 3 * blockDim.x] = data[idx3];

    __syncthreads();

    // Twiddle factors
    Complex tw = Complex(cos(2.0 * M_PI * k_offset_inFFTSize / outFFTSize), 
                        sign * sin(2.0 * M_PI * k_offset_inFFTSize / outFFTSize));

    Complex t0 = s_data[threadIdx.x];
    Complex t1 = s_data[threadIdx.x + blockDim.x] * tw;
    Complex t2 = s_data[threadIdx.x + 2 * blockDim.x] * tw * tw;
    Complex t3 = s_data[threadIdx.x + 3 * blockDim.x] * tw * tw * tw;

    // Compute FFT butterfly operations
    Complex a = t0 + t2;
    Complex b = t1 + t3;
    Complex c = t0 - t2;
    Complex d = (t1 - t3) * Complex(0, sign);

    // Store results back to shared memory
    s_data[threadIdx.x] = a + b;
    s_data[threadIdx.x + blockDim.x] = c + d;
    s_data[threadIdx.x + 2 * blockDim.x] = a - b;
    s_data[threadIdx.x + 3 * blockDim.x] = c - d;

    __syncthreads();

    // Write results back to global memory
    data[idx0] = s_data[threadIdx.x];
    data[idx1] = s_data[threadIdx.x + blockDim.x];
    data[idx2] = s_data[threadIdx.x + 2 * blockDim.x];
    data[idx3] = s_data[threadIdx.x + 3 * blockDim.x];
}


void fft_1d(Complex* data, int sizeFFT, int direction) {

    //(1) Host memory allocation for radix reversal and DMA 
    Complex * h_pinned_data;
    cudaMallocHost((void**)&h_pinned_data, sizeFFT*sizeof(Complex));
    // radix-reverse in the host memory
    int numBits = __builtin_ctz(sizeFFT);  
    #pragma omp parallel for  // Parallel processing for CPU speedup
    for (int i = 0; i < sizeFFT; i++) {
        unsigned int revIdx = mixedRadixReverse(i, numBits);
        h_pinned_data[revIdx] = data[i]; 
    }
    //(2) device memory allocaiton
    Complex* d_data;
    cudaError_t err;
    err = cudaMalloc((void**)&d_data, sizeFFT * sizeof(Complex));
    if (err != cudaSuccess) { std::cerr << "cudaMalloc failed: "  << std::endl;  return; }
    // Asynchronous copy using pinned memory for performance boost
    cudaMemcpyAsync(d_data, h_pinned_data, sizeFFT * sizeof(Complex), cudaMemcpyHostToDevice, 0);
  

    // (3) plan
    // dim3 blockSize(256);
    dim3 blockSize(BLOCK_SIZE);
    dim3 gridSize((sizeFFT + blockSize.x - 1) / blockSize.x);
    
    int numStagesRadix4 = 0;
    int temp_sizeFFT = sizeFFT;
    while (temp_sizeFFT % 4 == 0) {
        temp_sizeFFT /= 4;
        numStagesRadix4++;
    }
    // std::cout << "number of radix4 stages = " << numStagesRadix4 << std::endl;
    #if USE_SHARED_MEM
    // shared memory allocation
    size_t sharedMemSize = blockSize.x * sizeof(Complex) *SHARED_MEM_SIZE_FACTOR;//256;//4;
    #endif
    // (4) Call radix computations
    int inFFTSize = 1;
    for (int i = 0; i < numStagesRadix4; ++i) {   
        printf("radix-4 stage = %d\n", i);     
        
        #if USE_SHARED_MEM
        fft_radix4_shared<<<gridSize, blockSize, sharedMemSize>>>(d_data, sizeFFT, inFFTSize, direction);
        #else
        fft_radix4<<<gridSize, blockSize>>>(d_data, sizeFFT, inFFTSize, direction);
        #endif

        err = cudaDeviceSynchronize();
        if (err != cudaSuccess) { std::cerr << "fft_radix4 execution at stage " << i << " failed!" << std::endl; return; }

        inFFTSize *= 4;
        // display_device_array_cmplx(d_data, sizeFFT, "d_data =");
    }
    
    if (temp_sizeFFT == 2) {
        std::cout << "final radix2 stage: inFFTSize = " << inFFTSize  << std::endl;
        #if USE_SHARED_MEM
        fft_radix2_shared<<<gridSize, blockSize, sharedMemSize>>>(d_data, sizeFFT, inFFTSize, direction);
        #else
        fft_radix2<<<gridSize, blockSize>>>(d_data, sizeFFT, inFFTSize, direction);
        #endif

        err = cudaDeviceSynchronize();
        if (err != cudaSuccess) { std::cerr << "fft_radix2 execution failed!" << std::endl; return; }
    }
    // (5) scale if necessary
    if (direction == -1) scale<<<gridSize, blockSize>>>(d_data, 1.0/(float) sizeFFT, sizeFFT);
    // send the result back to host
    err = cudaMemcpy(data, d_data, sizeFFT * sizeof(Complex), cudaMemcpyDeviceToHost);
    if (err != cudaSuccess) { std::cerr << "cudaMemcpyDeviceToHost failed!" << std::endl; return; }
    
    cudaFree(d_data);
}

int main() {
    
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    float milliseconds = 0;

    // define input data ==================
    // int sizeFFT = 1<< 28; //Not working: cudaMalloc fails
    // int sizeFFT = 1<< 27; 
    // int sizeFFT = 1<< 26; 
    // int sizeFFT = 1<< 25; 
    // int sizeFFT = 1<< 23; 
    // int sizeFFT = 1<< 21; 
    // int sizeFFT = 1<< 20; 
    // int sizeFFT = 1<< 19; 
    // int sizeFFT = 1<< 18; 
    // int sizeFFT = 1<< 17; 
    int sizeFFT = 1<< 16; 

    // test vector
    std::vector<Complex> data(sizeFFT, Complex(1.0, 0.0));
    std::vector<Complex> data1(sizeFFT, Complex(0.0, 0.0));
    for (int i = 0; i < sizeFFT; ++i) {data[i] = Complex(i,0);}
    data1 = data;

    cudaEventRecord(start); // =============>
    // FFT
    int direction = 1; // 1: forward, -1: reverse
    fft_1d(data.data(), sizeFFT, direction);

    cudaEventRecord(stop); // <===============
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(&milliseconds, start, stop);
    double N = sizeFFT;
    double gflops = (5.0 * N * log2(N)) / (milliseconds * 1.0e6);
    std::cout << "Throughput: " << gflops << " GFLOPS, time:" << milliseconds<<  ", BLOCK_SIZE:" << BLOCK_SIZE << std::endl;


    ///////////////////////////////////////////
    // Comparison with cuFFT
    ///////////////////////////////////////////
    //(1) copy data
    cufftComplex* d_data;
    cudaError_t err;
    err = cudaMalloc((void**)&d_data, sizeFFT * sizeof(cufftComplex));
    if (err != cudaSuccess) { std::cerr << "cudaMalloc failed: "  << std::endl;  return; }
    // Asynchronous copy using pinned memory for performance boost
    cudaMemcpyAsync(d_data, data.data(), sizeFFT * sizeof(cufftComplex), cudaMemcpyHostToDevice, 0);
    //(2) Create cuFFT plan
    cudaEventRecord(start); // =============>
    cufftHandle plan;
    cufftPlan1d(&plan, sizeFFT, CUFFT_C2C, 1);
    size_t workspace_size;
    cufftGetSize(plan, &workspace_size);
    std::cout << "Workspace size: " << workspace_size/1024 << " K bytes" << std::endl;

    //(3) Execute FFT
    cufftExecC2C(plan, d_data, d_data, CUFFT_FORWARD);
    cudaEventRecord(stop); // <===============
    cufftDestroy(plan);
    cudaFree(d_data);
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(&milliseconds, start, stop);
    N = sizeFFT;
    gflops = (5.0 * N * log2(N)) / (milliseconds * 1.0e6);
    std::cout << "Throughput (cuFFT): " << gflops << " GFLOPS, time:" << milliseconds << " msec" <<std::endl;

    // IFFT
    direction = -1; // inverse
    fft_1d(data.data(), sizeFFT, direction);
    
    double maxDiff = 0;
    for (int i = 0; i < sizeFFT; ++i) {
        thrust::complex<double>  err = data[i] - data1[i];
        double errSquare = err.real()*err.real() + err.imag()*err.imag() ;
        maxDiff = fmax(maxDiff, sqrt(errSquare));
    }
    std::cout << "Max difference: " << maxDiff << std::endl;

    return 0;
}