Fast Fourier transforms are among the most heavily used operations in all of computational science, quietly underpinning everything from magnetic resonance imaging to turbulence simulations and radar signal processing. On a single graphics processor, the arithmetic itself is handled extremely well by vendor libraries such as NVIDIA’s cuFFT or AMD’s hipFFT. But the moment a transform grows too large for one device, or a simulation spans hundreds of GPUs across many nodes of a supercomputer, the real challenge shifts from arithmetic to communication. A new open-source library called FFTM, described in the journal SoftwareX by Nikolay M. Evstigneev, tackles exactly this problem, and it does so for a dimension that most GPU libraries have ignored: four.

The core difficulty in a distributed FFT is data ownership. Each MPI rank, typically mapped to one GPU, holds only a slice of the global tensor. A multidimensional transform requires that every axis be transformed in turn, but a rank can only transform an axis locally if it owns complete lines along that axis. So the computation alternates between bursts of local FFT work and full redistributions, called transpositions, in which ranks exchange enormous volumes of data so that the next axis becomes locally accessible. At scale, this communication, along with packing, unpacking, synchronization and temporary storage, can cost more than the local transforms themselves. Which decomposition and transport strategy wins depends delicately on tensor dimensions, message sizes, the process grid, device memory and the machine’s interconnect, which is why no single configuration is universally best.

Existing libraries cover important parts of this landscape. FFTW-MPI and PFFT serve CPU workloads, OpenFFT offers adaptive decomposition with communication autotuning on CPUs, and heFFTe, cuFFTMp, 2DECOMP&FFT and several others bring distributed 2D and 3D transforms to GPUs with varying transport options. Yet a striking gap remains: none of the reviewed dense GPU libraries exposes a coupled four-dimensional real-to-complex and complex-to-real plan spanning multiple GPUs and nodes, where a single distributed plan owns all four local-transform and redistribution stages. A fourth transformed axis is not equivalent to a batch of independent 3D FFTs; it demands its own layout and its own redistributions. Applications that need a joint spectrum in space and time, such as 3D magnetic resonance spectroscopic imaging, volumetric space-time filtering, light-field imaging, radar and turbulence analysis in combined wave-number and frequency space, have had to reconstruct this pipeline themselves, inheriting the burden of correctness, portability and machine-specific tuning.

FFTM addresses the gap with a header-only C++14 library for distributed 3D and 4D real-to-complex and complex-to-real transforms. Its architecture is deliberately conservative in its language requirements so that older hardware remains supported, and it is built around a clean abstraction boundary. A transform object owns the global dimensions, the processor grid, partition descriptors, local FFT plans, transpose executors and a reusable workspace. Multidimensional device arrays come from the companion SCFD header-only library, which provides backend-neutral tensors and transparent host-device operations. Crucially, all vendor-specific handles, stream types and complex types are confined to a thin wrapper layer, so neither the public API nor any transpose implementation contains cuFFT or hipFFT types. The same application code therefore runs unchanged on CUDA or HIP backends, and other local FFT libraries can be substituted provided suitable interfaces exist.

The library’s distributed model separates two concepts: partition, which answers which rank owns each global index, and layout, which answers how that rank stores its owned indices. In a slab decomposition only one global axis is divided among ranks, while a full pencil decomposition divides two axes in 3D and three axes in 4D over a process grid, leaving complete one-dimensional lines along the remaining transform axis. Redistributions then rotate ownership so each successive axis can be transformed. For 4D, ranks form a three-dimensional process grid, and the transform sequence walks through four local transforms interleaved with three redistributions. FFTM offers several transport modes, including point-to-point schemes that post nonblocking receives and sends and either wait for all requests or process completed payloads one at a time, overlapping communication with unpacking, plus collective alltoallv and alltoallw modes that delegate scheduling to MPI. Users can select whichever mode suits their cluster topology.

Because the fastest strategy varies by machine and problem, FFTM includes production presets that convert high-level choices into validated low-level settings, and a measured autotuning procedure for users who simply want the best performance on their hardware. Autotuning compares production-safe alternatives on the actual machine, scoring each candidate by the median over repetitions of the maximum rank wall time for a complete forward and backward pair, and caches the winner for later reuse. The results illustrate why measurement matters: on eight A100 GPUs with a 2048-cubed grid, a slab-pencil configuration beat a pencil-pencil one by 7.78 percent, while on two V100 GPUs with a 64-to-the-fourth-power 4D problem, the point-to-point wait-any mode was 2.79 times faster than the collective alltoallv mode for the same pencil decomposition.

The performance evaluation was conducted on homogeneous nodes each containing eight 80-gigabyte NVIDIA A100 GPUs and dual AMD EPYC 7742 processors, with one MPI rank per GPU, CUDA-aware MPI, and the HPC-X communication stack built with UCX and GDRCopy. Every reported time is the maximum rank wall time for a complete forward and backward pair, with plan creation excluded, and results were accepted only when the relative error stayed below ten to the minus eleventh power. Against matched open-source GPU competitors on the same hardware, FFTM was 1.03 to 1.55 times faster than the Egger implementation across memory-fitted 3D sizes on two to eight GPUs, and 1.04 to 1.68 times faster than GPU-FFT for a 1920-cubed grid on 8 to 32 GPUs. At equal node counts, it accelerated the 2048-cubed 3D case by 34.74, 16.83 and 12.39 times, and the 320-to-the-fourth 4D case by 51.71, 19.03 and 13.43 times, relative to tuned FFTW-MPI CPU configurations on one, two and four nodes respectively.

Strong scaling from 8 to 120 GPUs revealed the physics of the problem. For the fixed 2048-cubed workload, the pair time fell from 230.2 to 107.3 milliseconds, reaching 26.42 effective teraflops, but the efficiency relative to the eight-GPU baseline was only 14.3 percent, because MPI redistribution grew to consume 87 to 92 percent of total wall time at the largest scales while local FFT execution shrank to under 10 percent. Communication analysis showed that at two nodes the transfers reached 85 to 89 percent of the aggregate network line rate, but at 15 nodes, with smaller per-node payloads and greater peer fan-out, useful-payload efficiency dropped to 46 percent in 3D and 35 percent in 4D. Weak scaling told a similar story: aggregate throughput kept rising, but work-corrected efficiency fell to 13.4 and 9.3 percent at 120 GPUs, a limitation the author notes is largely imposed by the communication-volume lower bound itself rather than by FFTM alone, since independent libraries showed the same 3D trend.

The application demonstration ties the whole story together with the classic Taylor-Green vortex, an incompressible Navier-Stokes benchmark at Reynolds number 1600. A pseudo-spectral production run on a 2048-cubed grid used 64 A100 GPUs across eight nodes, requiring 25,752 time steps and 12.554 hours to reach the target time, with the maximum velocity divergence held to 7.51 times ten to the minus fifteenth. The computed kinetic energy and enstrophy histories matched published pseudospectral, finite-difference and discontinuous-Galerkin reference data with relative errors of 1.52 times ten to the minus sixth and 2.15 times ten to the minus fifth. Then came the showcase for the fourth dimension: 80 signed vorticity snapshots, each spectrally restricted to 512 cubed, were transformed as a single coupled 512-cubed-by-80 space-time array of over ten billion samples on eight GPUs, taking about 26.5 seconds each way. Filtering the lowest four temporal modes recovered the dominant large-scale dynamics, all within one distributed 4D plan and without any application-level global redistribution.

FFTM is released under the GNU General Public License version 3, with reproducible CUDA and HIP Docker capsules, checksum manifests and validation records, and it is already being used in an all-Fourier discretization project aimed at stabilizing relative periodic orbits in turbulent flows. Its limitations are honest ones: it handles only dense Cartesian real-to-complex and complex-to-real transforms, assumes one MPI rank per GPU by default, and does not replace sparse or nonuniform FFT packages. Future releases plan scheduling optimizations and NVSHMEM support for CUDA backends to attack the communication delays that now dominate at extreme scale. Even so, the library marks a meaningful step: a portable, vendor-agnostic, automatically tuned path to coupled 3D and 4D spectral computation on the machines where the largest simulations actually run, freeing application scientists from rebuilding the same distributed machinery for every new problem.

Subject of Research: Portable distributed three- and four-dimensional FFTs for multi-GPU systems

Article Title: FFTM: Portable distributed three- and four-dimensional FFTs for multi-GPU systems

Article References: Evstigneev, N. M. (2026). FFTM: Portable distributed three- and four-dimensional FFTs for multi-GPU systems. SoftwareX, 36, Article 103110. https://doi.org/10.1016/j.softx.2026.103110

Image Credits: AI Generated

DOI: 10.1016/j.softx.2026.103110

Keywords: FFTM, fast Fourier transform, multi-GPU computing, CUDA, HIP, MPI, autotuning, turbulence, spectral methods, high-performance computing, Taylor-Green vortex, SoftwareX