Mechanism for Calculating One-Dimensional Fast Fourier Transforms

In the last decade, many applications optimized for vector computers have moved to microprocessor based computers or MCs. The reason for the move is two-fold: microprocessors are cheap and the peak performance of the fastest microprocessor is converging to the peak performance of the fastest vector processor. However, obtaining a speed close to the peak speed for large applications on MCs can be difficult. Vector computer memory systems are designed to deliver one or more elements to the vector registers per clock period, while MCs are not designed to deliver as much bandwidth from their memory systems. Instead, the MCs rely on data caches to allow fast access to the most recently used data. Applications that have a very high reuse of the data in cache run well. For problems that have low data reuse, the performance is typically soley determined by the ability of the processor to load data from memory. Consequently, the goal for most algorithms is to access the memory system as infrequently as possible.

The fast Fourier transform (FFT) is an efficient way to calculate power-of-two one-dimensional discrete Fourier transforms (DFT). The FFT is described in the article by Cooley, J. W., and Tukey, J. W, entitled "An algorithm for the machine calculation of complex Fourier series", Math. Comp., 19: 297-301 1965, which is hereby incorporated by reference. Different algorithms are needed depending upon whether the datasets fit in cache. For datasets that are cache resident, the algorithms that have been developed for vector computers frequently work well in MC systems. This is because the cache on many microprocessors can deliver more bandwidth than their floating-point functional units can consume. This is similar to the case for vector computers. For example, the Stockham based algorithm, as described in the article by Cochrane et al., entitled "What is the fast Fourier transform?", IEEE Trans. Audio and Electroacoustics, AU-15: 45-55, 1967, which is hereby incorporated by reference, works well in this situation. Other variants of the Stockham algorithm exist, as well as other types of vector computer algorithms exist. These types of techniques work well for datasets that fit in cache.

However, when a problem exceeds cache size, the performance of Stockham based or other vector algorithms decreases dramatically in MCs. Consequently, one goal for large FFTs is to minimize the number of cache misses. Prior art techniques have attempted to reduce cache misses, such as the four-step and six-step approaches as described in the article by Bailey, D, entitled "FFTs in External or Hierarchical Memory", in The J. Supercomputing, 4: 23-35, 1990, which is hereby incorporated by reference. Other variants of these approaches exist. The basic four-step approach is composed of row simultaneous FFTs, a transpose, and a twiddle multiply. The basic six-step approach is composed of column simultaneous FFTs, three transposes, and a twiddle multiply. Although the six-step approach requires three data transposes, each column FFT may fit in the data cache.

One formulation of the six-step approach takes an input array X of size n=k.times.m, a work array Y of size n, and an array of roots of unity (twiddle factors) U of size n and comprises the following steps.

1. transpose X(k, m) to Y(m, k).

2. k simultaneous FFTs of length m using Y.

3. transpose Y(m, k) to X(k, m).

4. twiddle factor multiply U(k,m).times.X(k,m)=Y(k,m).

5. m simultaneous FFTs of length k using Y.

6. transpose Y(k, m) to X(m, k).

Although the individual FFTs may fit in cache, many misses may occur in the other steps, however the number of misses is lower than the Stockham based FFTs. Note that the FFT steps are short, one-dimensional contiguous FFTs. Steps 1 through 6 together are mathematically equivalent to performing one large FFT on the input data in array X of size n.

FIG. 1 graphically depicts the data movement on the six-step algorithm. The input is the vector X 11, which is represented as an array of size n which equals m.times.k. The first step 10 transposes the data in array X 11 into a work array Y 12 also of size n. Thus, the data in position i,j of array X 11, becomes the data of the position j,i of array Y 12. The transposition is performed so that the data will have better cache characteristics, i.e. the data in the cache is more likely to be reused. The second step 13 performs m short, one-dimensional contiguous FFTs, each of length k, on the data in array Y 12. Note that the size of the FFTs is k. Further note that short FFTs are usually about the square root in size of the original problem, and thus, more likely to fit into cache. The third step 14 transposes the data back from array Y 12 to array X 11. The fourth step 15 multiples the X array by the coefficients of array U 16, which has a size n, and stores the results into the work array Y 12. The array U 16 comprises the twiddle factors, which are previously calculated trigometric coefficients that are used for scaling. Note that in steps 4 through 6, the work array Y 12 is considered to be a k.times.m array, instead of a m.times.k array as in steps 1 through 3. This ensures that the step 5 FFTs are contiguous in memory. The fifth step 17 performs k FFTs, each of length m, similar to step 2, on the data stored in array Y. To complete the transform, the sixth step 18 transposes the data back into array X 11 from array Y 12. This implies that in step 6, the array X 11 is considered to be a m.times.k array, instead of a k.times.m array as in steps 1 through 3.

Since X 11, Y 12, and U 16 are all size n, and all will be stored into cache from the memory, then cache misses would only typically occur when n is greater than 1/3 of the cache size. FIG. 2 depicts typical size relationships of cache 21 and memory 22, where the memory 22 is shown to be much larger than the cache 21. Each block of data in the memory is the same size or larger than the cache, thus as additional blocks of data are loaded into cache, the prior data is overwritten. Each time a block of data is loaded into cache, cache misses occur. These cache misses are required as data must be moved from the memory and into cache. FIG. 3 depicts the size relationships between the cache 21 and the memory 22, however, the data blocks X 31 and Y 32 are larger than the cache. Therefore, as each portion of the block X 31 is loaded into memory, cache misses are accrued from the data movement. Note that cache misses occur when data is moved from memory to cache, as well as cache to memory. Typically, when data is stored, the data must first be loaded from memory to cache, and then stored back to the memory, thus accruing double the misses that a load accrues, as two data movements have occurred. For example, in storing data to Y 32, the cache first loads 34 the data into cache, and then stores 33 the data back into memory. Therefore, the mechanism of FIG. 1 will incur many cache misses if the size of n exceeds 1/3 the size of the cache, since the mechanism of FIG. 1 uses three arrays, X, Y, and U, all of size n.

Therefore, there is a need in the art for a FFT mechanism that reduces the number of cache misses that occur when the data block is larger 1/3 of the size of the cache.

These and other objects, features and technical advantages are achieved by a system and method which prior to calculating the FFT, decomposes the input data array in a plurality of squares, and then manipulates the data within the squares and within columns of the input array.

The inventive mechanism modifies the prior art six step approach into a seven step approach. The seven steps, taken as a whole, are mathematically equivalent to performing one large FFT on the input data. The input data X of size n is treated as s squares of size t.times.t, such that n=s.times.t.sup.2. In first step, each of the s squares of size t.times.t is treated separately, and their respective points are swapped across their main diagonals, i.e. the data at point i,j is exchanged for the data at the point j,i, for each square. In the second step, t FFTs, each of length t is calculated for each of the squares, for a total of s.times.t FFTs. The results are stored back into their respective locations. In the third step, the data is transposed in each of the squares as the first step. In fourth step, the s squares are multiplied by the twiddle coefficients of array U, which has a size n=s.times.t.sup.2. In the fifth step 75, t FFTs of length s.times.t are calculated. In steps 4 and 5, the s squares 61 and the array U 78 are viewed as comprising t one dimensional arrays of size s.times.t. The results of each of steps 4 and 5 are stored into a work array y, which has a size of s.times.t, and is small enough to remain in cache. In the sixth step, columns of y data are transposed and stored into columns of array X. In the seventh step, the data is transposed in each of the squares as the first and third steps.

Treating the data as squares allows the results of the transposition steps and the twiddle step, i.e. steps 1, 3-4, and 7, to be stored back onto the same memory location. This eliminates the need for a large work array which has a size equivalent to the input array. Since the large work array is not used, then cache misses occurring from loading and storing the work array do not occur. By using a small work array y in steps 4-6, which remains in cache, additional cache misses are saved.

The inventive mechanism allows for parallel processing of the different steps. Parallel threads can simultaneously process the different data structures. In steps 1-3 and 7, s threads can separately process each square. In steps 4-6, t threads can separately process each column or one dimensional array.

Therefore, it is a technical advantage of the present invention to calculate the FFT of a large input problem with fewer cache misses.

It is another technical advantage of the present invention that the memory requirements for calculating the FFT are reduced.

It is a further technical advantage of the present invention that the various operations of the invention can be parallel processed.

The foregoing has outlined rather broadly the features and technical advantages of the present invention in order that the detailed description of the invention that follows may be better understood. Additional features and advantages of the invention will be described hereinafter which form the subject of the claims of the invention. It should be appreciated by those skilled in the art that the conception and the specific embodiment disclosed may be readily utilized as a basis for modifying or designing other structures for carrying out the same purposes of the present invention. It should also be realized by those skilled in the art that such equivalent constructions do not depart from the spirit and scope of the invention as set forth in the appended claims.

For a more complete understanding of the present invention, and the advantages thereof, reference is now made to the following descriptions taken in conjunction with the accompanying drawings, in which:

FIG. 1 depicts the prior art FFT mechanism;

FIG. 2 depicts relationships of the cache and the memory;

FIG. 3 depicts the relationships of the cache and the memory when a data block is larger than the cache;

FIG. 4 depicts the normalized cache misses for the mechanism of FIG. 1;

FIG. 5 depicts the normalized cache misses for an aspect of the invention;

FIG. 6 depicts the decomposition of the input data into squares;

FIGS. 7A-7C depict the seven steps of the inventive FFT mechanism;

FIG. 8 depicts the normalized cache misses for the inventive mechanism of

FIGS. 7A-7C; and

FIG. 9 is a high level block diagram of a computer system adapted to implement the present invention.

As previously stated, the mechanism of FIG. 1 incurs many cache misses when n exceeds 1/3 of the cache size. The exact number of cache misses depends upon cache line size and other factors that are system dependant. A convenient mechanism to compare the number of misses is to normalize cache misses. Thus, a normalized cache miss is defined as follows:

normalized miss=n.times.bytes per point/cache line size, wherein a point is a point i,j in the input array.

Thus, when a load occurs and the data is not in cache, the cache misses and the data and must be loaded or read from memory, resulting in one normalized miss. When a store occurs and the data is not in cache, there are two normalized misses, one miss to read or load the data from memory into cache, and one miss to write or store the new data back to memory (this is also referred to as write-back). Thus, using 1 miss for loads and 2 misses for stores results in the six-step approach causing 17 misses, as shown in FIG. 4. Steps 1, 3, and 6 require moving data between memory locations, each of these steps causes three misses, one miss from a load of the data from X, and 2 misses from a store of the data to Y (Note that data cannot be stored until it has first been loaded). Steps 2 and 5 require the storing of data into Y, and thus each step incurs two misses. Step 4 incurs four misses, one from loading X, one from loading U, and two for storing into Y.

The first aspect of the invention is that if X can be transposed onto itself, then Y is not needed. This will cause fewer misses than the mechanism of FIG. 1. This will also reduce the memory space requirements by a third, as only two arrays will be stored, i.e. X and U. Moreover, this changes the cache requirements as only the two equal sized arrays, X and U, will need to fit into cache. Note that each transpose would only incur 2 misses instead of three misses, since the store of X has been preceded by a load of X. This can occur when X is a square, i.e. m.times.m. Thus, the transposition involves the exchange of two points, i,j with j,i. FIG. 5 depicts the data movement and cache misses when the Y array is eliminated. Note that step 4 saves a cache miss as the additional load is not required. Thus, by eliminating the Y array, a total of 13 cache misses occur, which is a savings of 4 cache misses over the mechanism of FIG. 1. However, note that the steps 1 through 6 together are mathematically equivalent to performing one large FFT on the input data in array X of size n.

The transposition of a square X or n, assuming a size of s.times.s, only requires swapping diagonals across the main diagonal, shown as follows.

complex X(s, s)

do j=1, s-1

do i=1, j

temp=X(n-j+i, i)

X(n-j+i, i)=X(i, n-j+i)

X(i, n-j+i)=temp

enddo

enddo

Note that the work array Y is not required. This transpose may be blocked for good cache line reuse. The six-step approach of FIG. 1 is performed with X.sub.i, j being overwritten in each step. This will achieve the cache miss savings as shown in FIG. 5.

However, not all inputs X or n are square, many are multiples of powers of 2, 3, and 5. However, all large numbers that are multiples of powers of 2, 3, and 5 can be decomposed into a plurality of squares. Thus, the elimination of the Y array can be generalized for all inputs having such powers. Note that if the powers of 2 ,3, and 5 are small, i.e. 0 or 1, then the resulting n is small, i.e. 1, 2, 3, 5, 6, 10, 15, or 30. These small n will typically fit into the cache, and thus will not incur misses from data movements to and from memory. The FFT can be calculated in a single step. All other n can be written as n=2.sup.2k+p .times.3.sup.21+q .times.5.sup.2m+r, wherein 0.ltoreq.p,q,r.ltoreq.1 and at least one of the powers is >1.

Assume for the remainder of this disclosure that the power of 2 is greater than 1, i.e. 2k+p>1. Note that either the power of 3 or 5 could be greater than 1, in which case that particular power would used to form the squares. If more than one of 2, 3, and 5 has a power greater than one, then the mechanism would select which of 2, 3, and 5 would be used to form squares.

As stated above, non-square rectangles can be decomposed into a plurality of squares. Thus, s squares, with each square being t.times.t can be formed from input n, with s=2.sup.p .times.3.sup.21+q .times.5.sup.2m+r and t=2.sup.k. Note that n=s.times.t.sup.2. This is depicted in FIG. 6 with s squares 61, each sized t.times.t. Note that the decomposition does not involve any data movement, and thus does not incur any cache misses.

FIGS. 7A-7C depict the seven steps of the inventive FFT mechanism. The seven steps together are mathematically equivalent to performing one large FFT on the input data in array X of size n. Note that since there is not a main diagonal as in the n=m .times.m case, points cannot be swapped as before. However, the transposes in steps 1 and 3 do not need to be complete transposes. The only requirement is that the data needed for the individual step 2 FFTs be in the same order as step 2 of the six step approach. Thus, in the first step 71, each of the squares is treated separately, and their respective points are swapped across their main diagonals, i.e. the data at point i,j is exchanged for the data at the point j,i, for each square. Each of the transposed squares is stored over itself. In the second step 72, t FFTs, each of length t is calculated for each of the s squares, for a total of s.times.t FFTs. The results are stored back into their respective locations. In the third step 73, the data is transposed in each of the squares as the first step 71. At the conclusion of the third step the data is back in proper order for the remaining steps.

In fourth step 74, the s squares are multiplied by the coefficients of array U 78, which has a size n=s.times.t.sup.2. The array U comprises twiddle factors, as discussed with regard to FIG. 1. In the fifth step 75, t FFTs of length s.times.t are calculated.

The sixth step 76 and the seventh step 77 together form a complete transposition. The goal is to transpose X(s.times.y.sub.1, y.sub.2) to X(y.sub.2, s.times.y.sub.1), wherein y.sub.1, y.sub.2 equal t. If the data is treated as three dimensional arrays, then the data X(y.sub.1, s, y.sub.2) must be transposed to X(y.sub.2, s, y.sub.1). This transpose may be achieved by:

complex X(y.sub.1, s, y.sub.2);

transpose X((y.sub.1, s), y.sub.2) to X((s, y.sub.1), y.sub.2); and

transpose X((s, y.sub.1), y.sub.2) to X(y.sub.2, (s, y.sub.1)).

The first transpose or step 6, operates on a column of data at a time. Once all the columns have been transposed, the second transpose or step 7 is performed. This part is the same operation as the partial transposes of steps one and three. This completes the transform.

In step 4, note that the s squares 61 and the array U 78 can be viewed as comprising t one dimensional arrays of size s.times.t. Thus, each of the multiplications of one array from the squares and one array of the array U would be stored into a work array y, which has a size of s.times.t. Note that for comparison, the large Y array used in the mechanism of FIG. 1 would have a size of s.times.t.sup.2. Thus, in the fourth step 74, a column from array U 78 is multiplied with a respective column from the array X, i.e. one column from all of the squares. This result is then stored in working array y 79. Note that array y 79 is usually smaller than cache, and should stay within the cache during further processing. In the fifth step 75, a FFT is performed on the data in the array y 79. The FFT has a length of s.times.t. Note that no cache misses will occur, because the data is already loaded into cache. Steps 4 and 5 will be performed t times, with t one dimensional arrays.

In the sixth step 76, columns of data are transposed. At the conclusion of step 5, data is in the y array. Each transposition reads the data from y and writes into a column of array X, such that upon completion of the last one dimensional array, the data is ready for step 7. Since the array y has been used in steps 4 and 5, then the array y is already in cache, then a cache miss does not occur for loading or storing y. However, a cache miss does occur for storing the result into X. Note that there is only one cache miss for the store because the load occurred in step 4. The transposition involves treating the vector y as an array y(y1,s) and transposing it to X((s, y1),y2), where it is stored in X.sub.j. In the seventh step 77, the data is transposed in each of the squares as the first step 71 and the third step 73.

FIG. 8 depicts the cache misses associated with the mechanism of FIGS. 7A-7C. As compared with FIG. 5, note that splitting the sixth step into two steps, has added a cache miss. However, by using the one dimensional array y, the number of cache misses has been reduced in the twiddle multiple step and the subsequent FFT step, for a overall reduction in cache misses by two, for a total of 11 cache misses.

The inventive mechanism of FIGS. 7A-7C also allows for parallelism or parallel processing of the different steps. As shown in FIG. 8, parallel threads can simultaneously process the different data structures. In steps 1-3 and 7, s threads can separately process each square. The processing of each square is independent of the others. In steps 4-6, t threads can separately process each column or one dimensional array. The processing of each column is independent of the others. Note that problems which fit into a cache should not be run in parallel as cache misses could occur and reduce performance.

Parallelism also depends upon the number of processors. For example, suppose there are z processors. The problem should be decomposed into at least z squares. If s and t are defined such that n=s.times.t.times.t, then the granularity of parallelism in steps 1-3 and step 7 is s, and the granularity of parallelism in steps 4-6 is t. For efficient parallelism z.ltoreq.s,t. This is met for large n. It can be shown that for z.ltoreq.16 and 31,000.ltoreq.n there are always s, t such that z.ltoreq.s,t. For example, if zk=16 and n=468,750=2.times.3.times.5.sup.7, the data can be partitioned as 30 squares of size 125. This gives 30 way parallelism in steps 1-3 and step 7, and 125 way parallelism in steps 4-6. To maximize the performance of the serial components of the parallel decomposition, the minimum s subject to the above constraints should be selected.

FIG. 9 illustrates a computer system 90 adapted to use the present invention. Central processing unit (CPU) 91 is coupled to bus 92. In addition, bus 92 is coupled to random access memory (RAM) 93, read only memory (ROM) 94, input/output (I/O) adapter 95, communications adapter 101, user interface adapter 98, and display adapter 99.

RAM 93 and ROM 94 hold user and system data and programs as is well known in the art. I/O adapter 95 connects storage devices, such as hard drive 96, to the computer system. Communications adaption 101 is adapted to couple the computer system to a local, wide-area, or Internet network 102. User interface adapter 98 couples user input devices, such as keyboard 103 and pointing device 97, to the computer system 90. Finally, display adapter 99 is driven by CPU 91 to control the display on display device 100.

CPU 91 may be any general purpose CPU, such as a HP PA-8200. However, the present invention is not restricted by the architecture of CPU 91 as long as CPU 91 supports the inventive operations as described herein. The HP PA-8200 has a 200 megahertz (MHZ) clock period with a peak speed of 800 million floating-point operations per second (MFLOPS). The processor has the ability to execute two memory operations and two floating-point instructions per clock cycle. The processor has a two megabyte (MB) data cache and a two MB instruction cache. The caches have 32 byte cache lines and are direct mapped. The maximum transfer rate from the data cache to the register file is 3200 MB per second (MB/sec). The processor supports multiple data cache misses from memory. Up to ten outstanding data cache misses are supported at one time. The maximum transfer rate from memory to cache is dependent upon the memory latency. For a system with a 500 nanosecond (ns) latency, 10 requests from memory can request a maximum of 640 MB/sec. Since the maximum bandwidth from memory to cache is much less than the cache to functional unit bandwidth, the number of cache misses should be minimized.

An example of a parallel processing system, which would be able to handle multiple threads, is the Exemplar V-Class server, which incorporates 16 of the PA8200 processors into a cache coherent, uniform memory access system. The processors are connected by a non-blocking crossbar to the memory system. The memory system has 32 banks and each bank is 32 bytes wide. The memory system latency is 509 ns and the memory system has been measured delivering data from memory to a PA8200 cache at 520 MB/s.

Although the present invention and its advantages have been described in detail, it should be understood that various changes, substitutions and alterations can be made herein without departing from the spirit and scope of the invention as defined by the appended claims