Scalability of Fast Fourier Transform (FFT)

Question

To use the Fast Fourier Transform (FFT) on uniformly sampled data, e.g. in connection with PDE solvers, it is well known that the FFT is an $\mathcal{O}(n\log(n)$) algorithm. How well do the FFT scale when processed in parallel for $n\to\infty$ (i.e. very large)?

Answer 1

This is more anecdotal evidence than demonstrated proof, but it appears that existing implementations for FFT's, such as FFTW, have a limit to their scaling ability.

When we started using LAMMPS's ${\bf k}$-space solvers in very large systems ($O(10^7)$ atoms), we found that the scaling continued, so long as we were able to keep the number of processors small enough that they could fit onto one rack. As soon as we tried to expand further (above about 4K processors, depending on the machine), the scaling broke down—apparently because the communication costs for pushing data between the processors became too large to maintain the scaling. [Recently, to circumvent this problem, they've introduced the ability to dedicate a certain partition of the processor allocation to FFT calculation.]

But the take-home message here is that FFT should scale up; however, there are sometimes unexpected limitations and interactions that come into play when one moves from the theoretical consideration of an algorithm's performance to its practical implementation on an actual HPC platform.

Answer 2

The fast multipole method (FMM) is $O(n)$ and has much lower communication requirements, so it provides a highly scalable discrete Fourier transform. Edelman, McCorquodale, and Toledo (1999) The Future Fast Fourier Transform? analyze this approach and conclude that FMM would be preferable to conventional FFT at large scale. Note that the FMM is only approximate, so the constants are worse if you need very high accuracy. Thanks to Jack Poulson for pointing this out in a discussion last week.

Answer 3

In the context of PDEs, it's important to recognize that the value of $n$ for a required 1D FFT will typically grow like the $d$th root of the total number of grid points, where the dimensionality $d$ is most often 3.

Searching for "parallel FFT" or "pseudospectral scalability" on Google Scholar yields a wealth of information that I'm unqualified to assess. But this seems like a nice recent example of what can be accomplished in practice:

A hybrid MPI-OpenMP scheme for scalable parallel pseudospectral computations for fluid turbulence

Abstract:

A hybrid scheme that utilizes MPI for distributed memory parallelism and OpenMP for shared memory parallelism is presented. The work is motivated by the desire to achieve exceptionally high Reynolds numbers in pseudospectral computations of fluid turbulence on emerging petascale, high core-count, massively parallel processing systems. The hybrid implementation derives from and augments a well-tested scalable MPI-parallelized pseudospectral code. The hybrid paradigm leads to a new picture for the domain decomposition of the pseudospectral grids, which is helpful in understanding, among other things, the 3D transpose of the global data that is necessary for the parallel fast Fourier transforms that are the central component of the numerical discretizations. Details of the hybrid implementation are provided, and performance tests illustrate the utility of the method. It is shown that the hybrid scheme achieves near ideal scalability up to ~20000 compute cores with a maximum mean efficiency of 83%. Data are presented that demonstrate how to choose the optimal number of MPI processes and OpenMP threads in order to optimize code performance on two different platforms.

Answer 4

If you have an infinite number of processors, the DFT can be determined in $O(n)$ time.

In the naive algorithm, you could put each point of output on a seperate node, and calculate that fourier tranformed point in $O(\log n)$ time. Any fast algorithm should be able to at least match this scaling.

However, you also need to collect all of the fourier transformed points in one array, which takes $O(n)$ time.