The Fast Fourier Transform on Hypercube Parallel Computers Clare Yung-lei Chu 87-882 November 1987 Department of Computer Science Cornell University Ithaca, New York 14853-7501 THE FAST FOURIER TRANSFORM ON HYPERCUBE PARALLEL COMPUTERS A Thesis Presented to the Faculty of the Graduate School of Cornell University in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy by Clare Yung-lei Chu January 1988 Qc Clare Yung-lei Chu 1987 ALL RIGHTS RESERVED The Fast Fourier Transform on Hypercube Parallel Computers Clare Yung-lei Chu, Ph.D. Cornell University 1988 The Fast Fourier Transform appears frequently in scientific computing. Therefore it is desirable to implement it efficiently on parallel computers. In this thesis, we investigate several different aspects of parallel Fast Fourier Transform imple- mentation techniques for distributed-memory message.passing systems such as hy- percube multiprocessors. We describe various Fast Fo'rrier Transform algorithms using a matrix notation. An error analysis is presented that considers the effect of different methods used in the computation of the Fo'rrier Transform coefficients as well as accumulated roundoff. New implementations of one and tw?dimensional Fast Fourier Transforms are presented along with comparisons with existing meth- ods. New algorithms for symmetric transforms are also developed and the results show excellent speedup when implemented on the Intel iPSC hypercube. Biographical Sketch The author was born on January 13, 1960 to Dr. Chong-Wei Chu and Mrs. Char- lotte C. Chu in Minneapolis, MN where her father finished up his Ph.D. in Fluid Dynamics and her mother an M.S. in Biostatistics. Her family then moved to Southern California when she was three and she has lived there since, first gradu- ating from Occidental College, Los Angeles with an A.B. degree in 1980 and then receiving a M.A. in Applied Mathematics from the University of California, San Diego, La Jolla. Her majors were sun-bathing and jogging, with a minor in wind- surfing. And her summers were radically fun since she worked at T.R.W. Defense Systems Group, Redondo Beach, as a Programmer/Analyst. The author's career at Cornell began when she decided to "experience" the East Coast. So she reviewed the Eastern universities that offer Applied Mathematics programs and chose Cornell both for its rural, scenic location and the quality of the faculty. Seasons changed, leaves fell, and wintry snow blew in a short spurt of spring, followed by sweltering summer heat. Still she vegged out at the Center for Applied Mathematics, studying matrix computations, partying at the International Living Center, jogging around the Plantations, and working out at the Ithaca Fitness Center. True she had lost her tan, but gained a world of knowledge. After graduation, she will work as a Senior Engineer at Northrop Aircraft in Hawthorne, California; regain her bronze complexion, continue to pursue Ultimate Fitness, and enjoy the Liberal Arts in sunny Southern California. Her East Coast experience was, as they say, quite character-building. Only thing, she was already a Mega-Character before she got here, so one can surmise that these five years have made her a bit of a Hyper-Character. SoCal Rocks! iii In the memory of my father, Dr. Chong-Wei Chu iv Acknowledgements My sincere gratitude to Professor Charles Van Loan, my thesis advisor, for his kindness and understanding during these trying years of dissertation creation. He was instrumental in directing me into this avenue of research and gave me many hours of advice. The personal attention he gave me these past few years means a lot to me. I am also grateful for the generous financial support provided by Professor Van Loan through the following grants, ONR Contract N00014-83-K- 640 and NSF Contract DCR 86-2310. I also wish to thank Professor Lars Wahlbin for his friendly help and interest in my well-being, and Professor Tom Coleman for serving on my committee and making himself available for an occasional chat. Finally I thank Professor Al Schatz for his willingness to substitute for Professor Wahlbin on my B-exam and his interest in my work. The Center for Applied Mathematics went through four Center Directors dur- ing my stay here at Cornell and I would like to thank Professor Phil Holmes, Professor Lars Wahlbin, Professor Mike Todd, and Professor Sid Liebovich for the advice and help they gave me. I also sincerely appreciate Mrs. Dolores Pendell for lending me a helping hand whenever I needed it and for cutting through some of the administrative matters that a graduate student has to face. My parents deserve appreciation for their love and care as well as providing me with an academic attitude, as do my brother Dale, and two sisters Constance and Charleen. I also appreciate the friendship of my former supervisor Mr. Richard Maitlen of T.R.W. for his interest in my career, as well as Professor Stanley "Gil" Williamson of U.C. San Diego for his confidence in my ability. My gratitude goes out to Professor Joan Moschovakis and Professor Benedict Freedman at Occidental College because if it weren't for their encouragement, I'd probably have been a Diplomacy and World Affairs major! I wish to thank all my friends here and in California for all the fun I had these last five years. I also appreciate my cronies here at Applied Mathematics for putting up with my talkative nature. And the fact that my physical condition survived Ithaca is due to the great staff and facilities at the Ithaca Fitness Center. I will also miss the ambience and people at the International Living Center where v I lived during my years at Cornell. My friends have all helped me immensely over the past five years and I would like to thank them: Daniel Medina, Tim Waag, Ann Monzon, Jada Valderrama, Amelia Wu, Irene Duffy, Richard Koo, Sherease Pratt, Jeter Iseley, John Dailey, Cindy Anfinson, Chris Hempel, Glen Swindle, Ernie DeMarie, Alex Espalin, Gabrielle Tayac, and Kym Fraser. And lastly I am grateful to my boyfriend, J??? A. Ayala for his love and concern. vi Table of Contents Biographical Sketch Acknowledgements 2 3 Preliminaries 1.1 Introduction 1.2 Parallel Scientific Computing 1.3 The Hypercube Architecture 1.4 Distributed Multiprocessor Computing Issues . 1.5 Definitions and Matrix Notation 1.6 Permutations 1.7 The Fast Fourier Transform . 1.8 Matrix Description of the FFT 1.9 The Binary Reflected Gray Code Multiplier Computation, Error Analysis 2.1 Introduction 2.2 Previous Work 2.3 Trigonometric Function Generation 2.4 Multipliers for the FFT 2.5 Error Analysis of the FFT 2.6 Practical Multiplier Generation iii v One Dimensional Parallel Fast Fourier Transforms 3.1 Introduction 3.2 One-Dimensional Complex FFT on the Hypercube . . 3.3 Implementation Results for the One-Dimensional FFT . 3.4 The Tw?Track FFT Implementation. . . 3.5 Distributed Multiplier Computation v?1 3 4 6 8 10 20 21 27 30 30 32 34 48 53 66 72 72 73 82 87 96 4 5 6 7 Comparison of Two-Dimensional FFT Methods on the Hyper cube 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 Introduction . . . . Two Dimensional FFT Algorithms Data Structures and Mapping Efficient Data Transposition on the Intel iPSC Comments on Data Ordering, Vectorization, Mixed Radix Imple- mentations, and Storage . . . Implementation and Timings Discussion and Model Conclusion Parallel Sine Transform for the Hypercube 5.1 Introduction 5.2 Symmetries and Properties of the DFT 5.3 Permutations and Communication on Distributed Processors 5.4 Development of the Parallel Sine Transform 5.5 The Sine Transform of a Sequence in Bit-Reversed Order 5.6 Implementation Results on the Intel Hypercube . 5.7 The DFT and IDFT of Real Sequences 5.8 The Cosine and Sine Transform . 5.9 Conclusion and Discussion Fast Symmetric Transforms 6.1 Introduction . . 6.2 Definitions and Matrix Notation 6.3 Description and Development 6.4 Error Analysis 6.5 Implementation on the Hypercube . . 6.6 Applications: Updating Shifted Sequences 6.7 Applications: Fast Poisson Solvers 6.8 Appendix: Proof for Type III and Type IV Symmetric Transforms 6.9 Fortran Code . . Distributed Mixed-Radix FFTs on the Hypercube 7.1 Introduction 7.2 Twiddle Factor Algorithm Splitting of the DFT 7.3 Procedure for Hypercube Implementation viii 103 103 105 108 110 113 117 127 132 134 134 136 140 149 153 155 157 161 167 169 169 171 173 184 186 189 191 192 199 204 204 205 206 References 211 List of Tables 2.1 2.2 2.3 2.4 2.5 3.1 3.2 3.3 3.4 3.5 3.6 3.7 4.1 4.2 4.3 4.4 4.5 4.6 5.1 5.2 5.3 6.1 6.2 6.3 6.4 6.5 6.6 Maximum Error Magnitude = IIcIIoo/u . Root Mean Square Error Magnitude = JICIIrms/U Double Precision Implementations . . Error Bounds for Weight Computational Methods Buneman's On-Line Generation of Trigonometric Functions Four Canonic Radix-2 FFT Algorithms One-Track FFT FORTRAN Code Forward FFT (Natural Implementation) Forward FFT (BRGC Implementation) Two-Track FFT FORTRAN Code Forward FFT (Tw?Track) Natural Forward FFT (Two-Track) BRGC Multiple Transforms and Vector Length. Total Execution Time. . . . . . Computational Time of 2D-FFT. Communication Time Including Blocking Overhead. Time Spent Blocked during Communication . Communication Times without Intervening Computation. Symmetric Sequences. . . Timing Results for the Parallel Sine Transform Timing Results for Conventional Complex Transform. Symmetric Transform Formulae Cooley et at. Sine Transform Errors Sine Transform (Type I) Errors . . . Algorithm Type12 n = 16384 Timings (Forward) Algorithm Type12 n = 16384 Timings (Inverse) . Updating of Shifted Sequences . 49 50 51 65 71 82 84 85 88 95 95 97 115 121 124 125 126 130 137 156 156 175 187 188 190 190 191 ix 7.1 Mixed Radix-2 Timings . x 209 List of Figures 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 5.1 5.2 5.3 5.4 n = 8, d = 2 mapping of data in natural order n = 8, d = 2 mapping of data in bit-reversed order n = 8, P = 4, signal flow graph (GS1, CT1) . Distributed Centleman-Sande butterfly. . . . . Distributed Cooley-Tukey Butterfly Walton's Implementation of the Distributed CT2 Butterfly n = 8, d = 2 mapping of data in natural order, cube in BRGC. Two-Track Forward Transform Data Swapping Pattern. Two-Track Inverse Transform Data Swapping Pattern Tw?Track Implementation of the Distributed GS Butterfly Two-Track Implementation of the Distributed CT Butterfly . . One and Two Dimensional Decimation Schemes . . 16 columns mapped into 4 processors 5-cube = 3-cube 0 2-cube.. 16 columns mapped into 16 processors Block Recursive Transpose. . . Blocks Exchanged by Processor # 2 Block Method 2D-FFT Communications for a 16-node hypercube Vector-Radix Exchanges for Processor p0i00 Permutation Pattern of Two-Track FFT Permutation Pattern of Two-Track En Permutation Pattern of Two-Track T? Permutation Pattern of Two-Track Gn xi 74 75 7(3 80 81 83 87 89 90 91 92 107 108 109 110 111 112 119 122 144 145 146 148 Chapter 1 Preliminaries 1.1 Introduction The parallel computation of numerical problems holds a tremendous amount of promise for speeding up large compute-intensive procedures. The hope is to be able to perform real-time computations on important applications such as image processing, the processing of seismic data, signal processing, and the solutions of partial differential equations. Buried inside a lot of scientific code is the discrete Fourier transform (DFT). The realization of a "fast" method for its computation in the 1960's [Cooley and Tukey (1965)] allows scientists to solve problems that were previously considered beyond the realm of practical computability. Many variations of the fast Fourier transform (FFT) were subsequently developed to cater to the parameters of specific problems and machine architectures. Having reached the limit with which a single processor can perform FFT computations, one is naturally led to consider the implementation of FFTs on parallel computers. This dissertation focuses on issues dealing with the design of new algorithms and procedures for the implementation of the FFT and fast symmetric transforms on distributed memory ensemble architectures, in particular the hypercube intercon- nected parallel machine. The complex FFT is one of the first algorithms to be implemented on the first experimental hypercubes where inefficiencies due to load imbalauce and com- munication overhead were noted. [Fox (1985)] While on the surface, the radix-2 FFT implementation appears straightforward and easy, further research reveals that there are many details with room for improvement. Furthermore the efficient implementation of symmetric transforms is more complicated and had not been 2 extensively studied before. The purpose here is to increase the understanding of parallel Fourier transforms by presenting new methods and techniques for improv- ing their execution rates on parallel message-passing computers. New heuristics for conceptualizing concurrent numerical computations are introduced in the form of a matrix language representing algorithm design, task allocation, data partitioning and inter-processor communication. The first chapter provides the reader with an introduction to the Fast Fourier Transform algorithm interpreted according to matrix factorizations. The use of matrix language is well-suited for the description of numerical algorithms especially for distributed computations. By getting away from the i-i level, replete with mul- tiple summations and backwards running indices, an overall view of the big picture can be obtained where blocks of data are partitioned among distributed processors and operated upon. The use of permutation matrices to describe communication is especially useful and provides a higher-level abstraction for algorithmic descrip- tion and development. Some basic permutations and notational conventions are introduced in Chapter 1. The fast Fourier transform is also conveniently denoted as a sequence of sparse matrix multiplications. Four canonical in-place FFT algo- rithms are presented along with some of their properties. Facility with the matrix notation allows us to transcend the tangled web of i4 indices and derive different, but equivalent, versions of the FFT by rearranging the matrix products. Chapter 2 presents an error analysis for the radix-2 FFT. A different approach is used where the FFT coefficients (or multipliers) are assumed to contain noise due to the method with which they are computed. This is a non-random, determinable error and dependent on the method that is used to calculate these trigonometric values. Several methods are highlighted, with their error properties examined. These errors are then incorporated into the error analysis of the FFT and are shown to contribute multiplicatively to the total FFT error. Chapter 3 describes a new method of implementing one-dimensional complex FFTs on the hypercube. Issues such as load balancing and the mapping of data to the nodes are examined. A new load balanced method is presented and shown to execute faster than the naive method. This method also requires no extra buffers for communication and therefore is more storage efficient. This method is used subsequent implementations of symmetric transforms. In Chapter 4, tw?dimensional FFTs are discussed and implemented using four different methods that vary in degrees of concurrency, communication costs, and data transposition requirements. A model is derived that accounts for vector computers at the nodes. Some conclusions are drawn as to the effect of architecture and hardware parameters on algorithm choice. Real symmetric transforms are the topic of both Chapters 5 and 6. The diffi- culties of implementing them by the usual sequential methods involving pre- and 3 post- processing of the data is highlighted. Two new methods are presented that remedy this. The Chapter 5 method involves a modified pre- and post- processing scheme that computes the sine transform on data both in natural and bit-reversed order so that all post- processing is processor local. The method of Chapter 6 computes the sine, cosine, quarter-wave sine and quarter-wave cosine transforms of a single sequence in one pass through the complex FFT subroutine. This new method is hence extremely simple to implement and fast, as four symmetric trans- forms are obtained in one FFT with additional processing limited to multiplication by diagonal matrices. It is also a numerically stable way of computing symmetric transforms since all of the multiplications applied to the data vector are roots of unity. Finally, Chapter 7 describes a particular implementation of a distributed mixed radix FFT for hypercube multiprocessors, demonstrating that hypercube FFTs are not limited solely to radix-2 methods. 1.2 Parallel Scientific Computing The area of scientific computing is being revolutionized by the development of su- percomputers and parallel machines. Large problems can now be speedily solved so that engineers can redesign and test problem parameters with faster turnaround times. More design parameters can be incorporated into a simulation to better model the physical properties of a particular problem. However as promising as these new developments are, the problem-solving methods must be carefully thought through to ensure numerical stability and accuracy as well as compu- tational efficiency. Therefore using parallel processing architectures requires the rethinking of solution algorithms. Efficient implementation of numerical algorithms on parallel architectures de- pends on a host of system parameters, such as memory organization, processor power, communication interconnection between the processors and processor load balancing. A shared-memory machine consists of a number of processors that can address a single global memory. Memory access confficts occur when two proces- sors try to read or write to the same memory address. Therefore some sort of protection is utilized. This overhead can degrade any potential parallel speedup. A message-passing system consists of independent processors each possessing its own private memory with limited memory access. Communication between the processors is relegated to message passing. Here an interconnection scheme is very important. The ideal situation is a "crossbar" interconnection where each processor possesses a direct link to every other processor. However if the goal is to create a system with a large number of processors, such a communication 4 pattern becomes prohibitively complex. Nearest neighbor meshes and generalized hypercube structures [Bhuyan and Agrawal (1984)] are usually implemented so that processors communicate directly with a few neighbors and indirectly through message forwarding. The interconnection schemes and degrees of system coupling will affect the implementation of numerical algorithms. For example, a fully implicit finite dif- ference algorithm that is efficiently parallelized on a shared memory machine mav be difficult to implement on message-passing architectures. [George et al. (1987)] Algorithms that work well on serial computers such as the preconditioned con- jugate gradient methods for solving sparse systems of linear equations become prohibitively slow on parallel computers because global information is required at each iteration to check for convergence. Consequently local relaxation methods whose convergence rates may be slower can become viable candidates on parallel architectures. [Kuo, Levy and Musicus (1987)] Thus, the fast computation of large-scale scientific problems requires investiga- tion into algorithmic design, architectural parameters, and the development and application of expert systems to integrate and aid in the interpretation of the re- sults obtained. In this dissertation we concentrate on the study of implementing a set of FFT algorithms onto a loosely coupled parallel system with the hypercube interconnection scheme. 1.3 The Hypercube Architecture The d-dimensional hypercube is a parallel architecture composed of P --H 2d node processors with local memories and an interconnection scheme where two nodes are joined if and only if their binary node labels differ at exactly one bit. Therefore nodes whose binary labels have Hamming distance exactly one are connected. This setting allows the number of nodes to grow linearly with the interconnection complexity only growing logarithmically. The hypercube connection scheme is a popular one because many other topologies such as meshes, trees, rings, and butterflies can be embedded in the binary d-cube graph. There are several parameters for evaluating an interconnection scheme for a network of computers. We can represent a computer network as a graph G = (V, E), a set of vertices V and edges E. The vertices represent the node processors. and the edges are defined between nodes where there exists a direct link between them. Each node should have a reasonably small degree, so that each processor has to only handle a few links. The maximum distance (diameter) from one node to another should be small. And there should be a large number of alternative paths from one node to another for fault tolerance. The binary d-cube possesses 5 the following properties that make it attractive. o+ Each and every vertex has degree d; o+ G is a connected graph of diameter d; o+ If A and B are any two nodes in a d-cube, then there are H(A, B) parallel paths of length H(A, B) between the nodes A and B. Parattet paths denote two paths from A to B such that there are no common nodes on the paths except for A and B. H(A, B) is the Hamming distance between A and B. These and other properties of the hypercube can be found in Saad and Schultz (June, 1985). Parallel computers configured with the hypercube topology are usually message- passing systems. The system is a multiple-instruction multiple-data (MIMD) ma- chine where communication between concurrent processes occur only in the context of messages rather than the use of shared variables. Each processor has a fast local memory without any shared memory among the processors. Several commercial machines based on the hypercube architecture have been built. The first hypercube computer, the Cosmic Cube, was constructed in 1983 at the California Institute of Technology. [Seitz (1985)] Since then, vendors such as Ametek, Inc., Floating Point Systems, Intel Scientific Computers, and NCUBE have developed hypercube computers, some with vector processor nodes and some with transputer communication co-processors. [Wiley (1987)] The Cornell Theory Center's 16 processor Intel iPSC/D4MX is the system used for the results gathered in this dissertation. The hardware consists of a Multibus- based System 286/310 microcomputer (cube manager) and computational units of sixteen high performance microcomputers with extended memory. The cube man- ager is connected to each node processor by an Ethernet communication link. Each node processor has its own numeric processing unit and local memory. Commu- nication between nodes is achieved entirely by queued message passing. Messages can be sent from any node to any other node. The message system automatically routes messages through the system, forwarding them if necessary. One can think of "virtual" communications channels that connect all the processors. All of the algorithms in this thesis are coded in FORTRAN, compiled by the Ryan McFarland compiler, and executed under the XENIX 286 release 3.4 Update 1 of the host operating system and the iPSC release 3.1.1 of the node operating system with the Exelan R3.3 networking software. Although the programs are written for the Intel iPSC System 286 hypercube and our timings results are de- rived from this particular system, other hypercube are discussed. By considering their configurations and properties, recommendations on algorithm characteristics relating to specific properties are given. The current hypercubes by Intel, Ametek, 6 XCUBE and Floating-Point Systems all have slow communication rates relative to computational speed. Therefore the implementations strive to both minimize the amount of necessary communication and also the distance that they have to travel. Some hypercubes are equipped with vector boards on the nodes and therefore al- gorithms are evaluated with respect to vector length. Other hypercubes, such as the Floating-Point T-Series, can interleave communication and computation simultaneously, making it possible for the programmer to "hide" the communica- tion overhead within the computational workload. In general there is no specific algorithm or method that is the best for all situations. Furthermore the state-of- the-art hypercubes do not come close to supercomputer performance. Therefore these machines are used primarily as research tools in the study of parallel al- gorithins with consideration of architectural constraints and hardware properties. The absolute running times of the results obtained should not be the criteria for judging the algorithms. Rather, the timings for one method should be considered only as relative to the timings obtained from another method. Models are created for the various algorithms to provide general guidelines concerning communication overhead, vector length, and storage requirements. The timing results are used roughly to check the model and also to identify the proportion of communication overhead relative to computational effort. 1.4 Distributed Multiprocessor Computing Is- sues Programming paradigms on message passing multiprocessor systems requires a whole new set of heuristics defining such factors as task allocation, message pass- ing protocols, computational complexity, communication requirements, processor mapping and topology, storage space, and synchronization. Conventional rules of measuring execution time primarily with computational complexity analysis are superseded by new standards encompassing all aspects of distributed concur- rent programming techniques. [Adams and Crockett (1984), Gannon and Van Rosendale (1984)] Given the proliferation of new parallel machines spanning a whole spectrum of architectural configurations [Haynes et aL (1982)], the task of creating and implementing algorithms to effectively exploit parallel processing capabilities is a complicated endeavor. Communication is a very important topic in distributed memory machines. A processor can access data in its local memory much faster than data resident in another processor. In fact general hypercubes are designed to compute at a rate that is 10 times faster than the rate of communication. [Wiley (1987)] Recall that each node can communicate directly with its nearest neighbors. Sending a message 7 to a more remote node requires intermediate processors to relay the message. Currently communication between neighboring processors is more efficient than communication between more distant nodes. Hence an important consideration in hypercube programming is to define the problem so as to minimize communication, Two message passing primitives are used: sendw and recvw. A processor ex- ecutes a sendw to initiate transmission of a message to another processor. The calling process is blocked until the message has been sent. A processor uses a recvw to receive information. The execution of recvw causes the calling process to be blocked until the message has arrived. Therefore, if the receiving process invokes recvw prior to message availability, execution of the routine is delayed. The processor is then blocked. Processors do not need to be directly connected in order to communicate. The networking software is designed to automatically forward messages from node to node until the message reaches its destination. Hence the particular path that a message will traverse is entirely transparent to the programmer. The invocation of a communication request involves a latency period, or communication start- up time that is an order of magnitude longer than the actual data transfer rate. Moving a vector of length m from one node to a neighbor takes the time 7' + rfltcomm where 7' is the communication start-up or latency and the constant tcomm is the elemental transfer time. It is usually the case that 7' » tcomm. [Saad and Schultz (Sept. 1985)] Therefore it is much more efficient to send one long message than several short ones. This is especially true on the Intel iPSC. A measure of how well a procedure is utilizing its resources on a multiproces- sor is the ratio between the cumulative active time and the total duration of the computation. A processor is active if it is either doing arithmetic, data loading or rearranging, and actively communicating. It is inactive, or idle, while it is blocked awaiting data that it must receive before it can proceed with further computation. The utilization percentage determines the effectiveness of a particular procedure on the multiprocessor. Several concepts directly affect processor utilization. First of all there is the allocation of tasks to processors. Tasks that can be performed in parallel are ideally suited to multiprocessor applications. However these tasks should be evenly distributed among the processors so that they all have approx- imately an equal amount of work to do. This is called Load balancin,g and is important in that processors with a higher workload may cause other processors to block and wait for critical portions of the total procedure. Processor utilization of the blocked processors decreases and the total time for the whole procedure `ncreases. If workload is approximately equal, blocking time would be reduced along with the total execution time. 8 Efficient use of the limited local memories of the node processors is another important consideration. Some computations require the data in the memories of more than one processor. Therefore message organization is important to minimize the need for buffering and additional storage space. Overwriting is possible for the straight exchange of two sets of data yHx However, if the data requirements were y f(x,y) x f(x,y) overwriting an existing array, say x by the incoming message y is clearly an algo- rithmic mistake. Hence an extra buffer is required to store the incoming message. Therefore, rearrangements in the way of algorithm implementation and data map- ping into the processor nodes may be called for. The mapping of node processors into various topologies such as meshes, rings, trees, etc. is a distributed computing issue as well as an algorithmic abstraction. Problems modeling tw? or three- dimensional physical phenomena such as heat transfer on plates or gas dynamics work well on mesh topologies. Other appli- cations which require global broadcast or search algorithms find better matches on tree structures. Much algorithmic implementation work is tnerefore concerned with finding the best match between topology and data movement and permuta- tion patterns inherent in the algorithm. Finally it should be noted that the hypercube, being a message-passing com- puter system, achieves synchronization primarily by the availability and type of message being sent. A processor cannot request a message, but must wait for mes- sages that it needs to arrive. Messages are tagged with an integer constant type that is used to synchronize the total computational effort. For example, iteration numbers can be used as a message t?e to keep processors using the correct data for the correct iteration. In any case, careful design of algorithms is required to make sure that the entire system does not become deadlocked due to processes waiting for nonexistent messages. In this dissertation we review all of the above issues and considerations with respect to the implementation of Fast Fourier Transforms. 1.5 Definitions and Matrix Notation The Kronecker Product notation is a way of expressing block matrices as tensor products (direct products) of lower order matrices. It gives a convenient way of describing matrices which have a certain repeating structure. 9 Definition 1.5.1 Kronecker Product: If A E ?Pxq and B e cmxn then the Kronecker product A 0 B is the p-by-q btock matrix a00B ao,q?iB AOB= . . . ??Pmxqm - a??1,0B ap?i,q?iB (A 0 B)(C 0 D) = (AC) 0 (BD) if the matrix muitiplications AC and BD are defined. Under certain circumstances, a commutativity relation is useful and the following lemma is a useful tool to have. Lemma 1.5.1 Let p, q, r, s be integers with the property that i) pq = rs ii) p divides r then if A E CP?P and B E C??? (A 0 Iq)(Ir 0 B) = (Ir 0 B)(A 0 Iq) Proof See Van Loan (1987) Lemma 6.3. _____*_____ An easy method for designating submatrices is useful for the description of matrix algorithms or procedures. The MATLAB "colon" method is described here. [Coleman and Van Loan, (1987), Chapter 5] There are two forms of colon arguments (starting valuej : (terminating valuei being shorthand for all the indices from the starting value, with an increment of unit value. For example 1 (1,2,3,4, 5?. When the increment if any other nonzero negative, the colon format is value to the terminating 5 designates the indices integer, either positive or (starting value) : (increment) : (terminating value) An example is 1: 2 : 5 denoting the indices (1,3,5), whereas 8: --H3:1 represents the set of indices (8,5,2). 10 A submatrix can now be designated by use of the colon notation for row and/or column indices. A couple of examples suffice to illustrate. ... aiq B = A(i : j,p: q) B --H ajp o+. ajq Letting a(k) denote the kth column of a matrix A the submatrix consisting of the odd columns of A is denoted as follows: B = A(:, 1 : 2 : n) = [a('), a???, a?5?, . . Finally the diag notation if used to designate a diagonal matrix. The argument is usually a vector or an algebraic expression so that and D=diag(x)?D= xo Xn?1 D=diag(?i), a0 a1 a??1 1.6 Permutations A permutation is a function defined on a finite set X that maps X onto itself in a one-to-one manner. Permutations occur in the discussion and modeling of switch- ing networks for loosely coupled distributed computing systems. Data transfers which are one-t?one are actually permutations that map data from one proces- sor to another processor. Permutations are also an important part of the FFT algorithm, manifested in the recursive splitting algorithm that describes and char- acterizes the FFT. Permutations can be described in several intuitive ways. function: If ir : fO, 1,... , n --H 1) fo, 1,... , n --H 11 is a permutation, we denote it by 7r(Oo) ir(\) . 7r(fl?ffi)) 11 The set of all available permutations P(X) forms a group (P(X), o) under composition. The theory of groups can be applied to their study. [Gilbert (1976)] base-p function: Let 7r be a function on a set of integers that can be represented by an arbitrary base p. Some permutations are more conveniently described in this manner. Let fbt?i, b??2, . , bo)p bt?lpt?1 + + ... + bip + b0, O 0) is precisely 2. The indices are taken modulo 2d Proof [See Johnsson (1985)] *_____ The first property shows that the BRGC is indeed a Gray code as claimed. Each element differs from its neighbors at exactly one bit. The second property is the exchange property of the BRGC where elements which differ only at the jth bit are exchanges of each other around the jth dimension. And the last is a global connectivity property since elements which differ by strides 2? (j > 0), no matter how large, remain at a constant Hamming distance of 2. This distance also does not increase with the dimension of the Gray code. Let Bd be the binary encoding of the numbers 0, 1,... , 2d --H 1. We now show that the binary encoding does not possess some of the important properties of the BRGC. Claim 1. The binary encoding Bd does not possess the nearest neighbor property. Example: For d = 4, 0111 differs from 1000 in all 4 bits. Claim 2. The binary encoding Bd has Hamming distance j for elements that are exchanges around the jth dimension. Example: For d = 4 the elements 0000 and 1111 that are exchanges of each other in all 4 dimensions are at a Hamming distance of 4. 29 Claim 3. The i?inary encoding D??possesses the property such that elements that differ by 2' when i and i + 2' are in the same (j + 1)-dimensional sub-cube have Hamming distance exactly one. This can be seen from the fact that flipping the jth bit results in a binary number that differs from the original one by a stride length of a power of 2. Example: For d = 4, 0100 and 0110 differ at only one bit. Therefore one can conclude that while the binary encoding is superior in terms of the particular powers-of-two connections needed for the FFT butteffly it falls short with regard to nearest neighbor connections and exchange connections, with a maximum Hamming distance that increases with the dimension of the code. The BRGC, however incorporates the best compromise between these three properties Chapter 2 Multiplier Computation, Error Analysis 2.1 Introduction There is no point in having an efficient computational procedure if the results produced are not to the precision required. Errors in finite precision arithmetic are unavoidable. The representation of numbers by a finite number of digits means that most quantities cannot be represented exactly on the computer. If a number x is represented by the nearest floating point number with precision t and machine base b, then fl(x) = where u is the unit roundoff defined by bl--Ht = 2 lel < u When arithmetic operations are done on floating point numbers, roundoff errors occur. These errors can then be magnified by a particular numerical method The amount of magnification of errors is related to the stability of the numerica method and its error propagation properties. Sometimes, two different ways of computing a quantity, while equivalent mathematically, can produce drastically different results when implemented in finite precision arithmetic. See Colub and Van Loan (1983, pp. 32--H38), for details on rounding errors. Mathematically, the FFT is a sequence of sparse unitary matrix-vector multi- plications based on the recursion x(k) LA?x(k?l) 30 31 where Ak is defined in Chapter 1. unitary, we have at each stage Here we have n = pt Since each ?#`pAk is IJx(k) II --H llx(k?1) II --H = jx(O) and II?II = lIFnXII = II xli letting x --H x(0) and y = x(t). In what follows we denote Ak = --H?Ak to simplify The use of finite precision arithmetic, however, causes the matrix-vector prod- uct to generate an error vector (1k such that the computed recursion is Akx(k?l) + (1(k) Letting e(k) --H x(k) --H x(k), we see that the error vector is also built up according to a recursion, e(k) A?e(k?l) + d(k), e(0) --H 0 or lle(k) ii < IlAkil IIe(k?1) Ii + 11(1(k) ii, Ak is unitary so it does not magnify errors, hence JJe(k) < IIe(k?l) ii + 11(1(k) ?, where lIe(0)ii = 0 lie(o)Ii --H 0 lld(k)ll = 0(u). llx(k)ll Because the FFT is performed in t = log2 n stages, the final roundoff error should be of the form lIe(t)iI ? ctE, II?II EI < h?cu [2(c%cj?1)(1 + ?j) --H c%?2](1 + Ej) Ic% --H (2c?c%?1 --H cj?2)j +I2c?c%?i --H --H ci I2c?c%?i(6ji + ej) --H c%??ejI +I2c?c%?i --H 2cicj?iI + Ic%?2 --H c??2J 2uIc?'I + 2Ic?fIc%?i --H c5?iI +I2ci?iIIc? --H ?1 + Ic%?2 --H c>?2I 2UIciI + 2IclIhftRiu + 2Ic??iIh?Cu + hf?R2u + 0(u2) ______*_____ The roundoff error for each step if 2u. The magnification of the error from the previous cosine Cj?1 is reflected in the term 2lcilhfAu. The error in the j --H 2 cosine Cj?2 is represented in the term hF.R2u Solving the recurrence gives us hfR --H = 2IciIx5?' + xj?2 hfR --H (IqI+7m2+iy' The same analysis can be carried out for the sines. 39 *______ This algorithm is unstable because it magnifies roundoff error at each stage by roughly a factor of IciI+7m2+i, a quantity greater than 21ci1. Furthermore. the roundoff error introduced at each step propagates exponentially into subsequent steps. Repeated Multiplication Another way to, generate FFT weights is to observe that exp(--Hi(j?)) = cos(j0) --H i sin(59). Let w? = exp(--Hi(5?)). Then w?+? = w Algoritbm 2.3.4 Repeated Multiplication c(O) .` 1; s(O) .` 0 c(1) COS(()); s(t) SIN(O) for j = 2 : n/8 --H 1 c(j)?c(1)-c(j --H1)+s(i).s(j --H1) s(j)& --Hs(1).c(j --H1)+c(1).s(j --H1) end Each step involves one complex multiplication which breaks down to two multi- plications and one addition for each trigonometric function. The multiplication w - ??;? is a Given's rotation; C1 __ cl si c??1 3> - ?1 Cl - - ?j--H1 Given's rotations are orthonormal and therefore do not magnify roundoff errors. Their properties are presented in Golub and Van Loan (1983, pp. 43-47). Theorem 2.3.4 IIc% and s% are computed via Algorithm 2.3.4 then with h3R,M --H 2j Ic%--Hc> < hfMu+o(u2)+o(u2) <--H h,RMu+Qu2)+o(u2) 40 Proof Since ci and si are computed by calls to the cosine and sine function w( have: < hCu Is?--HsiI < hCu Let yj = frj, sjjT and G be the Given's rotation. Then ft(Gy5?1) = GYj + E, where We have so that 11E112 < llYi--Hi 1120(u) ft[f1(cicj?i) +f((sisj--Hi)] kcici?i)(1 +?i) + (sisj--Hi)(l +?2)](1 + ei) --H c?c%?i +5?5%--Hil +lc?&5?i --H qc1?1 + ls?s%--Hi --H sisi--Hil 2ulc%l + lc?llc%?i --H ci?il + lci?illc? --H cii +ls?lls%--Hi --H 5i--Hil + lsj--Hills? --H sil 2ulc%l + 1c?1h,R.Afu +lQ?ilh?Cu + ls?lhjR.?]if u + lsi?ilh?Cu + 0(u2) Since cii < 1 and isiI < 1, the error from the previous cosine hjR?Mi is not magnified and the absolute roundoff error, 2Uicji, decreases as Cj .` 0. Therefore hR.M = this is quite a liberal upper bound for the error. The analogous computations give the same result for the sines except that the roundoff error 2Uisj I increases from near 0 to 2u as si 1. *_____ Logarithmic Recursion The methods in the previous two sections call COS and SIN only once namei?- to compute ci and ?i, and thus suffered degeneracy of accuracy as j increases. To garner more accuracy throughout the range of j, we could use the approach 41 of computing a logarithmic number of cosines and sines exactly and derive the remaining ones. A natural distribution is to compute (cj:j=2k, k=0,1,...? (sj : j --H 2k, k = 0,1,. . by calls to the COS and SIN routine. Using the recurrence formula: cos(A + B) = (2 cos(B)) cos(A) --H cos(A --H sin(A + B) = (2 cos(B)) sin(A) --H sin(A --H we derive a method of computing the rest of the Cj'S and Sj'5. Of course, this necessitates the storage of all the previously computed weights. In forward recur- sion, the angle B was always 0 and A ranged from 0, ..... . ,etc. In logarithmic recursion, the angle B is 0 at first, then 20, then 40, etc. Meanwhile for each of the B's = 2k0, A will range from 0 to (2k --H 1)0. The angle 1? represents the cosines that were computed by calls to COS. The following procedure generates Cj and 5>, i = 1,... , fl? Algoritbm 2.3.5 Logaritbmic Recursion n = 2?, 0 = ?r12 c(1) COS(0); s(1) SIN(0) for k = 1: t --H 4 j?2k c(j) COS(j0) s(j) `. SIN(j0) for m = 1 : j --H 1 c(j +m)?2c(j).?m)--Hqj --Hm); s(j +m)?2?j).s(m)+s(j --Hm) end end --H 2t--H3 c(j) .` COS(j0); s(j) SIN(j0) Notice that we stay in the first quadrant so that cos(A --H B) = cos(B --H A). This explains the subtraction by Cj?m in the main loop. For the sines, we have sin(A --H B) = --H sin(B --H A) explaining the addition by 5j-m. 42 If n = 2?, then j ranges from 0 to n/8. Let j=(bt?2... b1b0)2 and define t--H2 Pj = ?bj. s=O P? is the number of bits in j's binary expansion that are equal to 1. Since all fj : j --H 2k? were computed by direct calls to COS or SIN, c%? = C2k + e and = 52k + E, i.e. hLR = hDc = 1. 2k 2k Theorem 2.3.5 If c% and s? are comp?ted via Aigorithm 2.3.5 then c%--Hc?j ? hfRu Is%--Hsjl < hfRu with IhfRI = Hbi=i(IC2iI + 7m.12+1) Proof First we have and thus 1c%--Hcj1 < fl[2f1(C%kCm) --H [2C2kCm(1 + ?1) --H C2k?mi(1 + ci) --H 2c%?cm --H C%k?mI +12C%kCm --H 2C2kCmI + IC%k?m --H C2k?mI 2UIciI + 2IC?kIICm --H cm +2ICmIJc%k --H c2kI + IC%k?m --H C2k?mJ 2U1cj1 + 2Ic2kIhmLRu + 2IcmIhf?Smu + h2LkR_u + 0(u2) fl(c>) --H ci < 4uIciI + 2I?k IhmLRu + 2IcmIh??Cu + hfkR?mu + 0(u2) There are a logarithmic number of multiplications and additions for each j so that roundoff is of order P? However, this method also suffers from the instability analogous to that of forward recursion, although not to such a great extent. The error coefficient for Cm is multiplied by Ic2k I + j\Th12+1 where P2k = 1. ?`e have Pj = 1 + Pm due to this bk?1 = 1 term. Since the error coefficient for Cm includes all the other Ic2i I + j\Th.I2+1 where bi = 1, the error for Cj is simply Ic2k 1+ ?I2+1 times the error for Cm and is thus a product of the Icil + ?j2+1 for bi ines decrease as the argument goes from 0 to 1 and therefore error magnification is not as bad as forward recursion. There are only P5 steps to 43 contaminate. Letting the binary representation of j be ....... , b0, an approximate bound for the magnification factor for Cj would be fi I%I+7m12+1. i:b?=1 This says that the only factors Ic2i I + YFThi12+i's that enter into the error of j are those where j's ith bit is equal to 1. The same bound holds analogously for the sines, although sines increase from 0 to 1. Subvector Scaling ______*_____ This method also computes the Cj'5 and 3j'5 where j = 2k, k = 0, 1,... by directly calling the COS and SIN routines The rest of the weights are computed by noting = W?U?2 ji . [1OS2 ji At each stage k we have the weights v(O : 2k --H 1) and we can compute w(2k 2k+I --H 1) = wn2k w(O : 2k --H 1) by scaling our previous vector of weights by wn2k. The following algorithms com- putes the weight vector by this method. Algorithm 2.3.6 Subvector Scaling n = no n/2 O=(nTo) c(O) 1; s(o) 0 c COS(O); s SIN(O) c(I) c; s(1) s for q = 1 : t --H 4 L: 2?;L0: L/2 for j = 0 : L0 --H 1 cCL0+j)? c c(j) + s s(j); s(L0+j)? --Hs c(j) + c s(j) end no n0/2 c COS(O); s SIN(O) c(L) .` c, s(L) s end 44 Each w3 needs only a maximum of t --H 4 multiplications. The niultiplications are all Given's rotations and hence do not magnify errors. The following theorem proves that there are a logarithmic number of multiplications. Here we let the complex weight vector w(0 : n/8 --H 1) represent c(0 : n/8 --H 1) + is(0 : n/8 --H 1). Theorem 2.3.6 Let n --H 2? and j = (bt?3... bo)2, where bj is the ith bit in the binary expansion of j. Then comp'uting the weight vector w(O :n/4-1) by Al9o- rithm ?.3.6 gives w(j) as a resitit of --H1 complex muitiplications. Proof This is based on induction with q = number of bits involved. q=1 j = 0 w(O) = 1, no multiplications j = 1 w(1)=wn, direct computation, no multiplications Now assume this is true for q = k, i.e. j =(bk?1... b1b0)2. w(j) is a consequence of (bk?l +... + b0 --H 1) multiplications by the induction hypothesis. There are two cases. w(O?j) = w(j) result carries over since bj = 0 here, hence true for q = k + 1. w(1j) = wq??.w(j) one multiplication + what we had previously since bk = 1 here, we have a total of (bk + bk?l + ... b0 --H 1) multiplications, and hence the result follows for q = k + 1. *______ Theorem 2.3.7 If c% and % are comptted via Algorithm 2.3.6 then jQ--Hc5j < h,SSu Is%?sjI < h,S.Su with h,ss = 2Pj, with Pj the number of bits in j's binary expansion that is equal to 1. 45 Proof Each Ci C[log2jj Cj?[1og2jJ + C?og2jJ 5i--H[Iog2jJ Therefore the errors for each Cj are those of Cj??og2 ?j and 3i--HUog? ?j plus the roundoff for this step. Each Given's rotation adds roughly 2u to the roundoff, but does not magnify them. From the proof of Theorem 2.3.4, we see that Ic% cil < 2UICiI + C[I0g2jJ?hjS?S?0g25Ju+ 1Cj??0g2jjhD?0Cg2jju +s?og2 jj hjS?S?0g2 jj UISj??og2 jj hD?5g2 jj U + 0(u2) Recall that Pj is the number of ones in the binary representation of j, hence hjSS = 2Pi The sines are computed in the same manner so that a similar accounting of errors gives the same bound. _____*_____ The error properties of this method are extremely good. The stability is that of the Given's rotations combined with only a logarithmic number of multiplications. Recursive Bisection A stable method introduced by Hopgood and Litherland in the 1960`s involves interpolation. [Oliver(1975)] The underlying relations are specified by the trigono- metric identities: (2.3-7) (2.3-8) cosAcosB = 2--H1fcos(A--HB)+cos(A+B)i sin A cos B = 21 fsin(A --H B) + sin(A + B)J Notice that A is an angle that is half-way between the angles A--HB and A+B. This means that if we knew the cosines or sines for A --H B and A + B, we can then find the cosine or sine of A by one multiplication, one addition and the computation of cos Suppose n --H 2? and 0 = 7?W2 and we wish to compute the cosines Cj = cosit and sines Sj = sinj?r2, j = 1,..., ? --H 1. Algorithm 2.3.7 Recursive Bisection c(0) 1; s(0) 0 c(n/8) #2/2; s(n/8) #2/2 for L = 46 n/2L; 5 n/2L c(n/2L) & cos(2T/2L) s(n/2L) sIN(2t/2L) h & 1/2c(n/2L) for m?1:2L/8 --H 1 5 5 + 2? c(j) &h.(c(j --Hi)+c(j +i)); s(j) &h.(s(j --Hl)+c(j +1)) end end Again the powers-of-two multipliers are computed by direct calls to the COS and SIN functions. Therefore we have (cj, Sj : 5 --H 2kJ. To get the rest of the c?'s and Sj'5 we must compute them via Equations 2.3-7 and 2.3-8. How many steps it takes to get to a particular 5 depends on the number search steps S? it takes to "find" 5 when initially only the powers of two are available. For example if we had (1, 2, 4,8, 16, 32J, S3 = 1, since it comes right in between 2 and 4, which we have. S6 = 1 also, since it comes right in between 4 and 8. S14 = 2, since we have to first get 12 between 8 and 16, and then 14 between 12 and 16. If we write out these numbers in binary, we see that S1 is equal to the position of the highest order 1 minus the position of the lowest order 1. Let 5 = ..... . b1b0)2 Lemma 2.3.1 (2.3-9) Sj = max(k: bk = 1i --H min(k : bk = Proof Since we initially have all the powers of two, the number of search steps needed to find them is 0. Each binary representation of 2k has a 1 at the kth position and zeros everywhere else. Thus S2k = k --H k = 0. (5:5 = 2k + 2k--H11 can now be found by averaging between 2k and 2k+1 for each k. These numbers have a 1 at the kth position and another 1 at the k --H 1th position, thus S2k+2k?1 --H k --H --H 1) = 1. Next all the numbers (5:5 = 2k + 2?--H2,5 --H 2k + 2k--H1 + 2k--H2? can be found by interpolating between (5:5 --H 2k, 5 2k+2k--H1? These numbers all have a 1 at the kth position and another 1 at the k --H 2th position. Some of them also have a 1 at the k --H 1th position and others have a 0 there, but for our purposes, it does not matter. We can now see a pattern emerge where all numbers (5:5 --H 2k + (?tk=kl?s+l bj2?) + 2k--Hs? have search number Sj = 8, representing s bisection steps. Sj = [log2 m?J where m5 = 5/(5.AND. --H 5) *______ Lemma 2.3.2 47 Proof In two's complement arithmetic, j.AND. --H j determines the number which is represented by the trailing 1. For example, 1010.AND. --H 1010 = 1010.AND.0110 = 0010 Dividing j by j.AND. --H j has the effect of stripping away all the trailing zeros. Hence for our example, 1010/0010 = 0101. Now taking the floor of the logarithm base two of the result gives us the spread between the highest order 1 and the lowest order 1. ______*_____ Let ?? = log2(j.AND. --H j), then i? gives us the number of trailing zeros in j. Clj is the cosine of the angle which is half of the spread between Cj?Ij and and similarly ?1j is the sine of the angle which is halfway between ?j?1j and Theorem 2.3.8 I/c% and % are computed via Aigorithm 2.3?7 then Ic%--Hc51 < h,RBu $s%?sj! < h,RBu with h?B = 2Sj where Sj is described by Equation 2.3-9. Proof We have Ic% --H ci < Hence Ic% --H 2c% (C%?ij + c%+i?)I 1 1 1 2c?? c??4 ?+? +1 c1?j. --H 2cj1 2C?j j+lj --H 2cj1 c?+11 2 + 1hRB + 1hR.B 1c41 c, 2IcijI>?'i 21c41 ?+4 + 0(u2) --H cii < 2u Icil + max(S5?'j,S?+1j)u +U(u2? iciji cii and thus hjRB = 2Sj. The sines are computed by the same formula so analogous results hold. ______*______ 48 Experimental Results The results using single-precision arithmetic on a VAX shows that subvector scal- ing with the direct computation of the powers-of-two weights fares the best overall with recursive bisection a close second. Table 2.1 shows the magnification factor for the maximum component error of a vector c consisting of cosines. This is computed as the maximum error divided by the single precision machine epsilon, which on the VAX is 2.98023e --H 08. Table 2.2 shows the root mean square error of the vector c and thus takes an average error of all the components. Again, this error is divided by the single precision machine epsilon. When Forward Recursion and Repeated Multiplication is implemented in dou- ble precision, we do not get magnification of errors until n is quite large with respect to ??1 This is, of course, a machine dependent parameter. For Forward Re- cursion, this starts at n = 32768, and for Repeated Multiplication, n = 2097152. Table 2.3 compares the magnification factors of these two method with respect to single precision machine epsilon. The use of double precision for Forward Re- cursion delays the appearance of significant amounts of error until about n = io7 since the single precision unit roundoff for the VAX is about 10--H8. See Table 2.3. 2.4 Multipliers for the FFT We do not have to worry about the error coefficient and its ordering and placement in the FFT algorithm for the method of Direct Call to Library Functions since it is assumed to be constant for each Cj and s?. However for the other methods each multiplier and its error coefficient h,method must be linked explicitly to the components of x that it is applied to. We must therefore be able to pinpoint the exact multiplier that is applied to each component at each stage. The radix-two FFT algorithm for a vector x of length n = 2? runs through t stages (k .?...., t). Denote xj(k) to be the jth component of x at the completion of the kth stage. At the beginning of stage k, the multipliers are applied to x(k?1). The particular multiplier used in the CT2 butterfly involving Xj is wJ;mod2(r12) r=2k. The angle used at stage k is 0k --H 2ir/2k, hence Wr is directly computed. Since w?mod2(r12) is used for transforming x,(k?1) to x,???, the error coefficient associated with Xj 15 hjm?2(q?), r = 2k 49 Table 2.1: Maximum Error Magnitude = IIcIIoo/u n Forward Repeated?? Logarithinic 5ubvector Recursive RecursionMultiplicationRecursionScalingBisection 8 5.0 1.0 1.0 2.0 0.0 16 47.0 2.0 5.0 3.0 2.0 32 34.5 2.0 11.5 3.0 2.0 64 321.0 2.0 16.0 3.0 4.0 128 50.3 4.0 21.0 6.0 4.0 256 6353.4 14.4 57.0 4.0 4.0 512 14705.1 17.5 95.5 6.0 6.0 1024 18818.1 8.0 151.4 6.0 8.0 2048 373315.0 21.0 172.0 6.0 8.0 4096 682278.0 12.0 519.5 6.0 10.0 8192 1.87937e+06 35.0 717.5 6.0 12.0 16384 3.35480e+07 34.0 1920.9 6.0 12.0 32768 3.35512e+07 38.0 5137.7 8.0 12.0 65536 3.35528e+07 24.0 6746.0 8.0 12.0 131072 3.35536e+07 47.5 7550.2 8.0 12.0 262144 3.35540e+07 136.0 7952.3 8.0 14.0 524288 3.55542e+07 492.0 8153.3 8.0 14.0 1048576 3.35543e+07 5300.7 8253.8 12.0 14.0 50 Table 2.2: Root Mean Square Error Magnitude = IICIJrms/u n Forward Repeated Logarithmic Subvector Recursive Recursion MultiplicationRecursion ScalingBisection 8 4.S 1.5 0.5 1.5 0.0 16 31.3 1.5 1.9 2.1 1.0 32 19.2 1.1 5.3 1.9 1.1 64 146.7 1.0 6.5 2.0 1.4 128 31.0 1.6 7.3 2.0 1.4 256 3238.7 8.2 19.4 1.8 1.8 512 7981.9 10.0 31.5 2.1 2.1 1024 10514.6 3.2 50.6 2.0 2.1 2048 198593.0 11.9 56.4 2.2 3.3 4096 357597.0 6.2 181.3 2.2 4.2 8192 1.34489e+06 20.8 235.7 1.9 5.0 16384 1.59806e+07 20.3 661.4 1.9 4.9 32768 1.59795e+07 16.0 1793.2 1.9 5.0 65536 1.59789e+07 14.6 2291.6 1.9 5.2 131072 l.59787e+07 26.5 2536.0 1.9 5.4 262144 l.59785e+07 55.3 2657.3 1.9 5.6 524288 1.59785e+07 264.8 2717.8 1.9 5.8 1048576 1.59784e+07 2781.9 2748.0 1.9 6.1 51 Table 2.3: Double Precision Implementations n Forward Recursion Repeated Multiplication MAX RMSMAX RMS 32768 2.0 0.1 65536 2.0 0.2 131072 2.0 0.5 262144 2.0 0.6 524288 4.3 2.3 1048576 16.0 8.3 2097152 65.6 34.0 4194304 66.4 34.6 --H 8388608 1230.2 635.3 3.1e-02 1.5e-05 16777216 6388.3 3299.3 6.le-05 33554432 14912.5 7799.4 1.0 2.9e-04 67108864 18841.0 10058.6 2.0 5.2e-04 means it is < u.) This means that we are using the (j mod 2k--H1 )th cosine and sine generated by a particular method. Let m,(?)CT2 =jmod2k?l denote the order of creation that the multiplier applied to x(k?1) to produce x?) at the kth stage of the CT2 FFT algorithm. This notation serves as a shorthand in our multiplier-connected FFT error analysis. The multiplier used on Xj to get Xj at stage k in the Gentleman-Sande (GSl) algorithm is Xj at stage k is w3jmod(s/2), ? = 2t--Hk+i Therefore, the corresponding shorthand is, m;.k)GSl =jmod2t?k. The above indices for the weights are applicable to the two methods, Forward Recursion, and Repeated Multiplication, where at each stage, a different ? --H 97r/9k is used as the initial angle. This is because these two methods, when used efficiently and intelligently, do not store a table of cosines and sines, and thus to capture the 52 best possible error properties, i.e not let 5 get too large, we generate the multiplie':.? at each stage with the largest initial angle possible. The methods requiring the direct calculation of a logarithmic number of weigh?? are usually implemented in an FFT subroutine by initializing a table first. This table is of length 0(n) and hence includes all the cosines and sines of angles which are a multiple of ?2?W? ? --H 9? At each stage in the FFT algorithm one would perform a table lookup for the correct weight. The correct error coefficient would be corresponding to the index in the table that the weights are placed, since it has something to do with the order with which the weights are computed. In addition. recall that not all the weights are actually computed. The symmetries of the cosine and sine function (eqns. 2.3-1--H 2.3-6) are all that are needed to determine the rest. Hence for methods requiring table lookup (Logarithmic Recursion, Subvector Scaling and Recursive Bisection), m(k)CT2 J = G(2t--Hk(5 mod 2k?1)) mod (2t?3) = c(2k--H1(s mod 2t?k)) mod (2t?3) where C() is the Binary Reflected Gray Code of Chapter 1. The method of directly calling the Library Functions is most sensible for gen- eral purpose FFT computations where a large number of FFTs are done. Here a 0(n) table is created once at the initialization phase. Swarztrauber's FFTPACN pretabulates the weights during the initialization phase. The Digital Signal Pro- cessing Committee (1979) FFT subroutines also use direct call, but without storing the weights in table. They simply compute the cosines and sines as needed. The FFT subroutines in Press et ?1. (1985) use double precision Repeated Multiplication for computing the FFT weights, as does the FORTRAN code in Swarztrauber (1984a). Forward Recursion is not used by anyone known to the author, for obvious reasons. Logarithmic Recursion is not used in practice but was described primaril,? to illustrate that magnification errors can overwhelm roundoff errors. Repeated Multiplication consists of a stable set of computations, but has 0(5) steps for each c5. Meanwhile Logarithmic Recursion has only 0(log2 5) steps per Cj, yet its absolute error is much worse than that of Repeated Multiplication because of the exponential magnification of these 0(u log2 5) errors. Subvector Scaling is also an O(log2 5) procedure but with a stable algorithm. Hence its error bound is proportion to log2 5. Recursive Bisection is used by Buneman (1987) in his FFT codes. Buneman (1987) also describes a way of generating the %`S and Sj'5 in the natural order by requiring only 0(log2 n) of table storage. 53 2.5 Error Analysis of the FFT Floating Point Complex Arithmetic In this error analysis, we compute the componentwise error of the Cooley-Tukey radix-2 FFT algorithm. Since the arithmetic done is complex, we start by deriving a model for floating-point complex arithmetic. Lemma 2.5.1 (Floating-Point Real Arithmetic) If a and b are floating point numbers and 0 is a binary operator, then by the model offloating-point arithmetic in Golub and Van Loan (1983) we have fl(aOb) = (acb)(1 + e), lEl < u where 1f1(aob)--H(aOb)j < u, aOb # 0 IaCbI A complex number x is represented by the ordered pair x = (xR,xI). Let z be a floating point complex number and 2 be the computed quantity for z. Let u be a machine dependent constant, floating point machine precision. Then the difference between 2 and z can be quantifled by ?(z) where 12--H zi < ?z)u + 0(u2). We use this function ?(z) to keep an account of the amount of error accumulated during a series of floating point computations leading up to z. First we establish the model for floating-point complex arithmetic as it relates to addition and multiplication. Lemma 2.5.2 (Floating Point Complex Addition) If x = (xR, xi) and y = (YR, ui) are floating-point complex numbers and if 2= fl(x + y), then 12--Hzi < u6(z) z = fl(x + y) = fl(xR + YR) + ifl(x1 + yi) (xR + YR)(l + ER) + i(xi + yi)(1 + ei) z + eR(xR + YR) + iq(xi + yi) ERI, lEil < u where 6(z) = Izi. Proof 54 hence Iz--H zi2 < IeRI2IxR + YRI2 + leil2ixi + yiI2 < u2(Ix+y2) *_____ Lemma 2.5.3 (Floating-Point Complex Multiplication) If x = (xR, xi) and y = (YR, ?i) are floating-point complex numbers and if z = fl(xy), then where ?z) = 21z1. Proof Iz--HzI < u?z)+o(u2) z = fl(xy) = [xRyR(l+?) --Hxiyi(1+?)](i +eR) +i[xIyR(l + ?) + xRyI(l + ?)J(1 + El) --H XRyR --H xIyI + ?XRyR --H 62x1y1 + ER(xRyR --H xiyi) +i[XiyR + xRyI] + i?xiyR + i?xRyI + iEI(xIyR + xRyI) J?1I, I?2I, i?3I, I?4I, leRi, Ifil < u Recall that the norm of a complex number z = (zR, zi) is equal to the two-norm of the two-vector zR , hence fj?\)Jj --H Re[z] _ ul Im[z]o+ 2 2 ?1XR --H?xI YR + ?xR - 2 - + ER(xRyR --H xryi) Ei(xiyR+xRyi) 2 uIxIIyRI + uIxIIyiI + uIzI uIxI#21y1 + uIzI (1 + #2)uIzI or jz --H < (1 + #2)uIzI + 0(u2) and Iz --H < ?z)u + o(u2) with?z)=(1+#2)IzJ. 55 ______*_____ Having defined the bound ?(z), we now need to know how they accumulate from previous floating-point computations. For example, given computed quantities 1 and z? with discrepancies ?(zi) and 6(z2), respectively; the quantity z = ft(2?Cz%) has a new discrepancy ?(z) derived from discrepancies of the two terms and a contribution from the operation C performed. Lemma 2.5.4 (Compounding Errors?Addition) If z = zi +z2 and z = fl(2? + z%) where then where Proof I2?--Hzil < ?zi)u+0(u2) Iz%z--Hz2I < ?z?)u+0(u2) Iz--H zi < ?z)u + 0(u2) 6(z) = Izi + 6(zi) + 6(z?) Ii--H(zi+z2)I < Ii--H(2? +z?)I+l2?+z??(zi+z2)I we know from the lemma on complex addition that the first term on the right satisfies so Ii --H (2? + z?)I < 12? + i?Iu + 0(u2) = IzIu+O(u2) 12--H(z1+z2)I <--H Iz--H (2? + z?)I + k? --H zil + kTh --H z21 (Izi + 6(zi) + 6(z?))u + 0(u2) _____*_____ The corresponding result for multiplication is described in the next lemma. Lemma 2.5.5 (Compounding Errors--HMultiplication) Suppose z = z1z2 and z=fl(2?z?) where I2?--Hzil < 6(zi)u+O(u2) Iz?--Hz2I < 6(z?)u+0(u2) 56 then where Proof 12--HzI 0, an absolute error constant. Hence the error is of O(log2 n + ?l?%?21fl h_(k)). We next examine the radix-2 GS1 FFT algorithm and note a few facts. From the properties of the GS1 algorithm, we notice that each Xj has the same sequence of multipliers applied to it as in the CT2 algorithm, except in reverse order, and with one difference. If at step k, the (k --H 1)st bit of j's binary representation is equal to 0, then no multipliers were applied because Xj is the top member of the butterfly and hence created by a simple addition. Since we are interested in bounds on the maximum component error, we shall first look at the error analysis for Xj such that j's binary representation has nonzeros in every bit. Later we can say that if 5's binary representation has a zero in position k --H 1, then the contribution of that particular multiplier is deducted. The radix-2 GS1 algorithm has a butterfly which is the transpose of that for the radix-2 CT2 FFT. As in lemma 2.5.7, we have the following. Lemma 2.5.8 Let BT _ - 1 w -w and BT --H BT + ET. Let ET - 0 0 S(w)u --H?(w)u 1x4 --H xTI ?x?)u +0(u2) Ix? --H xBI S(r?)u +0(u2) 1w--Hwi _ ?w)u+0(u2) + XT x = + xB - the maximum absolute error in x+T and x+B - BTx = [x4, x?jT IIx+IIoo + I?(x?)I + I?(x?)I (1 + ?)IIx+IIoo + (1 + W)IIxIIooI?w)I +(1 + #)IIxIIoo + I?(x?)I + I6(x?)I 62 Proof The GSl butterfly is defined by x4+ fl(x%+xB) xB+ fl(wft(xT --H XII)) The XT+ term requires only a simple addition and hence = IXT+I + I?XT)I + I?(XB)I < IIx+IIoo + I?(XT)I + 16(XB)I The X+B term has 21X+BI + I?w)IIXT --H XII + [IXT --H XIII + I?XT)I + I6(XII)I?IwI (1 + ?)IIx+IIoo + (1 + #2)IIxII?I6(w)I +(1 + #2)IIxII? + I?XT)I + I6(XII)I _____*_____ Notice that the error bound for XII+ in lemma 2.5.8 is similar to the error bound for X+II in lemma 2.5.7. Therefore as far as maximum absolute error is concerned the radix-2 GSl and CT2 algorithms are almost equivalent. This was shown already in Tsao (unpublished). The Multiplier History of the GSl Algorithm is more complicated that that of the Cooley-Tukey Algorithm. Theorem 2.5.3 (Mttitiplier Histor?--HCS1) If X,(k) stands for the jth component of x at stage k of a radi?-2 FFT of iength n --H 2? then X = (k--Hi) (k--Hi) J Xj + Xji x(k) --H where j' = ? + 2k--H1 and to the multiptier used on Xj at stage k then m5 po:nts letting M,(k) be the set of all indices used in the process of creating X,(.k) we have Mj(k) = M3(.k?l) u Mfk--Hl) _ --H Mj(k?l) u Mfk--Hl) U 63 Proof The GS butterfly shows that computing x3??? requires no multiplications. Therefore the multiplier history at this step is simply the union of the multiplier history MJk--Hi) and MJ,k--Hi) The computation of xj(k?) requires one new multipli- cation, hence the result. _____*_____ The imbalance in the way multipliers are applied means that in the step in the error analysis where we wish to recurse on the I?(x,???'?)I and I?(x,(??i))j where = ? + 2k--Hi, we cannot assume that they are basically equivalent. This is because at each butterfly step, the top element is a pure and simple addition with no multiplier applied to it. Therefore we concentrate our error bound on the element with the most complicated set of multipliers, the element x> where j = 11... 12 To say something about the other components we use the following lemma. Lemma 2.5.9 In every GS butterfly between x,(k?1) and where the (k--H1)st bit of j is 0 and the (k --H 1)st bit of j' is 1, > I?(xj(k))I Proof If j' has a 1 at position (k --H 1) then Xjt is the xB for that butterfly step. Similarly Xj, with j having a 0 at position (k --H 1), is the xT. Now using Lemma 2.5.8 we see that I?x+B)I > 16(xT+)I _____*_____ The theorem on the Gentleman-Sande Error Analysis is therefore basically a theorem about the last component. Theorem 2.5.4 (CSJ: Error Analysis). Let x(0) be a complex vector of length n --H 2? with components xj?0? such that --H x,(O)I 2icil, however to simplify the order of magnitude, we can plug in 21c11 into = (2IciI) * *(2k?1 --H 1) into the error bound we get, lix --H xiioo < O(t2?12)u. Repeated Multiplication: Here = 2(2k--Hi --H 1), hence lix --H xiioo < O(tn)u. The following three methods all have the maximum index j at each step equal to --H 1, hence the maximum log2 j is approximately t --H 1. Logarithmic Recursion: Here ??X(hmj(k)) < (2icii) * *(t --H 1), so lix --H xiioo < O(tn)u. For Subvector Scaling and Recursive Bisection, ??X(hm(.h)) < 2(t --H 1), and thus Subvector Scaling: lix --H xiioo < O(t2)u. Recursive Bisection: lix --H xiioo < O(t2)u. Therefore, we see that although the ideal error bound for the FFT is O(t)u this is achieved only through the use of Direct Call. This is counterintuitive and shows that the order of the error bounds for the multipliers affects the total FFT error bounds multiplicatively instead of additively. This does not contradict the result of Ramos (1971) because in his case, the error bounds for the multipliers are constant and hence are of the same order as the roundoff error. 66 Double precision computation of the multipliers can be used for the two O(t2)u methods and the O(n)u method, whenever flu ? 1 where u is the double precision machine precision. This is why double precision Repeated Multiplication can be used for practical implementations requiring speed and low storage requirements. [Press et al. (pp. 394--H395, 1986)] We summarize our novel approach to FFT error analysis with respect to mul- tiplier computation and give some conclusions. o+ We verified that errors in computing the multipliers are represented in the same proportion as the roundoff error if the error bound for the multipliers is an absolute constant. The order of the error incurred in the computation of the multipliers is multiplied by t = log2 n in the final error bound of the FFT. Hence methods that bound the errors in the multipliers by O(t)u actually have a total error bound for the FFT of O(t2)u. Using double precision in the Repeated Multiplication method, is accurate as long as flu < 1. Therefore the double precision, covers up for the O(n)u absolute error in the multipliers, making them effectively 0. 2.6 Practical Multiplier Generation The six methods for the generation of FFT multipliers differ in both their numer- ical characteristics and their implementation properties. There are two ways in which FFT multipliers can be used. They can either be pre-computed and stored in a table as Algorithms 2.3.2--H 2.3.7 show, or the particular multiplier needed at each step can be computed right before it is used. We call the latter "just-in-time" computation. Not all of our methods are amenable to just-in-time computation. For example. Logarithmic Recursion and Subvector Scaling requires storage space of 0(n) for the n/2 multipliers needed for the FFT of length n = 2?. This is because at each stage. the next set of multipliers computed is a "scaling" of the ones previously computed. On the other hand, the method of Direct Call can compute multipliers on demand without need for any table storage. After all, when a cosine of a particular angle is needed, it is obtained by a direct call to the library function COS. The other three methods, Forward Recursion, Repeated Multiplication, and Recursive Bisection vary in the amount of storage required for just-in-time computation. The basic equations for Forward Recursion are cos(5?) = 2cos(?)cos([j --H ijo) --H cos([j --H 67 sin(j?) = 2 cos(?) sin([j --H --H sin([j --H 2j?) Hence six real numbers need to be stored for each iteration. Algoritbm 2.6.1 Forward Recursion (just-in-time) = 2r/n Cp & 1: sp & 0 Cc & COS(8); Sc & Ct & 2Cc for j = 2,.. t & Cc Cc & Ct Cc --H Cp Cp & t t & 5c 5c Ct 3c --H 5p end Repeated multiplication requires storage for only two complex numbers since & w is the basic iteration. Algorithm 2.6.2 Repeated Multiplication (just-in-time) 0 = 2r/n w & 2exp(--HiO) Wj & W forj = 2,.. Wj & W Wy end Recursive Bisection is a bit more tricky. The most obvious implementation requires a full table of 0(n). This is because the Recursive Bisection algorithm computes a full table of multipliers by first starting with eiw/4, eir/2 = i, and eiO = 1, then creating eir/8, ?3rqa, and then eir/l6, e3ir/16, e5ir/16, e7ir/16, and so on. Initially a short table of half-secants, .5 sec(ir/4), need to be supplied. For example, if n = 2? and 0 = 2?/n, then (t --H 2) half-secants or 0(log2 n) storage is required. However to go from one step to the next requires all of the previously computed multipliers. So to get the odd angles of 7r/2k, requires one to have 68 available, the multipliers for the angles r/2k?1 Obviously, memory space for the full table is needed. Furthermore, for the purposes of the FFT, specifically the Gentleman-Sande algorithms, the multipliers are not created in the correct order. Recall that the Gentleman-Sande algorithms start with multipliers of the smallest spacings, i.e. for step k, k = 1,. .. , t, we use powers of wL = exp(2?i/L) where L = 2t--Hk+1 Therefore the spacing of multipliers is exactly the reverse of their creation by Recursive Bisection, where at step k, k = 1t, the multipliers wL, L = 2k--HI are computed. This spacing, however, is consistent with the Cooley- Tukey algorithm where wL, L = 2k--H1 is utilized at step k. In any case, if spacing and storage considerations are a factor, Recursive Bisec- tion, whether pre-computed or not, requires a table of 0(n) to store the multipliers and another table of 0(log2 n) to store the half-secants. Buneman's On-Line Generation of Weights Buneman (1987) proposes a new method of on-line (just-in-time) creation of FFT multipliers with the use of two tables of order 0(log2 n). The table is dynamically 3 updated and the multipliers are created in the consecutive order, i.e. W?, w?, W?? etc. An initial table is supplied and holds the cosines and/or sines of the angles 4?, 2r, ?, r/2, 7r/4, ....... , 2r/2t requiring t + 2 spots. Another table of length t + 1 contains the half-secants .5 sec(2?), .5 sec(r)5 sec(2r/2t). The idea is that whenever an entry has been used, it is replaced by a new entry that is to be used at a later step. Just enough information is retained in the table so that the appropriate new entry can be completed by mid-point interpolation. To describe Buneman's method, we first define a "degree" to be equal to the smallest desired spacing of angles, i.e. cosj =--H cos2i?/2t so that even though the conventional degree is 7r/180, our "degree" is typically 2r/2t Our initial table now looks like this ccs2t+1,cos2t,cos2t?1,...,cos4,cos2,cos1 or simply, (symbolically), 2t+1, 2t, 2t--H1, . ,4,2,1 These initial values are all powers of two, with each slot in the table reserved for the odd multiplier of the power of two held originally. For example, the last slot holds progressively increasing odd degrees, cos 1 replaced by cos 3, replaced by cosS, etc. The second to the last slot has 2, 6, 10, etc. 69 A concrete example illustrates the order that multipliers appear in the table. Let n = 16 and 2r716 be the smallest increment, that is a "degree" is equivalent to r/8. This suffices for a 16?point FFT. Initially, the table holds cosines of or degrees 47r, 2r, ?, r/2, r/4, 7r/8 32,16,8,4,2,1 Writing everything in binary, we outline the sequence of events in Table 2.5. Keeping Table 2.5 in mind, we interpret the sequence of events by first noticing that all entries in the same column have the same number of trailing zeros. This is because they are all odd multiples of the same power of 2. If we index the columns in reverse order starting with zero for the last column, we see that the column index is equal to the number of trailing zeros in the binary representation of the angle degree. Each interpolation involves the entry preceeding the entry to be changed and one further back. Let K be the index of the angle that had just been called, residing in column p. ThenKistobereplacedbyK+2.2P. LettingI=2Pweseethat???+?p The two angles needed for interpolation are K +1 and K +31. Since K is an odd multiple of 1, K +1 and K +31 must be odd multiples of 2q1, (q > p). Therefore K = (2j+ 1)1 for some j, then K+I= (j+1)21 and K+31= (i+2)21. Case I: j is even. Then K + 1 is an odd multiple of 21 and is therefore found in the entry immediately preceeding our current entry. K + 31 is found further back in the table. Exactly where can be determined by how many trailing zero bits it has. Case 2: j is odd. Then K + I is the one further back in the table and K + 31 is the one immediately preceeding our current location. Where to find an entry can be determined by testing the number of trailing zeros it contains in its binary representation. A way to do this is to use two-s complement arithmetic. I=K.AND.-K for which log2 I is the number we're looking for. We already know where one of the entries to be used in the interpolation is. It is the one immediately preceeding our current one. To locate the other, we first form K --H I and test to see if the bit representing 21 is set or not by forming M = 2I.OR.(K --H I). If this bit was not set, we have just set it and know that we have M = K + I in our preceeding entry, i.e. j was even. If this bit was already set, we still have M = K --H I after the OR. operation and know that j was odd. In any case M is a number with p + 1 trailing zeros and is one of the interpolants. The other is found by adding 21 to M; M being either K + I or K --H I gives either K +31 or K + I. This number has 70 more than p+ 1 trailing zeros and has in fact L = log2(5vI + 2I.OR. --H (Al + 21)) of them and is thus the Lth entry from the end. This procedure is mathematically identical to the recursive bisection algorithin given in Oliver (1975) and differs only in the order with which the trigonometric functions are calculated. The minimal amount of information is kept which is needed to generate the subsequent weights. As soon as the weight is used, it is replaced by one which will be used in the future. The following procedure is adapted from a FORTRAN program of Buneman's. The array h( ) contains the half-secants of 2?, Ir, ?2, ?` etc. The array c( ) initially contains the cosines of 8n', 4r, 2ir, ?, etc. 8? and 47r are included so that we can get back the original table automatically after one complete period. The half-secant of ? is infinite and is irrelevant in the cosine routine. For a sine 2 routine, just remember that sin(2j + 1)?r2 is always 0. Algorithm 2.6.3 On-line Cosine Generation (Buneman) N = 2i realh(j+2), cCj+4) for K = 1 : 2 * N /* find L, the number of trailing zeros in K */ I = K.AND. - K L = log2(I) /* c(j+?.L) contains the Kth cosine */ use it and then replace it */ write c(j+4--HL-) /* find where the other entr? for use in interpolation is Thcated */ 1= 2*I+(2*I.OR.K--Hr) c(j+4-L) = h(j+2-L) *(c(j+3-L) ?cCj +4-log2 (LAND. I))) end Buneman has used this algorithm for several years in FFT subroutines running on VAXes and CRAXs. In Chapter 3, we demonstrate a modified version of Buneman's method for multiprocessor FFTs. 71 Table 2?5: Buneman's On-Line Generation of Trigonometric Functions 32 16 8 4 2 1 100000 010000 001000 000100 000010 000001 32 16 8 4 2 3 100000 010000 001000 000100 000010 000011 32 16 8 4 6 3 100000 010000 001000 000100 000110 000011 32 16 8 4 6 5 100000 010000 001000 000100 000110 000101 32 16 8 12 6 5 100000 010000 001000 001100 000110 000101 32 16 8 12 6 7 100000 010000 001000 001100 000110 000111 32 16 8 12 10 7 100000 010000 001000 001100 001010 000111 32 16 8 12 10 9 100000 010000 001000 001100 001010 001001 32 16 24 12 10 9 100000 010000 011000 001100 001010 001001 32 16 24 12 10 11 100000 010000 011000 001100 001010 001011 32 16 24 12 14 11 100000 010000 011000 001100 001110 001011 32 16 24 12 14 13 100000 010000 011000 001100 001110 001101 32 16 24 20 14 13 100000 010000 011000 010100 001110 001101 32 16 24 20 14 15 100000 010000 011000 010100 001110 001111 32 16 24 20 18 15 100000 010000 011000 010100 010010 001111 32 16 24 8 20 4 18 2 17 1 100000010000011000 010100010010010001 Chapter 3 One Dimensional Parallel Fast Fourier Transforms 3.1 Introduction The FFT has a lot of inherent parallelism. The signal flow graph shows clearly that at each step of the radix-p FFT, n/p independent butterfly operations are done on n/p mutually exclusive sets of p of data-points. These butterfly computations. being totally disjoint, can be done in parallel. The FFT is also a global operation. Each point in the DFT is computed as a weighted sum of every other point. This can be seen by the representation of the DFT as a matrix-vector multiplication Fnx. The FF'T achieves this result efficiently by using a divide-and-conquer split- ting procedure. This means that for a parallel distributed computing system, the nIp parallelism of the FFT can only be exploited profitably if the data-points that make up the disjoint pairs can be efficiently placed exactly where they are needed when they are needed. Implementation of a parallel one-dimensional FFT is primarily a data-routing problem. One of the first things to notice when studying the radix-2 hypercube archi- tecture and the radix-2 FFT algorithm is the similarity between the adjacency scheme of the hypercube and the signal flow graph of the FFT. The FFT but- terfly data flow graph can be mapped into the hypercube graph so that adjacent nodes in the radix-2 FFT signal flow graph are also adjacent in the hypercube. In this sense the hypercube is an ideal topology for the implementation of the FFT. This explains why it is easiest to implement a radix-2 FFT on the radix-2 hyper- cube. Several authors have either implemented or discussed the implementation of the radix-2 complex FFT in the hypercube, including Chamberlain (1986), Chan (1986), Swarztrauber (1986) and Walton(1986). FFTs of mixed radix can also be implemented on the hypercube by way of the 72 73 "twiddle factor" approach where a one-dimensional array is written in terms of a two-dimensional array. This method is discussed in Chapter 4 along with two- dimensional FFTs because the data-flow requirements are identical. A variation to this approach allows us to implement FFT's of length ? = 2k m on the radix-2 hypercube. We discuss this mixed radix-2 FFT in Chapter 7. In this chapter we study the implementation of a single one-dimensional com- plex FFT. Issues such as data distribution, load balancing, communication require- ments, processor topology and distributed multiplier computation are addressed. Comparisons for the four canonic in-place radix-2 algorithms of Chapter 1 are ex- amined to determine properties of the FFT algorithm with respect to distributed processing. Deficiencies in the conventional procedures are identified and a new mapping of data into two columns (the Two-Track mapping) is proposed to rem- edy these problems. The new method is load balanced, executes faster and allows in-place computation of the FFT without requiring extra buffer space for mes- sage handling. We also address the issue of topology and show that a mapping of the hypercube according to a Binary Reflected Gray Code (BRGC), which is not isomorphic to the FFT butterfly signal flow graph, is a perfectly viable way of implementing the FFT. This mapping is also more versatile and better suited for use on general numerical computational problems involving FFTs and other com- putations that require either nearest-neighbor connections and/or the exchange permutation. The last section of this chapter examines distributed multiplier computational issues and describes how each of the methods in Chapter 2 can or cannot be implemented efficiently on the hypercube. 3.2 One-Dimensional Complex FFT on the Hy- percube Let x be a vector of length n where n --H 2'. The vector x can be distributed in a hypercube of dimension d consisting of P nodes (P = 2d and d < t) by partitioning x into subvectors Xj each of length n5 = n/P, xo xl Xn?1 and placing Xj in processor?. This is called consecttive storage. Here it is helpful to imagine the hypercube as a column of P nodes. We call this mapping the natural ordered mapping of input in natural order. See Figure 3.1. 74 poo po? plo p1i Figure 3.1: n = 8, d = 2 mapping of data in natural order As noted before, in-place FFTs that take natural order input produce output in bit-reversed order. Specifically the procedures CT1 and GS1 ove?vrite x PnFnx. It is practical to also be able to implement the FFT algorithms that take bit-reversed input and produces natural ordered output. The CT2 and GS2 radix-2 procedures overwrite x FnPnX. The input is mapped so that Pnx is stored in the hypercube in consecutive order shown in Figure 3.2 The most direct method of implementing a distributed FFT on the hypercube is to map the input vector into processors consecutively in the natural order. [Cham- berlain (1986)] The hypercube graph is a compressed version of the radix-two FFT butterfly data flow graph. The signal flow graph for the original Cooley-Tukey al- gorithm (CTl) in Figure 3.3 shows how easy it is to map the FFT algorithm onto the hypercube architecture. For a transform of 2? points consecutively stored in the d-dimensional hypercube, the first d butterfly steps involve across processor operations and hence require communication. We refer to these steps as the in- termingling steps. The remaining t --H d butterfly steps are local and hence do not require communication. We call them the paratlel butterflies. The implementation of the intermingling butterflies requires communication. 75 pOO pOl plo p1i Figure 3.2: n = 8, d = 2 mapping of data in bit-reversed order Algorithrn 3.2.1 shows the direct implementation of the Cooley-Tukey butterfly I ? x x+Ay I --H? --H x--HAy and Algorithm 3.2.2 describes the Gentleman-Sande butterfly I I x = x+y ? --H? - ?Ve establish some general conventions for the description of the algorithms. ? A generic name referring to the node that the process is running on. On the Intel iPSC, this is the value of the function mynodeo). p A generic name referring to the processor that communication, whether a send or receive, is being done with. send(v, Lvp) A primitive containing three separate arguments. v is the address of the array or vector that is to be sent to processor p. Lv is the length of the vector or array sent. (Note: this is not the length in bytes, but rather a higher-level abstraction that we usually think of so that 1 unit in length is Figure 3.3: n = 8, P = 4, signal flow graph (GSl1 CTl) r J xo x2 x3 x4 x5 x6 xo x4 x2 x6 xl x3 76 X7 X7 77 equivalent to 1 element in the vector. If A is an m-by-n matrix written in column-major order, then sending A is represented by send(A, m n, p). recv(v, Lv, p) A primitive designating that the calling process is to receive an array v of length Lv from the process running on processor p. b A workspace buffer. n A constant referring to the length of the transform. t n --H P The number of processors. d The dimension of the hypercube p --H 2d, The index number that represents ?`s place in the node order that the hypercube is mapped in. For example B2 = fo, 1,2, 3J is the binary mapping of the nodes in a 2-cube in natural order, then id[2] = 2. However if G2 = fo, 1,3, 2? is the Binary Reflected Gray Code mapping of the nodes in a 2-cube, id[2J = 3. The basic Cooley-Tukey distributed butterfly involves two processors that send each other their data. One processor is responsible for the top portion of the butterfly x + Ay and the other one the bottom portion x --H ?y. The processor responsible for the bottom portion also computes ?y before the communication so it has a higher workload. Algoritbm 3.2.1 Cooley-Thkey Butterfly = processor id, n --H 2?, P = 2d, m = n/P */ processor ? has x(id[MJ m : (id[?] + 1). m --H 1) */ p = processor id of partner in this particular butterfly */ if (id[?] 0. The subscripts are to be taken modulo 2d This means that processors are at most a Hamming distance two from any processor which is ordered a power-of-two apart. Table 3.4 summarizes the implementation of the CT2 and GSl algorithm on a hypercube mapped in BRGC order. We choose to discuss the two variants where the multipliers are computed in natural order because they are simpler and more convenient to implement. The communication time ranges are roughly the same as those for the hypercube mapped in Natural order, see Table 3.3. 3.4 The Two-?ack FFT Implementation We now describe a new method that solves the problems of load imbalance and extra buffering needs with just minimal communication. We have already seen in Section 3.2 how an attempt to correct the load imbalance and buffering problems requires a processor to communicate two times during an intermingling step. By considering the second exchange in Walton's implementation we see that it is unnecessary. Recall that in efficient distributed in-place FFT implementations, no attempt is made to perform the bit-reversal permutation. Hence either the input 88 Table 3.4: Forward FFT (BRGC Implementation) total 9450 5420--H5435 3065--H3150 1720--H1790 955--H1030 n = 4096 CT2 comp comm/blocked 4595--H5385 50--H825/10--H790 2235--H3020 125--H830/10--H790 1090--H1680 105630/10--H595 725--H925 40--H390/15--H440 dim 0 2 3 4 dim 0 2 3 4 total 9205 4545--H5290 2555--H3305 1415--H1935 795--H1095 GS1 comp comm/blocked 4490--H5230 55/5--H10 2190--H2925 365--H375/70--H295 1075--H1630 300--H395/50--H270 530--H800 190--H300/90--H270 or the output is in bit-reversed scrambled order. Our new implementation method forgoes the second exchange and in doing so, introduces an additional scrambling of the data. This scrambling is shown to be manageable. The important point is that it is reversible. Hence our implementation achieves load balancing, in-place computations, no need for additional buffering space, as well as non-redundant communication requirements. The FET of a vector x of length n --H 2t) is broken into two pieces, the first half denoted by x(1) and the second half x(2) and distributed into P processors. ?# partition x(') and x(2) into subvectors of length n/2P, x?j'? and x?j2?, respectively. These are then placed into the processors such that the xs(k) are placed in processor i, with k = 1, 2. We define the left track to consist of the subvectors of x(1) and the right track the subvectors of x(2). At each intermingling step, each processor exchanges its portion of x?2? for the portion of x(1) that its exchange partner for that step has. The butterfly computation is done and the results are not sent back. Instead, we just accept that an additional permutation is put into the data. The inverse transform can be done by starting out with this permuted data and reversing the exchange steps until the final result is again in the correct order. This scheme is already load balanced and does not require extra buffers as it is similar to Walton's method. However in addition, by accepting an extra permutation, we have eliminated the 89 Step 1 Step2 Step 3 Step 4 pOo 0 8 0 4 0 2 0 2 19 1 5 1313 poi 2 10 2 6 4 6 4 6 3 11 3 7 5 7 5 7 pil 4 12 8 12 8 10 8 10 5 13 9 13 9 11 9 11 pi0 6 14 10 14 12 14 12 14 7 15 11 15 13 15 13 15 Figure 3.8: Two-Track Forward Transform Data Swapping Pattern need for extra communication. The permutation is illustrated with a simple example: (n = 16), (P = 4) in Figure 3.8, with elements for the next exchange at each step highlighted in bold. Notice that the last stage involves no intermingling and does not swap data between the columns, because parallel butterflies take place independently without need for communication. An inverse transform with data originating in the scrambled order is done with sub-vectors swapped in the reverse order as depicted in Figure 3.9. Figures 3.10 and 3.11 illustrate Two-Track Iiiiplementation of the Gentleman- Sande and Cooley-Tukey distributed butterflies and Algorithms 3.4.1 and 3.4.2 show how a distributed butterfly is done in the Two-Track method. Recall the two primitives send and recv used for the description of the algorithms. We also use a function flip(m, v) that takes the binary expansion of an integer m and complements the vth bit of that expansion. For example, flip(19, 4) = 27, since 19 = 100112 and 27 = 110112. The Gentleman-Sande GS1 FFT algorithm utilizing the Tw?Track data map- ping performs the distributed butterflies first. Each update step precedes the communication step during the distributed butterfly stages. And the last t --H d steps are done independently and locally in the processors. 90 Step 1 Step2 Step 3 Step 4 POO 0 2 0 2 0 4 0 8 138 3 1 519 pOi 4 6 4 6 2 6 2 10 5 7 5 7 3 7 3 11 pil 8 10 8 10 8 12 4 12 9 11 9 11 9 13 5 13 plo 12 14 12 14 10 14 6 14 13 15 13 15 11 15 7 15 Figure 3.9: Two-Track Inverse Transform Data Swapping Pattern Algorithm 3.4.1 GSl Two-Track FFT = processor id; n --H 2t, p --H 2d, m = (n/P)/2 */ /* processor M hoids: */ /* z(') --H x(rn . id[?] : m (id[ji] + 1) --H 1) */ /*z(2) =x(--H?2+rn.id[?: --H?2+??(id[?+1)--H1) */ 5 tells which m?ltiplier to initialize */ s n; for L = 1: d p = ii/m (mod s/2) w for j = 1 : m t z(1)(j) --H z(2)(j); z(1)(j) z(1)(j) + z(2)(j) z(2)(j) w t; w w? w end 5 s/2 /* comm?nicate */ v = d --H L; id[p] = ftip(?, v); p = id?'(id[p]) send(z(2), m, p); recv(z('), m, p) end /* do local FFT (CSlaigorithm */ call localGSi(z('), z2),m) The Cooley-Tukey CT2 FFT using the Two-Track data mapping commences 91 1 2 3 4 send x?2? recv y(l) update x?2? "??` to pi from pi y(l) P0J x(1) x(1) + y(l) Figure 3.10: Two-Track Implementation of the Distributed OS Butterfly y(2) y(2) y(2) I. __________ P1I y(l) x(2) y(l) y(l) recv x(2) from P0 update x(2) x(2) + y(2) y(2) x(2) y(l) ??? x(2) y(l) 92 P1I y(l) y(2) y(2) y(2) I __________ P0 x(1) 1 2 4 send x(2) recv y(l) x(2) "?" to pi from pi y(l) y(l) send y(l) recv x(2) to PC from PC x?2) 3 update y(l) x(1) + y(l) y(l) x(1) --H y(l) x(1) update y(2) x(2) y(l) ? x(2) y(l) x(2) x(2) + y(2) y(2) x(2) --H y(2) Figure 3.11: Two-Track Implementation of the Distributed CT Butterfly 93 with local independent FFTs and ends with distributed butterflies. Here the communication is done before the updating. Algorithm 3.4.2 CT2 Two-Track FFT = processor id; n --H 2t, p --H 2d, m = (n/P)/2 */ processor ? holds: */ /*z(1) --Hx(m?id[?]:m.(id[?]+1)--H1) */ /*z(2) x(--H?2+m.id[?]: --H?2+m?(id[ii]+1)--H1) */ /* do local FFT */ call ThcalCT2(z(1), z(2), m) s .` 2 for L = 1: d /* communicate */ v = L; id[p] = fhp(?, v); p --H id?'(id[p]) send(z(2), m, p); recv(z('), m, p) /* update */ WO p = ?/m (mod s/2) 4'? .. 4'?op for j = 1 m t??.z(2)(j);z(2)(j) ?z(1)(j)--Ht z(1)(j) .. z(1)(j) + t; 4'? w? 4'? end 5 ` 2 s end Data Permutation of the Two-Track Method The Two-Track method achieves its objective by utilizing an unconventional map- ping of data into the processors. The sequence of data transfers results in an output array that is permuted in another ordering. We describe this data permutation pattern in the language of permutation matrices and Kronecker products. Let x be a vector of length n = x(1 inside of P --H 2k processors. Every processor has a piece of x(') as well as a piece of x(2). During the communication of the Two-Track method, portions of x?'? are 94 swapped with portions of x(2), resulting in a new vector x which is also split into x(') and x(2) mapped into two columns. Algorithm 3.4.3 gives the data permutation pattern of the Forward Two-Track method. Recall from Chapter 1 that Bn is the matrix representing the Butterfly Permutation and Mn is the inverse perfect shuffle matrix of a sequence of length Algorithm 3.4.3 forj = 2 : k+ 1 r = x (Br ?InIr)X end From Theorem 1.6.6 we have (M2?2p? 0 1n12p)x The Inverse Two-Track data pattern starts from = (M2?2p? 0 Ln/2p)x and works backwards as in Algorithm 3.4.4 Aigorithm 3.4.4 for j = k + 1 : --H1:2 r = 2' (Br 0 Iqr)x end The transpose of Theorem 1.6.6 gives, x (112?2P? 0 1n12p)x the original x. If n = 2P, then the result of the forward Two-Track data permu- tation is that x --H M?2?x and similarly the Inverse Two-Track pattern gives back Timing Results The timing results in Table 3.6 were done with the processors mapped in natural order and those of Table 3.7 with processors mapped in BRGC order. The relevant lines of Table 3.3 are reproduced in Table 3.6 in italicized characters for easy comparison. Similarly, lines in Table 3.4 appear in italics below their corresponding entry in Table 3.7. 95 Table 3.5: Two-Track FFT FORTRAN Code Algorithm Code Number CT2 (3.12), (3.17) csl (3.13), (3.18) both (3.14), (3.15), (3.16) Table 3.6: Forward FFT (Tw?Track) Natural dim 0 0 2 2 3 3 4 4 0 0 2 2 3 3 4 4 total 8685 9915 4655 5655.5660 2485--H2490 3305-3320 1325--H1335 1870-1880 745--H765 1045-1055 n = 4096 CT2 comp 8685 4625--H4630 4830-5610 2445--H2460 2365-3270 1300--H1305 1165--H1840 690?95 565-1025 comm/blocked 30/0--H10 50-825/5-790 30/0--H10 50-940/10-900 25--H30/5--H10 40--H700/10-675 50--H70/0--H55 30-470/5-440 8365 9255 4490?30 4605-5280 2415--H2445 2290-2985 1295--H1320 1145-1670 700--H760 640-995 GS1 8365 4465?500 25--H30/5 4550-5225 55/5-10 2385--H2415 30/5 2235-2925 55-60/5-10 1265--H1285 25--H30/0--H10 1100-1620 45--H50/5-10 675?90 20?5/065 555-900 80-95/5-70 l7alicized timings are for the One-Track Method. 96 The implementation code (in FORTRAN) appears in the Appendix. Table 3.5 tells the subroutine number associated with each Two-Track FFT algorithm. As one can see, in each case these times for the Two-Track method are balanced and moreover faster than the slowest processor in Section 3.2. This is because a lot of the blocking inefficiencies have been eliminated. To summarize: o+ Load balancing is achieved by placing sub-vectors so that complete butterflies are done inside each processor during the intermingling steps. o+ Only n/2P elements are exchanged at each step rather than n/P elements. o+ No extra buffering is needed for the distributed butterflies. o+ Data is not sent "home", hence only one exchange per intermingling step is utilized. Timing results are faster than for one column method, because each processor only has to calculate n/2P multipliers rather than n/P multipliers (for the two algorithms with natural ordering of multipliers). The speedup of the Two-Track method over the one column method increases from 10% to 31% with increasing dimension for the GS1 algorithm, and from 14% to 38% with increasing dimension for the CT2 algorithm. 3.5 Distributed Multiplier Computation In distributed FFT's each individual processor requires a set of multipliers that is different from those required by other processors. This is because each processor is responsible for only a subset of the FFT butterflies. Let W? = e?2"?1? where n = 2? and p be a integer ranging from 0, . . ,n/2 --H 1 Then the powers of W?, p required by processor x, in natural order, are exactly those p whose binary representation satisfy the pattern, xi, j=O,l,?????2flp 1 (?denotes concatenation of binary numbers), and x is assumed to be in binarv form. (0 < x < P --H 1), with P --H 9d the numbers of nodes in a d-dimensional hypercube. The mapping the integers 0n/2 in consecutive order through a hypercube whose nodes are in natural order shows this. The highest order d bits denote the processor number. 97 Table 3.7: Forward FFT (Two-Track) BRGC n = 4096 total 8685 9450 4655 5420-5435 2565--H2650 3065-3150 1405--H1420 1720-1790 760--H870 955-1030 CT2 comp 8685 comm/blocked 4630 25--H30/5 4595-5385 50-825/10-790 2450--H2465 100--H195/0--H105 2235-3020 125-830/10-790 1295--H1305 100--H115/10--H95 1090-1680 105--H630/10-595 695--H700 65--H175/10--H155 725-925 40-390/15-440 dim 0 0 2 2 3 3 4 4 dim 0 0 2 2 3 3 4 4 total 8365 9205 4490--H4530 4545-5290 2525--H2560 2555-3305 1395--H1460 1415-1935 790--H865 795-1095 GS1 comp 8365 comm/blocked 4460--H4495 30/5 4490-5230 55/5-10 2380--H2410 125--H175/5--H110 2190-2925 365-375/70-295 1265--H1290 115--H175/50--H140 1075-1630 300-395/50-270 670695 110--H185/25--H165 530-800 190-300/90-270 *?Itahcized timings are for the One-Track Method. 98 Example n = 32, d = 2 pOO : foo?oo,oo?oi,oo?io,oo?ii? pOl : (o1?oo,oi?oi,oi?io,oi?ii? pi0 : (1o?oo,1o?oi,io?io,io?ii? pi1 : (ii?oo,ii?oi,ii?io,ii?ii1 Since the left-most d bits designate the processor number, we refer to these 4 bits as the processor field. Meanwhile, the right-most t --H d --H 1 bits run from 0, .., --H 1 and are designated the working field. We now show how each of the five algorithms from Chapter 2 can be distributed in the multiprocessor setting. We do not use the symmetry described on page 34 since each processor only has a selected subset of the trigonometric functions. Algorithm 3.5.1 Direct Call = processor id, n --H 2t, p --H 2d, m = (n/P)/2 */ ()= 2r/n --H 2t--Hd--Hi for 5 = 0 : c(j) COS(x +5)6; s(j) SIN(x +5)6 end Algorithm 3.5.2 Forward Recursion = processor id, n = 2t, p --H 2d, m = (n/P)/2 */ 6 = 2r/n = 2t--Hd--H1 c COS(6) tcl 2 c c(O) CO5(x6); s(O) SIN(x6) c(1) COS(x + 1)6; s(i) SIN(x + 1)6 for 5 = 2 : m ?j)?tc1 c(j --Hi)--Hc(j --H2) s(j)?tc1 s(j --H1)--Hs(j --H2) end 99 Algoritbm 3.5.3 Repeated Multiplication /*? = processor id, n --H 2t, p --H 9d m = (n/P)/2 */ = 2?/n = 2t--Hd--H1 . c COS(?); s SIN(6) c(O) COS(xO); s(O) SIN(x?) for j = 1 : m c(j) c c(j --H 1) + s s(j --H 1) s(j) --Hs.c(j --H 1)+c.s(j --H 1) end Logarithmic Recursion is not practical in the multiprocessor setting. A simple example illustrates this. Let n = 16, n/2 = 8, and d = 1. Let Cj denote cos(j()). Recall that Cj = c?5 if jO is in the first quadrant. PO Co = 1 Ci = cos(0) c2 = cos(20) C3 = 2C2 Cl --H C?l pi C4 = cos(4?) C5 = 2c4 --H C?3 C6 = 2c4 C2 --H C?2 C7 = 2C4 C3 --H C?1 One glance at the data that processor pi needs to compute C5, C6, and C7 shows that they are not local and hence require communication. Algorithm 3.5.4 Subvector Scaling = processor id, n --H 2t, p --H 2d, m = (n/P)/2 */ 8= 2?/n; m'. n/2 = 2t--Hd--H1 c COS(O); s SIN(8) c(O) COS(xO); s(O) SIN(xO) c(1) `. COS(x + 1)8; s(i) SIN(x + 1)8 forq= 1: t--Hd--H 1 L 2q; L0: L/2 for j = 0 : L0 --H 1 c(L0+j) ?c.?j)+s.s(j) s(Lo+j)?--Hs.?j)+c?s(j) 100 end m H m/2 = rIm c COS(?); S H SIN(?) c(L) c; s(L) 5 end Algorithm 3.5.5 Recursive Bisection = processor id, n --H 2?, P --H 2d, m = (n/P)/2 */ = 2r/n 2t--Hd--H1 . id[?]; y 2t--Hd--H1 (id[?] + 1) c(O) COS(xO); s(0) H SJN(x?) c(m) COS(y?); s(m) SIN(y9) for L = 1 : t --H d --H 1 m/2L; j m/2L c(j) COS(x + i)(); s(j) SIN(x + i)O h 1/2c(j) for k = 1 : 2L--H1 --H 1 j ?j+2i ?i)?h?(c(i?i)+c(j+i));s(j)?h.s(j?i)+c(j+j)? end end Storage Considerations The methods listed above generate a table of length (n/P)/2 complex numbers as the multipliers needed by each processor. This may or may not be desirable depending on whether there is sufficient storage and/or multiple transforms are done. If only one FFT is done and array storage space is precious, then methods which do not require table storage should be utilized. Such is the case with our implementations, where the multipliers are computed as needed. Recursive Bisection Recall that Buneman (1987) describes a procedure that implements the Recursive Bisection algorithm with tables of O(log2 n) rather than 0(n). 101 Now suppose we are in a multiprocessor setting where each processor requires only a subset of the weights. This situation is applicable to distributing FFT butterflies among many processors. Let n --H _ so that the multipliers needed have arguments 2rj/n ?=0,...,rn?1 2 Mapping the weights in the natural order among p --H 2d processors, d < t, we see that processor k needs to have the weights coskj, 5=0,???,2t--Hd--H1 (?denotes concatenation of two binary numbers). Example (t=5,d=2) p(OO): 0O?00O, O0?001, O0?O10, 00?011, 00?100, 00?101, 00?110, 00?111. p(01): 01?00O, 01?001, 01?O10, 01?011, 01?10O, 01?1O1, 01?11O, O1?111. p(10): 10?000, 1O?001, 10?010, 10?011, 10?100, 10?101, 1O?11Q 1O?111. p(11): 11?0O0, 11?001, 11?010, 11?011, 11?10O, 11?1O1, 11?11O, 11?111. Notice that the left-most two bits designate the processor number and the right-most ones run from 0 to 7. In general, we refer to the left-most d bits as the processor field, and the rightmost t--H d --H 1 bits, the working field. If we ignore the processor field, we notice that each processor needs to be able to generate 2t--Hd--H1 sines and/or cosines. This means that it must have a table of length t --H d + 1 to follow Buneman's scheme. In fact the procedure is exactly the same for each processor as if it were working on a 2t--Hd--H1 problem except for the initial values. r Instead of initially having the cosines of 2?, ir, 2' ?? etc., processor k would have the cosines of 2t--Hd--H1 + 2ir, 2t--Hd--H1 L:+r, 2t--Hd--H1 L?r+?r2, 2t?d?1.L: + ?r4, etc. In other words, the initial values are the cosines of the angles we would have had plus the quantity 2t--Hd--H1 k?r. The characteristics of the working bits are all that are used, even though the initial and subsequent values are for angles which have degree (+2t--Hd--H1 To illustrate, this is what processor 1 initially contains in our example: cos4O cos24 cosl6 cosl2 coslO cos9 cos 8 must also be computed and used, but does not enter into our discussion here. Pretend for a minute that we try to work with the binary representation of 40, 24, 16, 12, 10, and 9. 101000 011000 010000 001100 001010 001001 102 ?Ve quickly see that we run into trouble in that the first three columns do not have the proper number of trailing zeros for Buneman`s procedure. Instead, lookin? back at table 1, we see that if we add 8 to every number, the angles we need are generated in precisely the correct order. Therefore using the binary representations of the working field only, we see that despite what the initial values are, one can still get the correct order with which to generate the weights. Since the position and order is determined by the working bits and the values are determined by the initial values, we can generalize this to any other processor, let's say processor 3 whose initial values are cos56 cos4O cos32 cos28 cos26 cos25 and use the exact same algorithm. In fact Buneman's algorithm need only be changed in a few spots, the initialization of the cosine table i3 different and the loop runs through fewer values. Algorithm 3.5.6 Multiprocessor On-line Cosine Generation = proces?or id, n --H 2?, P --H 2d */ n --H real hct+1), cCt+3) /* generate the portion of the cosine table */ /* depending on processor number id[?] */ arg = 87r 2w ptn = ?nc --H 2t--Hd--H1 id[?] for i = 1: t+ 1 arg = arg/2.O ang = arg + inc *pin c(i) = COS(ang) end cCt+2) = COS(ang-pin) write c(t+2) for K = 1 : 2t--Hd--H1 --H 1 I = K.AND.-K L = log2(I) /* c(t+2-L) contains the Kth cosine; use it and then replace it */ write cCt+2-L) I --H 2*1 + (2*LOR.K-I) c(t+2-L) = h(t+2-L)*Cc(t+1?L)+c(t+2?log2(I.AND.?I))) end Chapter 4 Comparison FFT Methods 4.1 Introduction M?i1tiA?? 1 P??1?;? of Two-Dimensional on the Hypercube ????ensiona ?r trans??rms, as in the single dimensional case, can also be broken into pieces that can be done in parallel. Here the possibilities are even richer than in the one-dimensional case. This is because multidimensional transforms can be done either as a sequence of separable one-dimensional transforms, or by directly splitting them into blocks of smaller multidimensional transforms as in the vector-radix methods [Rivard (1977), Harris et ai. (1977)]. We study only the case of the two-dimensional Fourier transform because the discussion and algorithmic methods can be extended directly to computing higher diiii?nsionai Fourier transforms. Lower dimensional Fourier transforms can be computed by way of multidi- mensional transforms through the Twiddle Factor Algorithm (see Chapter 1). Therefore large one-dimensional Fourier transform problems can be approached from the multi-dimensional point of view. This multidimensional mapping of one- dimensional DFTs underlies the derivation of prime-factor FFTs, and the very efficient Winograd Fourier transform algorithm [Nussbaumer (1982)]. Work on two-dimensional FFTs on distributed processors has so far been restricted to the row-column approach. The strips method partitions the two- dimensional array, or matrix, into rows, mapping block rows into processors. The transform of each row is then found, and the matrix is transposed before a second row transform pass is done on rows that previously had been columns. We follow McBryan and Van de Velde (1985) in terming this approach the Transpose-Split (TS) 2D-FFT. Another row-column method partitions the matrix into blocks of submatrices. 103 104 assigning one block per node. The hypercube is then viewed as a two-dimensional cross-product of smaller dimensional hypercubes with distributed FFTs performed along both the rows and the columns. No transposing of data is needed here. ?? term this method the Block (B) 2D-FFT method. We present two new methods of implementing 2D-FFTs. The first one is a row- column approach that partitions data into strips much like the Transpose-Split method. The difference is that no transpose is done between the horizontal and vertical steps. Instead, the horizontal FFTs are done locally inside each processor and the vertical FFTs are distributed. We call this the Local-Distributed (LD) 2D-FFT method. Finally we implement a partial Vector-Radix (VR) 2D-FFT on the hypercube. What this means is that the individual 2D-FFTs that are done locally inside the processors are row-column 2D-FFTs, however, the intermingling steps use the vector radix update scheme. The Transpose-Split 2D-FFT is favored by several, including John Gustafson (1987), because all FFT computations are performed locally. The only communica- tion that takes place occurs within the transpose step. Gustafson has implemented the Transpose-Split FFT on a 1024-node NCUBE machine which has the pleas- ant property that each node can perform up to 9 simultaneous communications. thereby allowing the use of almost all the links of the hypercube during the trans- pose stage. He can reduce communication time by a factor of d, the hypercube dimension. Therefore all that is needed to effectively implement this method is a fast efficient matrix transpose procedure. See Ho and Johnsson (1986) for an in-depth analysis of the hypercube matrix transpose problem. The Block method was implemented by Miles et al. (1987) on the Floating Point Systems T-Series hypercube. By considering the signal flow graph of the radix-2 FFT algorithm, we see that this implementation requires communication during both the vertical and horizontal passes. At each step where the butterfly computation is split between two processors, each node exchanges with its neighbor exactly half of its data points. Each processor computes the butterfly updates for the points it possesses after which it contains updates for half of its own points and half of the points belonging to its communicating partner. Another exchange is then necessary to repatriate these updates. Therefore two exchanges are required for one butteffly step. This may seem inefficient unless the communications and computations are overlapped during the distributed butterfly calculations. This is indeed possible on the T-Series since each node possesses a transputer that allows a processor to send data to its neighbor in the next butterfly step even before it has totally completed the present step. Hence, two communication stages can be overlapped in one computational step. This method of implementation is referred to as the Block method since the matrix is mapped by sub-blocks into 105 the hypercube such that the (i, j)th block is mapped into the node whose label is the binary representation of i concatenated with the binary representation of j. [See Miles et a!. (1987)]. The powerful vector boards on the T-Series allows this method to be used advantageously. The Local-Distributed method and the Vector-Radix-2 method are imp?e- mented on the Intel iPSC. The Intel `iPSC, unlike the NCUBE and T-Series machines, does not allow simultaneous communication. This is detrimental for the Vector-Radix-2 method as the intermingling stages involve total exchanges between four processors instead of two. A total exchange within a subcube of processors means that each processor in the subcube exchanges data with every other processor in the subcube. This particular property of Intel communica- tion also means that the full cross-bar interconnection scheme cannot be simu- lated efficiently, and with Intel iPSC/System 286 capabilities, a transpose takes 2d = 2 log2 P steps to perform, as each node can only do one send followed by one receive in one direction at a time. Another drawback is that computation and communication cannot be overlapped and thus distributed FFTs will exhibit blocking during the intermingling butterfly steps caused by one processor waiting for data from another. The Local-Distributed method does not require a trans- pose and does distributed FETs along only one direction instead of two (the Block method). Meanwhile the Vector-Radix FFT performs local 2D-FFTs followed by intermingling stages requiring the summation and multiplication of the local FFTs by "twiddle" factors. The Vector-Radix FFT has a lot of potential that is not reflected in our implementation on the Intel iPSC precisely because of com- munication inefficiencies. But we think it is useful to offer it as an alternative to the row-column approach because of its rich parallelism. 4.2 Two Dimensional FFT Algorithms Row-Column The Discrete Fourier Transform (DFT) of a vector x of length n is defined as y: Fnx where Fn is the matrix consisting of powers of the nth root of unity W? = e-2ri/n [Fnijk = wn5k The two-dimensional (2-D) discrete Fourier transform (DFT) of a two-dimensional array X E Cnlxn2 is defined as Y H Fn1XFnT2 106 This matrix notation clearly demonstrates the row-column or column-row method of computing the 2-D transform, since matrix-multiplication is associative. If the fast Fourier transform (FFT) is used to evaluate the 1-D FFTs along both the rows and the columns, the number of complex multiplications required is n2 log2 n for n = n1 = ?2 In addition, a matrix transposition algorithm is required. A one-dimensional FFT of a long vector of length n = ?2 can be computed in a 2-D "fashion" by viewing it as a DFT of an array of size n1 x n2 [Nussbaumer (1982)], that is, by writing x as an array Xn1xn2. This notation is defined in Chapter 7. Using the Twiddle Factor Algorithm (Chapter 7), we can compute the fl-point DFT of x by an ni-point FFT of the rows, a point-wise multiplication of x by the twiddle factors, followed by another n?-point FFT on the columns. The matrix of twiddle factors T is defined [Tijk = wn?k, j =0,...,n1--HI; k =0,...,?2 --H1. and * denotes the point-wise multiplication of two matrices. Hence the DFT y of x is another two-dimensional array given by ynT2xni = Fni[(Tn1xn2) * [xnixn2Fn2j] The row-column or column-row method can be used to compute the horizontal and vertical DFTs. Vector-Radix The Vector-Radix FFT is a direct decomposition of the two-dimensional DFT into sums of smaller two-dimensional DFTs multiplied by "twiddle factors", (the matrix ?). Here a 2-D DFT is recursively broken down into four 2-D DFTs until onlv trivial 2D-DFTs need to be evaluated. The number of complex multiplications is 3 now -n2 log2 n, 25% lower than the row-column method. [Rivard (1977), Harris (1977?)] Moreover, no matrix transpose routine is required. The recursive block structure of the DFT matrix Fn is used in two-dimensions to derive the method. The matrix X E C??? is segregated into four sets; one over those samples X(j, k) for which j and k are both even, one for which j is even and k is odd, one for which j is odd and k is even and one for which both j and k are odd. The diagram in Figure 4.1 shows schematically the decimations of data in both the one and two dimensional decompositions. Theorem 4.2.1 Let X E C??? with n = 2?, then the 2D vector-radix splitting of the 2D-DFT of X is FnXFnT = (Fnfln)(MnXMTn)(llTnFTn) 107 Figure 4.1: One and Two Dimensional Decimation Schemes = In/2 ?n/2 - ?11 ?12 - - 1n/2 In/2 n/2 ?n/2 - - X2l X22 - n/2 --H n/2 where A?12 =djag(1,w?,...,w?fl/2?1) and xli = Fn/2X(O :2: n--H2,O :2: fl?2)FnT12 x12 = Fn/2X(O :2:n--H2,1 :2:n?1)FnT?2 x21 = Fn/2X(1 :2: n--H1,O :2 :n?2)F?T12 x22 = Fn12X(1 : 2 : n --H 1,1 : 2 : n --H l)F?T/2 Proof From Theorem 1.8.1 we have Fnfln= IA Fn/2 0 I --HA - 0 Fq2 Applying this to both sides of MnXMt gives the required decomposition. *______ Theorem 4.2.1 is the basic two-dimensional Cooley-Tukey (CT2) splitting of the Vector-Radix method for computing the 2-D FFT. The complete algorithm is obtained by recursive application of this basic decomposition. los 1111111111111111 1111111111111111 1111111111111111 1111111111111111 Figure 4.2: 16 columns mapped into 4 processors. Algorithm 4.2.1 Vector-Radix 2-D FFT --H 2t X PnXPTn for q = 1 : t L 2q AL/2 ?djag(1,w?,...,w?L/2?1) end AL/2 ? - 1L/2 n/L - AL/2 1n/L ? 1L/2 1L/2 - 1L/2 --HAL/2 - 4.3 Data Structures and Mapping The Local-Distributed method as well as the Transpose-Split method are based on the natural mapping of block rows into the hypercube, where block row i is stored in processor z. The particular way we map these block rows into the hypercube is equivalent to McBryan and Van de Velde's (1987) distributed co'umn mapping Figure 4.2 shows the mapping of a matrix into 4 processors. Local FFTs are done along the rows, with the vectors oriented perpendicular to the direction of the FFT computations. The Vector-Radix-2 method and the Block method both partition the matrix into submatrices, placing one submatrix in each processor. The Block method involves distributed FFTs along both rows and columns. The Vector-Radix-2 109 000?00 000?0i 000?10 000?1i 001?00 001?01 001?10 001?1i 010?00 010?01 010?10 010?1i 011?00 011?0i 011?10 011?ii 100?00 100?0i 100?10 100?ii 101?00 101?0i 101?10 101?1i 110?00 110?01 110?10 110?11 111?00 111?01 111?10 111?11 Figure 4.3: 5-cube = 3-cube 0 2-cube. intermingling steps involve total exchanges among two-dimensional subcubes. The adjacency structure of the communication in both methods requires viewing a dimensional hypercube as a "cross-product" of a k-dimensional by an i-dimensional hypercube where d = k + 1. To illustrate what we mean, Figure 4.3 gives a low dimension example. Suppose we have a 5-cube viewed as a "cross-product" of a 3-cube and a 2-cube. Each node has a binary label A with 5 bits that can be construed as the concatenation of a binary label b1b2b3 of 3 bits with a binary label c?c2 of 2 bits, i.e. A = b1b2b3c1c2 Now we place each node so that in the x-direction cic2 go from O0 to 11; and in the y-direction b1?b3 go from 000 to 111. This produces a two-dimensional array of nodes possessing "powers of two" adjacency for both rows and columns. Notice that we have exactly the adjacency structure necessary for a two-dimensional radix-2 FFT of a matrix mapped into the processors in a natural fashion. This scheme extends to arbitrary dimensions. Suppose we have a d-dimensional hypercube where d = k + 1+ m. We can write the binary label of each node A=bi?...bkc1c2...cifif2...fm and map our nodes in three dimensions. Denoting concatenation by exponentia- tion, we order the ....... Im from 0m to 1m along the x-direction, the c1c2... Cj from Ol to 11 along the y-direction and the ..... . bk from 0k to 1k along the z-direction. Therefore if d = d1 + d2 + ... + dj ( j= dimension of data struc- ture), a j-dimensional array can be mapped into the d-dirnensional hypercube so that in each direction k = 1, 2,. ,j, the interconnections are exactly those of an dk -dimensional hypercube. l10 1111111111111111 1111111111111111 1111111111111111 1111111111111111 Figure 4.4: 16 columns mapped into 16 processors. A block matrix mapping maps submatrices into the hypercube so that block Aij goes into processor i?j (where denotes concatenation). Figure 4.4 shows the mapping of a matrix into a hypercube of 16 nodes. In the Block method, the FFT computations must include distributed butterfly computations in both directions. This method can obviously be extended to higher dimensions. For the Vector-Radix method, higher dimensional FFTs require total exchanges among higher dimensional subcubes. 4.4 Efficient Data Transposition on the Intel iPSC An extensive study on matrix transposition on hypercube architectures was done by Ching-Tien Ho and S. Lennart Johnsson (1986). The algorithm we use in the implementations comes from a recursive block transpose algorithm described in Johnsson (1985a, 1985b) and McBryan and Van de Velde (1985) and coded by Ho. Basically a matrix of dimension n mapped into P processors can be fully transposed in a time of order (n2/P) log2 P. The following recursion fully describes the transpose algorithm. Algorithm 4.4.1 Transpose(A): A A00 A01 A10 A11 A01 H A10 Transpose(A00) Transpose(A10) Transpose(A01) Transpose(A11) 111 p Figure 4.5: Block Recursive Transpose. Figure 4.5 describes the recursive block transpose for an 8-by-8 block matri? mapped into a hypercube of 8 nodes. Each node receives one block row of the matrix. And the arrows denote the scope of the blocks that are transposed at each step. Figure 4.6 concentrates on processor 2 and highlights in black the portions of the matrix that are to be exchanged at each step. Since the matrix is mapped in the distributed columns formation, each dark block consists of contiguous elements and can be sent as a single message. The initial step involves an exchange between processors whose node numbers differ at the highest order bit, or bit number d--H 1. The messages are of length L/2, where L is the length of the portion of our matr'x resident in each processor. In our case, L is the number of data-points in a block row multiplied by the number of bytes per data-point. In the next step, each processor now has 2 messages of length L/4 and these are exchanged between processors wh6se node numbers differ at the 2nd-highest order bit, (bit number d --H 2). Counting the steps from k = 0 to k = d --H 1, we have the property that at the kth step, there are 2k messages of length L/2k+1 exchanged between processors whose node numbers differ at the (d --H k --H 1)th bit. A simple procedure describes where each chunk goes at each step. First we define a function `Thb't k --H 0 ifthe kthbit of? isO --H 1 if the kthbit of? is 1 For example: ithbit(20, 2) = 1 since 20 = 101002. A!gorithm 4.4.2 Recursive Block Transpose = processor id; n --H 2t; m = n/P */ 112 Figure 4.6: Blocks Exchanged by Processor # 2. /* processor ? has: */ /* Aid[? = A(id[j] rn : (id[ji] + 1) m --H 1,0 : n --H for k = 0 : d --H 1 /* partition block row Aid[i] into 2k+1 submatrices of width n/2k+l */ /* with the submatrices numbered from I = ... 2k+1 --H 1 */ w & n/2k+1 for j = 0 : 2k+1 1 A5 & Ai???(:,I w : (j + 1). w --H 1) end /* /* /* send 2k submatrices to neighbor in (d --H k --H 1)st direction */ receive 2k submatrices from neighbor in (d --H k --H 1)st direction */ determine processor p to exchange with */ id[pj = flip(id[?],d --H k --H 1); p = id?'[id[p]] if ithbit(?, d --H k --H 1) = 0 then send(fA5:j isoddj,m w,p) recv((A5:j isoddj,rn w,p) end if if ithbit(?, d --H k --H 1) = 1 then send(fA5:j iseven),m w,p) recv((A5:j isevent,m w,p) end if end 113 Notice that as the algorithm progresses, there are more and more messages of smaller and smaller sizes to exchange. This can cause a tremendous amount of overhead in terms of message start-ups. Since the start-up time is, in general. several orders of magnitude slower than the data transfer rate, we would prefer to buffer the smaller chunks into the largest possible array and transmit a fewer messages. Of course there are constraints to this approach also, taking into con- sideration the largest message size Bm and also the "copy" rate of the processors. Ho and Johnsson (1986) present an optimized version of this algorithm for the Intel iPSC taking into account the communication startup time, the copying time for buffering, and other peculiarities of the Intel cube such as its allowing commu- nication over only one port at time. They have discovered that the time it takes to copy 4K bytes is about 37 milliseconds and that the time for copying 256 bytes is approximated the same as one communication start-up. Hence their algorithm sends blocks of length at least 64 floating-point numbers without buffering. We have chosen to implement their algorithm for complex numbers. One should also note that if P = 2d, then after d steps in the block recursive transpose algorithm, the remaining transpose operations are local to each processor. 4.5 Comments on Data Ordering, Vectorization, Mixed Radix Implementations, and Storage One-di 1 fr fi - - -inensiona tran??jrms 0 ?ng vectors can also be implemented in procedures similar to the methods described here by utilizing the twiddte factor algorithm where the sequence is initially mapped in column order by distributed columns, i.e. block rows of columns. DFTs are taken along the rows, followed by a point- wise multiplication of twiddle factors. Then the whole array is transposed and the "column" DFTs are again done on the rows. Aigorithm 4.5.1 Twiddle Factor FFT of X = X?1 Xn2 X XFn2 X?T*X x XT X XFn1 Recall that * denotes the point-wise multiplication of two matrices. Since the results x(k1, k2) are in transposed order, the resulting array X, which actually contains xT is in the correct order. Therefore no extra transpose is needed to bring the data back to its original order. A two-dimensional DFT proceeds in exactly the same manner except for the omission of the twiddle factor multiplication. However, the resulting array is in transposed order. 114 The usual radix-two in-place FFT algorithms require a data permutation either at the start of the procedure or at the end. This permutation is the well-known bit- reversing permutation. In the Transpose-Split method, the FFTs are performed locally so that this permutation is also done independently and locally. ?`e can also avoid bit-reversing by using Stockham FFTs locally, albeit this necessitates an extra array of workspace. In any case, no distributed bit-reversing operation is necessary. The Local-Distributed methods, require up to an extra d communi- cations [Swarztrauber (1986a)], on a cube which can communicate simultaneously in all directions, to perform distributed bit-reversal. This is also discussed by Van Loan (1987). The Block and the Vector-Radix FET results in a need for distributed bit-reversal along both the horizontal and vertical directions, or an ex- tra 2#d communications if Swarztrauber's method is utilized. On the Intel cube, autosort FFTs incur too much communication overhead to be efficient. However if the data is to be transformed and then inverse transformed, there would be no need to unscramble in the transform domain since there are algorithms which take bit-reversed data on input and return the inverse transform in natural order. If the data array is mapped within each processor so that its vector orientation is perpendicular to the direction of the FFT, one can effectively vectorize the butterfly computations. The vector length is equal to the number of multiple transforms and hence one would prefer to do the FFT of m transforms in parallel rather than do one transform after the other. See [Hockney and Jesshope (1981) and Swarztrauber (1986a)]. The basic operation here is a vector SAXPY where V CLX + y v, x, and y are vectors with CL a scalar. Since the direction of the second FFT pass is perpendicular to that of the first pass, a transpose is needed between the two stages. For the Transpose-Split method this means that after the block transpose, one must complete the transpose by pefforming transposes on all the submatrices locally. For the two distributed methods, a local transpose of the array resident in each node is required to retain the correct orientation of the vectors. The Floating Point Systems implementation of the Block method consists of mapping the array X in sub-blocks into the nodes of the hypercube. Since the T-Series consists of vector boards, local transposes are done on each submatrix to keep the correct orientation for vectorizing. Vector length is also an important issue when doing multiple transforms on vector processors. Suppose we have an N-by-N array and p processors. The length of the vectors in the direction perpendicular to the FFT computations should be as long as possible (up to the length of the vector register) so as to take full advantage of vector operations. The Transpose-Split method requires each processor to do N/p FFTs of length N simultaneously and thus has an effective 115 Table 4.1: Multiple Transforms and Vector Length. horizontal vertical TS LD B VR TS LD B VR m N N N N N N N N nt N N N N N N N N p p p # of FFTs = vector length ?length of portion in each processor vector length of N/p. The Block method has each processor responsible for Niffip FFTs of which only N/? elements of each FFT are processor-local. The vector length here is N/# and is ? times longer than the Transpose-Split method. The Local-Distributed method has the same characteristics of the Transpose-Split method during its local phase, vectors of length N/p; however in its distributed phase, the effective vector lengths are N. These observations are summarized in Table 4.1 On the Floating Point Systems T-Series hypercube, the vector registers are of length 128, and thus the Block method is favored here since its average vector length is the longest. For example, if p = 64, and N = 1024, the Block method would give vectors of length 128, whereas the Transpose-Split produces vectors of length 16. Since the vector boards are so fast (5060 Mfiops for the 2D-FFT) [Miles et at. (1987)J, one would like to "fill up" the registers to minimize the loading and storage of data. In fact, this was a compelling reason for chosing the Block method for the FPS T-Series implementation. A further consideration is whether the problem at hand requires vectors whose lengths are not a power of two. As the powers of two get large, they are spaced further and further apart. A partial differential equation problem may take place on a grid which is 100-by-lOO. Most of the distributed methods in the literature presume the use of a radix-two FFT on vector length --H 2k (some k an integer). Therefore a point brought out by McBryan is that the use of the Transpose-Split method permits one to use any FFT subroutine that accepts general n. Given that there is a good transpose algorithm, such a method is more portable and simple for the end-user to implement. However, the Local-Distributed Method which we have proposed here can be used to handle transforms of length n = 2k? This is because the twiddle factor algorithm can be adapted so that the local in-house FFTs are of length m and the distributed FFTs running the other direction are of length 2k, where a radix-two algorithm is used. This method which we call 116 the Mixed-Radix-2 will be discussed in Chapter 7. The Vector-Radix method is generalized by Harris et al. (1977) to include arbitrary radices and non-square arrays. Finally, storage requirements differ for these methods. In general, the dis- tributed methods require an extra N2 /p points of storage per processor if one wants to do only one exchange per distributed butterfly step. This is because all information must be sent to the other processor before the updating is done. If one is allowed to incur two exchanges per distributed butterfly step, as in the FPS implementation, there is no extra buffer space necessary for storage. The Vector-Radix-2 method can also be done with two exchanges per intermingling step. Here each processor sends a quarter of its data to each of three members belonging to a specific two-dimensional subcube as itself. It also receives three quarter pieces from each member. Updates are done and another total exchange is necessary to repatriate the data. This means that no extra buffer is needed. The Transpose-Split method can be done totally in place, if unbuffered transpos- ing is used. However, as Ho and Johnsson (1986) have noticed, this can cause the transpose time to grow exponentially in the number of dimensions. Therefore, one extra buffer is needed to collect the noncontigous blocks and ship them all out at once. For the Intel iPSC, this buffer is not too large (only 64 floating point numbers). 4.6 Implementation and Timings The Transpose-Split, Local-Distributed, Block and Vector-Radix methods for the two dimensional FFT are implemented on the Intel iPSC/D4MX hypercube run- ning XENIX R3.4, iPSC Release 3.1 with Exelan R3.3 networking software. The code was written using Ryan-McFarland FORTRAN. Vector boards are not avail- able so that all computations within a node are done serially. The Transpose-Split method (Algorithm 4.6.1) uses local FFTs in both the vertical and horizontal stages and a recursive block transpose routine. The transpose code used was a modified version of Ching-Tien Ho's with the only change being the removal of the aforementioned buffering. The following algorithms give a general overview of what the whole system of processors are doing. For example, the construct do(in paraflel, i = 0,... ,P --H 1) means that all P processors are asynchronously performing the same operation. The construct do(for processor p(i)) (in parallel) 117 means that processor p(i) is performing that portion of the computations and that all processors are asynchronously performing their respective portions of work. with the necessary communications and updating. Algorithm 4.6.1 Transpose-Split 2-D FFT /* Find the 2D-FFT of X by using P processors partition X into block rows Xj i --H 0 .. , P --H 1 map Xj into processor p(i) do (in parallel, i = Xj XiF?T /* GSi FFT */ end do call Recursive Block Transpose (x) do (in parallel, i=O,...,P--H1) X?: XiFt /* Gsi FFT */ end do The method of implementation for distributed FFTs is that of Chapter 3 where the sequences are mapped in natural order into the hypercube. The isomorphism of the radix-2 FFT signal flow graph to the hypercube interconnection scheme is used to determine the inter-node communication pattern. The FFT subroutine implemented is the Gentleman-Sande (GS1) algorithm. The Local-Distributed method is illustrated by Algorithm 4.6.2 and the Block method by Algorithm 4.6.3. Figure 4.7 shows the communications necessary during the horizontal and vertical FFT passes for the Block method. Algorithm 4.6.2 Local-Distributed 2D-FFT /* Find the ?D-FFT of X by using P processors partition X into block rows Xj, i = .... . , P --H 1 map Xj into processor p(i) /* Local FFTs */ do (in parallel, i = 0,... ,P --H 1) Xj XiFTn /* csi FFT */ end do /* Distributed FFTs */ do (for processor p(i)) (in parallel) Xj [FnX]i /* distributed CSJ FFT */ end 118 Algorithm 4.6.3 Block 2D-FFT /* Find the 2D-FFT of X by ttsing */ an array of ? x processors */ partition X into blocks of snb-matrices Xij i=O?--H1,5=O map Xjj into processor p(i ?j) /* Horizontal Distribtted FFTs */ do (for processor p(iJ)) (in parallel) Xij [XF?Tjij end /* Vertical Distributed FFTs */ do (for processor p(i j)) (in parallel) Xij [FnX]ij end The Vector-Radix method (Algorithm 4.6.4) is implemented only partially. In other words, the local independent 2D-DFTs are done by a conventional row- column 2D-FFT subroutine found in Press et al. (1986), and only the intermingling steps involve updating by the Vector-Radix method. Since our hypercube has only 16 nodes, only two such intermingling steps are done. In order to perform in- place computation and do away with the need for extra buffers, the Vector-Radix updating was done so that each processor in the corner of its two-dimensional subcube was responsible for all the updating of that particular corner of the data for its three partners as well as for itself. For example, the processor in the south- west corner of the two-dimensional subcube for that iteration will send its north- west corner to the processor above it in exchange for the south-west corner of that processor. It will also send its north-east corner to the processor diagonally across from it in exchange for the south-west corner of that processor. And it will send its south-east corner to the processor to its right in exchange for the south-east corner of that processor. This processor will update all the submatrices it receives and then do a reverse exchange with all its partners to get its own updated corners back. The intermingling steps for processor pOlOO is illustrated schematically ill in Figure 4.8. Algorithm 4.6.4 Vector-Radix 2D-FFT /* Find the 2D-FFT of X */ an array of ? x processors */ partition X into blocks of sub-matrices Xij i = 0,..., ?--H 1, 5 = 0,..., ?--H 1 map Xij into processor p(i?) 119 hi h2 vi v2 Figure 4.7: Block Method 2D-FFT Communications for a 16-node hypercube 120 Figure 4.8: Vector-Radix Exchanges for Processor p0i00 /* Local ?D-FFTs of Xij do (in parallel, i = 0,... ,?--H 1, j = 0,.. .,?--H 1) .FT Xij F?/?XiJ n/?/? /* 2D-FFT */ end do /* intermingling steps for Vector-Radix method */ forq=O,...,d Xjj [AqXATq]jj end Here Aq is defined as in Chapter 1, (L = A = 1L ? 1L/2 ?L/2 q n - 1L/2 ?L/2 where ?L/2 =diag(1,w?,... ,wLL/2?1). Computation time as well as communication time is displayed for all four meth- ods. A portion of the communication time is reflected in a processor blocking while awaiting data that it is to receive. The total execution time is shown in Table 4.2. Table 4.3 shows the computation time required by each method while Table 4.4 displays the communication time whereas Table 4.5 displays the amount of time a processor spends blocked. Times are given in milliseconds and range from the fastest processor to the slowest processor. 121 Table 4.2: Total Execution Time. size Total Execution Time dim TS LD B VR 16 x 16 1 155 145--H160 2 90 80--H95 90--H110 160 3 75--H125 55--H70 65--H85 --H 4 140--H145 105--H150 90--H110 125--H190 32 x 32 1 710 665--H715 2 385--H390 340--H395 365?30 465--H470 3 230--H275 180--H230 235--H290 --H 2 1700 1550--H1750 1675--H1895 1850--H1895 3 910--H945 785--H935 1045--H1230 --H 4 525--H570 415--H585 425--H555 685--H840 128x 128 1 15130 15025--H15230 2 7825--H7840 7165--H7930 78208630 7735--H7820 3 3965 3575?145 4835--H5480 --H 4 2155 1800--H2225 1895--H2350 2345--H2490 256 x 256 3 17650--H17665 16870--H18605 22220--H24660 --H 4 9040--H9215 8135--H9615 8695--H10370 9040--H9215 122 We also mention that our implementation of the Block method differs from Floating Point System's in that only one exchange is incurred during each dis- tributed butterfly step. Since the Intel cannot overlap communication and compu- tation, the added cost of performing two exchanges would degrade the performance of the distributed methods without giving any basis for comparison. The results show timings that are roughly within 10% of each other for the four different methods. However if we look at the break-down of communication versus computational time, some interesting differences surface. The computational times of the Transpose-Split method and the Vector-Radix method are the most load balanced. The Transpose-Split computations exhibit the most parallelism as all of the actual FFT steps are independent and done in parallel. For the Vector-Radix method the individual 2D-FFTs done locally within each processor are also done entirely independently and in parallel. As expected, the computational times in the distributed methods show a small amount of im- balance, with greater gaps between the faster and slowest processor than shown by the Transpose-Split method and the Vector-Radix method. While the Local- Distributed method stays about even with the Transpose-Split method as far as computation is concerned, the Block method on the average takes a bit longer. Since the same FFT algorithm was used for the three row-column methods, we con- jecture that this might be due to the fact that processors in the Block method get interrupted during computation at two stages, during both the vertical and hori- zontal FETs, whereas the Local-Distributed method only gets interrupted during its vertical FET. The Transpose-Split method is only interrupted once, between the horizontal and vertical stages. Of course the effects of these interruptions can be minimized if the processors are able to simultaneously communicate and compute. The computational time for the Vector-Radix method is on the average faster than any of the row-column methods, especially on large problems. Keep in mind that we have not taken advantage of the 25% reduction of multiplica- tions since our Vector-Radix implementation does not recurse all the way down to the 2 x 2 trivial 2-D transform. Instead, only the intermingling steps are done via the Vector-Radix splitting and the processor local 2D-FFTs are done by the regular row-column approach. Hence even at this limited level, we see that the Vector-Radix shows potential in speeding up computation. The analysis of the communication times show a different story. Here the three row-column methods exhibit roughly the same range of times for commu- nication. As already mentioned, our implementation of the Vector-Radix 2-cube total exchange is very primitive since the Intel iPSC cannot communicate across more than one link at the same time. Hence a huge amount of blocking is seen to be responsible for the slowness of the Vector-Radix communication. The blocking time for the three row-column methods is about the same. This demonstrates that 123 Table 4.3: Computational Time of 2D-FFT. size Computation Time dim TS LD B VR 16 x 16 1 150/10 140--H155 2 90/10 70--H85 80--H100 90--H100 3 55/5 40?5 50--H70 --H 4 40/5 15?40 25?0 30?5 32 x 32 1 705/25 645?95 2 375/15 320--H380 345M10 380--H385 3 210/10 155--H215 210--H270 --H 1 15015/285 14805--H15010 72357250 2 7575/155 6935--H7705 7565--H8405 3 3870/90 3400--H3975 4660--H5315 _ 4 2015/60 1680--H2075 1775--H2225 1895--H1915 256 x 256 3 17325/315 16150--H17875 21515--H23975 --H 4 8815/185 7650--H9145 8205--H9910 8300--H8320 total computation time and time for internal transpose 124 the load imbalance of the distributed methods is not really much of a problem as far as blocking between send and recv are concerned. One interesting point is that the Transpose-Split communication results actually show an `ncrease in time for the X = 16 x 16 problem with increasing number of processors. This is due primarily to the added complexity of having to send smaller and smaller messages or extra buffering costs, Even though we have implemented the unbuffered trans- pose, the results from the use of buffered transposing still show this increase, See Section 4.7. In the next section we shall see that while the distributed methods require 0(log2 P) start-ups for communication, the transpose method could possibly re- quire up to 0(P) startups if not done carefully. We also consider the effective vector length of the different methods and hypothesize what would happen if the node processors have vector boards. 4.7 Discussion and Model Models of computational and communicational complexity are often useful in giv- ing general guidelines to the benefits of various methods of implementation. Since FFT implementations are usually communication bound, we first consider the analysis of simply transferring data among the processors as specified by the re- cursive block transform procedure and the distributed methods. The vectorization of multiple transforms are dealt with next. Finally we give an estimate of the total time required. Suppose we have an n-by-n array and P --H 2d processors. Throughout this discussion we shall assume that P is an even power of two. Assuming that P divides n, each processor would have n2/P points. The recursive block transpose algorithm requires d steps where n2/2P points are exchanged per step. Meanwhile both distributed FFT algorithms have d steps involving trans-processor butterflies. Each step requires the exchange of n2/P points. One can see immediately that twice as much data points are exchanged at each step by the distributed methods, Notice however, that these points are all contiguous, so that there is no overhead of sending multiple messages or need to copy into a buffer. However, due to the algebraic structure of the butterflies, an extra buffer array of n2/P points is needed for each processor since it cannot overwrite its array until after the butterfly computation. The extra buffer is not required in the Floating Point System implementation, however an extra exchange per trans-processor butterfly step is needed. Three extra buffers would be needed for the Vector-Radix method if we did not use this trick of exchanging twice per intermingling step. Therefore in our implementation, we incur the cost of the extra exchange and hence no extra 125 Table 4.4: Communication Time Including Blocking Overhead. size Communication Time/Blocked Time dim TS LD B VR 16 x 16 1 5 5 --H 2 5--H10 5--H10 5--H10 60A35 3 15--H70 60?5 5--H10 4 95--H110 65--H120 55-60 80--H155 32 x 32 1 5 15 --H 2 10--H15 15--H20 15--H20 85--H90 3 25?5 15--H20 15 2 30--H35 60 50--H55 220--H265 3 25--H70 45--H85 40?5 --H 4 40--H90 35?0 30?0 230--H385 128 x 128 1 110 220 --H 2 245--H265 220 220--H255 495--H570 3 80--H95 165--H170 160--H165 --H 4 115--H145 105--H120 105--H120 445--H585 256 x 256 3 310--H335 650--H710 660?95 - 4 230400 445--H510 435-495 860--H1355 126 Table 4.5: Time Spent Blocked during Communication size Blocked Time dim TS LD B VR 16 x 16 1 0 0 --H 2 0--H5 0--H5 0--H5 10--H45 3 560 060 0--H5 --H 4 15--H85 0--H110 5--H50 15--H90 32 x 32 1 0 0 --H 2 0--H5 0--H5 5 3065 3 0--H50 0--H10 0--H5 --H 4 565 0--H15 0-60 30--H175 64 x 64 1 5 10 --H 2 5 0 10 10--H175 3 5--H20 5--H55 0--H15 - 4 560 0--H15 5--H10 115--H355 128>c 128 1 15--H40 35 --H 2 20--H155 30--H35 3560 60--H340 3 10--H25 20--H30 20--H25 --H 4 45--H220 65--H170 60--H135 95--H715 127 buffers are necessary. Let r be the data startup time or communication latency time, Bm the maxi- mum packet size that can be transferred at a time, and tc the per element transfer time. Denote time by To?peCrthatOtdon the time required for a certain operation by a certain method. The total communication overhead is measured by tc (number of elements to be sent) + T (number of start-ups) First let's look at unbuffered transpose communication. Ho and Johnsson [1986] show that the complexity for unbuffered communication is TTs --H fl ____ min(d, log2 ) --H ___ Comm d--H)tc + 2BmP B$P The complexity for startups is 0(P) and grows exponentially with the dimension d = log2 P of the cube. This can be seen easily where, ignoring Bm, the complexity becomes TcToSmm = dffip2tc + (d?2i)2T. i=O When 7? > Bm, we must take into account these extra start-ups and the number of start-ups is O(P + log2 P [n2I2PBml). Buffered communication makes sense only when n2/P remains small and the complexity is approximately 0(log2 P) start-ups growing linearly with cube dimen- sion. Here one must also take into account the extra time required for buffering as well as the fact that the effective buffer is small, so that on large problems the transpose is essentially unbuffered. The Local-Distributed method has communication complexity as does the Block method. TcL0D?? = 2d??tc + 2drmax(1, T?0?? =7?%D?m Here the complexity only grows linearly with the number of cube dimension?. however, there is twice as much data to transfer and when the problem gets large, we get a measurement proportional to [n2/BmPl times the start-up costs. Thus 128 the complexity of the distributed methods is 0(log2 P) start-ups when n2/P < Bm and 0(log2 P [n2/BmP?) when n2/P > Bm. The communication complexity of our in-place Vector-Radix FFT necessitates two exchanges per intermingling step. Recall that this scheme is similar to that of the Floating Point Systems implementation and also Walton's one-dimensional FFT implementation where portions of the matrix or vector are exchanged and then sent "home". Here each portion is of size -34n2/P and each processor commu- nicates with three processors. Hence the communication complexity is TcVornRm = 4log2 WP(3ffi1tc # 4P If the three exchanges per intermingling step can be done simultaneously, the startup term becomes just 4log2 ?r In summary, TctSmm 2(P --H 1)T ? (log2 P)?ptc TDmm `?mm rrCVomRm (2 log2 P) ?p2tc+T (2 log2 P) %2tc + (2 log2 P) ??tc Table 4.6 illustrates data transfer time (in milliseconds) without any computa- tion. Two different transposes are compared for the Transpose-Split method, the unbuffered and buffered methods described in Ho and Johnsson (1986). One can see that for small problems, the buffered recursive transpose and the distributed methods take about the same time communicating, even though the distributed methods send twice as much data. This is because the data lie in contiguous loca- tions and thus require only one message. By comparing buffered and unbuffered transposing we can find the cut-off point beyond which it makes no difference whether to buffer or not to buffer. Here we see that buffering just does not matter when the size is larger than 64 x 64. The communication times of the distributed methods are comparable to buffered transpose times until the problem size gets larger enough so that the maximum packet size Bm is reached. Here we see that for the 128 x 128 problem, communication for the distributed methods is roughly twice the time for transposing. Since the number of data points moved is also twice 129 that of transposing, this is consistent. For d = 4, the timings are about equal since the message packets are small enough so that they can be all transferred in one step. Finally as expected the Vector-Radix communication was the slowest be- cause of the inefficient implementation on the Intel iPSC. It is expected, however that given an efficient "total exchange" capability where processors can commu- nicate simultaneously, that the communication times should speed up making the Vector-Radix method viable and feasible. Next we model computational time for the three row-column methods on vector nodes. Hockney and Jesshope (1981) give a model for performing lvi independent transforms of length N by using the best serial algorithm and vectorizing the arithmetic. Let ot?1 be the per flop computation time, and N112 be the vector length required to achieve half of the asymptotic performance. Then Tm?Itft = 5?N(N1/2 + M) log2 N This model cannot be used for the Vector-Radix method since the "multiple trans- forms" are done simultaneously in two-dimensions. However we know that the Vector-Radix method, when recursed down to the trivial 2 x 2 2D-FFT has 25% less multiplications than the row-column approach. Using this model we see that for the Transpose-Split method each processor does n/P transforms of length n twice, hence TToSmp = 2?(5fl(N1/2 + pfl) log2 n) For the Block method, each processor is responsible for 1/?/?th of the work for n/? transforms of length n and this occurs twice, so `fcBomp=2a(5?(N1i2+ ?_)log2n) Finally for the Local-Distributed, each processor does n/P transforms of length n during the horizontal phase, and (1/P)th of the work for n transforms of length n during the vertical stage. TcLoDmp = 5Qn(N1/2 + pfl) log2 n + 5Qp?(Nl/2 + n) log2 n The difference in computation time among the three methods comes from the N112 term, with TvB --H (n/ vF"'P)N112 TTs --H nN1/2 and TvLD --H 1/2(n + fl/P)N1/2. Comparing coefficients, we see that TcBomp < 7%oDmp < TTs when P> 1 comp 130 Table 4.6: Communication Times without Intervening Computation. Analysis of Communication w/o Computation size TS LD B VR dim unbuff. buff. 16 x 16 1 5 5 5 --H 2 10 10 10 10--H15 50--H90 3 20?0 65--H70 10?5 15--H20 --H 4 105--H110 75--H140 70--H110 65--H70 10--H115 32 x 32 1 10 10 20 --H 2 15 10 20 20--H25 45--H80 3 20?5 20 15--H20 25?0 --H 4 45--H140 20--H25 20--H25 25--H75 15--H190 64 x 64 1 30 25--H30 55 --H 2 35--H70 30?5 60 60--H65 55--H175 3 35--H75 30--H75 50--H85 50--H90 --H 4 105--H110 35--H90 35--H85 45--H105 85--H285 128 x 128 1 115 110 220--H225 --H 2 110--H115 110--H115 220--H225 225--H230 110--H235 3 140 90--H95 170 175--H180 --H 4 80--H125 80--H125 120--H150 175--H180 95--H385 256 x 256 3 335--H340 330--H335 660690 670--H760 --H 4 235--H240 235--H345 445--H450 450--H455 295--H740 131 Hence in terms of vectorization, the Block method is at an advantage which in- creases with increasing parallelism as measured by N112. All of the vectorized FFT implementations require an internal transform of the data in order to set up the vectors in the correct orientation. Since each processor has n2/P points the overhead here is n2/P. tcop?. Our final model combines the communication and computational portions of the row-column methods in a straight-forward manner. TTs --H 2?(5fl(N1/2 + pfl) log2 n) + 2(P --H 1)r + (log2 P)?tc P n TB?2?5 ?(N112+ )log2 )+(9 n2 ? log2P)(?tc+T) TLD --H 5Qn(N1/2 + ?pfl) log2 n + 5a--H??(N112 + n) log2 ?pfl + (2 log2 P)(?2tc + r) P An analysis of the communication time shows that the coefficient of the tc term is of the same order, but that of the T term is clearly against the Transpose-Split method. Since the T term represents the latency time or startup time for each communication and is several orders of magnitude larger than tc, it is obvious that as P increases, the overhead for transpose communication will become significant. Of course buffering can reduce this overhead. However the minimum number of communications needed for recursive block transpose is still of order log2 P. Therefore, the communication needs of the distributed methods present a lower bound for the transpose communication. 4.8 Conclusion After consideration of the model and timing results we can draw the following conclusions. The distributed methods are hurt by the interruption of computation dur- ing the trans-processor butterfly stages, hence their performance should be enhanced on systems which interleave computation and communication. o+ The complexity of communication for distributed methods is log2 P, there- fore they could show promise on systems where P is large. In the presence of vector processor nodes, the Block method exhibits better vectorization than any of the other methods since its working vector length Is approximate ? times longer. 132 The Transpose-Split method is superior on systems which support an efficient transpose algorithm since its major deficiency is the asymptotic complexity of recursive block transpose communication growing exponentially to the number of dimensions. Presently on the Intel iPSC/System 286 where each node is connected by an Ethernet communication channel, the penalties for traversing several nodes to get a message across increases linearly with the distance. Therefore on the current system, Transpose-Split is likely to be superior when the number of processors is small. o+ The Vector-Radix method is promising both in the fact that it is highly parallel and also that it requires 25% fewer complex multiplications than the row-column approaches. However a d-dimensional Vector-Radix FFT requires total exchanges between processors in a d-dimensional subcube dur- ing each intermingling step. Therefore it is possible to obtain good results for this method only if the hypercube in question is able to perform efficient "total exchange" communications. Future generations of hypercubes will likely support more efficient transpose algorithms and routing hardware, making the Transpose-Split methodology more efficient. We have already seen the difference between transpose meth- ods in Table 4.6. The Local-Distributed method is a compromise between the first two meth- ods. In our implementation it is competitive with the Transpose-Split method and can be used when one wishes to avoid the transpose. Our implementation results show that on the Intel iPSC/System 286 without vector boards on the nodes, there is essentially no difference between the three row-column methods. Efficient implementation of the Vector-Radix method de- pends on efficient total exchange communication and is therefore this method is promising given its faster computational potential. The final analysis is that this problem is highly system dependent, and one should be aware of the advantages and disadvantages of these different methods in order to best utilize the parameters of a particular system. Our discussion of higher-dimensional transform methods naturally leads to in- sights on further areas of research. A higher-dimensional DFT can be approached from any combination of the following: row-column, transposing along any di- mension, or distributed block manipulations. For example, one way of looking at a 3-dimensional transform is as a 2-dimensional transform and a i-dimensional transform mapped in planes within nodes of a linear array of processors. Another way of viewing it is to map the data into solid chunks within a three-dimensional 133 array of processors. Finally, the use of the vector-radix approach is promising in the sense that no transposing is involved and maximal parallelism is exhib- ited during the independent stages. In the distributed stages, there is a total exchange between subcubes of the same dimension as the transform. While our implementation of this total exchange on the Intel is extremely primitive, we be- lieve that this method can be viable in systems containing optimal total exchange algorithms invisible to the user. Transposing distributed data in three dimensions is extremely awkward, hence we are especially looking at vector-radix approaches to 3-dimensional FFTs as an avenue for further work. Chapter 5 Parallel Sine Transform for the Hypercube 5.1 Introduction In some situations, the DFT of a sequence with some sort of symmetry is required. This sequence might be real, or complex even or odd. The DFT of such a sequence usually has special structure, e.g., the DFT of a real sequence is conjugate even. Hence both storage space and computational time can be saved by utilizing the redundancy of information presented by symmetry. The implementation of symmetric FFTs on distributed systems is nowhere as straightforward as that of the complex FFT. While most symmetric FFT algo- rithms utilize the complex FFT as the core computation, they are implemented with a series of pre- and post- processing steps designed to take advantage of as- sociated symmetries to minimize the operation count. In the conventional single processor setting the savings of both storage and computational effort is a rea- sonable goal to pursue. We present a new parallel sine transform algorithm that minimizes the amount of communication required to pefform its limited pre- and post- processing steps. An important detail in multiprocessor FFTs is the fact that the output of the Cooley-Tukey complex FFT algorithm is generated in bit-reversed order. In the single processor case this array can be unscrambled easily with index permuta- tions. However in the case of the hypercube, the bit-reverse unscrambling takes d parallel transmission steps, where d is the dimension of the particular hyper- cube. [Swarztrauber (1986a)] If the data is not unscrambled, post- processing operations for the sine transform require a significant amount of communication between processors. Since the complex FFT itself is extremely efficient and also efficiently parallelizable, it would be a shame if the cost of the pre- and post- pro- 134 135 cessing operations and communications were to outweigh the benefits of parallel computation. Our new sine transform algorithm attempts to keep pre- and post- processing communications to a minimum at the expense of a relatively small number of redundant calculations. We believe that with the high overhead associated with communication on the Intel hypercube, the total return of this algorithm will be superior to the savings of some computational effort at the expense of a severe tie-up in communications. To summarize this chapter, we first describe the symmetric transforms and their properties. The reflection and exchange permutations are introduced to facilitate the matrix description of the transform procedure. We also show how to work these two permutations into bit-reverse permuted data vectors. A new permutation, the rectr3ive exchange is introduced to manipulate the reflection permutation of sequences in bit-reversed form. The implementation of permutations on a distributed processor system is done by communication between processors as well as local permutation of the data vector. We describe each of the relevant permutations and how to efficiently implement them on multiprocessor systems. The parallel sine transform of this chapter algorithm pares the necessary per- mutation/communication needs down to a bare minimum without doubling the length of the core complex FFT computation. This is done by partially exploiting symmetry and at the same time performing some redundant computations in par- allel. Our procedure is practical and can be utilized on data that is both input in natural order as well as in bit-reversed order. This provides a complete forward and inverse sine transform package. We finish by discussing possible implementations of the existing symmetric transform algorithms on the hypercube and show how the communication com- plexity of these methods can potentially hinder performance. Symmetric trans- forms are a subset of real transforms and as such, we discuss procedures for finding the FFT of real sequences and their possible implementation on hypercubes via the complex FFTs of Chapter 3. We then consider the conventional sine and co- sine transform methods of Cooley et ai. (1970) and Swarztrauber (1982). These methods pre- process the input sequence into either conjugate even sequences or real sequences and then utilize efficient algorithms for finding real FFTs. The symmetric transform is then extracted from the real FFT by post- processing computations. We show that the associated permutations involved with these data manipulations, when compounded to that of the real FFTs and bit-reversal permutation result in unwieldy communication requirements. This high degree of communication traffic means that the savings in communication time consis- tent with exploiting symmetry is negated by the extra communication overhead incurred. 136 5.2 Symmetries and Properties of the DFT We start with some definitions of the different kinds of symmetries that a given sequence may have. These symmetries can be described with the permutation matrices T, the reflection matrix for even and odd sequences, and E, the exchange matrix for quarter-wave even and odd sequences. The exchange permutation En is described by = [en?i,en??,... ,ei,eo] If y = EnX, then y consists of the elements of x flipped upside-down. The reflec- tion permutation matrix Tn is defined as follows Tn = 1 0 0 En?i If y = TnX, then y is best described by fixing the top element, and flipping the rest of the n --H 1 elements upside down. The reflection permutation can also be described via a circular upshift permutation followed by the exchange permutation. Tn = EnRnT where RnT is the circular upshift matrix defined by xo xi Xl x2 RTn f2 = Xn?2 Xn?1 Xn?1 xo The symmetric properties can now be described in matrix-vector notation in a compact format. Table 5.1 shows the different symmetries that a sequence might have along with an eight-point example. The following lemmas provide a link between the symmetries T and E and the DFT matrix Fn. These relationships allow us to describe some of the symmetry properties that Fn possesses in a simple matrix-oriented language. Lemma 5.2.1 FnTn = TnFn = Fn 137 Table 5.1: Symmetric Sequences Property Definition Example even (E) x= Tx [xo,xl,x2,x3,x4,x3,x2,xi] odd (0) x= --HTx [O,xl,x2,x3,O,--Hx3,--Hx2,--Hxij quarter-wave even (QE) x= Ex [xo,xl,x2,x3,x3,x2,xi,xO] quarter-wave odd (QO) x = --HEx [xo, xl, X2, X3, --HX3, --HX2, --Hxl, --Hx0] conjugate even(CE) = T?x [xo,x1,x2,x3,x4,x3,x?2,x1] conjugate odd (CO) x= --HTx [O,x1,x2,x3,O,--H?x3,--Hx2,--Hx1j Lemma 5.2.2 FnE = ?nFn where ?=diag(1,w?,... ,W???1), W? = exp(--H2ri/n). The next two lemmas are due to Swarztrauber (1986b) and provide a method of representing the DFT of quarter-even (QE-symmetric) and quarter-odd (QO- symmetric) sequences by purely real or purely imaginary sequences. This allows the programmer to save on array storage space as a complex number occupies twice as much space as a real number. These lemmas also simplify the splitting equations for the cosine and sine transforms by eliminating the need for complex arithmetic. Lemma 5.2.3 The DFT y of a QE-symmetric sequence x satisfies y = where ?? =diag(1,w?,... ,W???1) (wn = exp(--H2iri/n)) and can be represented in terms of a purely real sequence y. Proof Swarztrauber (1986b). ______*______ Lemma 5.2.4 The DFT y of a QO-symmetric sequence x satisfies = ?- 1--H y where ?? = diag(1,w?,... ,W???1) and can be represented in terms of a purety imaginary sequence y. 138 Proof Swarztrauber (1986b). Swarztrauber (1986b) devises quences. These equations provide We state them here in our matrix *_____ splitting equations for various symmetric se- the basic recurrence for symmetric transforms. notation. Splitting Equations for Real Sequences Let z be a real sequence of length 2n and let zi z(0 : 2: n --H 1) Z2 . z(1 : 2 : n --H 1) Since z, z1 and Z2 are real sequences their transforms are conjugate even symmet- ric. Let y = F2?z, Yi = FnZ1, and Y2 = FnZ2. This means that only n + 1 points of y need to be computed. Suppose yi(O : n/2) and y2(O : n/2) are available. The radix-2 splitting algorithm gives and y(O: n/2) = yi(O: n/2) + ?2ny?(O: n/2) y(n : --H1: n/2 + 1) = Tn/2Yl(O : n/2 --H 1) + Tn12?2ny2(O : n/2 --H 1) where ?2n = (1,W2n,... ,w2???2?1) (W2n = exp2riI2n). Tn/2?2n = ?2n since e--Hi(n--Hk)2r/2n --H eik2r/2n so that the second equation becomes y(n : --H1 : n/2 + 1) = yi(O : n/2 --H 1) --H ?2nY2(O : n/2 --H 1) Splitting Equations for Even Sequences Let z be a real even sequence of length 2n with zi comprised of the even points of z and z2 the odd points. If y, Yi, and Y2 are their transforms, respectively, we have the following syrnmetries. Yi is also E-syrnmetric, Y2 is QE-symmetric, and y is E-symmetric. Both y and Yi are real and Y2 = ?2?n1Y2 where Y2 is also real. Substituting this into the splitting equations for real sequences gives y(O : n/2) = yi(O : n/2) + y?(O : n/2) y(n : --H1: n/2 + 1) = yi(O : n/2 --H 1) --H y?(O : n/2 --H 1) Since all quantities are real, only real arithmetic is needed. Y2 is computed instead of Y2 so as to avoid complex arithmetic. 139 Splitting Equations for Odd Sequences Now let z be a real odd sequence of length 2n and zi and z2 defined as before, Yi is also 0-symmetric and Y2 is QO-symmetric, where y and Yi are strictly imaginary. Y2 = ?2?n1Y2 and Y2 is strictly imaginary. Substituting the whole lot into the splitting equations for real sequences gives, y(1 : n/2) = yi(1 : n/2) + y"'2(1 n/2) --Hy(n --H 1: --H1: n/2 + 1) = yi(1 : n/2 --H 1) --H y"'2(1 : n/2 --H 1) All quantities are purely imaginary so that we can multiply through by i and avoid complex arithmetic. Block Structure of the DFT Matrix The DFT matrix Fn has an interesting block structure when we consider real and imaginary parts separately. First we define some special vectors and matrices to simplify the following discussion. Wm = (1,1,...,1)TE1?m Vm = (--H1,1,--Hi,... ,(?1)m)T E 1?? Let En be the exchange matrix of order n defined as E? = I?(:, n --H1:1). The n-by-n real cosine and sine matrices are defined as cos(?/n + 1) ... cos(n?/n + 1) - - cos(nr/n + 1) ... cos(n2?/n + 1) and sin(?/n + 1) ... sin(n?/n + 1) - - sin(nir/n + 1) .. sin(n2ir/n + 1) Theorem 5.2.1 (Matrix Stritcture of Fn): Let n = 2m. Ifw = Wm?1, V = Vm?1, E = E??1, 0 = 0m-1, and S = Sm?i, then F2m = Re[F2mj - ilm[F2m] 1 1 Re[F2mJ = m-1 w 1 1 m-1 w m--H1 1 m--H1 wT 1 wT O v CE vT (?1)m vTE EC Ev ECE where 140 Im[F2rnj = rn -1 rn-1 Proof See Van Loan (1987). 1 m--H1 1 m--H1 o 0 0 0 o s o --HSE o 0 0 0 o -ES 0 ESE * ______ This matrix structure proves handy in the development of symmetric FFT algo- rithms and the description of cosine and sine transforms. 5.3 Permutations and Communication on Dis- tributed Processors A permutation of data mapped in a distributed fashion in separate nodes of a parallel computer requires communication between the processors. Sometimes permutations are unavoidable and on loosely coupled systems they can present quite a problem. Several permutations come into play when implementing real and symmetric FFT procedures. We have already seen the exchange permutation E, the bit- reversal permutation P, and the reflection permutation T. Because distributed FFTs are more efficiently implemented without distributed bit-reversal, the data in the transform domain is usually permuted by P. Hence a facility must be devel- oped for working with the combination of the exchange and reflection permutations with the bit-reversal permutation. To motivate, we first present an easy example concerning the exchange per- mutation. Suppose that during the course of our computations we are required to find x+Ex where x is distributed among processors in the consecutive order. This can be done by inter-processor communication that "flips" x upside down. We will suppress. for the moment, how this is done. Now suppose we really had Px available and wanted to find Px + P(Ex) What would be an easy way to implement the communications necessary? It would be extremely inefficient to first unscramble x Px, then find x Ex and finally scramble by P again. Since we have Px available, we would like to work 141 directly with Px. Luckily in the case of the exchange matrix E, PE = EP and therefore we find that Px + E(Px) can be implemented by a simple exchange permutation on the vector Px. Lemma 5.3.1 If n --H 2?, then PnEn = EnPn Proof We give an intuitive proof for the case n --H 2?. Recall from Chapter 1 that an exchange permutation is defined on the binary representation of the index j = (bt?1bt?2 . . b1b0)2 by complementing all of j's bits. Hence E(j) =(bt?1bt?2...?1?b0)2 Meanwhile the bit?reversal permutation is defined on the binary representation of the index j = (bt?1bt?2... b1b0)2 by reversing all of j's bits, i.e. p(j) = (b0b1.. bi?2bt?1)2 The commutativity comes from observing that if matters not a bit whether the bits of j are first complemented and then reversed in order, or vice-versa. The same observation works for the general radix n --H pt by defining complementation as the additive inverse, i.e. b + b = 0 mod p. ______*_____ We now turn to the reflection permutation T and show that things are not so easy here primarily because PT ? TP. Now supposing we want to find Px + P(Tx) and we only had Px to work with. What we now need is a permutation G such that PT = GP. It turns out that a special permutation G, that we call Recursive Exchange satisfies this criterion. A look at the structure of Gn n --H 2? shows that it is the direct sum of a sequence of smaller exchange matrices. Definition 5.3.1 Let Gn be a per'n?tation matrix which is a direct stm of smaller exchange matrices n --H E1 E1 E2 E4 E?12 142 G8 = can be defined recur?ivety as Example n/2 - 0 0 1 1 2 3 3 2 4 7 5 6 6 5 7 4 The following lemma provides a decomposition of the reflection matrix Tn into a direct sum of two permutation matrices of half the size coupled with a perfect shuffle of both the rows and columns. Lemma 5.3.2 n/2 Proof We apply the permutations on the right hand side to the vector (0 n --H 1)T 1. Mn(O : n --H 1)T --H (0: 2 : n --H 2)T (1: 2 : n --H 1)T - 2. Tn/2 0 - 2)T - - (0, n --H 2: --H2: 2)T 0 Eq2 - ((01 2 ?fl 1)T - = - --H 1: --H2: 1)T - 3. H? (0(n?fl??12: :?f?: 2)T - 1)T - 0 n--H1 n--H2 n--H3 2 The result of applying the permutations on the right hand side is exactly Tn(0 n --H 1)T 143 ______*______ Now we are ready to prove a result tieing together the reflection matrix Tn and the recursive exchange matrix Gn through the bit-reversal permutation Pn. Recall that Pn is always involved somewhere in distributed FFT procedures where unscrambling is not done. This next lemma is of great use in the analysis of bit- reverse permutation data that also exhibit symmetry. Lemma 5.3.3 PnTn = GnP? Proof We use induction. This is obviously true for n = 2. Now Tn = PnG?P? Hn Pn/2 0 G?12 0 ?n/2 o P?12 0 En/2 - - 0 P?12 --H Hn T?12 ? Mn o En/2 using the induction assumption Pn/2Tn12 = Gn12 Pn12 and the fact that EP = PE. Armed with this fact, we can find by the equivalent sum Px+PTx Px + G(Px) ______*______ Permutations on distributed vectors are equivalent to data movement between and within the processors of a loosely coupled system. In fact after having estab- lished the equivalency of the two concepts, we can use the permutation matrix as a shorthand for a sequence of sometimes complicated communication specifications. In this section we assume that the nodes of a hypercube are mapped in the Binary Reflected Grey Code (BRGC) order. The input vector is mapped in the tw?track fashion of Chapter 3. The following theorems tell how data mapped in the Tw?Track method of Chapter 3 is rearranged inside the processors when permuted by the various permutations. Keeping in mind that our main goal is to provide a setting where the corresponding components of x and let's say Tx are in the same processor so that they can be combined. For example, the Two-Track 144 o 8 0 4 0 2 0 1 19 1 5 1 3 23 2 10 2 6 4 6 4 5 3 11 3 7 5 7 6 7 4 12 8 12 8 10 8 9 5 13 9 13 9 11 10 11 6 14 10 14 12 14 12 13 7 15 11 15 13 15 14 15 H Figure 5.1: Permutation Pattern of Two-Track FFT mapping matches each element with its corresponding element that is half of the vector length away. Hence the computation x(n/2 : n --H 1) ? + - x(O : n/2 --H 1) - can be done locally without communication. Theorem 5.3.1 The permutation resttlting from the implementation of the FFT of the Two-Track method is equivalent to MnPn. (See Figure 5.1). Proof The in-place FFT results in a sequence permuted by Pn while Theo- rem 1.6.6 tells us that repeatedly pefforming super-butterfly permutations giving us Mn. Algorithm 5.3.1 Permutation Mn processors mapped in BRCC order */ = processor id; n --H 2t; p --H 2d; m = (n/P)/2 */ /* processor ? holds: */ /*z(1) --Hx(m-id[?]:m (id[i]+1)--H1) */ /*z(2)--Hx(--H?2+m-id[?: n2+m.(id[ii]+1)?1) */ for L = 1: d v = d --H L; id[p] = fhp(id[?], v) p = id?'(id[p]) send(z(2), m, p); recv(z('), m, p) end z(1) - - z(2) - H M2m z(2) *______ 145 o 8 0 15 1 9 1 14 2 10 2 13 3 11 3 12 4 12 4 11 5 13 5 10 6 14 6 9 7 15 7 8 Figure 5.2: Permutation Pattern of Two-Track En Theorem 5.3.2 Exchange. Suppose x is mapped into the hypercube by the Two- Track mapping. Then getting matching components of x and Ex into the same processor is equivalent t9 performing the exchange permutation on the second half of x only, i.e. applying In/2 0 ? (See Figure 5.?). 0 E?12 Proof The two-track mapping means that x(j) and x(j + ?2) are in the same processor. To get x and Ex together we want x(j) and x(n --H j --H 1) together. Therefore we need to replace But doing this is precisely the same as applying x(j+?)?x(n?j?1), n 2 In/2 0 0 En/2 The permutation pattern for an example n = 16 is given in Figure 5.2. Algorithm 5.3.2 Permutation En processors mapped in BRCC order */ = processor i&, n --H 2t; p --H 2d; rn = (n/P)/2 */ processor ? holds: */ /*z(') x(m.id[?]:m.(id[?j+1)--H1) */ /* z(2) --H x(?2 +m.id[?: --H?2+m.(id[?+1)--H 1) */ v = dl; p= flip(?,v) send(z(2), m, p); recv(z(2), m, p) z(2) Emz(2? *______ 146 o 8 0 8 1 9 1 15 2 10 2 14 3 11 3 13 4 12 4 12 5 13 5 11 6 14 6 10 7 15 7 9 Figure 5.3: Permutation Pattern of Two-Track Tn Theorem 5.3.3 Reflection. Suppose x is mapped into the hypercube 6y the two-track mapping. Then getting matching components of x and Tx into the same processor is equivajent to performing the reflection permutation on the second half of x only, i. e. applying `n/2 0 - 0 Tn12 - x. (See Figure 5.3). Proof In the two-track setting x(j) and x(j + ??2) are in the same processor. To get x and Tx together we want x(j) and x(n --H j) together. Therefore we need to replace x(j+--H?2)?x(n--Hj), j=Q...,--H?--H1 2 Doing this is precisely equivalent to applying `n/2 0 x 0 Tn/2 - Algorithm 5.3.3 Permutation Tn processors mapped in BRCC order */ 1 = processor id; n = 2? P --H 2d; m = (n/P)/2 */ processor ? holds: */ /*z(1) --Hx(m.id[?:m.(id[?]+1)?1) */ /*z(2) --Hx(--H?2+m.id[?: ?n2+m.(id[!J+1)?1) */ /* send top element to previous processor */ /* upshifl second column */ p = id?'[(id[?] --H 1) inod P] v = id?'[(id[?] + 1) mod P] ______*_____ 147 send(z(2)(O), 1, p); recv(t, 1, v) z(2? RmTz(2) z(2)(rn --H 1) t call exchange(z(2)) Theorem 5.3.4 Recursive Ezchange. Suppose x is mapped into the hypercube by the two-track mapping. Then getting matching components of x and Gx in the same processor is equivalent to first performing Mn on the sequence and then applying Gn/2 to the second half of x only. In other words, mapping In/2 0 0 Gn/2 - Mnx into the processors in the two-track manner. (See Figure 5.?). Proof In the two-track mapping, x(j) and x(j + n2) are in the same processor. Having matching components of x and Gx in the same processor means that x(j) and x(2['og??1 --H j mod ?Uog2iJ --H 1), j = 0,.. --H 1 are directly matched. The effect of Mn is x(j) 2 x(2j) x(2j+1) n J ?...., 2 --H 1 We want x(2j) together with x(2F'o?22>1 --H 2j mod 2[1og22?j --H 1) x(2j + 1) x(2??o?22?1 --H 2j mod 2Uog22jJ --H 1) This is equivalent to applying Gn/2 to x(j') where j' 2j + 1 is the change of variables. or *_____ 148 o 8 0 1 0 1 19 2 3 23 2 10 4 5 4 7 3 11 6 7 6 5 4 12 8 9 8 15 5 13 10 11 10 13 6 14 12 13 12 11 7 15 14 15 14 9 Figure 5.4: Permutation Pattern of Two-Track Gn Algoritbm 5.3.4 Permutation Gn /* processors mapped in BRCC order */ = processor id; n --H 2t; p --H 2?; m = (n/P)/2 */ processor ? holds: */ /* z(1) --H x(m . id[? : m (?d[?j + 1) --H 1) */ /*z(2) --Hx(--H?2+rn.id[?: --H?2+m.(id[?]+1)--H1) */ /* appL? M? only if Two- Track FFT was not done */ if(Two- Track FFT not done) then call Inverse Shuffle(z(1), z(2)) end if if (id[ji] # O.and.id?]#1) then id[p] 2?og2id[??id[i] mod ?Uog? id[?jJ --H 1) p = id?1(id[pj) send(z(2), m, p); recv(z(2),m, p) z(2) Emz(2) end if if (id[?] = 0) then z(2) Gmz(2) end if if (id[?] = 1) then z(2) Emz(2) end if Armed with this fact we can find Px+PTx Px+ G(Px). by the equivalent sum 149 Now that we have related permutations to hypercube communication we can simplify our ensueing discussion of algorithms. Using the "language" of a permu- tation matrix, let us say T?, we can talk about data movements in distributed systems by referring to Algorithm 5.3.3. This means that in the course of our discussion of symmetric transforms, we can simply say call Permutatibn(T?) and know that the communications detail in Algorithm 5.3.3 are performed. 5.4 Development of the Parallel Sine Transform In this section we present our new sine transform algorithm designed to be im- plementable on the hypercube with a simplified pre- and post- processing scheme. In fact, all of the post-processing is local to a processor and we require only one upshift (Rn) communication (or shift-vector operation of McBryan and Van de Velde (1987)) for the pre-processing, i.e. every processor sends one real number to its preceeding node and receives one real number from its succeeding node. In the case that the original sequence is generated by a formula instead of entered in, we can circumvent this communication by having each processor generate one extra real number. This method results in the sine transform being laid out so that just two exchanges are required for pre-processing the inverse sine transform. This simplified procedure is accomplished by doing some redundant computations and computing a complex FFT of length n, instead of a real FFT of length n as was the case for the conventional methods. Since the sine transform of an (n--H1)-point sequence is derived from the imag- inary part of the complex transform of the extended odd (2n)-point sequence of real numbers, we first consider the use of an fl-point complex FFT to find the DFT of this extended sequence. Note that the analogous cosine transform method can be derived by using the real part of the complex transform of the extended even (2n)-point sequence of real numbers. Let z be the real odd extension for our original (n --H 1)-point sequence x, 0 x(1 : 0n --H1) --HEx(1 : n --H 1) We can split this sequence into two sequences, one comprised of the even points and the other the odd points. Let zi z(O : 2 2n --H 1) 150 z2 z(1 :2: 2n --H 1) Swarztrauber (1986) has shown that these sequences exhibit symmetry also, that is TnZ1 = z(2n --H 2: --H2:0) = --Hz(0 : 2 : 2n --H 2) = --Hzi so that z is an odd sequence, and Enz2 = z(2n --H 1:--Hi: 1) = --Hz(1 : 2: 2n --H 1) = --Hz2 and thus z2 is quarter-wave odd. The result is that we have two real sequences each with special symmetry. We proceed in a similar vein as Procedure 2 except that we have two additional symmetries to exploit. Form the complex sequence S = zi + iZ2 and let (5.4-1) W = FnS = Yi + iy2 Yi and Y2 possess properties in addition to the fact that they are conjugate even Yi = Ty1 Yi = -Y1 Y2 = T?2 Y2 = -??1Y2 Y2 = ?2n1Y2 Y2 = -? Recall that Y2 is pure imaginary. Equation 5.4-1 gives (5.4-2) W = Yi + ??2??1?2 and w = y1 --H ??2nY2 (5.4-3) = --HYi + ??2n?2 Solving equations 5.4-2 and 5.4-3 gives Y2 = $diag k)/n 151 and Yi = ilm[w] --H idiag (tan k?/n) Re[w]. Because of the Yi and Y2 are conjugate even, we only have to solve for the first terms 0,.. , n/2. Finally w(n/2) = yi(n/2) + ?exp(ir/2)y2(n/2) --H y2(n/2) and yi(n/2) = 0 y2(n/2) = --Hw(n/2) = --Hilm[w(n/2)j The splitting equation for odd sequences gives y(1 : n/2) = y(n --H 1 : --H1 n/2 + 1) = --H21[yi(1 : n/2) + y2(1 : --H2[Y2(1 : n/2) --H yi(1 : From the definition of a sine transform, we know that b = 2iy so that we can multiply through by i and obtain the following equations b(1 : n/2) = b(n --H 1 : --H1: n/2 + 1) = yi(1 : n/2) + y2(1 : n/2) y2(1 : n/2 --H 1) --H y2(1 : n/2 --H 1) And we never had to compute y in the first place. None of our arrays need to be larger than n and the extended sequence z was used for illustration and derivation purposes only. In fact, Yi and Y2 can be real arrays of length n/2 + 1. Procedure 1. Given: A real sequence x of length n. Find: The sine transform Pn b by the parallelizable sine transform algorithm. o+ Form the complex sequence s as follows. __ x(O : 2 : n --H 1) + ix(1 : 2 : m --H 1) s `. --HTq2x(O 2 : n --H 1) --H ?En/2X(1 : 2 : n --H 1) This requires one linear-shift communication and one exchange permu- tation En/2 o+ Find PnW = PnFnS. This requires the communication pattern for dis- tributed FFTs using the Two-Track method. 152 o+ Set Pny1 = --HIm[P?w] + Pndiag(tan ?)Re[Pnwi Pn?2 = Re[Pnw]d?a?(1/cos?) with Pnyi(n/2) = 0 and PnY2(n/2) = Im[P?w(n/2)]. This step is local to a processor. 0 Let Local to processor. Pnb = Pny1 + Pn?2 Algorithm 5.4.1 Forward Sine Transform = processor id; n --H 2t,o+ p --H 2d; m = (n/P)/2 */ /* Initially, processor ? holds: */ /*x(2m.id[?:2m.(id[?+1)--H1) */ /* Form the two-tracks, ?(1) and ?(2) */ x(O : 2 : m --H 1) + ix(1 : 2 : rn --H 1) /*form n/2s(2) R?T12x(o:2: n --H 1) --H x(1 :2: n --H 1) */ p --H id?'[(id[? --H 1) mod P] v --H id?'[(id[?] + 1) mod P] send(Re[s(')(O)], 1, p); recv(temp, 1, v) s(2)(m --H 1) --Htemp --H ?(Im[s(1)(m --H 1)]) s(2)(o : rn --H 2) --HRe[s(1)(i : rn --H 1)] --H i(Irn[s(2)(O : rn --H 2)]) call Exchange (s(2)) /* code (5.3) */ w PnFns */ call Two-Track FFT(s(1),s(2)) /* code (5.2) */ /* post-processing code (5.5) */ /* k r?ns throigh the indices */ /*v(2m.id[?j :2m.(id[?]+1)--H1) */ /* where v = Pn(O: n --H 1) */ yi .` --HIm[w] + diag(tan L?r)Re[w] Y2 Re[w]diag(1/cos L:) /* b is the sine transform of x Pnb Yi + Y2 153 s Notice that we have chosen to compute all n points of Yi and Y2 when only n/2 +1 are needed. This is where the redundant computations come in. Therefore everything after the complex FET is local to a processor with the introduction of n extra multiplications and n/2 extra additions. However these extra operations are not noticeable in the parallel setting since the processors are working indepen- dently. If the redundant computations were not introduced, some of the processors would be sitting idle anyhow. Therefore we have achieved load balancing as well as saved 0(d) communication steps for a total of 2d --H 1 exchanges and one linear- shift. Notice the 2d --H 1 exchanges are solely used for the computation of the complex FFT on a hypercube mapped on a BRGC, therefore the pre-processing resulted in only one (optional) communication. b(O : 2 : n --H 1) + ib(1 : 2 : n --H 1) Tn/2b(O : 2 : n --H 1) --H ?En/2b(1 : 2 : n --H 1) 5.5 The Sine Transform of a Sequence in Bit- Reversed Order Even though the sine transform is its own inverse (scaled by a constant) our output from the forward sine transform is in bit-reversed order. Therefore we require an inverse sine transform which takes data originally in bit-reversed order and finishes with it in the natural order. We have Pnb. Recall that Pnb= Pn12b(O : 2 : n --H 1) Pn/2b(1 : 2 : n --H 1) We would like to construct MnPns so that we can use the inverse Two-Track method to obtain w .` FnS. Recall that From Theorem 1.6.7 we have Pn = 11n(12 ? Pn/2) hence Pns ? Hn Pn/2b(O : 2 : n --H 1) + ?Pn12b(1 : 2 : n --H 1) - Pn12Tn12b(O : 2 : n --H 1) --H ?Pn/2En/2b(1 : 2 : n --H 1) Since the permutations PT=GP and PE=EP P?12b(O: 2: n --H 1) + iPn12b(1 : 2: n --H 1) Gn12Pn12b(O : 2 : n --H 1) --H ?En/2Pn/2b(1 : 2 : n --H 1) - we can form MnPns 154 ??rking On MnPns the the Inverse Two-Track algorithm we get W H FnS The post-processing continues locally. Therefore the rest of Procedure la is identical to that of Procedure 1 except that we have sequences in the natural order. Since the sine transforms is its own inverse we have just outlined two procedures which can take both the sine transform and its inverse on data mapped either in the natural order or bit-reversed order. Procedure la. Given: A real sequence Pnb of length n. Find: The sine transform x by the parallelizable sine transform algorithm. o+ Form the complex sequence MnPns Pn/2b(O : 2 : n --H 1) + iPn12b(1 : 2: n --H 1) G?12Pn/2b(O : 2 : n --H 1) --H ?E?12 Pn/2b(1 : 2 : n --H 1) The communication patterns for G?12 and En/2 are needed. o+ Find w = FnPn(PnS). Need communication pattern for distributed FFTs via the Inverse Two-Track method. o+ Set kir xi = --HIm[w] + diag(tan T)Re[wi x2 = Re[w]diag(1/ cos ffi$) No communication is needed here. X = xi + x2 No communication is needed here. A!gorithm 5.5.1 Inverse Sine Transform = processor id; n = 2t; p --H 2d; m = (n/P)/2 */ /* Initially, processor ? holds: */ /*g(i) --Hg(rn.id[?] :m.(id[?j+1)--Hi) */ 155 /*g(2)?g(n2+m.id[?:?n2+m.(jd[?+1)??) */ /* where g = Pnb */ /* b is the sine transform of x */ /* Form the two-tracks, s(1) and ?(2) */ g(l) + ig(2) call Exchange(g(2)) call Recursive Exchange(g(2)) g(2) ?g(l) --H ig(2) /* w FnS */ call Two-Track FFT(s('),s(2)) xi --Hrrn[w] + diag(tan L:)Re[w] x2 Re[Wjdiag(1/cosL:) X .` xi + X2 Procedures 1 and 1a are efficiently parallelizable because in both cases the pre- processing requires minimal communication and the post- processing is entirely local. 5.6 Implementation Results on the Intel Hyper- cube The complete procedure was implemented on a four-dimensional Intel iPSC hy- percube with extended memory in the nodes. Since our goal is to compute a forward sine transform, scale it and then do an inverse sine transform, our timings are for both the forward and inverse transform done back to back. We have also implemented the Cooley et ai. (1970) algorithm on a single node and obtained a time of 121685 milliseconds for a sine and inverse transform of length 16384. The results show almost linear speedup and reflect the efficiency with which the complex FFT can be parallelized. The speedup* against the sequential algorithm is also very good. Recall that the Cooley et al. sequential algorithm does a complex FFT of length N/2 or (8192 for our example). Therefore it should be roughly twice as fast as our algorithm in the sequential case. But our algorithm using only one processor is only about 4% slower than the conventional algorithm. So we see that the layers of pre- and post- processing are quite wasteful even in the totally sequential case. One may wonder how this method would compare to doing a 2N-point complex transform of the extended odd sequence (the only other feasible method on the hypercube). Results for the N = 32768 sequence are given below. Comparing 156 Table 5.2: Timing Results for the Parallel Sine Transform Time (milliseconds) for Sine and Inverse Transform N = 16384 # proc time speedup speedup* 1 127515 1.00 0.95 2 67515?7535 1.89 1.80 4 36870--H36915 3.45 3.30 8 19315--H19375 6.58 6.28 16 10430--H10535 12.10 11.55 (* speedup against conventional algorithm) Table 5.3: Timing Results for Conventional Complex Transform Time (milliseconds) for Complex Forward and Inverse Transform N = 32768 # proc time speedup 1 197840 1.00 2 97785--H97825 2.02 4 52450--H52605 3.76 8 28800--H28865 6.85 16 15055--H15585 12.69 157 results shows that our method is faster than the direct computation of the 2N- point sequence by 45%, 43%, 49% and 48% for hypercubes of dimension 1, 2, 3, and 4 respectively. Sequentially, it is also 55% faster. 5.7 The DFT and IDFT of Real Sequences Real sequences exhibit special symmetry in the transform domain, that is their transforms are conjugate even. Thus an efficient algorithm computes the FFT of a real sequence of a fixed length n with the same amount of effort in terms of storage space and operations count as that of the complex FFT of length n/2. The economies come from the symmetry. This is indeed the case in the procedures described in Cooley et at. (1970). We repeat some of the procedures in matrix notation. In each case we first derive the procedure and then give a high-level description, mentioning communication patterns with our permutation matrix language. We also give the equivalent procedure in bit-reversed order and output in consecutive order for completeness. Procedure 2. The DFT of Two Real Sequences done as one Suppose xl and X2 are two real sequences. We wish to find Yi = FnX1 and Y2 = FnX2. Since the transforms of real sequences possess symmetry, and xi and x2 have zero imaginary parts, we can exploit this by forming the complex sequence z = x1 + ix2 and taking the complex FFT of z. Let w = FnZ = Yi + ?Y2. Yi and Y2 are conjugate even, so that we have w = y1+jy2 T'w = Ty?-iT? = Yi --H iY2 Yi = Y2 = -21[Tv+ w] wJ 2 Procedure 2. Given: Two real sequences xi and x2. Find: the DFT Yi and Y2 respectively of the two real sequences by using a radix-two in-place complex FFT which computes PnFnz (either OTi or GS1). hence 158 o+ Form the complex vector z x1 + ix2. No commumcation is needed if both xi and x2 are distributed among the processors in the same manner, o+ Compute Pnw = PnFnz. This requires the communication pattern of dis- tributed FFTs discussed in Chapter 3. o+ Find Pnyi = ?21[Gn(Pnw) + PnW] PnY2 = ?2?[Gn(PnW) --H PnWJ The communication pattern required for G? is needed. Procedure 2a. Given: Two real sequences P?xi and Pnx?. Find: The DFT Yi and Y2 respectively of the two real sequences using an in-place complex FFT which computes FnPnz (either CT2 or GS2). o+ Form the complex vector Pnz P?xi + iPnX2. No communication is needed. o+ Compute w = FnPn(Pnz). Requires the communication pattern for dis- tributed FFTs. o+ Find Yi = -[Tv+wJ 2, Y2 = L[TW?w] 2 The communication pattern for permutation Tn is required. Procedure 3. The DFT of a real sequence done as a complex FFT of half length Let x be a real vector of length n, with x --H M?2?x Denote xe = x(0: --H 1) and x0 = x"'(--H?2: n --H 1). Form the complex sequence z = xe + ix0 of length ?fl2 Notice that xe and x0 are also real sequences so we can use the above procedure to find = Fn/2xC and y0 = Fn12X0. Let w = Fn/2Z, then = -[Tw+w] 2 yO =?2Z[Tw?w] 159 The Radix-2 Splitting algorithm gives us --H In/2 In/2 ?Fn/2 xe --H ?Fn/2 - X0 --H? yo where ? = diag(l,????and W? = exp(--H2iri/n). Since y = Ty, we only need to find y(O : n/2). Hence y(O:n/2--H1) = y(n/2) = ye(o)?yo(o) This procedure is entirely reversible and can be used to find the IDFT or DFT of a conjugate even complex sequence. Procedure 3. Given: A real sequence x(O : n --H 1). Find: The DFT of y(O : n --H 1) by using a radix-two complex FFT of half the length and post-sorting of data in bit-reversed order. o+ Form two real sequences xe of the even points of x and x0 of the odd points of x. No communication is needed if x is mapped consecutively and n = 2t > p = 2k Let z = Xe + ixo. o+ Implement Procedure 2. o+ Form Pny = P??e +(PnA)Pny0. No communication is needed here. Note that Procedure 3 requires the same data pattern as Procedure 2. Procedure 3a. Given: A real vector Pnx. Find: The DFT of x by using a complex FFT of half the length. o+ Let x = PnX. Form P?12Z = x(O n/2 --H 1) + ix(n/2 : n --H 1). Processor local. o+ Find w = Fn/2Pn/2(Pn/2Z) This requires distributed FFT communication pattern. o+ Form yf and y0. The communication pattern for T? is needed. o+ Form ? = yC + Ay0. No communication is needed. 160 Procedure 4. IDFT or DFT of a Conjugate Even Complex Sequence Suppose we have a conjugate even complex sequence y = Ty. This means that y(o) y(1 : n/2 --H 1) y(n/2) Ey(i n/2 --H 1) Going backwards, we can find ye(0) = --H[y(O) + y(n/2)] 2 ye(1 : n/2 --H 1) = : n/2 --H 1) + Ey(1 : n/2 --H 1)] 2 Y0(o) = 2-1[y(O) --H y(n/2)] y0(1 : n/2 --H 1) = 21??1(1 : n/2 --H l)[y(i n/2 --H 1) --H Ey(1 : n/2 --H 1)] Form ? = yC + jy0 and find z = F?Hw Now the real sequence x is found xe = Re[zJ x0 = Im[z] Procedure 4. Given: A complex conjugate even sequence y(O: n/2) (we only need to store the first n/2 + 1 elements because of the sym- metry). Find: The IDFT or DFT of y. o+ Form ye and y0. Communication pattern for permutation Tn/2 is required. Let w = ye + iy0. o+ Find Pn/2Z = Pn/2Fn/2 w. The communication pattern for distributed FFTs is required here. o+ Set Pn/2x? = Pn/2X0 = No communication is needed. Re[P?12 z] Im[Pn/2zj Procedure 4a. Given: A complex conjugate even sequence Pn12y(0 n/2 --H1). Find: The DFT x by using a complex FFT of half the length. 161 o+ Form Pn/2YC and P?12y0. This requires communication for Tn permutation. Let Pn/2 W = + o+ Find z = Fn/2 Pn/2(Pn/2W). This requires communication for the dis- tributed FFT. o+ Let xe = Re[z] and x0 = Im[z]. Processor local. 5.8 The Cosine and Sine `I?ansform The DFT of a real even sequence corresponds to a special transform called the cosine transform. Similarly the DFT of a real odd sequence corresponds to the sine transform. These transforms come up in special applications ranging from image processing to the solution of elliptic partial differential equations with various boundary conditions. In this section we introduce these symmetric transforms in matrix notation. We then focus on the sine transform and discuss the methods of Dollimore (1973) and Cooley et at. (1970) to compute the sine transform efficiently by exploiting all possible symmetries. Finally we talk our way through a hypothetical implemen- tation of these methods on a parallel distributed system such as the hypercube and show how the plethora of symmetry exploitation causes a very complicated sequence of permutations (communication). Matrix Description Definition 5.8.1 The cosine series of a sequence x(j), j = .... . , n is defined to be - L wT 1 - a --H 2 i2 i2 0n--H1 V X ? 2?W vT ?21(L1)n - -2 Notice that the cosine series is only defined on n + 1 points x(O: Theorem 5.8.1 Cosine Series: Let z be the even eztension of x, z = x(O) x(1 : n --H 1) x(n) Ex(1 : n --H 1) a(O : n) = 2Re[y(O : then if y = F2nZ, then 162 Proof z is real even, implying that y is purely real. Then 1 wT 1 wT w C v CE 1 vT (--H1)? vTE w EC Ev ECE x(0) F?z = 1 x(1 : n --H 1) 2n x(n) Ex(1 : n --H 1) y(O) = 2T[x(o) + 2wTx(l : n --H 1) + x(n)] y(1 : n) = 2L?[X(O)W+ 2Cx(1: n --H 1) $ x(n)v] y(n) = ?21n?40? +2vTx(l : n --H 1) +(--H1)?x(n)] y(n + 1 : 2n --H 1) = 2?[X(O)W +2ECx(1 : n --H 1) +Evx(n)J Notice that y(n + 1 : 2n --H 1) = Ey(1 n), which is what we expect. Definition 5.8.2 The Inverse Co?ine Transform is defined as x = 1 wT 1 - Cn-1 4 a 21 vT 21(21)n - *______ Hence the cosine transform is its own inverse up to a multiplicative constant. Definition 5.8.3 The DFT ofa reat odd sequence corresponds to the s:ne serzes defined on a sequence x(j), j = 1,... 2 b(1 n) = ??Sfl?1X(1 : Notice that the sine series is defined on only n --H 1 points x( 1 : Theorem 5.8.2 Let z be the odd extension ofx. 0 x(1 : n --H 1) 0 --HEx(1 : n --H 1) then ify = F2nz, b(1 n) = 2z'y(1 : 163 Proof Since z is real and odd, the y = F2nZ is purely imaginary, therefore o 0 0 0 0 F2nz = ? ? 5 0 --HSE x(1 : n --H 1) 2n 0 0 0 0 0 o --HES 0 ESE --HEx(1:n--H1) y(1 : = ThZ[2Sx(l : n --H 1)] 2n --Hz y(n + 1: 2n --H 1) = --H[--H2ESx(1 : n --H 1)] 2n Here y(n + 1: 2n --H 1) = --HEy(1 : n) as expected. _____*_____ Definition 5.8.4 The Inverse Sine Transform is defined as X = Sn?1b so that the sine transform is also its own inverse up to a muttiplicative constant. The Sine Transform Aigorithm The most straightforward sine transform algorithm would be to follow the deriva- tion in Theorem 5.8.2 and extend the given sequence of length n into its odd extension. The FFT of this 2n-point sequence is then computed and the sine transform is found in the imaginary part of the transform of the extended se- quence. This method has a number of drawbacks in the single processor case. For one, it requires four times the array storage than necessary. This is because the extended sequence is both odd and real. But if a complex FFT were to be done, we would need 2n points of complex storage. It is also obvious that we are doing computationally more work than is necessary. However given the difficul- ties with which the conventional sine transform procedures of Dollimore (1973) and Cooley et at. (1970) are implemented on the hypercube, this naive method is actually recommended by some researchers for parallel implementation. This is because given an efficient complex FFT and parallel computing, one may not care too much about whether a 2n-point transform is done instead of an fl-point or n/2-point transform. It should be noted that it requires a non-trivial set of communication to convert an fl-point vector mapped in consecutive order into a 2n-point odd extension of itself. Try it with paper and pencil! 164 The sine transform can be calculated efficiently via the complex FFT by taking advantage of the symmetry in the extended sequence. These algorithms usually involve the use of a real FFT of length n coupled with some pre- and post- process- ing of the sequence to exploit the additional symmetry. There are two different ways of doing this. Each is the inverse of the other. The first method is usually attributed to Dollimore (1973) and is presented by Swarztrauber (1982) and Press et at. (1986). Procedure 5. Dollimore's Sine Transform Suppose we have a sequence x(O : n) where x(O) = 0. An auxilliary array is constructed so that the first term is even and the second term is odd. Let A3 --H diag(sin ?), j =O,1?... , n --H 1. The constructed sequence is (5.8-1) d = A8(x + Tx) + -2(x --H Tx) Taking the FFT of this sequence yields Y = Fnd = FnA5(x + Tx) + 2Fn(x --H Tx) with Re[y] = Re[F?]A8(x + Tx) Im[y] = -Im[F?](x --H Tx) 2 First we look at Re[y]. We see that TA3 = A8T = A8 so that Re[yJ = 2Re[F?]A8x Now Hence [Re[F?]A8]? = cos 2rpq sin 1 . +fl1 sin 2 2 fsrn(2??p7r) . ?ffl1??)? Re[yJ = RnRe[F?Jx --H RnTRe[FnJx where R?T is the down-shift matrix. For Im[y], we have Im[F?]T = --HIm[F?] and thus Im[y] = -2[Im[F?]](x --H Tx) = Im[F?]x 165 Im[y] gives us the even terms of the sine transform, and a recursion on Re[yj gives us the odd terms. So b1 2y(o) b2k Im[y?J b2k+l H b2k?l + Re[y?] rn2? 1 This algorithm is not efficiently vectorized because of this last recurrence. Tem- perton (1980) has shown however, that the numerical properties of Dollimore's method is superior to the vectorizable method we present next. The Sine Transform of Cooley et at. (1970) The inverse of Dollimore's algorithm is the one presented by Cooley et aL (1970). A complex conjugate even sequence is formed out of the original real sequence. Let b(O : n) be a sequence of real numbers with b(O) = 0. Form the complex conjugate even sequence z as follows z(O) & (2b(1),o) Re[z(j)] b(2j + 1) --H b(2j --H 1) Im[z(j)] b(2j) j=1,...,n/2--H1 z(n/2) (--H2b(n --H Procedure 4 is used to find d = FnZ. d is a real sequence and the inverse of equation 5.8-1 gives x --H -41?3?'(Td + d) + -1(Td --H 2 Implementation of the Conventional Sine Transform Methods on Dis- tributed Systems We outline how the sine transform can be implemented via Dollimore's algorithm and the algorithm of Cooley et at. and describe the permutations necessary. Here we see that the implementation of either of these two sine transform algorithms would be inefficient due to the complicated maneuvering of the permutations. Procedure 5. Given: A real sequence x. Find the sine transform b of x via Dollimore's algorithm. 166 o+ Form the auxilliary sequence d = ?3(x + Tx) + -2(x --H Tx) The communication pattern for Tn is needed. o+ Find Pny = PnFnd via Procedure 3. Let the even points of Pnb = Im[Pny]. Solve the recurrence for the odd points. This is difficult because the entries of y are in bit-reversed order and we require adjacent entries to be in adjacent processors. Hence a distributed bit-reversal must be done first before the linear shift nearest-neighbor per- mutations R?T and Rn are performed. Procedure 6. Given: A real sequence x. Find the sine transform b of x via the algorithm of Cooley et ai. o+ Set up the conjugate even sequence z. Nearest neighbor connections or RT and Rn permutations are required. o+ Use Procedure 4 to find the real sequence d = PnFnz. o+ Find P?x = 1??1(GnPnd + Pnd) + 1(Gnpnd --H Pnd) 43 2 The communication pattern for recursive exchange Gn is required. Since the sine transform is its own inverse, any sine transform algorithm can be used to invert a transform. Hence we outline Procedure 5 and 6 with bit-reversed input. Procedure 5a. Given: A real sequence PnX. Find the sine transform b via Dollimore's algorithm. o+ Form the auxilliary sequence Pnd = Pn?3(P?x + Gn(PnX)) + ??1(Pnx --H Gn(Pnx)) Here we need the communication pattern for Gn. o+ Find y = FnPn(Pnd) via Procedure 3a. o+ Let the even points of b = Im[yj, and the odd points are solved by recurrence. This is inefficiently parallelized. 167 Procedure 6a. Given: A real sequence Pnb. Fnd the sine transform x via the algorithm of Cooley et al. Set up the complex sequence Pnz. This is difficult because the real part of PnZ requires nearest neighbor connections, or shift permutations. but everything we have is scrambled up. Now we have Pnb, but there is no easv way to commute R?T and Pn. o+ Use Procedure 4a to find the real sequence d = FnPn(Pnz). o+ Find 1 1 x = ?s?'(Td + d) + -(Td --H 4 2 The communication pattern for Tn is required. A quick analysis of the communication needed to implement these two sine transforms show that the recurrence for Dollimore's algorithm is inefficiently par- allelized. Setting up PnZ for the algorithm of Cooley et at. is also difficult and requires unnecessary communications to unscramble data. So basically, to im- plement these conventional sine transform algorithms on distributed systems, we need to be able to do distributed bit-reversal permutations efficiently. Since manv algorithms require a transform followed by an inverse transform of the data, with- out need for natural ordered data in the transform domain, unscrambling of the data is unnecessary. Therefore the cost of pre- and post- processing operations and the communication needed to carry out these algorithms outweighs the benefits of their design to minimize operations counts. This is a typical problem in par- allel computation where problems which are optimized for a single processor can actually be far more difficult to implement and carry much more overhead to par- allelize than an algorithm which is less efficient in a single processor environment. This is because with parallel computing, we desire work to be as independent of each other to maximize parallelism. Sometimes this means that it is cheaper to do some redundant computations to defer the need for communication until abso- lutely necessary. This is how we approach our new sine transform, an efficientlv parallelizable version. 5.9 Conclusion and Discussion Although we have mainly concerned ourselves with the sine transform in this chapter, the exact same type of reasoning behind our parallel sine transform can be used for a parallel cosine transform algorithm. 168 The sine transform is used extensively in the Fourier method for the solution of boundary value problems in partial differential equations with Dirichlet-Dirichlet boundary conditions. [Swarztrauber (1977)] In fact, the tensor product method for finding the solution to the Poisson equation with Dirichiet boundary condi- tions involves a two-dimensional sine transform [see Lynch, Rice and Thomas (1964)]. ?Vhile in the single processor case, one would choose to use the matrix decomposition method [Buzbee (1973)], involving a fast transform in one direction and tridiagonal solves in the other, with the hypercube, it might be interesting to substitute the tridiagonal solves with fast transforms that are more efficiently parallelizable than tridiagonal solvers usually are. Sine and cosine transforms are also used extensively in image processing which may require transforms in two or more dimensions. Therefore fast efficiently implemented algorithms for vector computers and parallel systems are important and should be considered a basic building block in many other more complicated algorithms. Chapter 6 Fast Symmetric Transforms 6.1 Introduction Sine and cosine transforms arise naturally from the Fourier transform of symmet- ric sequences. They are used in a variety of problems in applied mathematics ranging from the digital coding of images [Clarke (1985)] to the fast solution of Poisson's partial differential equation with various boundary conditions [Swarz- trauber (1984b)J. Since they exhibit certain special symmetries, the computation of symmetric transforms should be up to four times more efficient, in terms of storage and op- erations count, than the complex FFT. This is because the transform of a real sequence is conjugate even, so there is a two-fold redundancy of information, Furthermore, the sine transform comes from the DFT of an extended odd real sequence and is itself odd. Therefore another two-fold redundancy of information is present. Hence a smart efficient Fast Symmetric Transform algorithm should require a quarter of the array storage and roughly a quarter of the computational effort of a complex FFT algorithm. There are several sequential algorithms that achieve this objective (discussed in Chapter 5 for the sine transform). [Cooley et ai. (1970), Swarztrauber (1982), Press et al. (1986)] These algorithms require pre- and post- processing and are very difficult to implement efficiently on a parallel memory-passing system. Therefore in Chapter 5 we present a parallel algorithm specifically for the sine transform. The technique of Chapter 5 can be used on other symmetric transforms but each case must be worked out individually. In this chapter we present a new algorithm that can simultaneously find the sine, cosine, quarter-wave sine and cosine transform of a single real sequence by computing only one complex FFT of the same length. Furthermore, all multipli- cation in our method is by roots of unity and hence the method is numerically stable. Portability is no problem because after an initial permutation of the data 169 170 any standard FFT subroutine can be used to find the results. Other methods used to find symmetric transforms include the method of trans- forming extended sequences. A sequence is extended so that its symmetry becomes obvious. For example, given a data vector Xj, ? = 1,... , n --H 1, we can extend ir into an odd function (0, xl, . . . , Xn?I , 0, Xn?1 , . . . , --Hxl The DFT of this extended function is both odd and conjugate even, and thus purely imaginary. Furthermore, the sine transform of Xj (up to a factor of 2i). is found in the imaginary portion of the first n points of the DFT of the extended function. [See Press et at. (1986)]. This method is as numerically stable as the Fast Fourier Transform, a sequence of unitary transformations, but requires the use of an extended sequence of twice the length. Hence, instead of reducing the work by a factor of four, we have actually increased the work by a factor of two. For higher- dimensional transforms, the factor increases as 2d where d is the dimension of the transform. This also means that 2d times the storage space is needed to transform these points. Therefore this method is neither efficient in the computational sense nor in the use of memory. The advantage however is that this method is portable to any architecture or machine where efficient FFT subroutines have already been worked out, implemented and tested. Direct methods are discussed by Chen, Smith and Fralick (1977), Yip and Rao (1980), Wang (1981a?), (1984), the gist of which is best explained by Swarz- trauber (1986b). The basic idea is to emulate the factorization of the DFT ma- trix. By factoring the discrete sine transform (DST) and discrete cosine transform (DCT) matrices into sparse matrices, a fast algorithm for these symmetric trans- forms can be derived with reduced operations counts since redundant operations are not done. Furthermore since the operations done are a subset of the ones done by the FFT algorithm, the direct methods possess the numerical stability of the complex FFT algorithm. The major drawback of these direct methods is the complexity of the signal flow diagrams. This makes it difficult to implement in parallel, as significant time could be taken just to find out where something has to go. See Wang (1984) and references therein for a complete description of the algorithms as well as the theoretical background leading up to these factorizations. This chapter is summarized as follows. The four types of cosine and sine transforms are first defined in terms of their representation as matrix-vector mul- tiplication. The new method is presented and a demonstration is given on how the symmetries are exploitable. This method is then implemented on the hyper- cube and timing results are analyzed. A few applications are discussed and in the appendix details of the proof for the type III and IV transforms are presented. 171 Finally the FORTRAN code for this set of procedures is given to show how easv it is to implement. 6.2 Definitions and Matrix Notation The four types of symmetric transforms are described by Wang (1981a--Hc) as a real approach to harmonic analysis. If f(x) is a bounded real function it can be expressed in terms of the orthonormal set of the trigonometric functions 7rLsinLr4+(n+0)x], n = . .., --H2,--Hi, 0, 1,2,.. n has to be an integer, but ? can be any real number. In practice, however, only integer and half integer (Q = k and OL = k + ?21, for k integer) harmonics are used. This is because in these two cases, symmetry properties of a sequence can be used to decompose the real transform into discrete cosine and sine transforms. Another reason is that it is easier to find fast algorithms for these two types of harmonics. The matrices for the four versions of the discrete cosine transforms (DCTs) and the discrete sine transforms (DSTs) corresponding to the integer and half- integer harmonics are described by defining their pqth element. The transform of a particular sequence is then a simple matrix-vector multiplication. DCT-I: 4Zn2[kpkq cos(pq?/n)J, p, q = 0, 1,... , DCT-II: ?[kpcos(p(q + 21)7r/n)j, p,q=0,1,... ,n --H 1. DCT-III: ?[kq cos((p + --H21)q?/n)J, p, q = 0, 1,..., n --H 1. DCT-IV: ?kos(p + 21)(q + -21)ir/n)j, p, q = 0, 1,..., n --H 1. DST-I: 7?2[sin(pqir/n)J, p, q ?...., n --H 1 DST-II: ?[kp sin(p(q --H 21)ir/n)J, p, q = 1,..., n. DST-III: ?[kqsin((p --H --H21)q?/n)], p,q = 1,... DST-IV: ?[sin((p + -21)(q + -21)?/n)J, p, q = 0, 1,..., n --H 1. where if l#0 orn 172 It can be shown by simple verification that DCT-I is its own inverse, DCT-II and DCT-III are inverses of each other, DCT-IV is its own inverse, DST-I is its own inverse, DST-II and DST-III are inverses of each other, and DST-IV is its own inverse [Wang and Hunt (1985)]. These symmetric transforms are not artificial creations, but appear naturally as the eigenvectors of the finite difference operator of Poisson's equation on a square with various boundary conditions. [Swarztrauber (1977)] DCT-I and DST- I are the usual cosine and sine transforms, respectively. And DCT-II and DST-II are the quarter-wave cosine and sine transforms with DCT-III and DST-III their inverses, respectively. These transforms are also used in digital signal processing applications and in image processing. [Ahmed et aL (1974), Kekre and Solanki (1978), Jain (1976, 1979), Wang and Hunt (1985)] Matrix Notation Each of these transform matrices can be simplified into an n-by-n matrix consisting purely of cosines or sines. We get rid of the k1 scaling terms. The matrix-vector multiplication of this is almost the corresponding discrete transform with a few adjustments that appear later. C1? --H kos(pq?/n)], C1?1 --H [cos(p(q + ?21)7r/n)], CnHI --H [cos((p + --H2)qir/n)], CInV = [cos((p + 2--H1)(q + ?21)7r/n)], S1? --H [sin(pq?/n)], SnH --H [sin(p(q + 2 sU' --H [sin((p + 2 SInV p,q=O,1,...,n--H1. p,q=O,1,...,n--H1. p,q=O,1,...,n--H1. p,q=O,1,...,n--H1. p,q=O,1,...,n--H1. p,q=O,1,...,n--H1. p,q=O,1,...,n--H1. p,q=O,1,...,n--H1. = [sin((p+ L)(q 1 2 Notice that these definitions are not exactly the same as the corresponding def- initions of the discrete transforms, however, once we have figured out let's sav C?1+1 x, it is easy to find the DCT-I of x by subtracting fi?21xo from each com- ponent, adding ?X? to the even components, and subtracting ?X? from the odd components. Finally the whole thing is scaled by 7?2 A summary of these "retrie'val" procedures appears later. 173 6.3 Description and Development In this section, we present our new method of computing the symmetric transforms so that those of type I and II and computed together and those of type III and IV are computed together. First we introduce a permutation matrix that is used in this analysis. We then give an overview of the relationships between the symmetric transforms, after which we derive these relationships in detail. Quarter-wave Splitting Permutations The permutation that is used to simplify the description of sequences with quarter- wave symmetries is defined as follows. I[n/21 = E[n/2J A highly related "manipulation" of a data vector is I[n/21 --H E [n/2J Mn In English, Kn takes all the even points of a sequence first, and then the odd points going backwards. L? does the same thing except it takes the negative of the odd points going backwards. Example 0 0 0 0 1 2 1 2 2 4 2 4 Kn --H 6 Ln 6 4 7 4 --H --H7 5 5 5 --H5 6 3 6 --H3 7 1 7 --H1 The Four Types of Symmetric Transforms All four types of cosine and sine transforms are related to a certain complex matrix-vector multiplication. Recall that the discrete Fourier transform matrix of dimension n is a matrix consisting of powers of w --H e-2ri/n. Let F? be the DFT matrix of order n and F2n be the DFT matrix of order 2n. We consider the upper-left n-by-n corner of the matrix F2n denoted by F2n = F2n(0 : n --H 1,0: n --H 1) 174 and show that all of these transforms involve either the term F2nX or F2n( Ax) where A is a particular scaling of x by a diagonal matrix of twiddle factors. i.e. 4n-th roots of unity. The matrix F2n provides the link between the various transforms. The diagonal matrices are n-by-n matrices consisting of A = diag(e?rik/2n) A2 = djag(e?rik/n) Table 6.1 displays the formulas involved. One can see that the type I transform can be obtained from the type II transforms and the type II transforms from the type IV transforms. From Table 6.1 we see that the DFT of KnX or LnX (permuted versions of x are all that are required to be computed via the FFT algorithm. Further the FFT is of length n rather than 2n. Basically two DFTs are needed for each set of four symmetric transforms: Fn(Knx) and Fn(Lnx) for the type I and type II transforms, and Fn(A2Lnx) and Fn(A2Knx) for the type III and type IV transforms. The first two are real transforms and the second two are complex transforms. Makhoul (1980) describes a way of finding C1?1X by find the real part of the transform AFn(Knx). He establishes the equivalence between Re[AF?(K?x)] and Re[AF2?x]. What we do is establish the same sort of procedure for the type II sine transform S1?1X Once we have C1?1X and SnHx available, we show how it is an easy step to finding the type I transforms. We form a sequence by putting C?11x in the real part and S'Jx in the imaginary part. We then scale the sequence by A--H' and then find C?1X in the real part of the resulting sequence and S1?X in the imaginary part. A similar procedure is used for the types III and IV transforms. The algorithms that follow describe the techniques used. Summary of Procedures Algorithm 6.3.1 describes how to find the cosine and sine transforms of type I and type II simultaneously. Algorithm 6.3.1 Procedure Typel2: v . K?x; z Lnx v `. F?v; z .` FnZ v A?nv; z .` A4nz C1?1X = Re(v); S1?1X = --HIm(z) Re(a) Re(v); Im(a) = Im(z) a A?'a C1?X = Re(a). S1x --H --HIm(a) 175 Table 6.1: Symmetric Transform Formulae transform formulae DCT-I C?X = Re[F2?x] --H Re[??l(c?Hx + ?S??X)] = Re[?F?K?xJ DST?II x = --HIm[?F2?x] = --HIm[?F?L?x] DCT-III x = Re[F2??xJ = Re[e??4fl??l(CI?Vx + iSInVX)J DST-III s x = --HIrn[F2??x] = ?Im[eTi/4n??l(cI?Vx + ?SInVX)j DCT-IV x = Re[e?rt 4??F2n?xJ = Re[e?wi/4n?Fn?2Lnx] DST-IV x = --HIm[e-11 4??F2n?x] 176 Algorithm 6.3.2 describes the procedure for obtaining the type III and type IN sine and cosine transforms simultaneously. Algorithm 6.3.2 Procedure Type34: v K?x; z LnX v ?2v; z ?2z v F?v; z Fnz v e?r?4n?v; z c?IVx --H Re(z); SInVx = --HIm(v) Re(a) = Be(z); Im(a) = Im(v) a e?14?A?1a cU1x = Re(a); Sifli1x = --HIm(a) We next describe how Algorithm 6.3.1 is derived. The derivation of Algo- rithm 6.3.2 is similar, but more complicated and is hence relegated to the Ap- pendix. Procedure 1: C1?? S?1? cin', and sin' from two real FFTs of length n The method of computing the type II cosine transform Ciflix by equating cUx = Re[?F?K?x] is derived and demonstrated by Makhoul (1980). Theorem 6.3.1 (Makhoul (1980)) Let x be a real vector oflength n, then cnI'x=Re[?F2?x]=Re[?F?K?x] Proof Makhoul (1980). *______ The type II DST can be computed in a similar manner to that of the type II DCT. We first present a lemma equating the imaginary part of the 2n-point DFT matrix multiplied by diagonal matrix of 4n roots of unity with the type 11 sine transform matrix. Lemma 6.3.1 ?F2n --H ??1F2n = --HiS11 n 177 Proof [?F2n --H ??1F2n]pq = ?--Hrip(2q+l)/2n --H ?rip(2q+l)/2n = --H2?sin(irp(2q + 1)/2n) ______*______ The following theorem shows that the type II sine transform has two expressions as imaginary parts the complex sequence ?F2nx and ?FnLnx, respectively. Theorem 6.3.2 Let x be a real vector oflength n then (6.3-1) S1?1X = --HIm[?F2?xj = Irn[?FnLnx] Proof Let y be the quarter-odd extension of x _ x -Ex Partition the 2n-point DFT matrix F2n as follows: F2_ --H [FL FR] Then letting Y be the DFT of y we have Y=F2 y?FLx?FREx However it can be easily shown that FRE --H so Y = FLx --H A?2FLx Defining F2n = F2n(O : n --H 1,0 : n --H 1) we have for the first n terms of Y Y(0 : n --H 1) = F2nX--H ??2F2nx --H ??1[AF2n --H ??1F2n]x --H ??`2iIm[?F2?xJ At the same time we have --H ??1F2n = iS1?1 hence Y(O : n --H 1) = ??12?S1?1X. 178 Thus S?11X = --HIm[?F2?xj and the first part of the proof is complete. The second part of the proof is to establish the relationship sUx = --HIm[?F?(L?x)J As above, we let y be the quarter-even extension of x, and let z(j) = w(j) - y(2j) y(2j+1) J =O,...,n--H1 and note that z = LnX and w = --HEz. Again let Y = F2ny. Denoting = F2n(:,O:2:2n--H1) F20? = F2n(:, 1: 2 : 2n --H 1) be the 2n-by-n matrices containing the even and odd columns of F2n, respectively, we have Y = F2ny = F2Cnz + F20?W = F2e?z --H F20?EZ Denoting the top half or F2C? and F02?, respectively: = F2?n(O:n?1,:) F0 2n = F20?(O : n --H 1,:) then we already know that Y(O : n --H 1) = F2Cnz --H F20?EZ = Fnz???2Fnz --H --H ??1Fn]z Y(O : n --H 1) = --H??`2?SIn'x hence this establishes the relationship S?11X = --HIm[?F?L?x] 179 ______*_____ This shows that the type II sine transform can also be found by an FFT of a permuted real sequence z oA length n. Expressions involving F2nX can also be found for the type I cosine and sine transforms. Theorem 6.3.3 If x is a sequence of length n, then C1?X = Re[F2?xj Proof The even extension of x is Let Y = F2ny, and F2? = [FL,FR], then F2ny --H But C1?X 21[(F2nY)(O : n --H 1)j, hence the result C1?X = Re[F2?x] FLx + FRTx (FL+?L)x 2Re [FL x] Theorem 6.3.4 If x is a sequence of length n, then S?1X = --HIm[F2?x] Proof The odd extension of x is _ x --HTx Let Y = F2nY, and F2n = [FL, FRi, then F2ny --H FLx?FRrrx (FL??L)x = 2iIm[FLxi But S1?X = -2?(F2ny)(O : n --H 1), hence the result S?1X = --HIm[F2?xJ ______*_____ 180 *_____ Now that all the relevant relationships are present: C1?1X = Re[AF2?xJ = Re[?F?(K?x)] s'Jx = --HIm[?F2?xj = --HIrn[?F?(L?x)] C1?X --H Re[F2?xJ S?1X = --HIm[F2?xJ we can produce the procedure that computes the four type I and II symmetric transforms. Let Then and gives A = ?F'2nX ReA = Re?F?v ImA = Im?F?z B = C?1X = ReB S?1X --H --H1mB Thus we need to find two real DFTs of length n. It is a well-known fact that one can find these two DFTs in one pass through the complex FFT subroutine. Thus by computing only one complex FFT of length n, we have managed to finci four sine and cosine transforms of length n. Procedure 2: CnIH, S1?11? CnIV? and SInV from two complex FFTs of length The type III and type IV symmetric transforms have a very similar relationship whereby the type III transforms can be computed in linear time once the type IV transforms are present. The type IV transforms can be found efficiently by taking the complex FFT of vectors of length n, i.e. a scaled version of either KnX or L?x. 181 The derivation and proof is long-winded and is thus relegated to the appendix. The important results are theorems 6.8.1 and 6.8.2 giving: CnIVx = Re[e??/4??F2?Ax] --H Re[e???14????2L?x] and theorerns 6.8.3 and 6.8.4 giving: SInVx = Also neeed are the relationships CnHIx = Re[F2?Axj and ?rm[e?Ti/4n?F2n?xi Im????4fl?F??2K??] SnIHx = ?Irn[F2n?x] Now Algorithm 6.3.2 is derived as follows. Let A = e???4fl?4??2??4?? Then so letting gives Re[A] = Re[e??14??F??2zj Im[A] = Im[e???/4fl?Fn?2vJ B --H C?111X = Re[BJ sU1x = Irn[B] Therefore we need to find two complex DFT's of length n Fn(?2Lnx) Fn(?2Knx) to get the four sine and cosine transforms of length n. 182 How to get the Transforms In the procedures above we have always assumed that X was a vector of length n, with subscripts ranging from 0 to n --H 1. However, the exact definitions of son?e of these transforms are slightly different. For example, DST-III assumes that the subscripts run from 1 to n. Furthermore there is that constant k1 to deal with. ?Ve now show how easy it is to get the exact definitions from the output of the procedures. DCT-I of x(O : Given: Y --H C'?x(O: n --H 1). Overwrite Y with the DCT-I of x. Y(O) Y(k) Y(n) 1 1 2 --H xo) -2xo + 2X? Y(k)--H(i--H 2)xo+(?1)k 1 ?Track FFT described in Chapter 3. o+ Post-processing is multiplication by ? a diagonal matrix and other local operations. Algoritbm 6.5.1 Type 12 Transform 188 Table 6.3: Sine Transform (Type I) Errors Type I Sine Transform n = 2 y = x random y = sin(x) 4 = 2 1.49O12e-O8 1.49O12e-O8 5.96O46&O8 8 --H 2? 5.96O46e-O8 5.96O46e-O8 2.98O23e?O8 16 --H 2? 5.96O46e-O8 5.96O46e-O8 2.OO234e-O8 32 --H 2? 2.98O23e-O8 5.96O46e-O8 5.96O46e-O8 64 --H 26 2.98O23e-O8 4.47O35e-O8 5.96O46e-O8 128 --H 2? 1.192O9e-O7 5.O2914e-O8 2.98O23e-O8 256 --H 28 1.78814e-O7 1.86265e-O8 1.192O9e-O7 512 --H 2? 5.96O46e-O8 2.98O23e-O8 5.96O46e-O8 1024 --H 210 2.38419e?O7 5.96046e-08 2.98023e-08 2048 --H 211 1.78814e-07 5.96046e-08 5.96046e-08 4096 --H 212 5.96046e-08 2.98023e-08 5.96046e-08 8192 --H 213 5.96046e-08 2.98023e-08 5.96046e-08 16384 --H 214 1.19209e-07 2.49129e-08 5.96046e-08 32768 --H 215 5.96046e-08 2.21189e-08 5.96046e-08 65536 --H 216 5.96046e-OS 2.21189e-08 2.98023e-08 189 /* ? = processor id; n --H 2t; p --H 2d; m = (n/P)/2 */ /* Data is entered into the hypercube as MnV */ /* Initially, processor ? holds: */ /* v(') --H v(2m.id[?:2:2m.Yd[?+1)?1) */ /* v(2) = v(2m id[?] + 1 : 2 : 2m (id[ii] + 1) --H 1) */ call Exchange(v(2)) /* code (6.6) */ z(1) z(2) --Hv(2) call Two-Track FFT(v('),v(2),z('),z(2)) /* code (6.2) */ v ?v z Of course, for any hypercube procedure the facility to start with bit-reverse permuted data is necessary in order to inverse transform a sequence. Here we outline the reverse procedure. Starting with (M2?2p? 0 1n12p)Pnx o+ Apply the permutation (112?2P? 0 1n/2P) to get back just Pnx. o+ Now we need to have PnKnX so that we can apply the inverse FFT to Pnv. But we have Pnx. So we apply E P?x = PnKnx Need to do exchange permutation of second track. o+ Do Two-Track FFT and find FnPn(Pnx). o+ Multiply by the diagonal matrix ?, and other local post-processing. Timing results are given for Algorithm 6.3.1 in Table 6.4 for the forward trans- form with data in natural order and in Table 6.5 for the inverse transform with data in bit-reversed order. The hypercube is mapped by the BRGC mapping. The speedup here is excellent. 6.6 Applications: Updating Shifted Sequences The stability of this method is good for updating shifted sequences as most of the error of the traditional methods is concentrated in the first few terms and last few terms. 190 Table 6.4: Algorithm Type12 n = 16384 Timings (Forward) dim total comp comm blocked copy diag per speedup 0 95735 94780 0 0 950 24310 23935 --H 1 49430 48620 335 50 480 12155 12360 1.94 2 25985 24935 810 175 235 6080 6495 3.69 3 13515 12790 600 265 120 3040 3380 7.08 4 7220 6560 600 290 60 1520 1805 13.26 Table 6.5: Algorithm Type12 n = 16384 Timings (Inverse) dim total comp comm blocked copy diag per speedup 0 84445 75765 0 0 8665 11395 21110 --H 1 43730 39125 440 65 4165 5705 10930 1.93 2 24250 20185 2060 495 2000 2855 6060 3.48 3 12600 10410 1220 530 1205 1425 3150 6.70 4 6720 5375 880 410 460 715 1680 12.57 Suppose we have sequence x = [xo,. . . ,xn?i]T and we wish to find the sine and cosine transforms of its shift x+ = ....... , xn]T. First let us define the upshift matrix R??1 --H [en?i, eo, e1,... , e??iJ. Then =Rn?1X?[0,.. . , xo]T + [0,... , xn]T If we had C1?X and S?1X, we can find C?1X+ and S1?X+ easily without having to do any more transforms. This is because C1?X+ = --H[0,..., xo]T + [0,..., xnjT) S1x+ --H S?1(R?1x --H [0,.. , xoJT + [0,.. , xn]T) Lemma 6.6.1 C1?R?1 = ?cCn1+?sS1n S1?R??1 = ?cSn1 --H ?sC1n tetting ?c = diag[ccs(?p/n)] and ?? = diag[sin(irp/n)]. Proof [Cn1Rn?']pq = cos(irp(q --H 1)/n) = cos(rpq/n)cos(?p/n) + sin(?pq/n)sin(rp/n) 191 Table 6,6: Updating of Shifted Sequences transform update DCT-I C?X+ = ?cCnX + ?sS?X + (xn --H xo)C?(:, n --H 1) DST-I S'nx+ = ?cS1nX --H ?8C1?X + (xn --H xo)5?'(:, n --H 1) DCT-II C1?1X+ = ?cCnHx + ?8S1?1X + (xn --H xo)C?H(:, n --H 1) DST-II SUx+ = AcSUx--HAsC'n'x+(xn--Hxo)Sn"(:,n?l) DCT?III C1?11X+ = AcCU'x + A351?11X + (xn --H xo)CU'(:, n --H 1) DST-III S1?11X+ = A S1?11X --H AsCnHIx + (Xn --H xo)S'?"(:, n --H 1) DCT-IV CInVx+ = AcCInVx + AsSInVX + (xn --H xo)c?IV(:, n 1) DST-IV SInVx+ = AcSnIVx?AsCInVx+(xn?xo)snIV(:,n?l) and sin(?p(q --H 1)/n) = sin(:"pq/n)cos(:"p/n) --H cos(irpq/n)sin(?p/n) Next notice that A[o,... biT?bA(:,n?1) that is, equal to b times the last column of A. Therefore updating shifted sequences involves some scaling only. summarizes the updating procedure. See also Yip and Rao (1987). 6.7 Applications: Fast Poisson Solvers. ______*_____ Table 6.6 Following Van Loan (1987) and others, the eigenvector matrix of the discrete finite difference operator of Poisson's equation with various boundary conditions can be written in terms of DCT or DST matrices. Here, types 1--H111 are involved. 2 QDD(n --H1) = ?[sin(?r/n)], p,q =1,... ,n --H1. ???=[sin(pqir/n)], p,q=1,...,n--H1. Dirichiet-DiricAlet 192 Neumann-Neumann Dirichlet-Neumann and Neumann-Dirichlet 9 # 1) = --H[kq cos(pqr/n)] p, q = 0, 1,.., , n = [kq cos(pq?/n)] p, q = 0,1 = ?2[kqsin((P --H1)qir/n)] p,q = 1n 2 Q??'?=[sin(p(q--H--H)ir/n)j p,q=1,..?,n. 2 9 1 QND(n) = ??[kqcos((p+ --H2)qr/n)J p,q = O,1,...,n--H 1. where QN?'D = [cos(p(q+ --H21)ir/n)] p,q =O,1,...,n --H 1. 121 or n. The stability of this method is not an issue until h3 < ? since the discretization error of finite differencing is already h2, and the error growth rate of the Coo1e?- et ai. algorithm is e/h. 6.8 Appendix: Proof for Type III and Type IV Symmetric Transforms We derive a procedure that allows us to find the type III and type IV sine anci cosine transforms all at the same time. Since the symmetry is a bit more compli- cated, this requires two complex FFT's of length n instead of one. 193 More notation is needed: F4'? --H F43? - F44? --H F4? - F40?1 --H F02?1 - F?n(:,0 : n --H 1) F4n(:,n : 2n --H 1) F4n(:,2n : 3n --H 1) F4n(:,3n : 4n --H 1) [F', F2, F3, F4] F4n(1 :2: 2n --H 1,0: n --H 1) F2n(1 :2: 2n --H 1,0: n --H 1) diag(1,--H1,...) djag(e?2wik/4n) The following lemmas are used in the main result. Lemma 6.8.1 Let D =diag(1,--H1,...) and ? = diag[exp(--Hp2iri/4n)], then F2E = F43? = DF41? F44?E = Proof Lemma 6.8.2 [F42nE]pq = [exp(--H2irip(2n --H q --H 1)/4n)] = [D??1?F41n]pq [F34? Jpq = [exp(--H2ir?p(q + 2n)/4n)] = [DF41?]pq [F44nE]pq = kxp(--H2ir:p(4n --H q --H 1)/4n)] --H [??1?F41n]pq F04?1 = F4n(1 :2: 2n --H 1,0: n --H 1) Proof We look at the pqth element of each matrix. [F04n1jpq = exp(--H2iri(2p + 1)q/4n) exp(?27r?ZnP(I --H ___ --H [F2n?Jpq *______ *______ 194 Lemma 6.8.3 e???/4fl?F2nA + ???4fl??l?2???I = 2cIv Proof Look at the pqth element and see that e????4?e?2???/4??--H4Tipq/4ne--H2wiq/4n + ewi/4ne2????4??4riPq/4n?2T?q/4n = 2cos is the pqth element of 2CI?V *______ Now we relate the type IV cosine transform CInVx with the transform of an ex- tended sequence. Theorem 6.8.1 Let x be a real vector oflength n Then CIVx --H Re [e?ri/4n?F2n?x1 Proof Let y be x -Ex --Hx Ex Y = F4n? = F41?x --H F4%Ex - F34?x + F44?Ex Notice that because of the D matrices, the even components of Y are zero, thus we are only interested in Y0 = Y(1 : 2 : 2n --H 1). Yo = 2F40?'x + 2[diag(exp(ri(2p+ 2F04?'x + 2eTi/2flA?2VO4?1x 2F?n?x + 2eT?2fl??2?2???1? 2e;ri14n ??1 (2Re[e?ri/4n ?F2n ?x]) 2eTiI4n ??1 [2cnIVxi cIVx --H Re [e?'r?'4??F2n?x] or *______ 195 Fast CosineiV Transform (FcTIv) More relationships concerning the 4n-point DFT matrix follow in the next lemma. Lemma 6.8.4 Proof F4n(:,O :2: 2n --H 1) = F2n(:,O : n --H 1) F4?(:,2n :2: 4n --H 1) = DF2?(:,O : n --H 1) F4n(:, 1: 2 : 2n --H 1)E = D??'?F2n(:, 0 : n --H 1) F4n(:,2n + 1 : 2: 4n --H 1)E = ??`?F2n(:,0 : n --H 1) [F4n(:, 0 : 2 : 2n --H 1)]? = kxp(--H2?i(p2q)/4n)j = [F2n(:, 0 : n --H 1)]pq [F4n(:, 2n : 2 : 4n --H 1)]? = kxp(--H2?ip(2q + 2n)/4n)] = [DF2n(:, 0 : n --H 1)jpq [F4n(:, 1 : 2: 2n --H 1)EJ? = [exp(--H2irip(2n --H 2q --H 1)/4n)J = [D??'F2n(:,0 : n --H 1)]pq [F4?(:,2n + 1:2: 4n --H 1)EJ? = Lemma 6.8.5 Proof [exp(--H2irip(4n --H 2q --H 1)/4n)J [??`?F2n(:, 0 : n --H 1)jpq F2%' = [F02n'Jpq = F2n(1 :2: 2n --H 1,0. n --H 1) exp ?2??(2p?+l)q exp(?2$? 22$q) [Fn ?2j? * ______ ______*_____ The following theorem relates CInVx to another transform, but this time of a sequence of length n. 196 Theorern 6.8.2 Let x be a real vector oflength n, then CInVx = Re[e??2/4fl?Fn?2z] Proof Let y = [x,--HEx, ?x, ExiT and define y(2j) y(2j+1) j=O,...,2n--H1 Notice that w' = Ez' and ZI = [z, --Hz] where x(2j), ffOcq by the contigtons column mapping into the hypercube. Step 2: Apply CTl to the rows of Xmxq. Step 3: Compute and store the twiddle factors TmxqPq. Step 4: Apply the twiddle factors. Step 5: Call FFTPACK to forward transform the columns of Xmxq of length rn. a: Factor m and store in array vsave. b: Calculate the forward transform of each column of Xmxq. Step 6: Call FFTPACK to back transform the columns of Xmxq. a: Use the saved array vsave. b: Calculate the inverse transform of the columns of Xmxq. Step 7: Apply the conjugate of the twiddle factors stored, TmxqPq. Step 8: Apply GS2 to the rows of Xmxq. End: Get back what you started with. Algorithm 7.3.1 Mixed-Radix-2 Forward FFT = processor id; n = mq; P = 2d; r = n/P; s = q/P */ /* Initially, processor ? holds: */ /*x(r'id[?j :r.(id[?+1)--H1) */ call Distributed CTl(x) /* compute Xmxq TmxqPqxmxq */ for j = 0 : q --H 1 = b?trev(s id[?] + j) for i = 0: m --H 1 h = r + j' m + i 208 x(i+j.m) ?wnh4i+j.m) end end call FFTPA CK(x) Algorithm 7.3.2 Mixed-Radix-2 Inverse FFT = processor id; n = mq; P --H 2d; r = n/P; s = q/P */ /* Initially, processor ? holds: */ /* y(r . : r (id[?] + 1) --H 1) */ /* where y =--H (yqT??Pq) */ call FFTPACK(y) /* compute YmxqTmxqPq */ for j = 0 : q --H 1 = bitrev(s . ?d[?] + j) for i = 0: m --H 1 h = id[?J r + j' m + i x(i+j.m) ?w?x?+j.m) end end call distributed GS?y) This procedure has been coded up on the Intel iPSC using FFTPACK for the processor-local FFT and distributed Cooley-Tukey and Centleman-Sande algo- rithms documented in Chapter 3 and Appendix A. Times in milliseconds for a 1280-point forward and inverse transform pair is given in Table 7.1. Our examples show that the minimum total time occurs roughly where the local transform and the distributed transform times are close. The overhead of twiddle factors is unavoidable, so to maximally exploit the parallelism in there one should use the largest cube available. Finally the factorization time for m increases as m increases. However if many transforms are done with the same m. the factorization can be calculated beforehand and stored, likewise for the twiddie factors. In our implementation, the factorization is calculated once and stored. as are the twiddle factors. A model for determining how large P should be takes into account both corn- putational complexity and communication costs. The computational workload is directly divided between the numerous processors, hence favoring the use of more processors. Yet as we add more processors, we get into more complex corn- munication because the distributed FFTs require communication complexity of 209 Table 7.1: Mixed Radix-2 Timings m x q = 1280 point MR2-FFT (Forward?and Inverse) log2p total dist. local twiddle factor 1 3540--H3640 2015--H2115 670 355 5 1 3265--H3360 1670--H1755 720 865--H875 10 1 3005--H3100 1410--H1490 720 8607875 15 1 2945--H3040 1190--H1270 875?80 850--H865 30 1 2810 995 940 820 55 1 2735 810 1040 780 110 1 2715 680 1125 695 215 1 2930 630 1335 540 439 q m 256 5 128 10 64 20 32 40 16 80 8 160 4 320 2 640 2 256 5 2 128 10 2 64 20 2 32 40 2 16 80 2 8 160 2 ____4 320 3 2?56 5 3 128 10 3 64 20 3 32 40 3 16 80 38 160 4 256 5 4 128 10 4 64 20 4 32 40 4 16 80 1815--H1915 1050--H1150 335 425 1680--H1775 880--H965 365 440 1550--H1655 745--H840 365 425440 1555--H1655 670--H755 440 415435 1465--H1515 545--H570 470 395425 1465 480 530 855 1520450580 280 965--H1040 580?50 170 215 890--H945 490--H535 185 210--H220 815?90 405475 185 220 800875 350410 220--H260 200--H225 790-840 310--H335 240 180--H215 820-875265--H300255--H270 145--H220 5 10 15 30 55 110 215 -5 10 15 30 55 105 5 10 15 30 55 500--H555 300--H360 85 460--H510 255--H300 90 105--H115 430480 230--H260 95 115 435490 200--H240 115 85--H115 445--H525 185--H245 115--H125 110 210 2 log2 P. Each processor has ?L points when n = rnq and P --H 2d Assuming that m is prime, the local FFTs cost O(m2q/P) flops. Multiplication by twiddle factors cost another ? multiplications and the computation for the distributed FFTs measures O(mq log2 q/P) flops. Hence, letting tcomm be the per element data transfer rate, r the communication start-up or latency period, and tcomp the per flop computational rate, a model for the Mixed-Radix 2 FFT is T = m2(p?)tc0mp + + log2 q + l)tcomp + 2log2 P(? + P Usually, tcomm > tcomp and T is several orders of magnitude greater than tcomm, Hence up until tcomp 2 log2 Ptcomm > p it is always better to use more processors. This is because an increased number of processors increases parallelism. Looking at the entries in Table 7.1, we can see that the local FFT portion, the computation of twiddle factors, and the factoriza- tion of m reflects linear speedup with the doubling of P, i.e., the timings for 2P processors is almost exactly half of the timings for P processors. Holding n and P constant, one might ask the question of how m and q should be chosen. That is, should we always choose the smallest q? The results in Table 7.1 suggest that there should indeed be the tendency to do this, although the timings between q = 256 and q = 16 differ by only io% for P = 16. Our timings are done on a hypercube without vector nodes. However we can conjecture that a smaller q would result in longer vectors of length m to be used during the distributed FFTs. In our implementation however, no use of this characteristic is made. Hence our timings show a slight increase for the timings of the smallest q value because of the increase in overhead for the factorization of m, m becoming more composite as q decreases. If this factorization time is not accounted for, i.e. the factorization of m is known and stored, then timings should decrease as q decreases. References L. Adams & T. Crockett [1984], "Modeling algorithm execution time on pro- cessor arrays," Computer 17, 38?3. H. Ahmed, T. Natarajan & K. R. Rao [1974], "Discrete cosine transform," IEEE Trans. Cornput. C-23, 90--H93. It. Alt [1978], "Error propagation in Fourier transforms," Math. Comp. Simul. 20, 37?3. G. D. Bergiand [July 1969], "A guided tour of the fast Fourier transform," IEEE Spectrum. S. Bertram [Mar. 1970], "On the Derivation of the Fast Fourier Transform," IEEE Transactions on Audio and Electroacoustics AU-18, 55--H58. L. N. Bhuyan & D. P. Agrawal [Apr. 1984], "Generalized Hypercube and Hyperbus Structures for a Computer Network," lEEE Trans. Cornput. C- 33, 323--H333. It. E. Blahut [1985], Fast Algorithms for Digital Signal Processing, Addison- Weslev Readin?. Massach ti P. Bois & J. Vignes [1980], "Software for evaluating local accuracy in the Fourier transform," Math. Comp. SimuL 22, 141--H150. R. N. Bracewell [1986], The Hartley Transform, Oxford University Press, New York. W. L. Briggs [July 1987], "Further Symmetries of In-Place FFTs," SlAM J. Sci. Statist. Cornput. 8, 644654. E. 0. Brigham [1974], The Fast Fourier Transform, Prentice-Hall, Englewood Cliffs, New Jersey. E. 0. Brigham & It. E. Morrow [Dec. 1967], "The fast Fourier Transform," IEEE Spectrum. 0. Buneman [1987], "Stable On-Line Creation of Sines or Cosines of Successive Angles," Proc. IEEE, (to appear). 211 212 B. Buzbee [1973], "A fast Poisson solver amenable to parallel computation.? IEEE Trans. Comput.C-22, 793--H796. R. M. Chamberlain [May 1986], "Gray Codes, Fast Fourier Transforms and Hypercubes," Chr. Michelsen Institute, CCS 86/1, Bergen, Norway. 0. W. C. Chan & E. I. Jury [Jan. 1974], "Roundoff Error in Multidimensional Generalized Discrete Transforms," IEEE Trans. Circuits and Systems CAS 21, 100--H108. T. F. Chan [Sept. 9, 1986], "On Gray Code Mappings for Mesh-FFTs on Binary N-Cubes," NASA Ames Research Center, RIACS Technical Report 86.17. W. Chen, C. H. Smith & 5. Fralick [Sept. 1977], "A Fast Computational Algorithm for the Discrete Cosine Transform," lEEE Trans. Comm. COM- 25, 1004--H1009. R. J. Clarke [1985], Transform Coding of Images, Academic Press, New York. W. T. Cochran, J. W. Cooley, D. L. Favin, H. D. Helms, R. A. Kaenel, W. W. Lang, G. C. Maling, Jr., D. E. Nelson, C. M. Rader & P. D. Welch [Oct. 1967], "What Is the Fast Fourier Transform?," Proceedings of the IEEE55, 1664--H1677. T. Coleman & C. F. Van Loan [1987], A Matrix Computation Handbook, (manuscript). D. 5. P. Committee [1979], Programs for Digital Signal Processing, IEEE Press, New York. J. W. Cooley, P. A. W. Lewis & P. D. Welch [1970], "The fast Fourier transform algorithm: Programming considerations in the calculation of sine. cosine and Laplace transforms," J. Sound Vibration 12, 315--H337. [June 1969], "The Finite Fourier Transform," IEEE Transactions on Audio and Electroacoustics AU-17, 77--H85. J. W Cooley & J. W. Tukey [Apr. 1965], "An algorithm for the machine calculation of complex Fourier series," Math. Cornput. 19, 297--H301. J. W. Cooley, P. A. W. Lewis & P. D. Welch [June 1967], "Historical Notes on the Fast Fourier Transform," IEEE Transactions on Audio and Elec- troacoustics AU-15, 76--H79. [Mar. 1969], "The Fast Fourier Transform and Its Application," IEEE Transactions on Education 12, 27--H34. M. J. Corinthios [June 1971], "The design of a class of fast Fourier transform computers," IEEE Trans. Cornput. C-20, 617623. 213 J. Dollimore [1973], "Some algorithms for use with the fast Fourier transform.? J. Inst. Math. Appi. 12, 115--H117. M. Drubin [May 1971], "Kronecker Product Factorization of the FFT Matrix.? IEEE Trans. Cornput.. D. E. Dudgeon & R. M. Mersereau [1984], in Multidimensional Digital Signal Processing, Prentice-Hall, Inc., Englewood Cliffs, NJ. D. B. Gannon & J. Van Rosendale [Dec. 1984], "On the Impact of Communi- cation Complexity on the Design of Parallel Numerical Algorithms," IEEE Trans. Comput.C-33, 1180--H1194. M. Gentleman [1975], "Error Analysis of QR Decompositions by Givens Transformations," Linear Algebra Appl. 10, 189--H197. W. M. Gentleman & G. Sande [1966], "Fast Fourier Transforms--HFor Fun and Profit," in 1966 Fall Joint Computer Conference, AFIPS Con Proceedings #29, Spartan, Washington, D.C., 563--H578. M. W. George, R. T. Ling, J. F. Mangus & W. T. Tbompkins [Mar. 10-12, 1987], "Application of Computational Physics Within Northrop," Supercomputing in Aerospace. W. J. Gilbert [1976], Modern Algebra With Applications, John Wiley and Sons, Inc., New York. J. A. Glassman [Feb. 1970], "A generalization of the fast Fourier transform.? IEEE Trans. Comput.C-19, 105--H116. G. H. Golub & C. F. Van Loan [1983], Matrix Computations, The Johns Hopkins University Press, Baltimore, MD. I. J. Good [Mar. 1971], "The Relationship Between Two Fast Fourier Trans- forms," IEEE Trans. Cornput.. J. Gustafson [Mar. 22-25, 1987], "Ensemble FFT's as a Function of Compute Communication Ratios," in Fast Fourier Transforms for Vector and Par- allel Computers Workshop, Charles F. Van Loan, ed., The Mathematical Sciences Institute, Cornell University, Ithaca, NY. D. B. Harris, J. H. McClellan, D. 5. K. Chan & H. W. Schuessler[1977]. "Vector Radix Fast Fourier Transform," Rec. 1977 IEEE Internat. CoL Acoust., Speech, Signal Proc.. L. 5. Haynes, R. L. Lau, D. P. Siewiorek & D. W. Mizell [Jan. 1982], ?A Survey of Highly Parallel Computing," Computer. C-T. Ho & 5. L. Johnsson[Sept. 1986], "Matrix Transposition on Boolean n- cube Configured Ensemble Architectures," Yale University, Department of Computer Science, YALEU/DCS/TR-494, New Haven, CT. 214 R. Hockney & C. Jesshope [1981], in Parallel Computers: Architecture, Pro- gramming and Algorithms, Adam Hilger, Bristol. A. K. Jain [1976], "A fast Karhunen-Loeve transform for a class of stochastic processes," IEEE Trans. Comm. COM-24, 1023--H1029. [1979], "A sinusoidal family of unitary transforms," IEEE Trans. PAMI- 1, 356--H365. 5. L. Johnsson[Feb. 1985a], "Data Permutations and Basic Linear Algebra Computations on Ensemble Architectures," Yale University, Department of Computer Science, YALEU/CSD/RR-367, New Haven, CT. [Jan. 1985b], "Communication Efficient Basic Linear Algebra Compu- tations on Hypercube Architectures," Yale University, Department of Corn- puter Science, YALEU/CDS/RR-361, New Haven, CT. D. K. Kahaner [Dec. 1970], "Matrix Description of the Fast Fourier Transform," IEEE Transactions on Audio and Electroacoustics AU-18, 442450. T. Kaneko & B. Liu [Oct. 1970], "Accumulation of Round-Off Error in Fast Fourier Transforms," J. Assoc. Cornput. Mach. 17, 637654. H. B. Kekre & J. K. Solanki [1978], "Comparative performance of various trigonometric unitary transforms for transform image coding," Internat. J. Electron. 44, 305--H315. H. Kitajima [Apr. 1980], "A Symmetric Cosine Transform," lEEE Trans. Com- put. C-29, 317--H323. W. R. Knight & R. Kaiser [Dec. 1979], "A simple fixed point error bound for the fast Fourier transform," IEEE Trans. Acoust. Speech Signal Process. ASSP-27, 615620. C-C. J. Kuo, B. C. Levy & B. R. Musicus [July 1987], "A Local Relaxation Method for Solving Elliptic PDEs on Mesh-Connected Arrays," SlAM J. Sci. Statist. Cornput. 8, 550--H573. B. Lui & T. Kaneko [1975], "Roundoff error in fast Fourier transforms (decima tion in time)," Proc. IEEE 63, 991--H992. R. E. Lynch, J. R. Rice & D. H. Thomas [1964], "Direct solution of partial difference equations by tensor product methods," Nurnerische Mathematik 6, 185--H199. J. Makhoul [Feb. 1980], "A Fast Cosine Transform in One and Two Dimensions," IEEE Trans. Acoust. Speech Signal Process. ASSP-28, 27-34. 0. A. McBryan & E. F. Van de Velde [Mar. 1987], "Hypercube Algorithms and Implementations," SlAM J. Sci. Statist. Comput. 8, s277?287. 215 D. Miles, P Kinney, J. Groshong & R. Fazzari [1987?, "Specification and Performance Analysis of Six Benchmark Programs for the FPS T Series," Floating Point Systems, Inc, P.0. Box 23489, Portland, OR. D. c. Munson, Jr. & B. Liu [1981j, "Floating Point error bound in the prime factor FFT," IFEE Trans. Acoust. Speech Signal Process. ASSP-29, 8?'7- 882. H. J. Nussbaumer [1982j, in Fast Fourier Transform and Convolution Algo- rithins, Springer Verlag, Berlin, Heidelberg. J. Oliver [197sj, "Stable Methods for Evaluating the Points cos(ir/n)," J. Inst. Maths Applics 16, 247--H257. A. V. Oppenheirn & 0. J. Weinstein [Aug. 1972j, "Effects of Finite Register Length in Digital Filtering and the Fast Fourier Transform," Proceedings of the IEEE6O, 957--H976. M. c. Pease [Apr. 1968j, "An Adaptation of the Fast Fourier Transform for Parallel Processing," J. Assoc. Coinput. Mach. is, 252--H264. I. Pitas & M. G.Strintzis [1983j, "Floating point error analysis of two effi- cient algorithins used in the fast computatio? of tw?dimensional DFT5," in Signal Processing II: Theories and Applications, H. W. Schussler, ed., North.Holland, Amsterdam, 219--H222. 5. Prakash & V V Rao tOct. 1982j, "Vector Radix FFT Error Analysis," IEEE Trans Acoust. Speech Signal Process. ASSP-30, 808--H811. W. K. Pratt [1978j, Digital Image Processing, Wiley, New York. W. II Press, B. P. Flannery, 5. A. Teukolsky & W. T. Vetterling [1986j, Numerical Recipes: The Art ofScientific Computing, Cambridge University Press, Cambridge. L. R. Rabiner [1979aj, "FFT Subroutines for Sequences with Spcial Properties," in Program5 for Digital Signal Processing, Digital Signal Processing Corn- mittee, ed., IEEE Press, New York, 1.3--H1--H1.3?22. [June 1979bj, "On the Use of Symmetry in FFT Computation," IEEE Trans. Acoust. Speech Signal Process. ASSP-27, 233--H239. G. U Ramos [Oct. 1971J, "Roundoff Error Analysis of the Fast Fourier Trans- form," Math. Comp. 25, 757--H768. G. K. Rivard [June 1977j, "Direct Fast Fourier Transform of Bivariate functions," IEEE Trans. Acoust. Speech Signal Process. ASSP-25, 250-52. Y. Saad & M. H. Schultz [June 1985j, "Topological Properties of Hypercubes," Research Report YALEU/DCS/RR.389. 216 [Mar. 1986], "Data Communication in Parallel Architectures," Research Report YALEU/DCS/RR-461. [Sept. 1985], "Data Communication in Hypercubes," Research Report YALEU/DCS/RR-428. C. L. Seitz [Jan. 1985], "The Cosmic Cube," Comm. ACM28, 22--H33. R. C. Singleton [Oct. 1967], "On computing the fast Fourier transform," Comm. ACMlO, 647654. H. Sloate [Jan. 1974], "Matrix Representations for Sorting and the Fast Fourier Transform," IEEE Trans. Circuits and Systems CAS-21, 109--H116. II. 5. Stone [Feb. 1971], "Parallel processing with perfect shuffle," IEEE Trans. Comput. C-22, 153--H161. P. N. Swarztrauber [1977], "The methods of sysclic reduction, Fourier analysis and the FACR algorithm for the discrete solution of Poisson 5 equation on a rectangle," SIAM Rev. 19, 490--H501. [1982], "Vectorizing the FFTs," in Parallel Computations, G. Rodrigue, ed., Academic Press, New York, 51?1. [1984a], "FFT algorithms for vector computers," Parallel Comput. 1, 4563. [1984b], "Fast Poisson solvers," MAA Studies in Numerical Analysis 24, 319--H370. P. N. Swarztrauber [1986a], "Multiprocessor FFT's," National Center for At- mospheric Research, Boulder, CO, (to appear in Parallel Comput.). [July 1986b], "Symmetric FFTs," Math. Comp. 47, 323--H346. C. Temperton [1979], "Direct Methods for the Solution of the Discrete Poisson Equation: Some Comparisons," J. Comput. Phys. 31, 1-20. [1980], "On the FACR(l) Algorithm for the Discrete Poisson Equation," J. Comput. Phys.34, 314-329. [1983], "Fast mixed?radix real Fourier transforms," J. Comput. Phys. 52, 340--H350. F. Theilbeirner [June 1969], "A Matrix Version of the Fast Fourier Transform," l'EEE Transactions on Audio and Electroacoustics AU-17, 158--H161. T. Thong & B. Liu [1977], "Floating point fast Fourier transform computation using double precision floating point accumulators," ACM Trans. Math. Soft. 3, 5?59. 217 N-K. Tsao [(no date given)], "The "Equivalence" of Decimation in Time and Decimation in Frequency in FFT Computations," Department of Computer Science, Wayne State University, Detroit, MI. C. F. Van Loan [1987], FFT's From the Matrix Point of View, (manuscript). S. R. Walton [Sept. 1986], "Fast Fourier Transforms on the Hypercube," Ametek Computer Research Division, Arcadia, CA. Z-D. Wang [1981a], "Harmonic Analysis with a Real Frequency Function I. Aperiodic Case," Appi. Math. Cornput. 9, 53--H73. [1981b], "Harmonic Analysis with a Real Function of Frequency II. Periodic and Bounded Cases," Applied Mathematics and Computation 9. 153--H163. [1981c], "Harmonic Analysis with a Real Frequency Function. III. Data Sequence," Applied Mathematics and Computation 9' 245-255. [Aug. 1984], "Fast Algorithms for the Discrete W Transform and for the Discrete Fourier Transform," IEEE Trans. Acoust. Speech Signal Process. ASSP-32, 803?16. Z-D. Wang & B. R. Hunt [1985], "The Discrete W Transform," Appi. Math. Cornput. 16, 19?. c. j. Weinstein [Sept. 1969], "Roundoff Noise in Floating Point Fast Fourier Transform Computation," IEEE Transactions on Audio and Electroacous- tics AU-17, 209--H214. P. D. Welch [June 1969], "A Fixed-Point Fast Fourier Transform Error Analysis," IEEE Transactions on Audio and ElectroacousticsAU-17, 151-157. P. Wiley [June 1987], "A parallel architecture comes of age at last," IEEE Spec- trum 24, 46--H50. P. Yip & K. R. Rao [1980], "A fast computational algorithm for the discrete sine transform," IEEE Trans. Comm. COM-28, 3O?31O. (1987], "On the Shift Property of DCTs and DSTs," lEEE Trans. Acoust. Speech Signal Process. ASSP-35, 4o4?O6.