0.10 Implementing ffts in practice  (Page 12/21)

The problem to be solved

In early versions of FFTW, the only choice made by the planner was the sequence of radices [link] , and so each step of the plan took a DFT of a given size $n$ , possibly with discontiguous input/output, and reduced it (via a radix $r$ ) to DFTs of size $n/r$ , which were solved recursively. That is, each step solved the following problem: given a size $n$ , an input pointer $\mathbf{I}$ , an input stride $\iota$ , an output pointer $\mathbf{O}$ , and an output stride $o$ , it computed the DFT of $\mathbf{I}\left[\ell \iota \right]$ for $0\le \ell <n$ and stored the result in $\mathbf{O}\left[ko\right]$ for $0\le k<n$ . However, we soon found that we could not easily express many interesting algorithms within thisframework; for example, in-place ( $\mathbf{I}=\mathbf{O}$ ) FFTs that do not require a separate bit-reversalstage [link] , [link] , [link] , [link] . It became clear that the key issue was not the choice of algorithms, as we hadfirst supposed, but the definition of the problem to be solved. Because only problems that can be expressed can be solved, therepresentation of a problem determines an outer bound to the space of plans that the planner can explore, and therefore it ultimatelyconstrains FFTW's performance.

The difficulty with our initial $(n,\mathbf{I},\iota ,\mathbf{O},o)$ problem definition was that it forced each algorithmic step to address only asingle DFT. In fact, FFTs break down DFTs into multiple smaller DFTs, and it is the combination of these smaller transforms that is best addressed by many algorithmicchoices, especially to rearrange the order of memory accesses between the subtransforms. Therefore, we redefined our notion of a problemin FFTW to be not a single DFT, but rather a loop of DFTs, and in fact multiple nested loops of DFTs. The following sections describe some of the new algorithmic steps thatsuch a problem definition enables, but first we will define the problem more precisely.

DFT problems in FFTW are expressed in terms of structures called I/O tensors, I/O tensors are unrelated to the tensor-product notation used by some other authors to describe FFTalgorithms [link] , [link] . which in turn are described in terms of ancillary structures called I/O dimensions. An I/O dimension $d$ is a triple $d=(n,\iota ,o)$ , where $n$ is a non-negative integer called the length , $\iota$ is an integer called the input stride , and $o$ is an integer called the output stride . An I/O tensor $t=\left\{d_1,,,d_2,,,...,,,d_{\rho }\right\}$ is a set of I/O dimensions. The non-negative integer $\rho =\left|t\right|$ is called the rank of the I/O tensor. A DFT problem , denoted by $dft(\mathbf{N},\mathbf{V},\mathbf{I},\mathbf{O})$ , consists of two I/O tensors $\mathbf{N}$ and $\mathbf{V}$ , and of two pointers $\mathbf{I}$ and $\mathbf{O}$ . Informally, this describes $\left|\mathbf{V}\right|$ nested loops of $\left|\mathbf{N}\right|$ -dimensional DFTs with input data starting at memory location $\mathbf{I}$ and output data starting at $\mathbf{O}$ .

For simplicity, let us consider only one-dimensional DFTs, so that $\mathbf{N}=\left\{(,n,,,\iota ,,,o,)\right\}$ implies a $DFT$ of length $n$ on input data with stride $\iota$ and output data with stride $o$ , much like in the original FFTW as described above. The main newfeature is then the addition of zero or more “loops” $\mathbf{V}$ . More formally, $dft(\mathbf{N},\left\{(,n,,,\iota ,,,o,)\right\}\cup \mathbf{V},\mathbf{I},\mathbf{O})$ is recursively defined as a “loop” of $n$ problems: for all $0\le k<n$ , do all computations in $dft(\mathbf{N},\mathbf{V},\mathbf{I}+k·\iota ,\mathbf{O}+k·o)$ . The case of multi-dimensional DFTs is defined more preciselyelsewhere [link] , but essentially each I/O dimension in $\mathbf{N}$ gives one dimension of the transform.