0.10 Implementing ffts in practice  (Page 14/21)

Roughly speaking, to solve a general DFT problem, one must perform three tasks. First, one must reduce a problem of arbitrary vectorrank to a set of loops nested around a problem of vector rank 0, i.e., a single (possibly multi-dimensional) DFT. Second, one must reducethe multi-dimensional DFT to a sequence of of rank-1 problems, i.e., one-dimensional DFTs; for simplicity, however, we do notconsider multi-dimensional DFTs below. Third, one must solve the rank-1, vector rank-0 problem by means of some DFT algorithm suchas Cooley-Tukey. These three steps need not be executed in the stated order, however, and in fact, almost every permutation and interleavingof these three steps leads to a correct DFT plan. The choice of the set of plans explored by the planner is critical for the usability ofthe FFTW system: the set must be large enough to contain the fastest possible plans, but it must be small enough to keep theplanning time acceptable.

Rank-0 plans

The rank-0 problem $dft(\left\{\right\},\mathbf{V},\mathbf{I},\mathbf{O})$ denotes a permutation of the input array into the output array. FFTWdoes not solve arbitrary rank-0 problems, only the following two special cases that arise in practice.

• When $\left|\mathbf{V}\right|=1$ and $\mathbf{I}\ne \mathbf{O}$ , FFTW produces a plan that copies the input array into the output array.Depending on the strides, the plan consists of a loop or, possibly, of a call to the ANSI C function memcpy , which is specialized to copy contiguous regions of memory.
• When $\left|\mathbf{V}\right|=2$ , $\mathbf{I}=\mathbf{O}$ , and the strides denote a matrix-transposition problem, FFTW creates a planthat transposes the array in-place. FFTW implements the square transposition $dft(\left\{\right\},\left\{(,n,,,\iota ,,,o,),,,(,n,,,o,,,\iota ,)\right\},\mathbf{I},\mathbf{O})$ by means of the cache-oblivious algorithm from [link] , which is fast and, in theory, uses the cache optimally regardless of the cachesize (using principles similar to those described in the section "FFTs and the Memory Hierarchy" ). A generalization of this idea is employed for non-square transpositions with a large common factor or a small differencebetween the dimensions, adapting algorithms from [link] .

Rank-1 plans

Rank-1 DFT problems denote ordinary one-dimensional Fourier transforms. FFTW deals with most rank-1 problems as follows.

Direct plans

When the DFT rank-1 problem is “small enough” (usually, $n\le 64$ ), FFTW produces a direct plan that solves the problem directly. These plans operate by calling a fragment of C code(a codelet ) specialized to solve problems of one particular size, whose generation is described in "Generating Small FFT Kernels" . More precisely, the codelets compute a loop ( $\left|\mathbf{V}\right|\le 1$ ) of small DFTs.

Cooley-tukey plans

For problems of the form $dft(\left\{(,n,,,\iota ,,,o,)\right\},\mathbf{V},\mathbf{I},\mathbf{O})$ where $n=rm$ , FFTW generates a plan that implements a radix- $r$ Cooley-Tukey algorithm "Review of the Cooley-Tukey FFT" . Both decimation-in-time and decimation-in-frequency plans are supported, with both small fixed radices (usually, $r\le 64$ ) produced by the codelet generator "Generating Small FFT Kernels" and also arbitrary radices (e.g. radix- $\sqrt{n}$ ).

The most common case is a decimation in time ( DIT ) plan, corresponding to a radix $r=n_2$ (and thus $m=n_1$ ) in the notation of "Review of the Cooley-Tukey FFT" : it first solves $dft(\left\{(,m,,,r,·,\iota ,,,o,)\right\},\mathbf{V}\cup \left\{(,r,,,\iota ,,,m,·,o,)\right\},\mathbf{I},\mathbf{O})$ , then multiplies the output array $\mathbf{O}$ by the twiddle factors, and finally solves $dft(\left\{(,r,,,m,·,o,,,m,·,o,)\right\},\mathbf{V}\cup \left\{(,m,,,o,,,o,)\right\},\mathbf{O},\mathbf{O})$ . For performance, the last two steps are not planned independently, but arefused together in a single “twiddle” codelet—a fragment of C code that multiplies its input by the twiddle factors and performs a DFTof size $r$ , operating in-place on $\mathbf{O}$ .