fft-gpu-programming-models

Six implementations of the same FFT, benchmarked head to head on NVIDIA P100 and V100, to find out which version has what tradeoff and performance, and to decide when it's a good idea to use which implementation.

#	Implementation	Approach	Effort to write
1	CPU	Single-threaded C reference	—
2	CUDA V1	Global memory, one kernel launch per butterfly stage	High
3	CUDA V2	Shared memory, all stages in a single launch — launch invalid above N = 2048, see below	High
4	CUDA V3	Shared memory + padding to avoid bank conflicts — same problem	High
5	OpenACC	parallel loop independent, collapse(2) on the same loops	Low
6	cuFFT	NVIDIA's library, cufftExecZ2Z	Lowest

What the runs showed

• cuFFT is ~12× faster than the best verified hand-written kernel (global-memory CUDA). If a vendor library covers what you need, use it.
• OpenACC is ~1.8× slower than the global-memory CUDA kernel — for a fraction of the effort, and timed with clock(), so treat the ratio as approximate.
• The shared-memory timings are not valid. Both kernels launch a single block of N/2 threads with N·16 bytes of block-shared memory. CUDA allows at most 1024 threads per block, and shared memory is kilobytes, so the launch is only legal up to N = 2048. The launch error was never checked, so at the benchmark sizes those timings most likely measured the bit-reversal kernel alone. That is also why padding and threads-per-block appeared to change nothing.

Stack: C · CUDA · OpenACC · cuFFT · NVIDIA HPC SDK · SLURM

Algorithm: radix-2 Cooley–Tukey FFT, forward transform, single GPU scaling

Where it ran: the COKA cluster at INFN — P100 and V100 nodes, skyvolta partition for the V100 runs, one GPU per job, exclusive.

Problem size	Value
N	2²³ = 8,388,608 points
Type	double-precision complex (double2), 16 bytes/point
Array	128 MiB
Operation count	5·N·log₂N = 0.96 GFLOP per transform
Timing	device-side; array resident on GPU, host transfer excluded for all six equally

All six are in one binary (src/fft_gpu_bench.cu) and run back to back on the same input.

Conclusion

At N = 2²³ on a COKA V100 (7.8 TFLOP/s fp64 peak):

Implementation	Time	GFLOP/s	% of peak	vs CUDA V1	vs CPU
cuFFT	~2 ms	482	6.2%	~12× faster	~1000×
CUDA V1 (global memory)	~25 ms	39	0.5%	baseline	~80×
OpenACC	~45 ms	21	0.3%	~1.8× slower	~45×
CPU, single-threaded	~2000 ms	0.5	—	—	1×
CUDA V2 / V3 (shared memory)	~~~9 ms~~	—	—	invalid launch	—

_{The plots were made before the shared-memory launch problem was found; ignore the V2/V3 lines.}

• Even cuFFT reaches only 6% of fp64 peak
• The global-memory transform moves the 128 MiB array once per stage, 23 stages
• So the difference between the valid implementations is how well each uses the memory hierarchy, not flop rate

→ This is a memory-hierarchy benchmark wearing an FFT costume.

What OpenACC actually costs you

• Took a fraction of the time to write — a few pragmas on loops that already existed
• The directives parallelise the loops as written, including the per-butterfly cos/sin; they don't restructure the algorithm or the memory traffic for you
• The result was ~1.8× slower than a hand-written global-memory kernel

→ Use OpenACC to get onto the GPU quickly. Use a library when one exists.

The shared-memory kernels — what went wrong

On the COKA P100 at N = 2²⁴:

Version	Time	Status
V1 — global memory	33 ms	valid
V2 — shared memory	7.2 ms	invalid launch
V3 — shared + padding	7.2 ms	invalid launch

• fft_shared_kernel<<<1, N/2, N·16>>> asks for one block of N/2 threads: legal only up to N = 2048
• No cudaGetLastError() after the launch, so the failure was silent
• Tell-tale signs, in hindsight: V2 and V3 identical; padding and block size changed nothing; the P100 at 2²⁴ "beat" the V100 at 2²³
• A correct version runs many blocks of ≤1024 threads, uses shared memory only for the butterfly stages whose pairs stay inside one block, and does the larger stages in global memory

Corrections

Re-reading src/fft_gpu_bench.cu against the numbers it had produced, I found the shared-memory kernels never ran at the benchmark sizes at all.

fft_shared_kernel<<<1, N/2, N*16>>> asks for one block of N/2 threads. CUDA caps a block at 1024 threads, and block-shared memory at tens of kilobytes, so the launch is only legal up to N = 2048. There is no cudaGetLastError() after it, so the failure was silent and the timer measured the surrounding work. Every conclusion that rested on those timings — "shared memory is worth 4.6×", and "OpenACC is 5× slower than CUDA" — is retracted, and the comparison is re-based on CUDA V1, the fastest kernel here that verifiably runs.

What gives it away, in the data itself:

• V2 and V3 produced identical times, though V3 exists solely to pad away bank conflicts
• Changing threads-per-block changed nothing
• A P100 at N = 2²⁴ "beat" a V100 at 2²³ — more work, older card, faster result

The useful lesson is sharper than "check your error codes": a result that doesn't respond to the parameter you're varying is usually not measuring that parameter. Three knobs did nothing, and that is the signature of a kernel that isn't running.

Doing this properly means splitting the shared-memory stages into ≤1024-thread blocks — shared memory for the butterflies whose pairs stay inside a block, global memory for the wider stages — timing all six paths with CUDA events so the OpenACC clock() measurement becomes comparable, and validating at each benchmark size rather than only at small N.

Caveats

Caveat	Detail
Shared-memory variants invalid at benchmark sizes	Single-block launch exceeds CUDA's per-block limits above N = 2048; errors unchecked
Timing methods differ	CUDA and cuFFT paths use cudaEventRecord; the OpenACC path uses clock() around its data region — process CPU time, not wall time. Not equivalent measurements
Small-N results are an artefact	Below ~N = 2¹⁵ the OpenACC line comes out fastest in the plots, which can't be right for this kernel. Ignore that end
Fix not done yet	Split the shared-memory stages into ≤1024-thread blocks, check cudaGetLastError() after every launch, time all six with CUDA events, validate each size, re-run
Times read off plots	The raw logs were not retained

Correctness

validate.py compares each implementation's output against NumPy's FFT. The output writers in fft_gpu_bench.cu are commented out, so validation was only done at small N, where the single-block launch is legal. At that size all six agree with the reference to ~1e-14; the shared-memory variants have not been validated at the benchmark sizes.

Build and run

The benchmark has CUDA kernels, cuFFT and OpenACC in one file, so it needs the NVIDIA HPC SDK. nvcc on its own can't compile the OpenACC parts.

cd src
nvc++ -O3 -acc -gpu=cc70 -cuda -lcufft -lm fft_gpu_bench.cu -o fft_gpu_bench.x
./fft_gpu_bench.x 8388608          # N must be a power of two

CPU reference alone, no GPU toolchain needed:

gcc -O3 -fopenmp -o fft_cpu.x src/fft_cpu.c -lm && ./fft_cpu.x 1048576

On a cluster:

sbatch slurm/bench.sh

• slurm/bench.sh is written for COKA — skyvolta partition, one V100, nvhpc/24.5
• Change the partition and module for anywhere else
• Set -gpu=cc60 for P100, cc70 for V100, cc80 for A100

Layout

src/
  fft_gpu_bench.cu        all six implementations + benchmark harness
  fft_cpu.c               standalone CPU reference
  snapshots/              earlier development stages
    stage1_global_shared.cu
    stage2_padded.cu
benchmarks/               plots
slurm/bench.sh            size sweep, 2^11 to 2^23
validate.py               NumPy check
presentation.pdf          course presentation

Where this came from

Course	P2.2 — GPU Programming, MHPC, ICTP / SISSA Trieste
Written	February–June 2026
Public since	21 June 2026, in prabhkodes/low_level_optimisations
This repo	Same code, better filenames, result written down properly. Source unchanged
Course repository	Belongs to SISSA, private — can't be linked