FFTs in External or Hierarchical Memory
David H. Bailey
December 30, 1989
Ref: Journal of Supercomputing, vol. 4, no. 1 (March 1990), p. 23{35
Abstract

Conventional algorithms for computing large one-dimensional fast Fourier transforms(FFTs), even those algorithms recently developed for vector and parallel computers, arelargely unsuitable for systems with external or hierarchical memory. The principal reason
for this is the fact that most FFT algorithms require at least m complete passes throughthe data set to compute a 2m -point FFT.
This paper describes some advanced techniques for computing an ordered FFT on acomputer with external or hierarchical memory. These algorithms (1) require as few as twopasses through the external data set, (2) employ strictly unit stride, long vector transfers
between main memory and external storage, (3) require only a modest amount of scratch
space in main memory, and (4) are well suited for vector and parallel computation.
Performance gures are included for implementations of some of these algorithms on
Cray supercomputers. Of interest is the fact that a main memory version outperforms thecurrent Cray library FFT routines on the Cray-2, the Cray X-MP, and the Cray Y-MP
systems. Using all eight processors on the Cray Y-MP, this main memory routine runs atnearly two giga ops.

The author is with the Numerical Aerodynamic Simulation (NAS) Systems Division at
NASA Ames Research Center, Mo ett Field, CA 94035.
1

Introduction

The development of numerous advanced architecture computers has posed a considerable challenge to computer scientists. Many numeric algorithms that were completely
satisfactory for traditional serial computers are unsatisfactory for these advanced systems.
This phenomenon is particularly pronounced in the case of algorithms for evaluating one
dimensional fast Fourier transforms.
One reason for this diculty is the fact that many modern computers, particularly
those with interleaved main memories, do very poorly with data that is accessed with amemory stride that is a large power of two. By far the most popular sizes of data to betransformed using FFTs are powers of two, and traditional implementations of FFTs forsuch data sets involve heavy use of power of two memory strides. Fortunately, it is possible
to devise alternative FFT algorithms that do not rely on power of two strides. Indeed,some FFT algorithms can be performed using exclusively unit stride data access in innercomputational loops [3], [5], [6], [9], [10]. Even for systems with external or hierarchicalmemory systems, these unit stride algorithms are a de nite improvement over conventionalalgorithms, since unit strides improve the locality of accesses to and from external memory.
However, many FFT algorithms, both traditional and modern, still require roughly mpasses through the data set to compute a 2m -point FFT. The number of required passes
can be signi cantly reduced by using radix-4 or radix-8 variations of these algorithms,
but the number of passes remains proportional to m. Since such external data access is
usually a crucial bottleneck in such computations, it would be highly desirable to reduce
this number to a bare minimum.

The Basic \Four Step" FFT Algorithm

There is one algorithm in the FFT literature that is quite e ective in reducing thenumber of passes through the dataset. Recently variants of this algorithm were featured
in papers by Agarwal and Cooley [1, p. 150], Ashworth and Lyne, [4, p. 219], andSwarztrauber [10, p. 202 - 203]. Swarztrauber used this technique as a starting point for avery ecient hypercube FFT, and both [1] and [10] noted the suitability of this algorithm
for systems with nonlocal memory systems, including hierarchical and distributed memory
designs. However, as it turns out, this algorithm was actually rst presented over twenty
years ago in a paper by Gentleman and Sande [8, p. 569]. This early paper even described
the application of this algorithm to a system with hierarchical memory. Unfortunately,
this algorithm appears to have been largely forgotten in the interim, as a number of more
recent papers have suggested much less ecient methods.
This algorithm, which shall hereafter be referred to as the \four step" FFT algorithm,
can be stated very succinctly. Let n = n n be the size of the transform. Note that ndoes not necessarily need to be a power of two. On many systems, the implementationof this algorithm is most ecient when n and n are as close as possible to pn. In thefollowing and hereafter, matrices will be assumed to be stored in memory columnwise as
in the Fortran language. The FFT of n complex input data values can then be obtainedby performing the following four steps:
1

2

1

2

2

1. Perform n simultaneous n -point FFTs on the input data considered as a n  ncomplex matrix.
2. Multiply the resulting data, considered as a n  n matrix Ajk , by e ijk=n. The 
sign is the sign of the transform.
3. Transpose the resulting n  n complex matrix into a n  n matrix.
4. Perform n simultaneous n -point FFTs on the resulting n  n matrix.
Several important features of this algorithm should be noted: rst of all, note that bothof the simultaneous FFT steps can be performed using exclusively unit stride data access,which is optimal on virtually any computer system. Secondly, this algorithm produces anordered transform (provided the simultaneous FFTs are ordered) | it is not necessary
to perform a bit reversal permutation, which is inecient on many advanced computer
systems. Finally, note that only three passes through the external data set are required toperform this algorithm | the second step can be performed on a block of data after the
rst step, before it is returned to memory. This bounded number of passes is in accordancewith the I/O complexity results in [2].
1

2

1

1

1

2

2

2

2

2

1

2

1

2

1

Main Memory Performance Results using the Four Step FFT

Depending on implementation, the four step FFT algorithm may actually require aslightly larger number of oating-point arithmetic operations than conventional FFT algorithms. In spite of this slight handicap, it is remarkably ecient even for a single processor
vector computer transforming data in main memory. As can be seen in tables 1 and 2,
a straightforward implementation of this scheme is up to 10% faster than Cray's libraryroutine on the Cray-2 and up to 20% faster than Cray's library routine on the Cray Y-MP.
The percentage results on the Cray X-MP are very close to those on the Cray Y-MP, whichis to be expected since the CPU and memory designs of the X-MP and Y-MP systems arevery similar, di ering mainly in speed of operation. For these tests, the four step FFT
algorithm was implemented using a simple Fortran program; assembly code was employed
only within the Cray library simultaneous FFT routine (CFFTMLT), which is called bythis Fortran program to perform steps 1 and 4. The transpose step (step 3) was performed
without power of two strides by employing a diagonal technique, as mentioned in [6, p.
85]. The Cray-2 library 1-D FFT routine (CFFT2) used in table 1 is an assembly-coded
implementation of an algorithm described by the author in a previous paper [6]. The Cray
Y-MP library 1-D FFT routine (CFFT2) used in table 2 is essentially the same routinethat has been available for some time on the Cray X-MP systems.
The CPU times shown in both tables 1 and 2 are for forward 2m -point FFTs followed
by inverse FFTs, averaged over ten trials, in seconds. All mega ops performance guresin these tables are computed based on 10m2m oating-point operations, even though thefour step routine may perform slightly more than this gure. These tests were run in a
typical daytime environment, and so the results re ect a normal amount of memory bankcontention. The computers used for these tests belong to the Numerical Aerodynamic
3

Simulation (NAS) Systems Division at NASA Ames Research Center. This particular
Cray-2 system has a clock period of 4.1 nanoseconds (ns), and has 268 million words of 80
ns DRAM main memory. The Cray Y-MP system used for these tests has a clock periodof 6.3 ns and 33 million words of bipolar main memory. This Y-MP system was the rst
Y-MP delivered by Cray. Newer Y-MP systems have a faster clock (6 ns), and thus theseresults would be correspondingly better on the newer systems.
The results listed in tables 1 and 2 are single processor results | no attempt was madeto employ more than one processor. However, with the new \autotasking" feature now
available on Cray systems, it is possible to study the performance of a program using all
available processors, with only a minimum of changes to the source code. When autotasking
was invoked on the Fortran program mentioned above, performance levels very nearly eighttimes the single processor levels were achieved on the eight processor Y-MP. These results
are shown in table 3. This very high speedup underscores the suitability for the four stepFFT algorithm for parallel processing.

FFTs on Data in External or Hierarchical Memory Systems

The Cray-2 is noted for its very large main memory. Most Cray-2 systems include 268
million 64 bit words of main memory, although recently Cray has shipped a 536 million word
system. However, the performance of the Cray-2 on many codes in a normal production
environment is not outstanding, due to severe memory bank contention, a direct resultof the relatively slow operation speed of DRAM memory chips. Most Cray X-MP and
Y-MP systems utilize a faster technology (bipolar) in main memory, so that memory bankcontention is very much reduced. However, bipolar memory chips are not available innearly the density of equivalent generation DRAM chips, and so as a result the largest
main memory currently available for Y-MP systems is 33 million 64 bit words. Y-MP
systems typically have eight CPUs, so this means an average of only four million words
per processor. Systems that support interactive as well as batch users must be even more
restrictive in the amount of main memory that can be allocated to a single job.
As a result, users of the Cray X-MP and Y-MP systems who wish to perform large one
dimensional FFTs are led to consider utilizing the solid state disk (SSD) available on thesesystems. SSD systems with a capacity of up to 536 million words are now available on
the Y-MP. Users of the ETA-10 or the IBM 3090/VF systems have an analogous choice inutilizing the virtual memory system, which is a large semiconductor memory similar to the
Cray SSD, but which does not require explicit programmer input/output commands. Userson other systems can even consider utilizing disk drives, although the relative slowness ofsuch devices compared to main memory is a bottleneck even with the best of algorithms.
In addition to minimizing the number of data accesses to an external memory device, anobvious consideration in designing an ecient algorithm for such systems is to minimize
the amount of scratch space required in main memory. Clearly if an external memoryalgorithm requires a substantial scratch array in main memory, then the largest transformsize will again be limited by the available main memory. In addition, it will be assumed in
the following that the amount of external memory is also limited and must be conserved.
4

Size
m
8
9
10
11
12
13
14
15
16
17
18
19
20

Four Step FFTTime MFLOPS
0.000542.5
0.000860.9
0.001376.4
0.0021106.6
0.0036137.8
0.0074143.8
0.0145158.5
0.0300163.9
0.0559187.5
0.1248178.6
0.2426194.5
0.4971200.4
1.0260204.4

Cray Library FFTTime MFLOPS
0.000457.2
0.000681.8
0.0010106.0
0.0021109.4
0.0038130.6
0.0073145.2
0.0138165.7
0.0327150.2
0.0660159.0
0.1260176.8
0.2555184.7
0.5763172.9
1.1863
176.8

Table 1: The Four Step FFT vs. Cray's Library Routine on the Cray-2

Size
m
8
9
10
11
12
13
14
15
16
17
18
19
20

Four Step FFTTime MFLOPS
0.000368.68
0.0005102.240.0008128.160.0013168.910.0024201.900.0049215.880.0103222.570.0212231.670.0443236.780.0935238.390.1976238.810.4117241.960.8587244.23

Cray Library FFTTime MFLOPS
0.0001137.850.0003151.880.0006161.270.0013168.080.0028173.510.0060178.090.0126181.86
0.0265185.240.0557188.260.1167190.970.2439193.440.5090195.701.0635197.20

Table 2: The Four Step FFT vs. Cray's Library Routine on the Cray Y-MP
5

Size Time MFLOPS Speedup
12 0.00079
625.093.096
13 0.00138
771.493.574
14 0.00218 1053.61
4.734
15 0.00376 1308.07
5.646
16 0.00667 1571.73
6.638
17 0.01318 1690.03
7.089
18 0.02566 1838.72
7.700
19 0.05275 1888.53
7.805
20 0.10882 1927.12
7.891
Table 3: Cray Y-MP Performance of the Four Step FFT Using Eight Processors
It will also be assumed for the time being that the nal result in external memory mustbe physically ordered | index schemes or \virtual" orderings of external blocks will not
be allowed.

Reducing the Scratch Space Requirement in the Four Step FFT

As presented above, a straightforward implementation of the four step FFT algorithm
requires scratch space for several di erent purposes. These are as follows: 2n cells for the precomputed root of unity table.
 2n cells of scratch space for the simultaneous FFT steps.
 2n cells of scratch space for the transpose step.
The scratch space requirement for the simultaneous FFT steps can easily be reduced
by noting that the n simultaneous n -point FFTs (i.e. in step 1 of the four step FFT)
may be performed in batches of v rows, where v is the natural vector length of the systembeing used. If the simultaneous FFTs employ an algorithm, such as the Stockham FFT,
which requires a scratch array the same size as the input data array, then only 4vn scratchcells are required. This gure may be reduced by one half if an in-place algorithm can beeciently used for the simultaneous FFTs. Note that if the individual processors do not
rely on vector processing, then only one row need be fetched at a time, and these scratch
space gures drop to only 4n cells and 2n cells, respectively. For step 4 of the four step
FFT, the corresponding scratch space gures may be obtained by replacing n by n in
the above discussion.
However, the scratch space requirements for the simultaneous FFT steps in reality are
dependent more on the block size b of an ecient input/output (I/O) transfer between
main and external memory. In other words, if the natural I/O block length is 128, then
128 rows of the n  n complex matrix should be fetched into main memory, or else the
1

2

2

2

2

2

1

2

6

1

I/O operations will be highly inecient. Thus it follows that a main memory scratch spaceof size 2bn + 2vn is needed for the rst step of the four step algorithm. In an similarmanner, the last step of the four step FFT requires 2bn + 2vn scratch cells. The second
term of each of these expressions may be omitted if an in-place algorithm can be eciently
used for simultaneous FFTs in main memory.
The scratch space for the two FFT steps could be reduced to virtually zero if an FFTalgorithm somewhat more complicated than the four step FFT were used. This algorithm
is as follows:
2

2

1

1

1. Transpose the input data set, considered as a n  n complex matrix, into a n  nmatrix.
2. Perform n individual n -point one dimensional FFTs on the resulting n  n matrix.
3. Multiply the resulting n  n complex matrix Aij by e ijk=n .
4. Transpose the resulting n  n matrix into a n  n matrix.
5. Perform n individual n -point one dimensional FFTs on the resulting n  n matrix.
6. Transpose the resulting n  n complex matrix into a n  n matrix.
1

1

2

2

2

2

2

1

1

1

2

2

1

1

1

2

1

2

2

1

2

2

1

This algorithm, which could be termed by analogy the \six step" FFT algorithm, is very
well suited for distributed memory systems, as the individual one dimensional FFTs can be
performed in individual processors. Its main memory scratch requirement is only 4n cells
for step 2 and 4n cells for step 5 (per processor). As before, if an in-place FFT algorithm
can be eciently used in main memory, then these gures can be reduced by one half.
However, there are other ways of performing FFTs on systems such as MIMD hypercubes[10], and the six step FFT has the serious disadvantage of requiring an additional two
transpose steps, which typically are the chief bottlenecks on any system with a distributed
or external memory.
Reducing the size of the precalculated root of unity table used in step 2 of the four stepFFT algorithm is somewhat trickier. Nonetheless, it can be reduced in size to virtuallyzero with only a slight increase in overall run time, by using what may be termed the
dynamic block scheme for roots of unity. The full size n  n root of unity table can be
written as U (j; k) = jk where = e i=n. Let B (r; s) denote a block of dimensions a  bwithin the matrix U . Note that
U (j + ra; k + sb) = j ra k sb
= jk jsb rak rasb
= U (j; k)U (j; sb)U (ra; k)U (ra; sb)
2

1

1

2

( +

)( +

)

7

2

Thus an a  b block B (r; s) in the interior of U can be dynamically computed as follows:
B (r; s) = top left block [i:e: B (0; 0)] \spike" from top edge
 \spike" from left edge
 upper left corner element of B (r; s)This scheme can be visualized as in gure 1. It can be implemented with only one additional
complex multiplication in the innermost loop. The storage requirements for those subsetsof U array that must be precomputed are as follows:
Description
Upper left basic block\Spikes" from top edge
\Spikes" from left edge
Block corners

Size
2ab
2an =b2bn =a2n=ab
2

1

If we assume
p that n = n = a = b (which is an optimal choice), then the total
space is only 8 n cells, a suciently small amount that this data can be kept in main
memory. In tests of this scheme on Cray systems, the author merely selected a and b
to be 64, the natural vector length. With this choice, only a few thousand cells of main
memory are required even for multimillion point transforms. Performance tests of FFTs
using this scheme indicates that it adds only about ve percent to the total run time (forlarger transforms), and the accuracy of the dynamically calculated roots is excellent.
1

2

2

2

Transposing Arrays in External Memory

The transpose step (step 3 of the four step FFT) is perhaps the most challenging
to perform eciently on a data set residing in external memory. Before discussing thismatter in detail, it should be recalled that the array to be transposed consists of complex
data. In the following discussion it will be assumed that the real and imaginary parts ofthis data are stored in completely separate memory locations, not interleaved as is theFortran convention. In this way the problem of transposing a complex array reduces totransposing two real arrays. In fact, separate storage of real and imaginary data avoids
a signi cant performance degradation in computing with complex data on a number ofsystems, including the ETA-10 and the Cray-2, since the standard Fortran COMPLEX
data format requires stride two access.
Probably the most ecient algorithm currently known to transpose data in external
storage is due to Fraser [7]. A particularly attractive aspect of this algorithm is that it caneasily be tuned for maximum eciency on a given system. It is easier to exhibit an example
of Fraser's algorithm than to precisely state it. Suppose one wishes to transpose a 2  2matrix, which resides on an external random access dataset, into a 2  2 matrix. Suppose
also that the size of an ecient I/O block is 64 = 2 , that two main memory bu ers of size
512 = 2 are available, and that an external scratch dataset of size 2 is available. Let the
8

7

6

9

15

8

8

7

Figure 1: The Dynamic Block Scheme for Roots of Unity

9notation (0 1 2 ... 12 13 14) denote the binary digit positions in the reverse binary
expansion of an index in the 2 -long input array. Then the steps required to transpose
this array can be compactly presented as follows:
15

E1
M1
M2
E2
M1
M2
E1

0
0
8
8
8
8
8

1 2 3 4 5 6 7 8 9 10 11 12 13 14
1 2 3 4 5 8 9 10 11 12 13 14 6 7
9 10 0 1 2 3 4 5 11 12 13 14 6 7
9 10 0 1 2 11 12 13 14 3 4 5 6 7
9 10 0 1 2 11 12 13 14 3 4 5 6 7
9 10 11 12 13 0 1 2 14 3 4 5 6 7
9 10 11 12 13 14 0 1 2 3 4 5 6 7The notation at the beginning of each line indicates the source of the data in each operation:
E1 denotes external dataset number 1, M2 denotes main memory bu er number 2, etc.
Note that the transfers between external memory and main memory only alter locations
6 through 14, and leave locations 0 through 5 unchanged (i.e., 64-long contiguous blocks
are preserved), and that transfers between two main memory bu ers only alter locations
0 through 8 (i.e. only a ect data within a single 512-long main memory bu er).
The rst step, from external to main, involves fetching contiguous blocks of size 64
from disk with a block stride of four (i.e. fetch the rst 64-long block, skip three blocks,
fetch the fth 64-long block, etc.). The rst step is done in batches of 8 blocks, so that 512
words are fetched to one of the main memory bu ers before proceeding. The second step,which is performed between the two main memory bu ers, is to transpose the resulting
512-long array, considered as a 64  8 matrix, into a 8  64 matrix. In the third step, theeight 64-long blocks in the main memory bu er are stored out to external memory, thistime with a block stride of eight. This completes one pass through the external data set.
In the next pass, eight contiguous 64-long blocks are fetched into main memory, and theresulting 512-long array is transposed in a block fashion that preserves 8-long contiguous
sections. Finally, the resulting 64-long blocks are stored back to external memory, again ina manner that achieves a certain block permutation. The array has now been transposedin just two passes. With an adjustment of the parameters (for example, with a block sizeof 32 and a memory bu er size of 1024), the transposition could be achieved in a single
pass.
Even from the above example, the power and generality of Fraser's technique can beappreciated. Unfortunately, Fraser's algorithm cannot in general be performed in place(i.e. using only one external dataset), unless one relaxes the requirement for a physicallytransposed array (by utilizing pointers to index the external data blocks instead). However,in special cases typical of common FFT sizes, there are other methods that can be done inplace and still produce a physically transposed array.
Consider rst the case where n = n , so that the matrix is square. In that case a
block interchange technique can be used to transpose the array in a single pass, in place.
This can be done by simply considering the external n  n matrix to be decomposed into
square blocks of size b on a side, where b is the block size of an ecient I/O operation.
1

2

1

10

2

The square blocks down the diagonal can be transposed simply by fetching the blocks one
at a time into main memory, transposing them using any ecient main memory scheme,
and storing the resulting matrices back in the same locations. The o -diagonal square
blocks can be fetched in opposing pairs, transposed in main memory, and then stored backin opposite locations. One diculty in applying this scheme is when the main memory
block size b is a power of two (which it almost always is). Transposing matrices whose
dimensions are powers of two in main memory, using the straightforward scheme, results
in severe memory bank con icts on many vector supercomputers. However, such arrays canbe transposed completely without bank con icts by fetching and storing opposite diagonals,
as is described in [6, p. 85]. The main memory scratch space requirement for the entire
scheme is 2b cells.
For the common case of power of two FFTs, it can be assumed that either n = nor else n = 2n . In the second case, it does not appear possible to transpose the arrayin one pass, in place, using only full block I/O transfers. However, such arrays can betransposed in just two passes, in place, using only full block transfers, as follows. First,
consider the n  n external array as two blocks of size n  n , and transpose each ofthese two square blocks in place, as described in the previous paragraph. This completes
the rst pass. Now consider the resulting data array in external memory to be a n  nmatrix. Inspection of an example shows that the columns of the resulting array need to
be de-interleaved | column 2j; 0  j < n needs to be moved to column j , and column
2j +1 needs to be moved to column j + n (here the columns are numbered beginning withzero). This de-interleaving could be done b rows at a time using a main memory scratch
array of size bn , but this task can be done more eciently and without need of substantial
scratch space by moving the columns in permutation cycles. For example, suppose n = 8
and n = 4, so that after the block transpose operations we have a 4  8 matrix. Then the
rst cycle would consist of storing column 1 in main memory, moving column 2 to column
1, column 4 to column 2, and column 1 (from main memory) to column 4. The second
cycle would consist of storing column 3 in main memory, moving column 6 to column 3,
column 5 to column 6, and column 3 (from main memory) to column 5. Columns 0 and
7 do not need to be moved. Note that this column movement procedure requires only 2ncells of main memory scratch space. The dominant scratch space requirement for this caseis thus 2b (for each of the two square block transpositions), the same as the case n = n .
2

1

1

2

2

1

2

1

2

2

2

2

2

1

1

2

2

2

1

Performance Results Using the Minimal Scratch Space FFT

2

The above procedure has been implemented and tested on the Cray Y-MP, using one
processor and the SSD external memory device. The SSD I/O primitives SSREAD and
SSWRITE were called directly from the Fortran program. As before, the Cray librarysimultaneous FFT routine (CFFTMLT) was used in steps one and four of the four stepalgorithm. This routine is not an in-place FFT, so that a scratch array in addition to
the space for the data is required. Since the SSD is a rather limited resource, Fraser's
algorithm was not employed for the transpose steps | the in-place schemes described in
the previous section were employed instead. The block length b for ecient I/O transfers
11

Size
m
16
17
18
19
20
21
22

ScratchSpace
3382506682666683301328426
132868226491302650154

Using Memory
Time MFLOPS
0.0704149.000.1574141.590.2897162.900.6391155.861.2263171.022.7007163.075.2996174.12Using SSD
Time MFLOPS
0.116989.68
0.252988.10
0.3094152.520.6908144.211.3065160.522.8179156.295.6132
164.39

Table 4: Minimal Scratch Space FFT Performance Results
between main memory and SSD (or between main memory and disk) on the Cray Y-MP
system is 512.
Table 4 includes results not only for an actual external memory (SSD) implementationof the above scheme on the Cray Y-MP, but also for a modi ed version of the program
where the Fortran routines handling I/O actually just transfer data to a block of main
memory, insteading of referencing the Cray SSD primitives. With the latter gures one
can actually see how much of the performance degradation is due to the algorithm andhow much is due to ineciencies in the Cray I/O system routines. The total amount of
main memory scratch space for this algorithm, including space for precalculated roots ofunity, is also included in this table.

Performing an FFT with Only Two Passes

The schemes that have been described so far produce a physically ordered FFT on an
external dataset in three or four passes. If one is willing to relax the requirement thatthe nal result be physically ordered, or if one is willing to allow a scratch dataset in
external memory of the same size as the input dataset, then the entire FFT operation canbe performed in only two passes (subject to certain conditions). The author is indebted toPaul N. Swarztrauber for this observation.
As in the four step FFT above, it will be assumed in the following that n = n n andthat b is the block size for ecient I/O operation. Also, all references to matrices will, as
before, assume columnwise storage. For simplicity, it will be assumed for the time beingthat two bu ers of size 2bn cells each are available in main memory, although it will later
be seen that only one bu er this large is necessary. Similarly, it will be assumed for the time
being that a scratch dataset equal in size to the input dataset is available in the external
memory device, although it will be seen later than this scratch dataset is not necessary if
one does not mind using pointers. This algorithm can then be stated as follows.
12

1

1. Consider the data in external memory as a n  n complex matrix. Fetch the datab rows at a time into one of the main memory bu ers. For each batch of b rows,
1

12

2

perform b simultaneous n -point FFTs on the b  n array in main memory, using
the second main memory bu er as a scratch array.
2. Multiply the resulting data in each batch by appropriate roots of unity as in the four
step algorithm.
3. Transpose each of the resulting b  n complex matrices into a n  b matrix, using
the second main memory bu er as a scratch array, and store the resulting data on
the scratch dataset in contiguous order. Store successive batches of data in successive
contiguous sections on the scratch dataset.
4. Consider the resulting data in the scratch dataset as a n  n complex matrix. Fetch
the data b rows at a time into one of the main memory bu ers. For each batch ofb rows, perform b simultaneous n -point FFTs on the b  n array in main memory,using the second main memory bu er as a scratch array, and return the resulting b
rows to the same locations on external storage from which they were fetched.
2

2

2

2

2

1

1

1

As before, this FFT is an ordered transform | no bit reversal transposition is necessary.
The reduction of the number of passes from three to two is accomplished by combining thefour step FFT with Fraser's transposition algorithm.
Let r = max(n ; n ). Then at least one main memory bu er of size 2br is required in
the above to hold b rows of the fetched data. However, the second main memory bu er canbe sharply reduced in size in many cases of interest. The additional scratch requirement
for performing the simultaneous FFTs in steps 1 and 4 can be reduced to only 2vr byperforming the FFTs in batches of v rows, where v is the natural vector length of thesystem. If an in-place algorithm is used for the simultaneous FFTs, then this scratch
requirement can be completely eliminated.
Also, in the most common case of power of two transforms, the additional scratch spaceneeded for performing the main memory transpose in step 3 above can be reduced to only
2r cells by applying techniques similar to those mentioned above for transposing power of
two arrays in external memory. One di erence in this case is that the second dimension n
can be much larger than the rst dimension b. Nonetheless, the basic scheme of transposing
the square sub-blocks in place and then moving columns in permutation cycles can also be
applied for this application.
Main memory space to hold precomputed roots of unity can be reduced from 2n to
only 8r by using the dynamic block method described above. Thus the total main memorystorage requirement for power of two transforms can be reduced to only 2(b + 5)r cells
using this algorithm.
The requirement for a separate scratch dataset in external memory can be eliminatedby utilizing a block indexing scheme. At the end of step 3 above, the blocks of datathen in main memory can be returned to the same set of blocks in external memory fromwhich they were fetched, provided a table is maintained of where they are kept. Actually,a table is not even necessary | the permutation involved here is a simple index digit
1

2

2

13

Size
m
16
17
18
19
20
21
22

ScratchSpace
2071786682666683301328426
132868226491302650154

Using SSD
Time MFLOPS
0.0824127.26
0.1659134.340.2705174.440.5607177.67
1.1960175.352.4966176.405.1872177.89

Table 5: Two Pass FFT Performance Results
permutation. However, the ultimate user of the transformed data would also need to use
the same indexing mechanism to access the data.
This \two pass" FFT algorithm has been implemented on the Cray Y-MP using SSD.
A separate SSD scratch array was used instead of the virtual block scheme just mentioned.
This implementation also employed many of the same procedures discussed above to conserve main memory scratch space. The resulting performance gures are shown in table
5. As expected, these results are even higher than the SSD gures in table 4. In fact, theperformance gures listed in gure 5 are almost as high as those for Cray's main memory
FFT, which are listed in table 2.

Conclusion

The excellent data locality of the four step FFT algorithm and its derivatives clearly
is a signi cant advantage for a number of advanced computer systems. In addition, thefact that most of the computation in these schemes reduces to simultaneous FFTs permits
some rather high performance implementations. It has also been demonstrated that someapparent weaknesses of the basic algorithm, such as its large root of unity and scratch
space requirements and its reliance on array transpositions, can be largely eliminated by
employing some advanced techniques.
The performance gures in tables 4 and 5 show that very large FFTs can be eciently
computed using a Cray Y-MP with SSD. In fact, with 33 million words of main memoryand 268 million words of SSD, it should be possible to perform a FFT as large as 2 =134; 217; 728 complex points, provided the SSD device can hold precisely 2 data elements
and no fewer. Such favorable results might not be possible on other systems with slower I/O
to external memory, but the techniques that have been presented should greatly improve
the performance reduction that otherwise occurs.
Another important limiting factor in performing very large FFTs in external memory,which has not been mentioned yet, is the fact that there is often a signi cant wall clock delay
in performing I/O of any sort, even if the CPU time performance is acceptable. Wall clock
performance is particularly important when one is using almost all of main memory, so that
27

28

14

other jobs cannot be utilizing CPU resources when one's own job is waiting for I/O. Suchwall clock delays can be mitigated by overlapping computation and I/O where possible,
and by performing several I/O operations concurrently, provided the overall system I/O
bandwidth is not a limiting factor. Also, some systems have I/O primitives that allow data
to be accessed in external memory in sequences of contiguous blocks with a constant skip
distance between blocks. This is exactly the situation in all of the algorithms mentioned
above, and thus such an I/O function could be expected to substantially improve the wall
clock performance of this FFT scheme, and perhaps the CPU time performance as well.

References
1. Agarwal, R. C., and Cooley, J. W., \Fourier Transform and Convolution Subroutines
for the IBM 3090 Vector Facility", IBM Journal of Research and Development, vol.
30 (1986), p. 145 - 162.
2. Aggarwal, A., and Vitter, J. S., \The Input/Output Complexity of Sorting andRelated Problems", Communications of the ACM, vol. 31 (1988), p. 1116 - 1127.
3. Armstrong, J., \A Multi-Algorithm Approach to Very High Performance 1D FFTs",
Journal of Supercomputing, vol. 2 (1988), p. 415 - 433.
4. Ashworth, M., and Lyne, A. G., \A Segmented FFT Algorithm for Vector Computers", Parallel Computing, vol. 6 (1988), p. 217 - 224.
5. Bailey, D. H., \A High-Performance Fast Fourier Transform Algorithm for the Cray2", Journal of Supercomputing, vol. 1 (1987), p. 43 - 60.
6. Bailey, D. H., \A High-Performance FFT Algorithm for Vector Supercomputers",
International Journal of Supercomputer Applications vol. 2 (1988), p. 82 - 87.
7. Fraser, D., \Array Permutation by Index-Digit Permutation", Journal of the Association for Computing Machinery, vol. 23 (1976), p. 298 - 309.
8. Gentleman, W. M., and Sande, G., \Fast Fourier Transforms { For Fun and Pro t",
AFIPS Proceedings, vol. 29 (1966), p. 563 - 578.
9. Swarztrauber, P. N., \FFT Algorithms for Vector Computers", Parallel Computing,
1 (1984), p. 45 - 63.
10. Swarztrauber, P. N., \Multiprocessor FFTs", Parallel Computing, vol. 5 (1987), p.
197 - 210.

15