A Fast Fourier Transform Compiler
Matteo Frigo
MIT Laboratory for Computer Science
545 Technology Square NE43-203
Cambridge, MA 02139


     !"
February 16, 1999
Abstract

/ /

250

The FFTW library for computing the discrete Fourier transform (DFT) has gained a wide acceptance in both academia
and industry, because it provides excellent performance on
a variety of machines (even competitive with or faster than
equivalent libraries supplied by vendors). In FFTW, most of
the performance-critical code was generated automatically
by a special-purpose compiler, called #%$'&%()(+* , that outputs
C code. Written in Objective Caml, #$'&%()(,* can produce
DFT programs for any input length, and it can specialize
the DFT program for the common case where the input data
are real instead of complex. Unexpectedly, #%$-&%()(,* “discovered” algorithms that were previously unknown, and it was
able to reduce the arithmetic complexity of some other existing algorithms. This paper describes the internals of this
special-purpose compiler in some detail, and it argues that a
specialized compiler is a valuable tool.

/

200

/
0

150

100

50

0

/
/
0 0

0

0

0

/ /
0 0

/ /
/

0

/

FFTW
SUNPERF

0 0 / / /
0
0 0
0 /
0 /
0 / /
0 0 0/ 0/

transform size

Figure 1: Graph of the performance of FFTW versus Sun’s Performance Library on a 167 MHz UltraSPARC processor in single
precision. The graph plots the speed in “mflops” (higher is better)
versus the size of the transform. This figure shows sizes that are
powers of two, while Figure 2 shows other sizes that can be factored into powers of 2, 3, 5, and 7. This distinction is important
because the DFT algorithm depends on the factorization of the size,
and most implementations of the DFT are optimized for the case
of powers of two. See [FJ97] for additional experimental results.
FFTW was compiled with Sun’s C compiler (WorkShop Compilers
4.2 30 Oct 1996 C 4.2).

1 Introduction
Recently, Steven G. Johnson and I released Version 2.0 of
the FFTW library [FJ98, FJ], a comprehensive collection of
fast C routines for computing the discrete Fourier transform
(DFT) in one or more dimensions, of both real and complex
data, and of arbitrary input size. The DFT [DV90] is one
of the most important computational problems, and many
real-world applications require that the transform be computed as quickly as possible. FFTW is one of the fastest
DFT programs available (see Figures 1 and 2) because of two
unique features. First, FFTW automatically adapts the computation to the hardware. Second, the inner loop of FFTW

(which amounts to 95% of the total code) was generated automatically by a special-purpose compiler written in Objective Caml [Ler98]. This paper explains how this compiler
works.
FFTW does not implement a single DFT algorithm, but it
is structured as a library of codelets—sequences of C code
that can be composed in many ways. In FFTW’s lingo, a
composition of codelets is called a plan. You can imagine
the plan as a sort of bytecode that dictates which codelets
should be executed in what order. (In the current FFTW
implementation, however, a plan has a tree structure.) The
precise plan used by FFTW depends on the size of the input (where “the input” is an array of 1 complex numbers),
and on which codelets happen to run fast on the underly-

.

The author was supported in part by the Defense Advanced Research
Projects Agency (DARPA) under Grant F30602-97-1-0270, and by a Digital
Equipment Corporation fellowship.

To appear in Proceedings of the 1999 ACM SIGPLAN
Conference on Programming Language Design and Implementation (PLDI), Atlanta, Georgia, May 1999.

1

250

200

150

100

50

2
2

2 2
2

3

2

besides noticing common subexpressions, the simplifier
also attempts to create them. The simplifier is written in
monadic style [Wad97], which allowed me to deal with
the dag as if it were a tree, making the implementation
much easier.

FFTW
SUNPERF

2 2
3
32
2
3 3 2 2
3
2
2
2
2
3 3
32
3
3
2
3
3 2 32
3 3 3
3
3 3

3. In the scheduler, #$'&%()(,* produces a topological sort of
the dag (a “schedule”) that, for transforms of size 4+5 ,
provably minimizes the asymptotic number of register
spills, no matter how many registers the target machine
has. This truly remarkable fact can be derived from
the analysis of the red-blue pebbling game of Hong and
Kung [HK81], as we shall see in Section 6. For transforms of other sizes the scheduling strategy is no longer
provably good, but it still works well in practice. Again,
the scheduler depends heavily on the topological structure of DFT dags, and it would not be appropriate in a
general-purpose compiler.

0

transform size

Figure 2: See caption of Figure 1.

ing hardware. The user need not choose a plan by hand,
however, because FFTW chooses the fastest plan automatically. The machinery required for this choice is described
elsewhere [FJ97, FJ98].
Codelets form the computational kernel of FFTW. You
can imagine a codelet as a fragment of C code that computes
a Fourier transform of a fixed size (say, 16, or 19).1 FFTW’s
codelets are produced automatically by the FFTW codelet
generator, unimaginatively called #%$'&%()(+* . #$'&%()(,* is an
unusual special-purpose compiler. While a normal compiler
accepts C code (say) and outputs numbers, #%$'&%(,(,* inputs
the single integer 1 (the size of the transform) and outputs
C code. The codelet generator contains optimizations that
are advantageous for DFT programs but not appropriate for
a general compiler, and conversely, it does not contain optimizations that are not required for the DFT programs it generates (for example loop unrolling).
#$'&%()(,* operates in four phases.

4. Finally, the schedule is unparsed to C. (It would be easy
to produce FORTRAN or other languages by changing
the unparser.) The unparser is rather obvious and uninteresting, except for one subtlety discussed in Section 7.
Although the creation phase uses algorithms that have
been known for several years, the output of #%$'&%()(+* is at
times completely unexpected. For example, for a complex
transform of size 17698: , the generator employs an algorithm due to Rader, in the form presented by Tolimieri and
others [TAL97]. In its most sophisticated variant, this algorithm performs 214 real (floating-point) additions and 76
real multiplications. (See [TAL97, Page 161].) The generated code in FFTW for the same algorithm, however, contains only 176 real additions and 68 real multiplications, because #%$-&%()(,* found certain simplifications that the authors
of [TAL97] did not notice.2
The generator specializes the dag automatically for the
case where the input data are real, which occurs frequently
in applications. This specialization is nontrivial, and in the
past the design of an efficient real DFT algorithm required
a serious effort that was well worth a publication [SJHB87].
#%$-&%()(,* , however, automatically derives real DFT programs
from the complex algorithms, and the resulting programs
have the same arithmetic complexity as those discussed by
[SJHB87, Table II].3 The generator also produces real variants of the Rader’s algorithm mentioned above, which to my
knowledge do not appear anywhere in the literature.

1. In the creation phase, #%$'&%(,(,* produces a directed
acyclic graph (dag) of the codelet, according to some
well-known algorithm for the DFT [DV90]. The generator contains many such algorithms and it applies the
most appropriate.
2. In the simplifier, #%$'&%()(+* applies local rewriting rules
to each node of the dag, in order to simplify it. In traditional compiler terminology, this phase performs algebraic transformations and common-subexpression elimination, but it also performs other transformations that
are specific to the DFT. For example, it turns out that if
all floating point constants are made positive, the generated code runs faster. (See Section 5.) Another important transformation is dag transposition, which derives
from the theory of linear networks [CO75]. Moreover,

2 Buried somewhere in the computation dag generated by the algorithm
are statements of the form ;=<?> @BA , C<?>EDFA , GH<?;I@JC . The
generator simplifies these statements to GK<ML
> , provided ; and C are
not needed elsewhere. Incidentally, [SB96] reports an algorithm with 188
additions and 40 multiplications, using a more involved DFT algorithm that
I have not implemented yet. To my knowledge, the program generated by
N OQP RSRST performs the lowest known number of additions for this problem.
3 In fact, N OQP RSRST saves a few operations in certain cases, such as UV<JWX .

1 In the actual FFTW system, some codelets perform more tasks, however. For the purposes of this paper, we consider only the generation of
transforms of a fixed size.

2

not describing some advanced number-theoretical DFT algorithms used by the generator.) Readers unfamiliar with the
discrete Fourier transform, however, are encouraged to read
the good tutorial by Duhamel and Vetterli [DV90].
The rest of the paper is organized as follows. Section 2
gives some background on the discrete Fourier transform
and on algorithms for computing it. Section 3 overviews related work on automatic generation of DFT programs. The
remaining sections follow the evolution of a codelet within
#%$-&%()(,* . Section 4 describes what the codelet dag looks like
and how it is constructed. Section 5 presents the dag simplifier. Section 6 describes the scheduler and proves that it minimizes the number of transfers between memory and registers. Section 7 discusses some pragmatic aspects of #$'&%()(,* ,
such as running time and memory requirements, and it discusses the interaction of #%$'&()(,* ’s output with C compilers.

I found a special-purpose compiler such as the FFTW
codelet generator to be a valuable tool, for a variety of reasons that I now discuss briefly.

Y Performance was the main goal of the FFTW project,
and it could not have been achieved without #%$-&%()(,* .
For example, the codelet that performs a DFT of size 64
is used routinely by FFTW on the Alpha processor. The
codelet is about twice as fast as Digital’s DXML library
on the same machine. The codelet consists of about
2400 lines of code, including 912 additions and 248
multiplications. Writing such a program by hand would
be a formidable task for any programmer. At least for
the DFT problem, these long sequences of straight-line
code seem to be necessary in order to take full advantage of large CPU register sets and the scheduling capabilities of C compilers.

Y Achieving correctness has been surprisingly easy. The
DFT algorithms in #%$'&%(,(,* are encoded straightfor-

2 Background
This section defines the discrete Fourier transform (DFT),
and mentions the most common algorithms to compute it.
Let ] be an array of 1 complex numbers. The (forward)
discrete Fourier transform of ] is the array ^ given by

wardly using a high-level language. The simplification phase transforms this high-level algorithm into optimized code by applying simple algebraic rules that
are easy to verify. In the rare cases during development when the generator contained a bug, the output
was completely incorrect, making the bug manifest.

d)ef
g
e ph
^`_ acb6
] _ m'bon d q
l
hjik

Y Rapid turnaround was essential to achieve the perfor-

(1)

efjwjd

is a primitive 1 -th root of unity, and
where n d 6Zrtsju,v
x=y
az{1 . In case ] is a real vector, the transform ^ has
the hermitian symmetry

mance goals. For example, the scheduler described in
Section 6 is the result of a trial-and-error process in
search of the best result, since the schedule of a codelet
interacts with C compilers (in particular, with register
allocators) in nonobvious ways. I could usually implement a new idea and regenerate the whole system within
minutes, until I found the solution described in this paper.

^|_ 1F}Jacb 6M^~'_ acb q
where ^ ~ _ acb is the complex conjugate of ^ ~ _ ab .
The backward DFT flips the sign at the exponent of n d ,
and it is defined in the following equation.

Y The generator is effective because it can apply problem-

def
g
ph
^`_ acb6
] _ m'bn dƒ‚
€
hjik

specific code improvements. For example, the scheduler is effective only for DFT dags, and it would perform poorly for other computations. Moreover, the simplifier performs certain improvements that depend on
the DFT being a linear transformation.

(2)

The backward transform is the “scaled inverse” of the forward DFT, in the sense that computing the backward transform of the forward transform yields the original array multiplied by 1 .
If 1 can be factored into 1„6 1 f 1 , Equation (1) can be
s
f
f
rewritten as follows. Let m†6ƒm f 1
sE‡ m s , and aˆ6‰a ‡ a s 1 .
We then have,

Y Finally, #%$'&%(,(,* derived some new algorithms, as in
the example 1Z6[8\: discussed above. While this paper does not focus on these algorithms per se, they are
of independent theoretical and practical interest in the
Digital Signal Processing community.

^`_ a f ‡ a s 1
d-Š ef ‹ŒŽ
g
h Š ik

This paper, of necessity, brings together concepts from
both the Programming Languages and the Signal Processing
literatures. The paper, however, is biased towards a reader
familiar with compilers [ASU86], programming languages,
and monads [Wad97], while the required signal-processing
knowledge is kept to the bare minimum. (For example, I am

fb 6
(3)
d+‘ef
g
e p’ h 

e p• h Š
e poŠ h Š
—– n d-Š
‚
]l_ m f 1 sE‡ m s bn d+ ”“ n d
h  ik

This formula yields the Cooley-Tukey fast Fourier transform algorithm (FFT) [CT65]. The algorithm computes 1 s
3

transforms of size 1 f (the inner sum), h it multiplies the result
e p• Š
by the so-called twiddle factors n d
, and finally it computes 1 f transforms of size 1 (the outer sum).
s
If ˜+™\š›œ1 f q 1
s 6ž8 , the prime factor algorithm can be
applied, which avoids the multiplications by the twiddle factors at the expense of a more involved computation of indices. (See [OS89, page 619].) If 1 is a multiple of Ÿ ,
the split-radix algorithm [DV90] can save some operations
with respect to Cooley-Tukey. If 1 is prime, #%$'&%()(+* uses
either Equation (1) directly, or Rader’s algorithm [Rad68],
which converts the transform into a circular convolution of
size 1 }{8 . The circular convolution can be computed recursively using two Fourier transforms, or by means of a
clever technique due to Winograd [Win78] (#%$'&%()(+* does not
employ this technique yet, however). Other algorithms are
known for prime sizes, and this is still the subject of active
research. See [TAL97] for a recent compendium on the topic.
Any algorithm for the forward DFT can be readily adapted to
compute the backward DFT, the difference being that certain
complex constants become conjugate. For the purposes of
this paper, we do not distinguish between forward and backward transform, and we simply refer to both as the “complex
DFT”.
In the case when the input is purely real, the transform
can be computed with roughly half the number of operations
of the complex case, and the hermitian output requires half
the storage of a complex array of the same size. In general, keeping track of the hermitian symmetry throughout
the recursion is nontrivial, however. This bookkeeping is
relatively easy for the split-radix algorithm, and it becomes
particularly nasty for the prime factor and the Rader algorithms. The topic is discussed in detail in [SJHB87]. In the
real transform case, it becomes important to distinguish the
forward transform, which takes a real input and produces an
hermitian output, from the backward transform, whose input
is hermitian and whose output is real, requiring a different
algorithm. We refer to these cases as the “real to complex”
and “complex to real” DFT, respectively.
In the DFT literature, unlike in most of Computer Science, it is customary to report the exact number of arithmetic operations performed by the various algorithms, instead of their asymptotic complexity. Indeed, the time complexity of all DFT algorithms of interest is ¡¢›œ1V£¤'˜ˆ1 , and

a detailed count of the exact number of operation is usually
doable (which by no means implies that the analysis is easy
to carry out). It is no problem for me to follow this convention in this paper, since #$'&%()(,* produces an exact arithmetic
count anyway.
In the literature, the term FFT (“fast Fourier transform”) refers to either the Cooley-Tukey algorithm or to any
¡¢›œ1V£¤'˜E1  algorithm for the DFT, depending on the author.
In this paper, FFT denotes any ¡¢›¥1£¤'˜¦1 algorithm.

3 Related work
Researchers have been generating FFT programs for at least
twenty years, possibly to avoid the tedium of getting all the
implementation details right by hand. To my knowledge, the
first generator of FFT programs was FOURGEN, written by
J. A. Maruhn [Mar76]. It was written in PL/I and it generated
FORTRAN.4 FOURGEN is limited to transforms of size 4+5 .
Perez and Takaoka [PT87] present a generator of Pascal
programs implementing a prime factor FFT algorithm. This
program is limited to complex transforms of size 1 , where
1 § must be factorable into mutually prime factors in the set
4 q : q Ÿ q‘¨q‘©qjª«qj¬%q 8\­%® .
Johnson5 and Burrus [JB83] applied dynamic programming to the automatic design of DFT modules. Selesnick
and Burrus [SB96] used a program to generate MATLAB
subroutines for DFTs of certain prime sizes. In many cases,
these subroutines are the best known in terms of arithmetic
complexity.
The EXTENT system by Gupta and others [GHSJ96] generates FORTRAN code in response to an input expressed in
a tensor product language. Using the tensor product abstraction one can express concisely a variety of algorithms
that includes the FFT and matrix multiplication (including
Strassen’s algorithm).
Another program called #$'&%()(,* generating Haskell FFT
subroutines is part of the &«¯,(°\± benchmark for Haskell
[Par92]. Unlike my program, this #%$'&()(,* is limited to transforms of size 4+5 . The program in &«¯+(°± is not documented
at all, but apparently it can be traced back to [HV92].
Veldhuizen [Vel95] used a template metaprograms technique to generate ²)³)³ programs. The technique exploits the
template facility of ²)³,³ to force the ²)³)³ compiler to perform
computations at compile time.
All these systems are restricted to complex transforms,
and the FFT algorithm is known a priori. To my knowledge, the FFTW generator is the only one that produces real
algorithms, and in fact, which can derive real algorithms by
specializing a complex algorithm. Also, my generator is the
only one that addressed the problem of scheduling the program efficiently.
4 Creation of the expression dag
This section describes the creation of an expression dag. We
first define the &«¯,´$ data type, which encodes a directed
acyclic graph (dag) of a codelet. We then describe a few
4 Maruhn argues that PL/I is more suited than FORTRAN to this
program-generation task, and has the following curious remark.

One peculiar difficulty is that some FORTRAN systems produce an output format for floating-point numbers without the
exponent delimiter “E”, and this makes them illegal in FORTRAN statements.



5 Unrelated

4

to Steven G. Johnson, the other author of FFTW.

Y A split-radix algorithm [DV90], if 1 is a multiple of Ÿ .

*)µ,¶%¹H$·º,&%»-¼ ¯,´$ƒ¸ º,»'¼
»'¼
¯,(
±«$+½¿¾À& ±«$+½
¹`Á
¹|Ç ¯,Â,´M¯,(MÃÂ+½ °-Â'±«Ä,$¦¾ÆÅ%Â'½ °-Â'±Ä+$
¹HÊ *»¯+½%Ë $ȯ,(ÈÃ%Â+½ °-Â-±ËÄ,$¦¾ÆÅÂ+½ °-Â'±«Ä+$„É€&%¯,´$
¹HÌ Ä ¼ Ë ¯,(È&%¯,´$MÄ«° *
¹¢Î+¼ ° $ » Ë ¯,(ƒ&«¯,´$ÍÉl&«¯,´$
°&
¯,(·&«¯,´)$
Figure 3: Objective Caml code that defines the
which encodes an expression dag.

Ï'Ð\ÑtÒ

Otherwise,

Y A prime factor algorithm (as described in [OS89, page
619]), if 1 factors into 1 f 1 s , where 1 pß6 Þ 8 and
˜+™\š›œ1 f q 1 s  6{8 . Otherwise,
Y The Cooley-Tukey FFT algorithm (Equation (3)) if 1
Þ 8 . Otherwise,
factors into 1 f 1 , where 1 p6{
s
Y (1 is a prime number) Rader’s algorithm for transforms
of prime length [Rad68] if 1=6 ¨ or 1JàÍ8\: . Otherwise,

data type,

ancillary data structures and functions that provide complex
arithmetic. Finally, the bulk of the section describes the function ()(,*,#%$'& , which actually produces the expression dag.
We start by defining the &«¯,´$ data type, which encodes
an arithmetic expression dag. Each node of the dag represents an operator, and the node’s children represent the
operands. This is the same representation as the one generally used in compilers [ASU86, Section 5.2]. A node in
the dag can have more than one “parent”, in which case the
node represents a common subexpression. The definition
of &«¯,´)$ is given in Figure 3, and it is straightforward. A
nodeº,is»'¼ either a»'¼ real number (encoded by the abstract data
±«$'½¿¾À& ±«$+½ ), a load of an input variable, a store of
type
an expression into an output node, the sum of the children
nodes, the product of two nodes, or the sign negation of a
node. For example, the expressionÊ Ó» }l
ËÖÔ ,ÕÆÁ where Ó and
Î+¼ Ô are
» Ë
input
variables,
is
represented
by
Ä
)
¯
,
Â
ƒ
´
Ó
×
°&
ØQÁ
¯,Â,´MÔÚÙtÛ .
º,»-¼
The structure
±«$+½ maintains floating-point constants
with arbitrarily high precision (currently, 50 decimal digits),
in case the user wants to use the quadruple precisionº,»'floating¼
point unit on a processor such as the UltraSPARC.
±«$'½ is
implemented on top of Objective Caml’s arbitrary-precision
rationals. The structure Ã%Â+½ °-Â-±Ä,$ encodes the input/output
nodes of the dag, and the temporary variables of the generated C code. Variables can be considered an abstract data
type that is never used explicitly in this paper.
The &«¯,´$ data type encodes expressions over real numbers, since the final C output contains only real expressions.
For creating the expression dag of the codelet, however, it
is convenient to express the algorithms in terms of complex
¼ numbers. The generator contains a structure called
²¯ ¶Ä,$+Ü , which implements complex expressions on top of
the &«¯,´$ data type, in a straightforward way.6 The type
¼
²¯ ¶Ä,$+Ü¿¾Ý$+Ü,¶½ (not shown) is essentially a pair of &«¯,´$ s.
We now describe the function ()(+*)#%$'& , which creates a
dag for a DFT of size 1 . In the current implementation,
()(,*,#%$'& uses one of the following algorithms.

Y Direct application of the definition of DFT (Equation (1)).
We now look at the operation of ()(,*)#%$-& more closely.
The function has type

Ø
()(+*)#%$'&ßáJ°&*Zâ)ã Ø °&*Íâ)ヲ¯
°&*Zâ)ã ° &*Íâ)ヲ¯

¼

¼ ¶Ä,$+Ü¿¾Ý$+Ü,¶½ Ù â)ã
¶Ä,$+Ü¿¾Ý$+Ü,¶½ Ù

The first argument to ()(,*,#%$'& is the size & of
» the trans
°
)
&
¶
form. The second
argument
is
a
function
¼
Ø * » with type
°&)*Zâ,ãȲ¯ ¶Ä,$+ÜE¾$+Ü,¶½ . The application °&)¶ *{°)Ù returns a complexË expression that contains the ° -th input. The
third argument °t#,& is either 8 or }ä8 , and it determines the
direction of the transform.
Depending on the size 1 of the requested transform,
()(+*)#%$'& dispatches one of the algorithms mentioned above.
We now discuss how the Cooley-Tukey FFT algorithm is implemented. The implementation of the other algorithms is
similar, and it is not shown in this paper.
The Objective Caml code that implements the CooleyTukey algorithm can be found in Figure 4. In order to understand the code, recall Equation (3). This equation translates almost verbatim into Objective Caml.
With reference
¼
to Figure 4, the function application * ¶åçæ)è computes the
inner sum of Equation (3)
¼ for a given value of m s , and it returns a function of a f . (* ¶å is curried over a f , and therefore
Ø ¼
a f does not appear explicitly in the definition.) Next, * ¼ ¶å
æ)è{°%å+Ù is multiplied by the twiddle factors, yielding * ¶%è ,
that
¼ is, the expression in square braces in Equation (3). Next,
* ¶«é computes the outer
summation, which is itself a DFT
¼
of size 1 . (Again, * ¶«é is a function of a f and a , curried
s
s
over a .) In order to obtain the a -th element of the output of
s
f
the
Ø ¼ transform, the index a is finally mapped into a and a s and
* ¶«éZ°%å·°-èÙ is returned.
Observe that the code in Figure 4 does not actually perform any computation. Instead, it builds a symbolic expression dag that specifies the computation. The other DFT algorithms are implemented in a similar fashion.
At the top level, the generator
invokes ()(,*,#%$'& with
Ë
the size » & and the direction °t#,» & specified ¼ by the user.
°&)¶ * function
The
is set to ( &7°·â)ヲ¯ ¶Ä+$+Ü¿¾jÄ,¯,Â,´
Ø
»
ÃÂ+½ °-Â'±«Ä,$¦¾S°&,¶ *È°)Ù , i.e., a function that loads the a -th
input variable. Recall now that ()(,*,#%$'& returns a function

6 One subtlety is that a complex multiplication by a constant can be implemented with either 4 real multiplications and 2 real additions, or 3 real
multiplications and 3 real additions [Knu98, Exercise 4.6.4-41]. The current generator uses the former algorithm, since the operation count of the
dag is generally dominated by additions. On most CPUs, it is advantageous
to move work from the floating-point adder to the multiplier.

5

ȓ
» Ë
Ä,$+*Ƚ%$%ê‰
¼ ê'¯,¯Ä,$+µ«ë-* $+µç&å &«è{°&)¶ * °t#,&‰¸
Ä,$+*M
Ø »* ¶åíæ)èȸM()(,*)» #%$-&ƒ
Ø &å
Ë
( ¼&Íæ å„â,ã{°&,¶ * æ å·Él&«èM³‰æ)èÙ,Ù °t#,&7°&
Ä,$+*M* ¶«èÍØî°%Ë ålæ)èȸ
¼
$+Ü+¶‰¼ &
°t#+&{ÉZ°%å„É·æ)èÙ„Ø ï«¼ Ʉ* ¶åçæ)èZ
Ë °%僰&
Ä,$+*M* ¶«éÍ°%åç¸M()(,*)#%$-&ƒ&«è * ¶«èZ°%å+Ù °#,&
°& Ø »
¼
Ø ¼
Ø
( &7°·â)ã·* ¶«é ° ¯,´ƒ&å+Ù °·ð„&å+Ù)Ù

¼
Á Ø
»
* ¼ 
¶ åí¸ Á Ø °&,¶ »
* ¼ ¶ ƒ¸ Ø °&,¶ »
* ¼ ¶ ƒ¸ Ø °&,¶ »
* ¼ ¶«èƒ
¸ Á Ø °&,¶ »
* ¼ ¶«éƒ
¸  Á Ø °&,¶ »
* ¼ ¶ ƒ
¸  Ø °&,¶ »
* ¼ ¶ ƒ
¸  Ø °&,¶ »
* ¶ M
¸

Á Ø » »  Õ°&,¶

¯ *,¶ Ø,Ø * ÛÙ ¼
Ø ³» » å„ÕɄ*
 ¯ *,¶ *
ÛÙ ¼
Ø,Ø
³
„
å
„
É
Á Ø » » Õ *

¯ *,¶ Ø,Ø * åÛÙ
Ø ³» » âå„
Õ É„*
 ¯ *,¶ * åÛÙ
Ø,Ø
³
âå„É„*

Figure 4: Fragment of the FFTW codelet generator that implements the Cooley-Tukey FFT algorithm. The infix operator ñÚò
computes the complex product.
 The function ÒóôõÏJö computes
the constant exp ÷¥ø ùú,û üþý ÿ
 .

» »
Ø » »
¯ +* ¶ * , where ¯ *+¶ *{°)Ù is a complex expression that
computes the a -th element of the output array. The top level
Ç
Ø » »
builds a list of *%¯+½%$ expressions that store ¯ *,¶ *{°)Ù
xÈy
into the a -th output variable, for all
aFzß1 . This list
Ç
of *%¯+½$ s is the codelet dag that forms the input of subse-

Õ
* Õ 
Û Ù×
* Õ ÛÙ×
* Õ ÛÙ×
* Õ ÛÙ×
* Õ å\ÛÙ×
* Õ å\ÛÙ×
* Õ å\ÛÙ×
* å\Ø)ÛØ Ù×
¸
å·É€
Ø *
¶«éل
â
Ø)Ø
¸
å·É€
*
Ø
¶ Ùl
³
Ø)Ø
å·É€Ø *
¼¸
¶ Ø)ـ
Ø â
¸¼
å·É€Ø *
¶ Ùç³
¼

Ø
¶ å+¼ لâ

¼ ɀ* ¶Ù)ÙØ ×
¶%è¼ Ùl³
ɼ €* ¶«é«Ù)ÙØ ×
¶ ل

¼ â
„
É
*
¶ Ù)Ø Ù×
¼
¶ Ùl¼ ³

Ʉ* ¶ Ù)Ù×

¼
ɀ* «¶ èÙ,Ù
¼
ɀ* ¶å+Ù,Ù
¼
ɀ* ¶ Ù,Ù
¼
ɀ* ¶ Ù,Ù

Figure 5: C translation of a dag for a complex DFT of size 2,
as generated by ÚÒîÏ . Variable declarations have been omitted
from the figure. The code contains many common subexpression
(e.g., 
 ô and 
'ô! ), and redundant multiplications by " or ý .

quent phases of the generator.
We conclude this section with a few remarks. According
to the description given in this section, ()(,*,#%$'& contains no
special support for the case where the input is real. This
statement is not completely true. In the actual implementation, (,(,*)#%$'& maintains certain symmetries explicitly (for
example, if the input is real, then the output is known to
have hermitian symmetry). These additional constraints do
not change the final output, but they speed up the generation process, since they avoid computing and simplifying the
same expression twice. For the same reason, the actual im¼
plementation memoizes expressions such as * ¶僰Úè?°%å in
Figure 4, so that they are only computed once.7
At this stage, the generated dag contains many redundant
x
computations, such as multiplications by 8 or , additions
x
of , and so forth. ()(,*)#%$-& makes no attempt to eliminate
these redundancies. Figure 5 shows a possible C translation
of a codelet dag at this stage of the generation process.

applies a series of local improvements to every node. For explanation purposes, these improvements can be subdivided
into three categories: algebraic transformations, commonsubexpression elimination, and DFT-specific improvements.
Since the first two kinds are well-known [ASU86], I just discuss them briefly. We then consider the third kind in more
detail.
Algebraic transformations reduce the arithmetic complexity of the dag. Like a traditional compiler, the simplifier performs constant folding, and it simplifies multiplications by
x
x
, 8 , or }ä8 , and additions of . Moreover, the simplifier
applies the distributive property systematically. Expressions
of the form #$ ‡ #% are transformed into #›&$ ‡ % . In the

same way, expressions of the form # f $ ‡ # $ are transformed
s
into ›'# f ‡ # $ . In general, these two transformations have
st
the potential of destroying common subexpressions, and they
might increase the operation count. This does not appear to
be the case for all DFT dags I have studied, although I do not
fully understand the reason for this phenomenon.
Common-subexpression elimination is also applied systematically. Not only does the simplifier eliminate common
subexpressions, it also attempts to create new ones. For example, it is common for a DFT dag (especially in the case of
real input) to contain both $I}(% and %E}($ as subexpressions,
for some $ and % . The simplifier converts both expressions
to either $ })% and }K›&$ }*% , or }K›+%H})$ and %†}*$ , de

pending on which expression is encountered first during the
dag traversal.
The simplifier applies two kinds of DFT-specific improvements. First, all numeric constants are made positive, possibly propagating a minus sign to other nodes of the dag. This
curious transformation is effective because constants generally appear in pairs # and },# in a DFT dag. To my knowl-

5 The simplifier
In this section, we present FFTW’s simplifier, which transforms code such as the one in Figure 5 into simpler code.
The simplifier transforms a dag of &«¯+´$ s (see Section 4) into
another dag of &«¯+´$ s. We first discuss how the simplifier
transforms the dag, and then how the simplifier is actually
implemented. The implementation benefits heavily from the
use of monads.
5.1 What the simplifier does
We first illustrate the improvements applied by the simplifier
to the dag. The simplifier traverses the dag bottom-up, and it
7 These performance improvements were important for a user of FFTW
who needed a hard-coded transform of size 101, and had not obtained an
answer after the generator had run for three days. See Section 7 for more
details.

6

2
0

.
/

1

ø

3

2
0

-

.
/

1

ø

3

adds
muls
size (not transposed)
complex to complex
5
32
16
10
84
32
13
176
88
15
156
68
real to complex
5
12
8
10
34
16
13
76
44
15
64
31
complex to real
5
12
9
9
32
20
10
34
18
12
38
14
13
76
43
15
64
37
16
58
22
32
156
62
64
394
166
128 956
414

-

Figure 6: Illustration of “network” transposition. Each graph defines an algorithm for computing a linear function. These graphs
are called linear networks, and they can be interpreted as follows.
Data are flowing in the network, from input nodes to output nodes.
An edge multiplies data by some constant (possibly ý ), and each
node is understood to compute the sum of all incoming
edges. In
354
2 076 . /
this
example, the network on the left computes
and
4
1
ø 086 / . The network on the right is the “transposed” form
of the first network, obtained by reversing
The
net4
3 all edges.
4 new
3
2
. 6 1 work computes the linear function 0 4:9 6 ø and /
9 .
In general, if a network computes 0 4;9=< / for some matrix ,
0 . (See [CO75] for a
the transposed network computes /
proof.) These linear networks are similar to but not the same as ex Ò\Ï> , because
pression dags normally used in compilers and in t
in the latter case the nodes and not the edges perform computation.
A network can be easily transformed into an expression dag, however. The converse is not true in general, but it is true for DFT dags
where all multiplications are by constants.

edge, every C compiler would store both # and },# in the
program text, and it would load both constants into a register
at runtime. Making all constants positive reduces the number
of loads of constants by a factor of two, and this transformation alone speeds up the generated codelets by 10-15% on
most machines. This transformation has the additional effect
of converting subexpressions into a canonical form, which
helps common-subexpression elimination.
The second DFT-specific improvement is not local to
nodes, and is instead applied to the whole dag. The transformation is based on the fact that a dag computing a linear function can be “reversed” yielding a transposed dag
[CO75]. This transposition process is well-known in the Signal Processing literature [OS89, page 309], and it operates a
shown in Figure 6. It turns out that in certain cases the transposed dag exposes some simplifications that are not present
in the original dag. (An example will be shown later.) Accordingly, the simplifier performs three passes over the dag.
It first simplifies the original dag ? yielding a dag @ . Then,
it simplifies the transposed dag @7A yielding a dag BCA . Finally, it simplifies B (the transposed dag of B A ) yielding a
dag D .8 Figure 7 shows the savings in arithmetic complexity that derive from dag transposition for codelets of various
sizes. As it can be seen in the figure, transposition can reduce the number of multiplications, but it does not reduce
the number of additions.
Figure 8 shows a simple case where transposition is bene-

adds muls
(transposed)
32
84
176
156

12
24
68
56

12
34
76
64

6
12
34
25

12
32
34
38
76
64
58
156
394
956

7
18
14
10
35
31
18
54
146
374

Figure 7: Summary of the benefits of dag transposition. The table
shows the number of additions and multiplications for codelets of
various size, with and without dag transposition. Sizes for which
the transposition has no effect are not reported in this table.

ficial. The network in the figure computes EI6·ŸGF-›4'Ó ‡ :+Ô .

It is not safe to simplify this expression to EB6 ª Ó ‡ 8t4'Ô ,
since this transformation destroys the common subexpressions 4'Ó and :+Ô . (The transformation destroys one operation
and two common subexpressions, which might increase the
operation count by one.) Indeed, the whole point of most
FFT algorithms is to create common subexpressions. When
the network is transposed, however, it computes Ó 6 48FڟHE
and ÔK67:IFtŸHE . These transposed expressions can be safely
transformed into ӆ6 ª E and ÔI6784JE because each transformation saves one operation and destroys one common subexpression. Consequently, the operation count cannot increase.
In a sense, transposition provides a simple and elegant way
to detect which dag nodes have more than one parent, which
would be difficult to detect when the dag is being traversed.
5.2 Implementation of the simplifier
The simplifier is written in monadic style [Wad97]. The
monad performs two important functions. First, it allows
the simplifier to treat the expression dag as if it were a tree,
which makes the implementation considerably easier. Second, the monad performs common-subexpression elimination. We now discuss these two topics.
Treating dags as trees. Recall that the goal of the sim-

8 Although one might imagine iterating this process, three passes seem
to be sufficient in all cases.

7

Ë ¼
Ä,$'¼ *M¼½$«êƒÂ)Ä'# ° ¶·Ü‰¸
$ Ø ¯«» >° R°&#
( & º+êÚ»'*¼ °-¯'&
Ë »-¼
¹HÊ » Ë ÂÍâ)ã & ‰Â
Ä ¼ ÂÍâ,ã Ë ¼
Ë » Ë
¹HÌ ¼ Â'¶ Ë ÍØ Â)Ä'# ° ¶ȉã)ã+¸ ¶Ä 
° $ Ë ¼T
 S ±لâ)ã
»
ÂÄ'# ° Ë ¶ ¼ ƒÂÍã)ã'¸‰( &Z» ÂPM=â)ã
ÂÄ'# Ë ° ¼ ¶ „
± M â)ã
Ë ±?Ø ã)ã+¸È( &ÈO

*
°
$

Â
N
M!S|O
±
M‘Ù
¹¢Î'¼ » Ë
°& Ë ¼ ‰â)ã
Ët»-¼ » Ë

Â
°& 
¹|Ç Ä'# ° Ø ¶ ƒÂÍã)ã'¸
*%¯+½%$ Ë ¼U
Å S Âلâ)ã
»
ÂÄ'# ° » ¶ ƒÂÍؑã)Ç ã'¸‰( &Z
 M=â)ã
Ø P
%
½
+
$
*
+
½
&


%
*
'
¯
%
½
$
Å
P
S=N
 M‘Ù,Ù
¹
»
ÜZâ)ã·½%$+* ½,& ·Ü Ù
Ü

;
L
C

K

>

L

Figure 8: A linear network where which dag transposition exposes
some optimization possibilities. See the text for an explanation.

plifier is to simplify an expression dag. The simplifier, however, is written as if it were simplifying an expression tree.
The map from trees to dags is accomplished by memoization,
which is performed implicitly by a monad. The monad maintains a table of all previously simplified dag nodes, along
with their simplified versions. Whenever a node is visited
for the second time, the monad returns the value in the table.
In order to continue reading this section, you really should
be familiar with monads [Wad97]. In any case, here is a
very brief summary on monads. The idea of a monadic-style
program is to convert all expressions of the form

Ä,$+*ÈÜ͸ÍÂ{°&

Ø

Figure 9: The top-level simplifier function VWX!YZ['ô\ , written in
monadic style. See the text for an explanation.

construction from [ASU86, page 592]. The monad maintains a table of all nodes produced during the traversal of
the dag.
Each time a new node is constructed and returned,
»
in the
½%$'* ½,& checks whether the node appears elsewhere
»
dag. If so, the new node is discarded and ½$+* ½,& returns
the old node. (Two nodes are considered the same if they
compute equivalent expressions. For example, Ó
‡ Ô is the
same as Ô
.)
Ó
‡
It is worth remarking that the simplifier interleaves
common-subexpression elimination with algebraic transformations. To see why interleaving is important, consider for
example the expression ÓH}„Ó] , where Ó and Ó] are distinct
nodes of the dag that compute the same subexpression. CSE
x
rewrites the expression to Ó }KÓ , which is then simplified to .
This pattern occurs frequently in DFT dags.

±MÜ Ù

into something that looks like

»
»
Ø
‰ã)ã+¸M( M
& Ü{â,ヽ%$'* ½,& ±ÍÜÙ

The code should
Ø be read “call ( , and then name the result
Ü and return ±‰Ü Ù .” The advantage of this transformation
is that the meanings of “then” (the infix operator ã)ã+¸ ) and
»
“return” (the function ½%$'* ½,& ) can be defined so that they
perform all sorts of interesting activities, such as carrying
state around, perform I/O, act nondeterministically, etc. In
the specific case of the FFTW simplifier, ã)ã+¸ is defined so
as to» keep track of a few tables used for memoization, and
½%$+* ½,& performs common-subexpression elimination.
Ë ¼
The core of the simplifier
Ë ¼ is the function ÂÄ'# ° ¶ , as
shown in Figure 9. ÂÄ'# ° ¶ dispatches on the argument
Ü (of type &«¯,´$ ), and it calls a simplifier function for the
appropriate case. If the node has subnodes, the subnodes
are
Ì ¼ Ë
simplified
first.
For
example,
suppose
Ü
is
a
°
$
node.
Ì ¼ Ë
Since
Ë a¼ ° $ node has two subnodes  and ± , the function
ÂÄ'# ° ¶ first calls itself recursively
Ë ¼ on  , yielding ÂNM , and
then on ± Ë , yielding
OM . Then, ÂÄ'# ° ¶ passes control to the
±
¼ Ë
Ë ¼ Ë
function *° $  . If both ÂPM and ±OM are constants, * ° $ 
Ë ¼ Ë
computes the product directly. In the same way, * ° $ 
x
takes care of the case
either ÂNM or ±QM is or 8 , and so
Ë ¼ where
Ë
on. The code for *° $  is shown in Figure 10.
The neat trick of using memoization for graph traversal
was invented by Joanna Kulik in her master’s thesis [Kul95],
as far as I can tell.
Common-subexpression elimination (CSE) is » performed behind the scenes by the monadic operator ½%$+* ½,& .
The CSE algorithm is essentially the classical bottom-up

6 The scheduler
In this section we discuss the #%$'&()(,* scheduler, which produces a topological sort of the dag in an attempt to maximize
register usage. For transforms whose size is a power of 4 , we
prove that a schedule exists that is asymptotically optimal in
this respect, even though the schedule is independent of the
number of registers. This fact is derived from the red-blue
pebbling game of Hong and Kung [HK81].
Even after simplification, a codelet dag of a large transform still contains hundreds or even thousands of nodes, and
there is no way to execute it fully within the register set of
any existing processor. The scheduler attempts to reorder the
dag in such a way that register allocators commonly used in
compilers [Muc97, Section 16] can minimize the number of
register spills. Note that the FFTW codelet generator does
not address the instruction scheduling problem; that is, the

8

¼ Ë
»
Ä,$+*È
¹€½%ØÝÎ'$%¼ ê *» °Ë $ ·¸Í( &êÚØ * °-¯'&
Ø
 S Ø ±Ù·â)ã Ƀâ,ÂZ
Ë °& ¼ Ë T
Ët»'¼ ɀ± » Ë ¸)¸%ãÈâ ÂZɀ±كÉÙ
°& 
¹€Ø * Î+° ¼ $ » Ë Â^SF±لã)Ø ã+¸
Ø
ÂTË S ¼ °& Ë Ø ±Ù·â)ã É·Â{Ët»'Ƀ¼ â-± » Ë ¸)¸%ãÈâ ÂZɀ±كÉÙ
¹€Øjº+»'¼ * ° $ º+ »'¼ Â^SF±لã)ã+Ø ¸ ¼«» °& 
»'¼
Ë
 S غ,»'¼ ±كâ,¼«ã » É
Ä-* °¶Ä'µ„*_«¯„& ±%$+½ ÉÙ
Ë »'¼^
¹€Øjº+»'¼&  Ì ¼ ±«Ë $'½¿Øjº,¾ »-¼ Älƒ± Ù
 S غ,° »'¼ $
± S ê)Ù)كâ,ã »
Ë »'¼^
¼«» N
& Ë  ¼ Ë ±«Ø $'½¿¾ Älƒ± كã)ã'¸M( &ÍÜZâ)ã
Ü Sõº,ê»'Ù ¼
¹€Øjº+»'¼ *° $  P
Ë
Â
^
S ± )
Ù
_`%$'&
±«$'½¿Ø ¾S° !ë R$+½%¯lÂÍâ,ã
Ë »'¼ º+»'¼
±«$+½E¾ Rº,$+»'½%¼ ¯
É Ë É€± ¸)¸ã
ÉÙ
¹€Øjº+»'¼& 
» S ± )
^
Ù _`%$'&
±«$'½¿Ø ¾S° ë+¯'&«$lÂÍâ)ã
É{Ë å·¼ ɀ± ¸)¸ㄱ ÉÙ
¹€Øjº+»'½%¼ $+* ½+& „±
º,»'¼
Â
^
S ± )
Ù
_`%$'&
«
±
'
$
¿
½
S
¾
Ët»'¼ » Ë
Ø ° ë ¯'&«$lÂÍâ,ã

°
&
„±
¹€Ø
Øjº+»'¼
Ë
Ë Éȼ âå·Ë Él± Ø ¸,¸%ãÈâ-± É)Ù
Â
T
S
M
ë
Â
±
O
M‘Ù,ـâ)ã
± M!S=ÂÙ
¹€Ø
»
ØjÌ ¼ * Ë ° Ø $  O
 S ±لâ)ã·½%$+* ½,& 
T
° $
 SF± Ù)Ù
T

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

ih

Ë

Figure 11: Illustration of the scheduling problem. The butterfly
graph (in black) represents an abstraction of the data flow of the
fast Fourier transform algorithm on 8 inputs. (In practice, the graph
is more complicated because data are complex, and the real and
imaginary part interact in nontrivial ways.) The boxes denote two
different execution orders that are explained in the text.

 Ò>Y
\ , which simplifies the
Figure 10: Code for the function Ya!Z,
product of two expressions. The comments (delimited with bÝò ò>c )
briefly discuss the various simplifications. Even if it operates on a
dag, this is exactly the code one would write to simplify a tree.

written to some temporary memory location, and the same
process is repeated for the four inputs in the bottom part of
the dag. Finally, the schedule executes the rightmost column
of the dag. In general, the algorithm proceeds by partitioning
the dag into blocks that have d input nodes and j›œ£¤'˜Pd

depth. Consequently, there are j›œ1V£o¤+˜ˆ1gf›'dJ£o¤+˜Nd
 such
blocks, and each block can be executed with ¡¢›'d transfers

between registers and memory. The total number of transfers
is thus at most ¡¢›œ1V£¤'˜E1gfŽ£o¤+˜Pd , and the algorithm matches

the lower bound.
Unfortunately, Aggarwal and Vitter’s algorithm depends
on d , and a schedule for a given value of d does not work
well for other values. The aim of the FFTW generator is
to produce portable code, however, and the generator cannot
make any assumption about d . It is perhaps surprising that
a schedule exists that matches the asymptotic lower bound
for all values of d . In other words, a single sequential order
of execution of an FFT dag exists that, for all d , requires
¡¢›œ1V£¤'˜¦1gf£o¤+˜Pd  register spills on a machine with d registers. We say that such a schedule is cache-oblivious.10
We now show that the Cooley-Tukey FFT becomes cacheoblivious when the factors of 1 are chosen appropriately (as
in [VS94a, VS94b]). We first formulate a recursive algorithm, which is easier to understand and analyze than a dag.
Then, we examine how the computation dag of the algorithm
should be scheduled in order to mimic the register/cache behavior of the cache-oblivious algorithm.
Consider the Cooley-Tukey algorithm applied to a trans-

maximization of pipeline usage is left to the C compiler.
Figure 11 illustrates the scheduling problem. Suppose a
processor has 4 registers, and consider a “column major” execution order that first executes all nodes in the diagonallystriped box (say, top-down), and then proceeds to the next
column of nodes. Since there are 8 values to propagate from
column to column, and the machine has 4 registers, at least
four registers must be spilled if this strategy is adopted. A
different strategy would be to execute all operations in the
gray box before executing any other node. These operations
can be performed fully within registers once the input nodes
have been loaded. It is clear that different schedules lead to
different behaviors with respect to register spills.
A lower bound on the number of register spills incurred by
any execution of the FFT graph was first proved by Hong and
Kung [HK81] in the context of the so-called “red-blue pebbling game”. Paraphrased in compiler terminology, Theorem
2.1 from [HK81] states that the execution of the FFT graph
y
of size 1B6È4 5 on a machine with d registers (where d
1 )
requires at least eV›œ1V£¤'˜ˆ1gf£o¤+˜Pd register spills.9 Aggarwal

and Vitter [AV88] generalize this result to disk I/O, where a
single I/O operation can transfer a block of elements. In addition, Aggarwal and Vitter give a schedule that matches the
lower bound. Their schedule is constructed as in the example
that follows. With reference to Figure 11, assume again that
the machine has d96 Ÿ registers. The schedule loads the
four topmost input nodes of the dag, and then executes all
nodes in the gray box, which can be done completely using
the 4 registers. Then, the four outputs of the gray box are

10 We say “cache-” and not “register-oblivious” since this notion first
arose from the analysis of the caching behavior of Cilk [FLR98] programs
using shared memory. Work is still in progress to understand and define
cache-obliviousness formally, and this concept does not yet appear in the
literature. Simple divide-and-conquer cache-oblivious algorithms for matrix multiplication and LU decomposition are described in [BFJ k 96].

9 The same result holds for any two-level memory, such as L1 cache vs.
L2, or physical memory vs. disk.

9

form of size 1J6Z4'5 . Assume for simplicity that # is itself a
power of two, although the result holds for any positive integer # . At every stage of the recursion, we have a choice of
the factors 1 f and 1 of 1 . Choose 1 f 6 1
s
s 6ml 1 . The
algorithm computes l 1 transforms of size l 1 , followed by
¡¢›œ1  multiplications by the twiddle factors, followed by l 1
more transforms of size l 1 . When 1?zmj›'d , the trans
form can be computed fully within registers. Thus, the number n ݴ1 of transfers between memory and registers when

computing a transform of size 1 satisfies this recurrence.
n

›¥1  6;o

4Hl 1pn ›ql 1  ‡ ¡¢›œ1 
j ›ud



7 Pragmatic aspects of #%$'&%(,(,*
This section discusses briefly the running time and the memory requirements of #$'&%()(,* , and also some problems that
arise in the interaction of the #%$'&%()(+* scheduler with C compilers.
The FFTW codelet generator is not optimized for speed,
since it is intended to be run only once. Indeed, users of
FFTW can download a distribution of generated C code and
never run #%$'&%(,(,* at all. Nevertheless, the resources needed
by #%$-&%()(,* are quite modest. Generation of C code for a
transform of size 64 (the biggest used in FFTW) takes about
75 seconds on a 200MHz Pentium Pro running Linux 2.2 and
the native-code compiler of Objective Caml 2.01. #$'&%()(,*
needs less than 3 MB of memory to complete the generation.
The resulting codelet contains 912 additions, 248 multiplications. On the same machine, the whole FFTW system can be
regenerated in about 15 minutes. The system contains about
55,000 lines of code in 120 files, consisting of various kinds
of codelets for forward, backward, real to complex, and complex to real transforms. The sizes of these transforms in the
standard FFTW distribution include all integers up to 16 and
all powers of two up to 64.
A few FFTW users needed fast hard-coded transforms of
uncommon sizes (such as 19 and 23), and they were able
to run the generator to produce a system tailored to their
needs. The biggest program generated so far was for a complex transform of size 101, which required slightly less than
two hours of CPU time on the Pentium Pro machine, and
about 10 MB of memory. Again, a user had a special need
for such a transform, which would be formidable to code by
hand. In order to achieve this running time, I was forced to
replace a linked-list implementation of associative tables by
hashing, and to avoid generating “obvious” common subexpressions more than once when the dag is created. The naive
generator was somewhat more elegant, but had not produced
an answer after three days.
The long sequences of straight-line code produced by
#%$-&%()(,* can push C compilers (in particular, register allocators) to their limits. The combined effect of #%$'&()(,* and of
the C compiler can lead to performance problems. The following discussion presents two particular cases that I found
particularly surprising, and is not intended to blame any particular compiler or vendor. Ë
The optimizer of the $+#ê âå¾ å¾îå compiler performs an
instruction scheduling pass, followed by register allocation,
followed by another instruction scheduling pass. On some
architectures, including the SPARC and PowerPC procesË
sors, $+#ê employs the so-called “Haifa scheduler”, which
Ë
usually produces better code than the normal $+#ê /#ê)ê
scheduler. The first pass of the Haifa scheduler, however,
has the unfortunate effect of destroying #%$'&()(,* ’s
Ë schedule
(computed as explained in Section 6). In $+#ê , the first
instruction scheduling pass can be disabled with the option

when 1srtj›ud
Pv
otherwise ‚

The recurrence has solution n ݴ1
 6 ¡¢›œ1V£¤'˜ˆ1gfŽ£¤'˜wd  ,
which matches the lower bound.
We now reexamine the operation of the cache-oblivious
FFT algorithm in terms of the FFT dag as the one in Figure 11. Partitioning a problem of size 1 into l 1 problems of
size l 1 is equivalent to cutting the dag with a vertical line
that partitions the dag into two halves of (roughly) equal size.
Every node in the first half is executed before any node in the
second half. Each half consists of j›'l 1 connected compo
nents, which are scheduled recursively in the same way.
The #%$'&()(,* scheduler uses this recursive partitioning
technique for transforms of all sizes (not just powers of 2).
The scheduler cuts the dag roughly into two halves. “Half
a dag” is not well defined, however, except for the power
of 2 case, and therefore the #%$'&()(,* scheduler uses a simple
heuristic (described below) to compute the two halves for the
general case. The cut induces a set of connected components
that are scheduled recursively. The scheduler guarantees that
all components in the first half of the dag (the one containing
the inputs) are executed before the second half is scheduled.
For the special case 1ç674'5 , because of the previous analysis, we know that this schedule of the dag allows the register
allocator of the C compiler to minimize the number of register spills (up to some constant factor). Little is known about
the optimality of this scheduling strategy for general 1 , for
which neither the lower-bound nor the upper-bound analyses hold.
Finally, we discuss the heuristic used to cut the dag into
two halves. The heuristic consists of “burning the candle at
both ends”. Initially, the scheduler colors the input nodes
red, the output nodes blue, and all other nodes black. After
this initial step, the scheduler alternates between a red and a
blue coloring phase. In a red phase, any node whose predecessors are all red becomes red. In a blue phase, all nodes
whose successors are blue are colored blue. This alternation
continues while black nodes exist. When coloring is done,
red nodes form the first “half” of the dag, and blue nodes
the second. When 1 is a power of two, the FFT dag has a
regular structure like the one shown in Figure 11, and this
process has the effect of cutting the dag in the middle with a
vertical line, yielding the desired optimal behavior.
10

Ø
Å%¯°Ú´È(¯,¯ %Å ¯°t´Ù

Ø

Å%¯°t´M()¯)¯ Å%¯«°Ú´Ù
x
»
´¯ » ±Ä,$ƒÂ ×
´¯ ±Ä,$€±E×
¼
¾)¾ Ä«°t($+* ° ¼ $ƒ¯,(MÂ
¾)¾ Ä«°t($+* ° $ƒ¯,(·±
y

structure is independent of the input. Even with this restriction, the field of applicability of #%$'&%(,(,* is potentially huge.
For example, signal processing FIR and IIR filters fall into
this category, as well as other kinds of transforms used in
image processing (for example, the discrete cosine transform
used in JPEG). I am confident that the techniques described
in this paper will prove valuable in this sort of application.
Recently, I modified #%$'&%()(+* to generate crystallographic
Fourier transforms [ACT90]. In this particular application,
the input consists of 2D or 3D data with certain symmetries.
For example, the input data set might be invariant with respect to rotations of 60 degrees, and it is desirable to have
a special-purpose FFT algorithm that does not execute redundant computations. Preliminary investigation shows that
#%$-&%()(,* is able to exploit most symmetries. I am currently
working on this problem.
In its present form, #%$'&()(,* is somewhat unsatisfactory
because it intermixes programming and metaprogramming.
At the programming level, one specifies a DFT algorithm,
as in Figure 4. At the metaprogramming level, one specifies
how the program should be simplified and scheduled. In the
current implementation, the two levels are confused together
in a single binary program. It would be nice to have a clean
separation of these two levels, but I currently do not know
how to do it.

x

x
y

»
´)¯ ±Ä,$ƒÂ× ¼
¾,¾BÄ«°t($+* ° $·¯,(‰Âß¾)¾

¾)¾
)¾ ¾
x
y

»
´)¯ ±Ä,$€±¦× ¼
¾,¾BÄ«°t($+* ° $·¯,(ƒ±

¾)¾

y

Figure 12: Two possible declarations of local variables in C. On
the left side, variables are declared in the topmost lexical scope. On
the right side, variables are declared in a private lexical scope that
encompasses the lifetime of the variable.

Ë
»
Ë Ë
â+(+&%¯â ê `%$,´ Ä,$)â«°& & , and on a 167 MHz UltraSPARC I,

the compiled code is between 50% and 100% faster and
about half the size when this optionË is used. Inspection of
the assembly code produced by $+#ê reveals that the difference consists entirely of register spills and reloads.
Digital’s C compiler for Alpha (DEC C V5.6-071 on Digital UNIX V4.0 (Rev. 878)) seems to be particularly sensitive to the way local variables are declared. For example, Figure 12 illustrates two ways to declare temporary
variables in a C program. Let’s call them the “left” and
the “right” style. #%$'&%(,(,* can be programmed to produce
code in either way, and for most compilers I have tried
there is no appreciable performance difference between the
two styles. Digital’s C compiler, however, appears to produce better code with the right style (the right side of Figure 12). For a transform of size 64,Ë for example,
and
Ë
Ë com-Ë
piler flags â-&«$!_» êƒâ>_ Ë âz‰âË ,Â'& ° ÂÄ«°-Â
â,Â'& ° Â+½)#
â+(+¶ ½%$)¯+½´)$+½‰â'* &%$*`«¯ *Íâ *)´å , a 467MHz Alpha
achieves about 450 MFLOPS with the left style, and 600
MFLOPS with the right style. (Different sizes lead to similar results.) I could not determine the exact source of this
difference.

9 Acknowledgments
I am grateful to Compaq for awarding me the Digital Equipment Corporation fellowship. Thanks to Compaq, Intel, and
Sun Microsystems for donations of hardware that was used
to develop #%$'&()(,* . Prof. Charles E. Leiserson of MIT provided continuous support and encouragement. FFTW would
not exist without him. Charles also proposed the name
“codelets” for the basic FFT blocks. Thanks to Steven G.
Johnson for sharing with me the development of FFTW.
Thanks to Steven and the anonymous referees for helpful
comments on this paper. The notion of cache obliviousness
arose in discussions with Charles Leiserson, Harald Prokop,
Sridhar “Cheema” Ramachandran, and Keith Randall.

8 Conclusion
References

In my opinion, the main contribution of this paper is to
present a real-world application of domain-specific compilers and of advanced programming techniques, such as monads. In this respect, the FFTW experience has been very
successful: the current release FFTW-2.0.1 is being downloaded by more than 100 people every week, and a few users
have been motivated to learn ML after their experience with
FFTW. In the rest of this concluding section, I offer some
ideas about future work and possible developments of the
FFTW system.
The current #%$'&%(,(,* program is somewhat specialized to
computing linear functions, using algorithms whose control

[ACT90] Myoung An, James W. Cooley, and Richard Tolimieri.
Factorization method for crystallographic Fourier transforms. Advances in Applied Mathematics, 11:358–371,
1990.
[ASU86] Alfred V. Aho, Ravi Sethi, and Jeffrey D. Ullman.
Compilers, principles, techniques, and tools. AddisonWesley, March 1986.
[AV88]
Alok Aggarwal and Jeffrey Scott Vitter. The input/output complexity of sorting and related problems. Communications of the ACM, 31(9):1116–1127,
September 1988.

11

[BFJ k 96] Robert D. Blumofe, Matteo Frigo, Chrisopher F. Joerg,
Charles E. Leiserson, and Keith H. Randall. An analysis of dag-consistent distributed shared-memory algorithms. In Proceedings of the Eighth Annual ACM
Symposium on Parallel Algorithms and Architectures
(SPAA), pages 297–308, Padua, Italy, June 1996.
[CO75] R. E. Crochiere and A. V. Oppenheim. Analysis of linear
digital networks. Proceedings of the IEEE, 63:581–595,
April 1975.
[CT65]
J. W. Cooley and J. W. Tukey. An algorithm for the machine computation of the complex Fourier series. Mathematics of Computation, 19:297–301, April 1965.
[DV90] P. Duhamel and M. Vetterli. Fast Fourier transforms: a
tutorial review and a state of the art. Signal Processing,
19:259–299, April 1990.
[FJ]
Matteo Frigo and Steven G. Johnson. The FFTW web
 ôp}q~~a|{ ÒÚÐX|€‚qW>ƒY„ Za‚ ÒÑ|†!~!‡ˆ .
page. {>|
[FJ97]
Matteo Frigo and Steven G. Johnson. The fastest Fourier
transform in the West. Technical Report MIT-LCS-TR728, MIT Lab for Computer Science, September 1997.
The description of the codelet generator given in this
report is no longer current.
[FJ98]
Matteo Frigo and Steven G. Johnson. FFTW: An adaptive software architecture for the FFT. In Proceedings of the IEEE International Conference on Acoustics,
Speech, and Signal Processing, volume 3, pages 1381–
1384, Seattle, WA, May 1998.
[FLR98] Matteo Frigo, Charles E. Leiserson, and Keith H.
Randall. The implementation of the Cilk-5 multithreaded language. In Proceedings of the ACM SIGPLAN ’98 Conference on Programming Language Design and Implementation (PLDI), pages 212–223, Montreal, Canada, June 1998. ACM.
[GHSJ96] S. K. S. Gupta, C.-H. Huang, P. Sadayappan, and
R. W. Johnson. A framework for generating distributedmemory parallel programs for block recursive algorithms. Journal of Parallel and Distributed Computing,
34(2):137–153, 1 May 1996.
[HK81] Jia-Wei Hong and H. T. Kung. I/O complexity: the
red-blue pebbling game. In Proceedings of the Thirteenth Annual ACM Symposium on Theory of Computing, pages 326–333, Milwaukee, 1981.
[HV92] P. H. Hartel and W. G. Vree. Arrays in a lazy functional
language—a case study: the fast Fourier transform. In
G. Hains and L. M. R. Mullin, editors, Arrays, functional languages, and parallel systems (ATABLE), pages
52–66, June 1992.
[JB83]
H. W. Johnson and C. S. Burrus. The design of optimal
DFT algorithms using dynamic programming. IEEE
Transactions on Acoustics, Speech and Signal Processing, 31:378–387, April 1983.
[Knu98] Donald E. Knuth. The Art of Computer Programming, volume 2 (Seminumerical Algorithms). AddisonWesley, 3rd edition, 1998.
[Kul95] Joanna L. Kulik. Implementing compiler optimizations
using parallel graph reduction. Master’s thesis, Massachussets Institute of Technology, February 1995.

[Ler98]

Xavier Leroy. The Objective Caml system release 2.00.
Institut National de Recherche en Informatique at Automatique (INRIA), August 1998.
[Mar76] J. A. Maruhn. FOURGEN: a fast Fourier transform program generator. Computer Physics Communications,
12:147–162, 1976.
[Muc97] Steven S. Muchnick. Advanced Compiler Design Implementation. Morgan Kaufmann, 1997.
[OS89]
A. V. Oppenheim and R. W. Schafer. Discrete-time Signal Processing. Prentice-Hall, Englewood Cliffs, NJ
07632, 1989.
[Par92]
Will Partain. The Ï'ÐX!Z
‰ benchmark suite of Haskell
programs. In J. Launchbury and P. M. Sansom, editors, Functional Programming, Workshops in Computing, pages 195–202. Springer Verlag, 1992.
[PT87]
F. Perez and T. Takaoka. A prime factor FFT algorithm implementation using a program generation technique. IEEE Transactions on Acoustics, Speech and Signal Processing, 35(8):1221–1223, August 1987.
[Rad68] C. M. Rader. Discrete Fourier transforms when the
number of data samples is prime. Proc. of the IEEE,
56:1107–1108, June 1968.
[SB96]
I. Selesnick and C. S. Burrus. Automatic generation of
prime length FFT programs. IEEE Transactions on Signal Processing, pages 14–24, January 1996.
[SJHB87] H. V. Sorensen, D. L. Jones, M. T. Heideman, and C. S.
Burrus. Real-valued fast Fourier transform algorithms.
IEEE Transactions on Acoustics, Speech, and Signal
Processing, ASSP-35(6):849–863, June 1987.
[TAL97] Richard Tolimieri, Myoung An, and Chao Lu. Algorithms for Discrete Fourier Transform and Convolution.
Springer Verlag, 1997.
[Vel95] Todd Veldhuizen. Using C++ template metaprograms.
C++ Report, 7(4):36–43, May 1995. Reprinted in C++
Gems, ed. Stanley Lippman.
[VS94a] J. S. Vitter and E. A. M. Shriver. Optimal algorithms for
parallel memory I: Two-level memories. Algorithmica,
12(2–3):110–147, 1994. double special issue on LargeScale Memories.
[VS94b] J. S. Vitter and E. A. M. Shriver. Optimal algorithms for
parallel memory II: Hierarchical multilevel memories.
Algorithmica, 12(2–3):148–169, 1994. double special
issue on Large-Scale Memories.
[Wad97] Philip Wadler. How to declare an imperative. ACM
Computing Surveys, 29(3):240–263, September 1997.
[Win78] S. Winograd. On computing the discrete Fourier transform. Mathematics of Computation, 32(1):175–199,
January 1978.

12