Sec. 8-7 Extensions and Scientific Applications 

Extending the k-sums to infinity we find: 
T calc -n (1 + 112 + 114 + ...) = 2n 

and 

151 

(8-13a) 

T comm -n 112(1 + 1/--12 + 1/2 + ...) = 3.4n 1/2 (8-14a) 
In other words, we obtain a value for the communication overhead f c that is increased by a factor of 1.7 as compared with Eq. (8-10), but is still proportional to n-1/2. Here n is the grain size on the finest grid. The presence of the coarse grids with less concurrency only gives rise to the factor 1.7 and does not qualitatively change the efficiency behavior. The point is the same as that discussed in Section 3-6; the stages involving coarse grids (small nk) have high overheads but they don't take very long, so that the net performance is dominated by the high efficiency stages on the fine grid. An improved multigrid for- malism for the hypercube with less ovefhead on the coarse grids is presented in [Freder- ickson 87]. 
Our discussion of finite element problems in this chapter has tended to emphasize two-dimensional grids with a comparatively regular and simply connected structure. Such simple grids aid in the explanation and visualization of interprocessor communica- tion paths. A critical issue, however, in the applicability of these ideas to finite element problems of practical importance is their generalization to higher-dimensional and irregu- lar grids. 
Consider first the extension to three dimensions. As shown in Fig. 8-6, whereas the two-dimensional protocol requires keeping track of eight distinct "types" of communi- cated degrees of freedom, in three dimensions there are 26 types of face, edge, and comer nodal points which must be shuttled between processors. This increase in com- plexity alone might be enough to suggest that simple mesh communication schemes are cumbersome in three dimensions. In fact, such schemes are not difficult to implement until irregular grids are considered. 
As long as the modeled domain can be mapped by a topologically rectangular decomposition with no holes or niches, all shared degrees of freedom in a three- dimensional grid are guaranteed to reach their proper destinations after six communica- tion steps (one in each spatial coordinate with each sign). This elegant method of routing edges, which is robust in a simply structured domain, forms a central part of the particle dynamics algorithm discussed in Chap. 16. Complications occur, however, when we consider more general kinds of finite element domains. 
As we pointed out at the beginning of the chapter, a principal advantage of the finite element method is its facility in handling irregular domains. Thus a practically use- ful concurrent finite element program should not be limited to simple grids. It is a difficult programming task, however, to design a neighbor-passing procedure sufficiently general to handle every possible contingency of routing around "holes" in a complicated grid. General message routing systems provide the solution to this problem. 
Non-nearest neighbor communication as provided by the crystal-router utility (Chap. 22) or by interrupt-driven systems (Sec. 14-8) can provide the robust 

152 

Elliptic Problems in Two Dimensions ll: The Finite Element Method Chap. 8 

communications needed for general finite element problems, though at some additional communication cost. In such a scheme, the original distinction between "accumulated" and "sent" degrees of freedom is retained, but explicit specification of routing directions with relative addresses such as "right" and "up" is abandoned. Instead, communicated degrees of freedom are accompanied by absolute addresses corresponding to their desti- nation processors(s) and the routing system takes care of the rest. 
Of course, the decomposition and routing utilities in such a system should be designed in order to take best advantage of the direct communication channels of a par- ticular machine. As long as this is true, the performance of the finite element program with routed communications should not be much worse than optimal, since the number of nodal points requiring long-distance forwarding willbe small. 

Distinct Types of Communicated Nodal Points in 2 and 3 Dimensions 

~ 

..- 

Two dimensions: 
8 communicated types 

+Et=:f=:=f? 

Three dimensions: 
26 communicated types 

l'lgure 8-6 Comparisoo of the communicated degrees of freedan in two- and three- dirnensiooal finite element irnplernentations. 
Another problem of importance in the concurrent finite element field is the automa- tion of the domain decomposition task. While this issue does not deal directly with con- currency in the finite element solution process, it is important because it determines how much more difficult for the end user the concurrent program will be to use than its sequential equivalent In formulating automated domain decomposition procedures, the communication protocol and data organization should be designed in such a manner as to minimize the user's awareness of the particular concurrent processor's communication restrictions. Significant progress has been made in the context of a two-dimensional 

Sec. 8-7 Extensions and Scientific Applications 

153 

example problem [Flower 86, Nour-Omid 87b] and powerful techniques are reviewed in Sec. 3-8 [Fox 86g]. We expect that the tasks of domain decomposition and mesh genera- tion may well be combined in future systems. The emerging generation of computers will be capable of handling the extremely large finite element problems for which automatic decomposition will be a necessity. 

9 

9-1 Introduction 
In this chapter we consider algorithms in which a nontrivial computation must be performed on all pairs of objects in a data base. Many situations arise in which this type of calculation is necessary .In molecular dynamics calculations, such an algorithm is needed when the intermolecular force law has a Coulombic or other long-range com- ponent A similar situation arises in Newtonian gravitational n-body simulations, in which the force law also has infinite range. Computing the product of two large polyno- mials is also a long-range problem because every pair of monomials must be combined to produce a term in the result By definition, a long-range algorithm is characterized by a running time proportional to the number of pairs, or roughly, the square of the number of objects. Thus, long range problems with more than a few thousand objects pose significant challenges for even the fastest computers. 
In order to present the important components of the long-range force algorithm without obscuring the ideas with unnecessary detail, this chapter is divided into several sections. Each of these sections presents a successively more complicated version of the algorithm. In this way, the most fundamental ideas are developed relatively cleanly, and the progression to more complex problems is accomplished by adding features to the fun- damental framework. In Sec. 9-2, we consider the problem of computing the potential energy of a collection of particles which interact by a pairwise potential. This subprob- lem disregards Newton's third law. In Sec. 9-3, we take advantage of Newton's third law to reduce the necessary work by half, as is normally done in the sequential case. In Sec. 9-4, we compute a unique force for every particle in the system, and in Sec. 9-5, we comment on the application of the algorithm to the calculation of correlation functions. 

9-2 Total Potential Energy 
In order to inb"oduce the basic communication Sb"ategy and to outline the approach, we begin with a computation of the total potential energy of a system of particles that interact via a pair potential. This is very much simpler than computing the forces on all the particles because the final result will be a single number rather than a different force for every particle. For the rest of this chapter the number of processors is denoted by N , and the number of particles is N p .A particle is an aggregate data structure containing all the information necessary to describe the state of a physical particle in the simulation. For example, mass, charge, position, velocity, acceleration, etc. 

155 

156 

Long Range Interactions Chap. 9 

The total potential energy of a system of such particles is given by 
NN 1 p p 
U = 2LL<!>(Xi,Xj) (9-1) i=lj=l 
In Code 9-1, the function pot-seq computes the double summation in Eq. (9-1) for an array of particles. We use the auxiliary function pot-2arr, which uses two arrays of particles, in anticipation of the requirements of the parallel algorithm. The function pair yotential contains the physics of the computation, and returns the potential energy contribution, <!>, of its two arguments. For example, pair yotential might compute the Newtonian potential of two massive objects. 

<I>(.t; , ~ ) = - 

Gm.m. I J 
IX'; -~ I 

Code 9.1 Sequential Evaluation of Potential Energy 

9-2.1 Parallel Decomposition 
First we observe that the positions of all particles must eventually be communi- cated to every processor. If there are Np particles, and Nlocal are stored in some proces- sor, then that processor must somehow obtain the attributes of N p -N local particles and compute N local N p pairwise interactions. The time spent in computation is thus: 
Tcomp :: max(Nlocal) Np tpair (9-3) 
where tpair is the time required to compute the interaction potential of a pair of particles. We take an optimistic view and assume that we can do the communication without 

Sec. 9-2 Total Potential Energy 

157 

bottlenecks and that all data transfers are equally fast Then the number of attributes read by any given processor provides an estimate of the communication time 
Tcomm :::: max(Np-Nlocal) texchg (9-4) 
where texchg is the time required to communicate the attributes of one particle to a neigh- boring processor. In Eqs. (9-3) and (9-4) max(x) refers to the largest value of x in any processor in the system. If texchg :::: tpair' then From Eqs. (9-3) and (9-4) we see that com- putation is more time consuming than communication, by a factor which grows approxi- mately linearly max(N local ). Thus, the optimum decomposition will minimize T comp and will be such that max(N local) is minimized. In other words, an equal number of particles is stored in every processor. By definition, there is no locality in long-range problems, so the decomposition need not be concerned with the proximity of particles in space. All decompositions which balance the load are equally good. Assigning particles to proces- sors at random, with the constraint that an equal number is assigned to each processor, is a perfectly good way to decompose long range problems. 
Unfortunately, the number of particles is often fixed by considerations other than algorithmic convenience, and this number might not be divisible by the number of pro- cessors. In any case, the difference in Nloca/ between any two processors can be made to be either one or zero. Our subsequent discussion assumes that all processors have exactly the same number of particles. The more general case is essentially the same, and can be treated with minimal cost by placing dummy particles in some processors. 

9-2.2 Communication Structure: The Ring Topology 
The discussion in the previous section relied on the optimistic assumption that the communication could be done without bottlenecks. In this section we discuss a very sim- ple interconnection topology which verifies that assumption. 
In order to evaluate all of the pair potentials, every pair of particles must meet in some processor at some time in the calculation. In other words, every particle must visit every processor, and while there, its interactions with all the particles in that processor must be tallied. We propose to connect the processors in a ring, as in Fig. 9-1. Every particle can visit all the processors in the ring by taking N -1 steps in either the clockwise or counterclockwise direction. Crucial to the practicality of this approach is that the communications can be done concurrently. Therefore, we can think of the ring as a pipe- line. In a single time interval of duration, texchg , it is possible to transfer N particles one step around the ring, with each processor transmitting and receiving exactly one particle. 
Now we ask: How much time does it take, in communication alone, to have every particle visit every processor? Since the communication is done concurrently, this inter- val is equal to the time required for a single processor to transmit and receive all N p par- ticles. This time is simply N p texchg .This estimate, however, is based upon the additional assumption that no processor will finish its own work and fall idle waiting for its neigh- bor to send more particles. If all processors contain the same number of particles at all times, then such a load imbalance will not occur. With equal numbers of particles in each processor, the computational load between communication steps is exactly the same 

158 

(b) 

Long Range Interactions Chap. 9 

'- 

u 

-. 

:0 

/ 

~ 

, 

~ 

; 

~ XJ 

I'lgure 9.1 (a) A single object, objl' can visit each of four processors by following the path denoted by the arrows, (b) By perfonning the communicatioo in parallel, four ob- jects, objl through obj4' can each visit the same four processors in the same amount of time. In each step, each processor writes its object in the directioo of the outgoing arrow and reads a new object fran the directioo of the incoming arrow. After only three steps (three units of time), all four objects have visited all four processors. 
in all processors, so dlere need be no waiting for communication, As soon as one proces- sor is ready to read, its neighbor should be ready to write. 
In considering dle parallel algoridlm, we have already concluded dlat we must move particles from processor to processor, computing dle interactions between two arrays of particles, dlose which started in a given processor, and dlose which are passing through. We can use pot-2a" from Code 9-1 to compute dle interactions between dlese two arrays. In Code 9-2, pot yar is a procedure dlat computes dle potential in parallel. 
The procedure pot yar contains three procedure calls and a loop not found in pot-seq. The first procedure is ring-setup which maps dle abstract ring topology, in which each processor has a clockwise and counterclockwise neighbor, onto dle virtual machine. Since a ring is really a one-dimensional grid, ring-setup can obtain dlis infor- mation from dle decomposition routines described in Sec. 4-6. The ring is closed by con- necting dle processors at locations I and N . 
The for loop in dle function pot yar , which runs over each processor in dle ring, contains a call to rotate, which is dle basic communication operation for dlis algoridlm. Rotate moves all of dle particles in dle array, Q, around dle ring by an amount of one 

Sec. 9-2 Total Potential Energy 

159 

processor in me clockwise direction. The contents of me array are replaced by me parti- cles obtained from me counterclockwise neighbor. Rotate is easily implemented using vm-shift of Sec. 4-4. 

Code 9.2 Parallel Evaluation of Potential Energy 
Finally, we find the subroutine global-add. Before the call to global-add, each processor has accumulated only a fraction of the total of all pair interactions. To obtain the total, we must compute an overall sum of all the partial sums stored in each proces- sor. The function vm-combine , also introduced in Sec. 4-4, is used in global-add to sum the contributions from each processor. 

9-2.3 Efficiency of the Algorithm 
Problems involving long range forces seem to be ill-suited to parallel processing. After all, isn't the communication overhead extremely high, since every particle must visit every processor? The answer is that algorithms for long-range problems do, indeed, require considerable communication. The computational time, however, is so much larger that the communication overhead is relatively unimportant. The quantity of interest is not the total time spent communicating, but rather the speedup obtained by use of the concurrent processor. 
Let T N (N p ) be the time to compute the potential energy in an N p -body problem on N processors. The overhead, f , of the implementation has been defined in Eq. (5-5) as: 
NT 
(N 
)-T(N 
)f = N p 1 p (9-5) 

T l(Np) 
The parallel program pot yar has a running time given by: 

160 

Long Range Interactions Chap. 9 

TN(Np) ::: N (tpot-2ano + trot) (9-6) 
where tpot-2ano is dIe time to complete a call to pot-2arr and trot is dIe time to complete a callto rotate. These are expressed in terms of texc/18 and tpoir as 
2 

2 
[Np 
Jtpot-2arr =Nlocaltpair = N tpair 

and 

Np 
trot :::Nlocaltexch8 = Ntexch8 

The overhead approaches zero if 

(9-8) 

Nlocal tpair 
Thus, the algorithm win perform efficiently as long as Nlocal is large, i.e., as long as the amount of computation assigned to each processor is large. Increasing the number of processors with a fixed problem size reduces the efficiency, but the perhaps more likely scenario in which problem size grows with machine size leaves the efficiency approxi- mately constant. 

1 

~«1 

(9-10) 

9-3 Symmetry of the Pair Potential 
The discussion. until now. has overlooked one crucial fact that the pair potential function <I> often satisfies the symmetry relation: 
<I>(Xl.XV=<1>(X2.XJ (9-11) 

In dynamics problems, Eq. (9-11) is a consequence of Newton's Third Law. Even in nondynamical problems, such as the computation of correlation functions, a similar sym- metry can be used to reduce the required amount of computation by half. Unless we can use Eq. (9-11) in the parallel algorithm, it will be impossible to achieve a true efficiency greater than 50%, because the parallel algorithm will always do twice as much work as the optimal sequential algorithm. It is easy to incorporate the symmetry of Eq. (9-11 ) into the sequential program. The subroutine pot- sym -2arr illustrated in Code 9-3 makes use of the symmetry , and may replace pot-2arr in Code 9-1. 

Sec. 9-3 Symmetry of the Pair Potential 

161 

Code 9.3 Sequential Evaluation of Potential Energy Exploiting SymmetlY of Potential Functioo 
The subroutinepot-sym-2arr is called by a sequential program in exactly the same way as pot-2arr .The idea here is to take each particle, p i, in array p and compute its interaction with those particles in Q starting at i + 1. All pairs interact exactly once in this way, and self-interactions are never computed. 
In the parallel program, it is harder to avoid computing self-interactions. The ith particle in the moving array must interact with the i th particle in the stationary array at least some of the time. If the interaction between two such congruent particles is com- puted every time, however, it will be computed independently in two processors and incorrectly counted twice in the total energy. We must introduce a function, should-self -interact, which sets the flag self -interact to 1 (YES) when the interaction should be computed between two congruent particles and to 0 (NO) otherwise. The func- tion pot-sym yar2 in Code 9-4 is a parallel version of pot- sym-2arr , i.e., it computes the interaction between two arrays of particles in a parallel context, properly accounting for the interactions of congruent pairs. It may replace the call to pot-2arr in Code 9-2. 
Now we turn to the details of should-self-interact, in Code 9-5. Consider two pro- cessors with ring locations Pa and Pb , and assume Pb > Pa. We will compute the interac- tion between congruent pairs in P a and, hence, not in Pb , if (Pb -P a) < N 12. Conversely, if (Pb -P a) > N 12, we compute the interaction in Pb .If (Pb -P a) is equal to N 12, we must use some criterion other than their separation to decide when to compute the self - interaction. Furthermore, to avoid the possibility of idle processors when (Pb -Pa) = N /2, we would like each processor to compute exactly half of the interac- tions. One possibility is to compute the interactions between even-numbered pairs in Pa and to compute those between odd-numbered pairs in Pb . 
Note that this version of should-self -interact does not allow for computation of bona fide self-interactions. It always returns NO for the interaction of a particle with itself, i.e., when step is 0. This is most often desirable, especially in dynamics problems in which the self-energy is irrelevant, and often divergent. Nevertheless, it is sometimes useful to compute the self-interactions. This is accomplished by exchanging YES and NO in the first two if clauses of Code 9-5. 
The discussion of the efficiency in Sec. 9-2.3 applies to the improved algorithm with little modification. We have reduced the computation by half, both in the sequential and parallel cases. Nevertheless, the communication costs have remained constant. 

Long Range Interactions Chap. 9 

Code 9.4 Parallel Evaluation of Potential Energy Exploiting Symmetry of Potential Function 

Code 9.5 The Routine should-se/f -imeract 
Thus, the bulk of the program runs almost twice as fast, but the speedup is slightly reduced. The new overhead is 

f :: 2-N texch!- Np tpair 

(9-12) 

Sec. 9-4 Force Calculations 

163 

9-4 Force Calculations 
We now generalize the algorithm of the previous section to compute and record the force felt by every particle rather than simply to calculate the total potential energy. We assume the physics of the problem is taken care of in a subroutine, pair Jorce, which computes the force between two particles, and adds it to the total force stored in each. 
The parallel program force ~r in Code 9-6 differs from pot yar in Code 9-2, in two ways. First, it calls force- sym yar2 , which is the same as pot-sym yar2 , except that pair Jorce is used to compute the interactions between particles. Second, the function global-add is replaced by the function recombine. In the force calculation, the moving particles have felt the effects of some of the particles, while the stationary array has felt the effects of the others. Clearly, we must bring these two arrays back together in the original processor and combine the results. 
Bringing them back together is easy. At the end of the calculation, after N calls to rotate, the moving array has traveled all the way around the ring and is back where it started. The ring topology has returned the moving data without any extra effort. The combining operation is also simple. It is only necessary to add the force accumulated by each particle in the moving array to that accumulated by the corresponding particle in the stationary array. The procedure recombine achieves this. 

Code 9.6 Parallel Evaluation of the Mutual Force Exerted by a System of Panicles 

9-5 Calculation of Correlation Functions 
One practical use of the ideas in this chapter was motivated by the study of experi- mental data taken from studies of turbulent plasmas in a Tokamak [Zweben 83]. An algorithm proposed by Grassberger and Procaccia for computing the correlation dimen- sion associated with a dynamical system [Grassberger 83] was adapted for use in a paral- lei processor [Theiler 86]. This dimension is given by the power-law behavior of what is called the correlation integral. 

164 

Long Range Interactions Chap. 9 

For a dynamical system with time dependence, x (t ), the correlation integral, C (r ), measures the probability that a given pairwise distance rjj = I x(tj ) -x(t j ) I will satisfy rjj < r. To estimate the correlation integral from experimental data consisting of N o discrete measurements of x(t) requires the construction of a histogram of all N; values of rjj . 
We construct the histogram in parallel as follows. Two copies, one stationary and one moving of N o IN = N local values of ~ are stored at each node of the hypercube. Also, at each node is a separate histogram. As the moving copies travel around the loop, the pairwise distances between the moving and stationary vectors at each node are computed, and the local histograms are updated. After the moving copy of the data has traveled all the way around the ring, the separate histograms in each node must be combined This can be performed either by the sophisticated technique that will be described in Chap. 19 or by using the function vm-combine function discussed in Sec. 4-4. In any case, this last step takes a negligible amount of time, and a high performance concurrent calculation of the correlation integral is achieved. 
Recently the basic sequential correlation algorithm has been improved, reducing the N 6 time complexity to N o log N o [Theiler 87]; this new algorithm has not yet been implemented on a concurrent machine. It is interesting to note that it uses a similar clus- tering technique to that introduced by Barnes and Hut [Barnes 86] for the N p body prob- lem discussed in Secs. 9-1 to 9-4. As discussed in Sec. 23-4.5, this fast Np bodyalgo- rithm is currently being implemented on the hypercube. 

9.6 Conclusion 
The techniques described in this chapter apply to any problem that requires the examination of every pair of objects in a data base. The general outline for treating such problems in parallel is shown in Code 9- 7. 

Code 9.7 General Algoritlun for Examining Every Pair of Items in a Database 
Large problems of this type can be treated very efficiently with parallel computers. Large only means that the number of objects in the data base is much larger, say by a 

Sec. 9-6 Conclusion 

165 

factor of 50 for a typical machine, than the number of processors. The high efficiency of long-range problems is quite surprising, since they are among the most communication intensive of regular problems. Although the amount of communication is indeed large, it still represents only a small fraction of the very costly calculation, and, thus, long-range algorithms perform very efficiently on parallel computers. An additional result is that the necessary communication topology is extremely simple; a ring of processors in which each processor has only two neighbors is all that is needed for an efficient implementa- tion. Objects can be communicated around a ring of processors concurrently, making additional connectivity superfluous. 

10 

10-1 IntrOOuction 
The numerical solution of many problems reduces in part or fully to various matrix operations. Typical of these operations are multiplication of matrices and vectors, inver- sion or LU decomposition, and the determination of eigenvalues and vectors. These operations can be further classified by the characteristics of the matrix. The matrix can be "full" (i.e., with an insignificant number of zero matrix elements), banded (nonzero elements clustered about the diagonal), or sparse with a more complicated structure in the location of the nonzero elements. We find that matrix algorithms generally run well on concurrent computers and that the hypercube, in particular, is well suited to essentially all of these algorithms. In this volume, we will discuss three cases in detail: multiplication of full matrices in this chapter; the LU decomposition of banded matrices in Chap. 20; and matrix-vector multiplication in Chap. 21. Chapter 21 will also survey the key issues in some of the other matrix algorithms. The aim of presenting these three particular topics is to illustrate the typical issues in the use of concurrent computers to solve matrix problems. 
The algorithms in this chapter are intended to be applicable to arbitrary concurrent computers, but as in previous chapters, the explicit performance measurements given in Sec. 10-5 correspond to a specific hypercube implementation. The Software Supplement to this volume contains detailed listings of the hypercube code for several matrix algo- rithms including those covered here and in Chaps. 20 and 21. 
After these preliminary remarks, we turn our attention to the problem at hand: full matrix multiplication which constitutes one of the simplest matrix problems. In Chap. 21, we will discuss the case of matrix-vector multiplication, which might have been expected to be simpler but, in fact, for the most interesting sparse matrix cases, can require sophisticated techniques. 
This chapter treats the following issues related to concurrent matrix multiplication: 10-2, parallel decomposition of the problem; 10-3, the specific case of square subblock decomposition; 10-4, rectangular and column decompositions; 10-5, results and practical considerations; and 10-6, a summary. Section 10-7 concludes with a comment on the use of concurrent processors to calculate spreadsheets. The work covered here is based on the discussion in [Fox 85b] and related research is covered in [Kung 80], [Sameh 85] and [Johnsson 87]. 

167 

168 

Malrix Algoritluns I: Malrix Multiplication Chap. 10 

10-2 Decomposition or the Matrix Multiplication Problem 
As in earlier chapters, the design of our concurrent algorithm is based upon the assumption of each node running a single process, with those processes communicating via messages. Decomposition involves breaking up the underlying data-in this case ini- tial and final matrices-into parts and associating a part of each matrix with each proces- sor. We will find that the algorithm will require substantially more movem~nt of data than the previous examples in this book. 
For the main discussion, we will need only to make use of processor interconnec- tions which define a two-dimensional periodic mesh. This contrasts with the choice made in Chap. 7 for the case of the iterative sparse matrix solution of systems arising from partial differential equations. In that case, the necessary processor topology was governed by the spatial dimension of the underlying physical system. In the present case, the two-dimensional"world" of matrix rows and columns forms the basis for the decom- position. We will examine in Sec. 10-3.2 a part of the algorithm in which the perfor- mance of the concurrent algorithm could be improved by a richer node connectivity than that provided by a two-dimensional mesh. 
We are interested in performing the multiplication: 
C =A .B (10-1) 
where C, A , and B are full M x M matrices. For simplicity, we will only discuss the multiplication of square matrices. 
We begin by assuming that matrices A and B are identically decomposed into sub- blocks as shown in Fig. 10-1. We will demand that the resultant product, C, has the same decomposition so that it could serve as input for further processing steps. The algo- rithm will first be described for the case of square subblocks. The extension of the algo- rithm to rectangular subblocks and their extreme limit, a pure-row or column decomposi- tion, will follow in Sec. 10-4. Among our eventual goals in this chapter will be to show that the choice illustrated in Fig. 10-1 is preferable to rectangular decompositions. 

ADO I ~l 

A1O I All 

A20 I A21 

A30 I A3l 

~2 

Al2 

~3 

A13 

A221 A23 

A32 

A33 

Sec. 10-2 Decomposition of the Malrix Multiplication Problem 

169 

We will assume mat me multiplicand matrices, A and B, already reside in me memory of me ensemble processors, distributed piecewise in me two-dimensional fashion described. Our subsequent performance analyses do not consider me costs asso- ciated wiili me loading of A or B .This is not unreasonable if A and B came from previ- ous processing steps of some larger algorithm within me concurrent processor. On me oilier hand, if one is interested in using me ensemble as a "flow-through" matrix multi- plier, attention must be paid to me potentially significant costs of ensemble I/O. This is a general problem of essentially all parallel computers and can in principle be solved by a parallel I/O disk system or oilier special purpose devices for me rapid parallel loading and dumping of data. The cost of me matrix calculation grows in proportion to M3, while me 1/0 cost is proportional to M2; in principle, boili can be speeded up by a factor that approaches N, me number of processors. The following discussion, however, will analyze me specific problem of direct matrix operations, wiiliout further consideration of me application and hardware dependent issues of data pre- and postprocessing. 

10-3 Square Subblock Decomposition for Matrix Multiplication 

10-3.1 The Basic Algorithm for Square Subblocks 

The fundamental steps of the algorithm are illustrated schematically in Fig. 10-2 for the case of a 4 x 4 mesh of processors. The square submatrices are denoted by hats: e.g., C~ denotes the matrix subblock of Cij for which 0 ~ i < M /4, 0 ~ j < M /4. In gen- eral, CIm is comprised of the elements of Cij for which IM /4 ~ i < (l+l)M /4, kM/4~j «k+1)M/4. Square submatrices can be manipulated as if they were single elements of the matrix. That is 

c/k = ~A111 .link (10-2) II 
where t!te multiplication in the sum is actually matrix multiplication of the square subma- trices, A III and B nk . 
The concurrent matrix multiplication proceeds as follows: Step I: (Fig. 10-2a) 
The diagonal subblocks of A are successively communicated in a horizontal direc- tion so that all processors in the first "row" receive a copy of A 00, all processors in the second "row" receive a copy of All, and so on. In the terminology of Chap. 4, the diago- nal subblocks of A are broadcast along each row of processors. 
Step 2: (Fig. 10-2b) 
The broadcast A subblocks are multiplied by the B subblocks currently residing in each processor and stored in C subblocks. 
Step 3: (Fig. 10-2c) 
The B subblocks of each processor are shifted to the "upper" neighboring proces- sor and replaced by the corresponding B subblock from the "lower" neighbor. The periodic mesh connection is used to "roll" the topmost row into the bottommost row. 

170 

(a) 

If{){)= If 01= Ifo2= Ifo3= 

l~lO= I~ll= 1f.12= 1f.13= 

i f20= I ~21= I ~22= I ~23= 

It3°= it31= 

~32= 1~33= 

c 

Matrix Algoritluns I: Matrix Multiplication Chap. 10 

= 

fjOO I fjOl 

JilO I Jill 

1}20 I 1}21 

1l3O I 1l31 

jJO2 I jJO3 

~121 ~13 

1l22 

jj231 

1}32 

~33 

A 

B 

FIgure 10-2(8) Step ate of the square subblock algorithm for a 4 x 4 machine. See the text for an explanation of each step. 

(b) 

:t;OOtJOO I ~1JOl ~-p;OOf1O21 ~1JO3 

:t;llglOI:t;llgJJ I:t;JJg121:t;JJg13 

~22't20 11'22~211 ..-22..22 I ~22't23 

I ~33n30 I ~3n3J I ~3n32 ! ~33n33 

c 

= 

ADO 

All 

A22 

A33 

ADO i AoO I AoO 

All 

All 

A!J 

A22 

A22 

A22 

A331 A33 

A33 

T 

BoO 

li1O 

B2Q 

1i3O 

fjOl 

All B 

1}21 

1331 

JjD2 

1112 

§22 

1132 

liO3 

fj]3 

B23 

1133 

B 

Figure lO.2(b) Step two of the square subblock algorithm for a 4 x 4 machine. See the text for an explanation of each step. The matrix T holds the block of A utilized at this stage of the algorithm. 
Step 4: (Fig. 10-2d) 
The A subblocks located in the processors one position to the right of the diagonal are now broadcast throughout the corresponding row, as in Step 1. 
Step 5: (Fig. 10-2e) 
The copied A subblocks are again multiplied by the B subblocks, with the results added to the partial results in the C subblocks. The above pattern of rolling blocks and performing block multiplications is continued until the B subblocks have returned to their original processors. At the completion of this cycle, the C subblocks contain the complete product matrix, distributed in the same manner as A and B . 

Sec. 10-3 Square Subblock Decomposition for Matrix Multiplication 

(c) 

= 

171 

c B 
FIgure lO-2(c) Step three of the square subblock algorithm for a 4 x 4 machine. See the text for an explanation of each step. 

A 

(d) 

c 

A 

B 

Figure lO-2(d) Step four of the square subblock algorithm for a 4 x 4 machine. See the text for an explanation of each step. 

172 

(e) 

c 

Matrix Algorithms I: Matrix Multiplication Chap. 10 

= 

AO1 

AI2 

AO1 

AJ2 

A231 A23 

AO1 

.4 

A23 

AO1 

AI2 

A23 

A30 I A30 I A30 I A30 

1110 I 1111 

1}20 I 1}21 

B30 I B3/ 

Boo I BOl 

§12 

1l22 

832 

BO2 

1113 

1}23 i 

833 

1}03 

T 

B 

F1gure lO-2(e) Step five of the square subblock algorithm for a 4 x 4 machine. See the text for an explanation of each step. The matrix T holds the block of A utilized at this stage of the algorithm. 

Code 10-1 Matrix Multiplicatioo Algorithm 
The pseudocode for this algorithm is contained in Code 10-1 

10-3.2 Performance Analysis for Square Subblocks 
Suppose that we multiply an M x M malrix on an {jif x {jif array of processors in the manner described in the previous section. The size of the subblock in each processor is m x m. where m = M /{jif .The "broadcast-multiply-roll" cycle of the algorithm is repeated {jif times. For each cycle we need to broadcast an A submatrix across each 

Sec. 10-3 Square Subblock Decomposition for Matrix Multiplication 

173 

row of processors. It is implemented using the broadcast collective communication rou- tine introduced in Chap. 4. This implementation is illustrated in Code 10-2. 

Code 10-2 Naive Broadcast Algorithm 
On some actual machines, including Mark II hypercubes, this approach is not the most effective since such a broadcast requires time of order m1{fl-l) tcomm , and the initiating node in each row constitutes a sequential bottleneck. This estimate assumes that although the broadcast subblocks are identical, they must be separately sent to each destination node. On more advanced machines (e.g., the Caltech/JPL Mark III hyper- cube) provision may exist for hardware assist in such a broadcast to a subcube of the hypercube ensemble. However, we shall here consider a simple algorithm for hardware which has no broadcast enhancement to improve upon the above theoretical performance. We will instead set up a pipeline as shown in Fig. 10-2(d) with the originating node on each row just sending the subblock to one other node. Here we need one and not {fl -1 copies of each message to be made in the originating node. This subblock is then succes- sively forwarded {fl -1 steps along the entire row. Such an algorithm requires a time for each cycle given by: 

m2 tcomm + ({fi/ -2) tstart (10-3) 
Here, as usual, tcomm is the time required to pass a floating-point word between pro- cessors and tstart is the startup time of the pipeline per step of the pipe. Typically, tstart is comparable in value to tcomm and the first term in Eq. (10-3) dominates as usually m2 » ({fi/ -2). In Chap. 20, we will study an example in which the tstart term is rela- tively more important as one is communicating fewer (m not m1 numbers, and the pipe- line startup time will be apparent in the measured hypercube performances. We emphasize that the exact form of Eq. (10-3) depends on the details of the hardware, and it is written for the special case of the Caltech Mark n hypercube. This transmits messages in small (64-bit) packets with low latency. Further it takes approximately equal times (1/2 tcomm) to read or write a message. In this case, tstart -2tcomm (32-bit) = tcomm (64-bit) is the time taken to read and write a single packet 
This modified algorithm for subblock pipelining is illustrated in Code 10-3. 

174 

Matrix Algorithms I: Matrix Multiplication Chap.10 

Code 10.3 Pipe Broadcast Algoridun 
Optimal broadcast algorithms on particular concurrent machines are important and surprisingly nontrivial. Two papers studying them for dIe hypercube architecture are [Fox 86b] and [Ho 86a]. One can reduce the coefficients of both tcomm and tSlart but these refinements are not important unless m is very small or tcomm very large. This issue will be discussed in more detail in Chap. 21. 
The time taken to roll B is: 

m2 tcomm (10-4) 
Note that in both Eqs. (10-3) and (10-4), we have assumed a negligible latency for initiat- ing a message between two processors, as is actually the case for such machines as the Mark II Caltech hypercube. We will comment on this assumption at the end of Sec. 10- 6. 

Finally, the time required to compute the submatrix product for contribution to C 

is: 

2m3 tcalc (10-5) 
where, as usual, tcalc is the time required for a floating-point multiply or add. Therefore, the total computation time of the matrix multiplication is given by: 
T =..IN"" (2m3tcalc +2m2tcomm+(..JN"" -2)tstart) (10-6) 

or 

-liT +-/N" (-/N"-2)tstart \J.v-,} N -./N 
The first term of (10-7) corresponds to a concurrent speedup of magnitude N while the last two terms represent communication overhead. Before presenting actual timing results for this algorithm, we will generalize the discussion to the case of rectangular sub- blocks and its extreme limit of column decomposition. 

T= 

2M3tcalc 

2 2M tcomm 

+ 

Sec. 10-4 Rectangular Subblock Decomposition for Matrix Multiplication 

175 

10-4 Rectangular Subblock Decomposition for Matrix Multiplication 

10-4.1 The Rectangular Block Algorithm 

We now extend dte previous algoridtm to dte case in which dte submatrices in each processor are no longer square. The situation is illustrated in Fig. 10-3. An M x M matrix is decomposed onto an N o x N 1 processor array where, for simplicity of presenta- tion, we restrict ourselves to dte case N o ~ N 1. The submatrix in each processor is nio x ni 1 where nio = M IN o and ni 1 = M IN 1. Furthermore, we assume dtat dte aspect ratio, r =niolnil' is an integer, so that an integral number of square subblocks of size ni 1 x ni 1 fits in each processor. The issue of padding a matrix which may not perfectly fit widtin dtis decomposition is discussed in Sec. 10-5. 

r. 

l'lgure 10.3 An M x M matrix decomposed onto a 2 x 4 processor mesh (No=2,N1=4). mo= MINo= MI2;ml= MINi= M14. 
The rectangular algorithm proceeds by a strategy similar to the square case. The only difference is that now each processing node splits its time among several square subblocks. 

10-4.2 Performance Analysis for Rectangular Subblocks 

The broadcast-multiply-roll cycle of dle algoridlm is repeated N 1 times. For each such cycle, dle time required for dle horizontal broadcast is: 
ni t tcomm + (N 0- 2) tstart (10-8) The time to multiply subblocks is: 

A3 
2A2A r2ml tcalc = mlmOtcalc 

and dte time to roll vertically is: 

mlmotcomm 

(10-9) 

(10-10) 

171\ 

Malrix Algorithms I: Malrix Multiplication Chap. 10 

Therefore, the total computation time for the rectangular subblock algorithm is: 
T = N 1 [2mrmotcalc + tcomm (mr+ mlmO) +(N o-2)tstart] (10-11) 

or 

+N1 [tcomm(ml + m1mo> + (No-2)ts/art~ N 
We would naturally like to determine the optimal choice of the aspect ratio, r , holding the total number of processors N o N I = N and the matrix size M fixed. The communication overhead from (10-12) is contained in the term: 
N I tcomm (nil + m1mo) + N I(N 0- 2) tstart (10-13) 
If 2 is negligible compared with N 0' the pipeline startup time only depends on N I N o=N and is independent of r. Under this assumption, our task is to minimize the remaining term: 

T= 

3 2 M tcalc 

00-12) 

tcomm N 1 (ml + mlmO) = tcomm M (mi + mo) (10-14) 
Since M is held fixed, the minimum is found at mi = mo, or r = 1. In other words, matrix multiplication with square subblocks is most efficient This conclusion does not extend at this point to the case of a pure-row storage algorithm which would utilize one- dimensional ring decomposition since in that case the timing formula Eq. (10-11) does not hold. For pure-row or column-storage schemes, there is no broadcast step in the algorithm, so that Eq. ( 10-11 ) takes the modified form: 
2M3 M2 
Trow.co!umll = N tcalc +Ntcomm N (10-15) 

10-4.3 Comparison of Square Subblocks with Row/Column 

To complete the theoretical analysis, we will write down the concurrent efficiencies as defined in Chap. 5 for the square subblock and the pure-row algorithms. For the subblock algorithm, the efficiency obtained from Eq. (5-5) is: 
1 1 (10-16) 

1+- 
..J1i"" tca/c 
where, as usual, n denotes the grain size which is in this case the number of matrix ele- ments per node (n = M2/N). In (10-16), S is typically a small term proportional to tstarl and is due to the pipeline startup time. In the same notation, the efficiency for the pure- row scheme obtained from (10-15) is: 
1 1 

- Eblock = 

1 + f block 

trow = 

1 ~+5 

l+frow 

I 1+- ..[ii" 

tcomm .!...IN" 2 

tcalc 

(10-17) 

Sec. 10-4 Rectangular Subblock Decomposition for Matrix Multiplication 

177 

The overhead f row associated with the pure-row algorithm is larger by a factor m- /2, where N is the number of processing nodes. 

10-5 Timing Results for Matrix Multiplication on the Caltech Hypercube 

10-5.1 Load-Balanced Decompositions 

To verify the validity of the preceding theoretical analysis of matrix multiplication algorithms, we examine the performance of matrix multiplication via the subblock and row techniques when run on the Caltech/JPL Mark II hypercubes. To serve as an efficiency baseline, a third program was run on a single-node processor multiplying matrices via the usual sequential algorithm. This baseline measurement yields the time on a single processor T l(M) needed in Eq. (5-4). Similar programming optimization steps were taken in the design of the three test programs. For example, indexing over- heads and related activity were reduced to a minimum in each case so as to yield a fair comparison of the algorithms. 
Timings for the sequential algorithm running on a single-node processor are shown in Fig. 10-4. The times are fit very well by the functional dependence: 
T l(M) = Tsillgle IIOde = 119~ M2+ 123~ec M3 (10-18) 
for M x M matrix multiplication. This result was obtained using single precision, 32-bit arithmetic throughout the algorithm. The quadratic term (which is manifested in Fig. 10-4 as the very small but nonzero y -intercept) reflects the times taken to initialize the matrix C to zero at the beginning of the multiplication and to do some simple index arith- metic. As much index arithmetic as possible was removed from the innermost loops of the program. 
Figure 10-5 shows timing results for a hypercube mapped to a 4 x 4 mesh and run- ning the square subblock algorithm. Plotted in Fig. 10-5 is TIM2 versus M, so that the asymptotic floating-point speed (the cubic term of Eq. (10-7» plots as a linearly increas- ing term; communication overhead appears as a constant (nonzero y -intercept). The pipeline startup time of the algorithm appears as a I1M2 term at small M. Examination of the plot shows that this implementation of concurrent matrix multiplication is quite efficient. Even for modest matrix sizes such as 16 x 16 or 32 x 32 (total matrix size) on a 4 x 4 machine, the pipeline startup time and communication overhead are small effects. A linear fit to the points representing timings for the larger matrix sizes is also shown in Fig. 10-5. The slope derived from this fit is 7.64 1.1 sec, almost exactly one sixteenth of 1231.1 sec, the coefficient of M3 in Eq. (10-18). 
An alternative presentation of the 4 x 4 timing data is given in Fig. 10-6 in which we plot lIe -1 versus ll{ii .As in the discussion of Chap. 5, this is done so as to reveal the overhead terms of Eq. (10-16) .The linear behavior for smallll{ii (large matrices) shows that the functional form ofEq. (10-16) holds. A linear fit yields the coefficient of ll{ii and the result: 

178 

Time/M2 (msec) 

Matrix Algorithms I: Matrix Multiplication Chap. 10 

20 

16 

12 

4 

0 

ii 

20 

40 

60 80 Matrix size, M 

100 

120 

140 

l'lgure 10-4 Matrix multiplication timings for a single node processor running the stan- dard, sequential algorithm. The times normalized by M2 are plotted for clarity. Also shown is the fit given in Eq. (10-18). 

tca/c 
This is a reasonable figure for the crystalline operating system on the Mark II hypercube machine discussed in Chaps. 4 and 14. One must bear in mind. however. that some of what we have lumped into the category of "communication overhead" (quadratic. M2. terms in the timing equation (10- 7» are actually attributable to other costs such as index initialization. which are dependent upon the details of user-written code. The lesson to be learned from this is that if one tries to write fast code by removing as much index arithmetic as possible from the innermost loop of the program. the sequential and con- current implementation are both speeded up. but the coefficient of 1/~ in the hypercube overhead is increased. Thus the strategy lowers the coefficient of the cubic term (giving a faster code) but increases the quadratic terms which have dependencies on n and N identical to that of "real" communication overhead. This will be further illustrated below. 

tcomm 

= 1.39 

(10-19) 

To show that the result of Eq. (10-16) is universal, i.e., independent of total machine size and dependent only upon the grain size, n, we examine different ensemble sizes. Fig. 10-7 is similar to the plot of Fig. 10-6, but with 8 x 8 machine results added. For n > 100 (a 10 x 10 submatrix in each node), the behavior is almost exactly universal. The deviation from the universality at small n is due to the larger pipeline startup time on the larger machine. 

Sec. 10-5 Timing Results for Matrix Multiplication on the Caltech Hypercube 

179 

1.25 

1.00 
Time/M2 (msec) 
0.75 

0.50 

0.25 

0.00 

Figure 10.5 Timings for the square subb1ock algorithm on a 4 x 4 machine. Times are normalized by M2, as in Fig. 10-4. 

.l-1 £ 

1/,fn 

l'lgure 10.6 The same timings as in Fig. (10-5), but plotted as overhead ( 11£-1) versus 11{;i" (11 is the grain size or the number of matrix elements per node). 

180 

~ 

Matrix Algorithms I: Matrix Multiplication Chap. 10 

1.0 

0.8 

0.6 

0.4 

0.2 

0.0 

0.0 

0.1 

0.2 

0.3 

0.4 

0.5 

I/In 

F1gure 10.7 The same as Fig. (10-6), but with timings fr<Xn a 8 x 8 machine added, showing universal behavior. 

We next turn our attention to the purely row-based algorithm which requires a one-dimensional (ring) decomposition. This program was timed using the same hyper- cube machine configurations as above. As expected, the communication overhead for this algorithm is greater than for the subblock case. The overheads are shown in Fig. 10-8 along with the previous results of Fig. 10-7 for comparison. A breakdown in the univer- sal scaling behavior is quite evident for this one-dimensional decomposition. 
From these timings, the coefficients of 1/..fn are 3.61 and 7.49 for the 16- and 64- node cases, respectively. The ratio of these slopes, 2.07, is in excellent agreement with the theoretical value of 2, predicted by Eq. (10-17). The mtio of 3.61 to 1.39 (Eq. (10- 19)), however, is 2.6, significantly different from the factor of 2 expected from Eqs. (10-16) and (10-17). The reason for this disagreement is the previously discussed effect of indexing costs masquemding as communication overhead. Due to the large aspect ratio, r, of the row decomposition, much more index initialization work is required than in the square subblock case. This hidden cost causes the row decomposition to have even more overhead than predicted by Eq. (10-17). 
The added cost associated with such index operations is a common occurrence. Frequently, a pamllel algorithm executes the same total number of floating-point opera- tions as the equivalent sequential algorithm. This is the case for all of the algorithms dis- cussed in this chapter. Communication overhead and load imbalance give rise to overt performance degradations, which reduce the concurrent efficiency. This other type of degradation, however, is rooted in the integer arithmetic associated with the indexing of 

Sec. 10-5 Timing Results for Matrix Multiplication on the Caltech Hypercube 

181 

.1 E 

F1gure 10.8 The same as Fig. (10-7), but with timings iran the one-dimensional column algorithm addOO. The large canmunicatioo overltead and non-universal behavior ob- tainOO with this decomposition are significanL 
the various arrays involved. A typical matrix algorithm incorporates inner loops con- structed as follows: 
for (i running over ma rows){ 
Do some index arithmetic dependent upon decomposition to find associated columnsfor the processor. 
for (j running over mi columns) { 
Inner loop identical to sequential code so that asymptotic efficiency is 1 }} 
The parallel algorithm requires more indexing operations than the sequential algo- rithm, leading to a new type of overhead to which we will refer as software overhead . As indicated above, the software and communication overheads have similar dependen- cies on n and N. Though typically a small effect, it. should be kept in mind in any attempt to rigorously account for concurrent efficiencies. We will comment on this again in Chap. 20 for another matrix algorithm. 
These experimental timing results confirm the theoretical assertion that the row decomposition scheme is less efficient than the subblock decomposition and that the size of the overhead increases with machine size as {f[ .in accordance with Eq. (10-17). 

182 

Matrix Algorithms I: Matrix Multiplication Chap. 10 

10-5.2 Padding for Arbitrary Matrix Size-Load Imbalance 

Any practical matrix multiplication algorithm must be able to multiply matrices of arbitrary size-not just those whose order allows an exact mapping onto the concurrent processor with each processor having equal numbers of matrix elements. The simplest remedy for this problem is to pad the multiplicand matrices with a unit matrix block, in order to fill out the mapping. 
Figures 10-9 and 10-10 present the results of the padding algorithm, examining a continuous range of matrix sizes M. The square subblock and column cases were run on a Mark II 64-node hypercube. Both algorithms display sawtooth perfonnance curves whose upper envelopes are described by Eqs. (10-16) and (10-17). It is noteworthy that the severity of the sawtooth degradation is much larger for the one-dimensional decom- position. The reason for this difference is readily understood. If we set out to multiply 65 x 65 matrices using a row decomposition on a 64-node ring, it requires padding by 63 rows up to a matrix rank of 128. If the same multiplication is done blockwise on an 8 x 8 machine, however, the amount of padding necessary is only 7 rows, resulting in a 72 x 72 multiplication. 

1.0 

.8 

.6 

.4 

.2 

+ + +++ 
+++++ +++ + + + 
++ +++ + +++ + ++ + +++ + ++ +++ + + + ++ ++ + ++ ++ +++ + + + + + + 
+ + + + Efficiency On A 
++ +++ 64NodeHypercube 
With A Square Submatrix Decomposition 

+ + 
+ 
+ ++ 

+ + ,. + + T + + .t .+ T 

J t 

++ 

+, + 

0 50 100 150 200 250 300 Matrix size, M 
FIgure 10.9 Efficiency versus matrix size for a continuous range of matrix sius, illus- trating "sawtooth" inefficiency effecL The square subblock algorithm on a 8 x 8 machine is illustrated. 

0 

350 

The sawtooth effect may in some cases be a much more important practical con- sideration governing the choice between the two algorithms than Eqs. (10-16) and (10- 17), although the amplitude of the effect diminishes for large matrices. 

Sec. 10-5 Timing Results for Matrix Multiplication on the Caltech Hypercube 

183 

E 

1.0 

.8 

.6 

.4 

.2 

.0 

0 50 100 150 200 250 300 350 Matrix size, M 
FIgure 10-10 Same as Fig, (10-9), but for the ooe-dimensional colwnn algori!Jlln on a 64 node machine, 
10-6 Summary of Matrix Multiplication 
In the previous sections, we have demonsb"ated a simple but efficient method for implementing matrix multiplication based upon a generalized block decomposition scheme. Both theoretical and experimental performance analyses suggest that low-aspect ratio block decompositions offer the most favorable concurrent performance characteris- tics, 

It is significant to note that purely row (or column) based algorithms may be more easily adapted from sequential programs, due to the simpler communication topology. Whether the additional costs incurred by the selection of such a nonoptimal method are justified may be application-dependent or may even depend upon the manner in which matrices are stored and handled by the programming language. For example, it may be expedient to make use of FOR1RAN's natural treatment of column-based algorithms into which substantial optimization effort has been placed over many years. It is unfor- tunate that such algorithms are nonoptimal for concurrent decompositions. 
In the performance analysis of Sec. 10-3 and 10-4, we have ignored any latency in the communication system, i.e., we have assumed that tcomm is independent of message size. As seen in Fig. 4-2, this is a good approximation on the Caltech hypercubes but not on the INTEL iPSC. Note that row/column storage involves the transmission of fewer messages (by a factor of 2) but in total more bytes (by a factor of -.[f{ /2) than the square subblock decomposition. Thus, on machines with substantial latency it can, in fact, be best to use row or column storage [Geist 86, Moler 86]. 

184 

Matrix Algoritluns I: Matrix Multiplication Chap. 10 

10.7 Comments on the Decomposition of Spreadsheets 
The calculation of spreadsheets is a major use of personal computers for both financial and database applications. Although spreadsheet calculation is not generally regarded as the epitome of large-scale computation, it is of interest to consider the appli- cability of concurrent processing to such problems. 
In Fig. 10-11, we show a simple spreadsheet; the typical structure is that of a rec- tangular array of cells. Each cell can be assumed either to be fixed (i.e., entered as data) or to be a function of the contents of other cells. In general, such a functional depen- dence results in a highly irregular pattern of connections between cells. The proper implementation of such a problem on a concurrent processor requires sophisticated message-routing and load-balancing techniques which will be discussed further in Chaps. 21, 22, and 17, respectively. 

A SAMPLE SPREADSHEET DECOMPOSED 
ON A 16 NODE CONCURRENT PROCESSOR abcdefgh k mnop 
1 
~~ICI~~'~~m~ I 'I 
2a~OZ ~ 'b. 

5 

6 

h~ 

---GLOBAL CONSTANTS 

TYPICAL PROCESSOR BOUWARIES 

A ~l 

UNUSED 

GRAND TDTAL 

l'lgure 10-11 A sample 16 x 16 spreadsheet decomposed on a generic 16 node con- current processor with 16 cells in each node. 

Sec. 10- 7 Comments on the Decomposition of Spreadsheets 

185 

Here we will be satisfied to observe that as in the matrix case, the spreadsheet encompasses a domain which is two-dimensional; however, it does not necessarily incor- porate nearest-neighbor interactions. In general, we might expect the same optimal decomposition schemes as for the case of matrices. In this situation, the scattered square decompositions discussed in Chap. 20 are likely to be best. For simplicity, however, Fig. 10-11 only shows the simpler contiguous block decomposition which we have used in this chapter for matrix multiplication. Our example considers the basic cell operation to be: 

M..=c..+ad..M 
l.+Re..M.l(10- 
20)'I 'I 'I I p 'I , 
where M lj and Mil store overall properties of the given columns and rows, respectively. This calculation may be performed concurrently after global broadcast of the constants a and ~. The row and column paramete~s Mlj and Mil need to be piped along columns and rows just as we saw for the matrix A in the multiplication algorithm of Sec. 10-3.2. 
Often in spreadsheet calculation one needs to form partial or full sums of all cells in a row or column. If, for example, the rows denote people, the columns years, and the cells annual expenditures, then the subtotals represent the total expenditure in a given year or an individual's expenditure summed over all years. Operations of this kind are straightforward in their adaption and perform well on a hypercube concurrent processor; the best such special algorithm for this task will be described in Chap. 19 and Sec. 21- 3.4. 

The above discussion leads us to expect the typically efficient concurrent imple- mentation, with overhead proportional to 1/-rn , where n cells reside in each node. The only real difficulty with practical concurrent spreadsheet calculation is in the load imbal- ance caused by inhomogeneity (differences) of the cell definitions. Such load-balancing problems will be treated in Chap. 17. 
It would appear that, in principle, concurrent processors are well suited to spreadsheet applications, and recognition of this point could even motivate the generation of real code. With this possibility in mind, it may not be unreasonable to expect a future generation of personal computers or workstations built around a concurrent architecture. 

11 

The Fast Fourier Transform 

11-1 Introduction 
In 1965 Cooley and Tukey [Cooley 65] described an algorithm for computing the discrete Fourier transform (FT) of a series in a number of operations of order M logM , where M is the number of points in the series. This algorithm provided a tremendous improvement over the order M2 algorithms previously used. This algorithm and its vari- ants are known as Fast Fourier Transforms (FFTs). The FFT is one of the principal algo- rithmic tools used for signal processing, image processing, and solving differential equa- tions. Fourier transforms can be used to compute convolutions, analyze spectra, and smooth and identify features in images. In many cases, the reduction in computation from order M2 to order M logM has changed the Fourier transform from a theoretical abstraction to a powerful numerical tool. 
In this chapter we will consider finite, complex Fourier transforms. This is not an unduly restricted class, since the bulk of numerical FFTs are of this type. Infinite transforms are extremely important, but they must invariably be approximated by finite transforms before being treated numerically. We will also restrict ourselves to the binary FFT , in which the number of points is a power of two. Other algorithms, which actually perform better than the ones discussed here exist to treat the more general case. Such algorithms have been shown to work extremely well on parallel processors [Aloisio 87], and are discussed briefly in Sec. 22-3.2, in the context of the crystal Jouter . 

11-2 The FFf Algorithm 
The FFT can be cast as either a recursive or an iterative algorithm. Most scientific applications use the iterative algorithm because of its slightly enhanced speed. not to mention the practical fact that FOR1RAN does not support recursion. The recursive algorithm is easier to understand and program. but is usually slower because of the over- head for the additional subroutine calls. Since simplicity is more important in this illus- trative discussion than speed. we will discuss a parallel implementation of the recursive FFT algorithm in this chapter. The issues involved are essentially the same in both cases. and a clear understanding of this chapter should allow one to adapt an iterative version of the algorithm to a parallel machine. The Software Supplement contains implementations of both the iterative and recursive algorithms. 
We will start by deriving the recursive FFT algorithm. Consider a sequence of M complex values. (Xo. Xi. XM-J. and define their Fourier transform as 
M-l .k Xk = LXj exp(21tiL) (11-1) j=O M 

187 

188 

The Fast Fourier Transfonn Chap. 11 

For reasons that will become clear, M is taken to be a power of twO. It is possible to compute this sum precisely as it is written; compute x k for each k from zero through M -1 as a sum of M terms, in a total time proportional to M2, However, if we split the summation into two parts, with one containing the even elements of x and the other con- taining the odd, a dramatic reduction in time is realized. To effect this improvement, we begin by writing: 

-.M.-1 
-2 .2.'k 2 .(2.'+1)k Xk= LX2j'exp(21tl~)+ LX2j'+leXp(21tl~) 
j'=O M j'=O M 
We now make the following definitions: M'=M 

.M.-I 

2 

ei = X2i 

and 

(the even elements of x ) 

Next, restrict consideration to k < T, and write Eq. (11-6) as: 

(11-5) 

which can be expressed in tenns of e and if, the transfonns of e and 0. -- 
(2.k)- Xk == ek + exp 1tl- °k M 
x M == e k -exp(21ti~ )if k (11-8) k+- M 
2 
Thus, to compute an M point Ff , we compute two M /2 point Ffs, e and if, and combine them. The combination requires M /2 complex multiplications and M complex additions. If we let T (M ) be the time required to perfonn an M point Ff, then: 

T(M)=2T(M/2)+M(t:t+~) (11-9) 2 
where t* is the time required for a complex multiplication, and t:t is the time required for 

Sec. 11-3 Sequentiallmplementation 

189 

a complex addition. The M /2 point FTs can be computed recursively, terminating when M is equal to 1. Since the one-point FT requires no work at all, T(I) = 0. The recursion relation, Eq. (11-9), can then be solved for T(M), giving 

t T(M)= (t~+f)MlogzM 

(11-10) 

11-3 Sequential Implementation 
The program fft-l in Code II-I implements the recursive algorithm just described, on a sequential computer. To avoid unnecessary data movement, the routine acts on "arrays" that are not stored in contiguous memory locations. The elements of the array x are understood to be stored at x[O], x[P], x[2p], ..., x[(M -I)p], where p is initially I at the zeroth level of recursion. Thus, the Ff of the even elements of x is computed by dou- bling the increment, p , and halving M, while that of the odd elements is obtained by starting at x[p], and doing the same. Upon completion x holds the Ff of the input array. The original values are lost. The sequence of procedure calls for an 8-point recursive FFf is shown, diagramatically in Fig.II-I. 

Code 11.1 The Sequential Recursive FFI' 
As shown in Code 11-2,fft-combine combines the twoM 12 point FTs, e and {f , produc- ing the M point FT. It assumes that the input vector, x, contains the elements of e and elements of {f at alternate locations. However, the order of the elements of e and {f is permuted. Schematically, we can write the vector x as: 
x= 
(e.{f.e.{f e. 0. 
)(11-11 
)I. 1. I, I, IMI2-1 IMI2-1 
where (i o i 1 ...iM/2 -J is a permutation of (0 1 ...M 12 -1 ). In order to properly evaluate Eq. (11-8), we require the function, permute which maps 0 into io, I into it, etc., and which is discussed further in Sec. 11- 7. 

190 

The Fast Fourier Transfonn Chap. II 

-cFFf(I.8..+1) ~ 
{FFf(2.4..+1) ~COmbine(2.4..+1) I FFf(I.8..+5) 4 FFf(4.2..+1) 
-CFFf(I.8..+3) ~ 
FFf(2.4..+3) ~ Combine(2.4..+3) - FFf(I.8..+7) 
E ~E ~E ~E ~E ~E ~E ~ RO Rl R2 R3 R2 Rl RO 

r FFT(4,2,x) -l 

FFT(8.1.x)-1 

-cFFf(I,8,x) L, ~ FFf(2,4,x) ~ Combine(2,4,x) 
I FFf(I,8,x+4) 

I 
-CFFf(I,8,x+2) ~ 4 FFf(2,4,x+2) -I' Combine(2,4,x+2) 
FFf(I,8,x+6) 

-7 Combine(4,2,x) 

--7 Combine(4,2,x+l) - 

~ Combine(8,I,x) 

Figure II-I Schematic representation of an 8-point recursive FFr. Rn refers to the nth level of recursion. The notation FFr(N,p,x+m) indicates that the N data points x[m], x[p+m], x[2p+m],..., x[(N-l)p+m] are to be transformed. Thus m is the offset from the start of the data array, and p is the increment between the data points. The notation for the canbine routine is similar. 

Code 11-2 Routine to Combine the Even and Odd Partial FFrs in the Sequential Algorithm 

11-4 Parallel Decomposition 
Now we turn to the parallel decomposition of the FFr algorithm. The natural space over which to decompose is that of the elements of x. In this approach, we divide the M elements of x equally among N processors. While not strictly necessary , 

Sec. 11-4 Parallel Decomposition 

191 

simplicity dictates that N, the number of processors, be a power of 2 so that each proces- sor holds an equal number of elements. The simplest such decomposition involves storing the first M IN elements in processor 0, the next M IN elements in processor 1, and so on. As the algorithm proceeds, the computed FT involves progressively fewer points. Even- tually, there remains only one point in each processor. At even deeper levels of recursion, the number of points in each transform is less than the number of processors. At this point, several transforms are computed independently, and simultaneously, in disjoint subsets of the processors. In such cases, if the first element of any given sub-problem is stored in processor p, then the next corresponding element is found in processor p +NIM, and so on. 
With this data decomposition, the parallel algorithm is straightforward. It must take account of three distinct cases: If M is 1, stop the recursion. If M is less than or equal to the number of processors, then only one point of the array to be transformed is stored locally. The processor assists in computing either e or 5 , according to whether the single element is even or odd. Then, interproces- sor communication brings e and 5 together in the same processor, and Eq. (11-8) is evaluated. 
If M is greater than the number of processors, the computation is the same as in the sequential case, but extends only over the values stored locally. 
Code 11-3 displaysfft-2, a parallel version offft-l. The subroutinefft-combine-2 in Code 1 1-4 is almost identical tofft-combine, which was used in the sequential imple- mentation. The only difference is that the inner loop spans only those elements of the array stored on a single processor. 
Interprocessor communication enters only in the subroutine fft-combine-comm. shown in Code 11-5. This routine performs the combination operation for an M point FFT when the number of processors is equal to or less than the number of points in the transform. It must decide whether it holds an even or odd element and which processor holds the complementary element. Then it must exchange a complex value with that pro- cessor, carry out the complex arithmetic, and store the result. 
With the addition of routines for complex arithmetic and index permutation, Codes 11-3, 11-4, and 11-5 form a complete module for computing FFTs in parallel. Neither the sequential routines, fft-l and fft-combine , nor the parallel ones have been optimized. Precisely the same optimizations can be applied to both. These include avoiding the subroutine overhead associated with complex arithmetic by in-line expan- sion, avoiding repeated calls to permute by using a look-up table for the complex roots of unity which are keyed to the global index, avoiding the subroutine overhead of the recur- sive algorithm by switching to the iterative algorithm, and making use of the regular pro- perties of the 2nd , 4'11 , and 8'11 roots of unity to avoid redundant floating-point operations using a radix-8 algorithm. These optimizations neither restrict nor hinder the paralleliza- tion or the decomposition developed so far. 

192 

The Fast Fourier Transfonn Chap. II 

Code 11.3 The Concurrent Recursive FFf 

Code 11.4 Routine to Combine the Even and Odd Partial FFfs 

Sec. 11-5 Communication Structure 

193 

Code 11.5 Routine to Exchange and Combine Elements of x in the Concurrent Algorithm 

11.5 Communication Structure 
The parallel routine exchanges data between processors. In order to implement this operation in a real program, the virtual machine function vm-shift would be used. To determine how well the algorithm will work on a given real machine, we must know whether the virtual connections used by vm-shift are present in hardware. If they are present, we can expect the algorithm to perform very well. Thus, it is important to know exactly which processors are connected in the algorithm and what hardware topology is necessary to provide these links directly. 
When fft-combine-comm is first called, the number of elements involved in the transform is equal to the number of processors. Thus, there is an element in each proces- sor, and even and odd elements are stored in processors with adjacent processor numbers. At the next level of recursion, a point is stored in every second processor, and even and odd elements are stored in processors whose processor numbers differ by two. Proceed- ing deeper, we find that even and odd elements are stored in processors whose processor numbers differ by four, eight, sixteen, etc. In other words, a processor requires commun- ication with those processors whose processor numbers differ by exactly one, two, four, eight, etc. This connectivity is precisely that of the hypercube. Thus, in order to avoid message forwarding in the FFT algorithm, the hardware communication graph should contain a hypercube as a subset The use of the hypercube for FFTs was explored quite early by Allakhverdiyev and Sarafaliyeva [Allakhverdiyev 81]. 

194 

The Fast Fourier Transfonn Chap. II 

11-6 Performance on the Hypercube 
To simplify the perfonnance analysis, we will assume that the machine intercon- nection topology contains a hypercube. Under this assumption, an exchange requires a fixed amount of time in each of two processors, and no other processors are affected, as they would be if messages required forwarding. Following the discussion in Chap. 5, we define the communication overhead, f , according to: 
NT 
(M 
)-T(M 
)f = N 1 /1 1 1 ", 

T l(M) ,~~ ~-I 
where TN(M) is the time to complete an M -point FFf on N processors. The parallel algorithm splits neatly into two stages, M > N and M ~ N. The following recursion rela- tion holds for M > N : 

M M 
TN(M)=2TN(2)+ W't> (11-13) 
where ~ is the time spent in one iteration of the loop in fft-combine-2 .Note that this relation is similar to Eq. (11-9). The solution is slightly different because the recursion terminates when the number of points equals the number of processors, i.e. when N = M . 

ForM :$;N we have 

M 
TN(M)= TN(2) +'t< 
where 't< is the time spent in a single ca11 to jft-combine-comm. Eq. (11-15) is 

(11-15) 

The solution of 

TN(M) = '[< log2M 
In particular, TN(N) = '[<logzN , which can be substituted into Eq. (11-14) giving M 
TN(M) = -('[> log2 M + (2'[< -'[» log2N) 2N 

Thus, the overhead is given by: 

(11-16) 

(11-17) 

(11-18) 't> log2 M 
The effects due to details of the implementation and the particular machine's archi- tecture are contained in 't> and 't< .On most parallel computers, floating-point arithmetic is much more costly than integer arithmetic, so we make the following generalizations about't> and 't< : 

f= 

2't< -~ log2 N 

't> = t. + 2t:t ::: 10 tcalc 

(11-19a) 

Sec. 11-6 Perfonnance on the Hypercube 

195 

't< = t. + t~ + tshift :::: 8tcalc + 2tcomm (11-19b) 
where t. and t:l: are the typical times for complex multiplication and addition operations, respectively, and tshift is the time required to shift a complex value, equal to twice the value of tcomm as defined in Sec. 3-5. Note that when N=l, Eq. 11-19 and Eq.11-17 are equivalent to Eq.l1-10. With the approximations of Eq.11-19, the overhead is 
2tcomm + 3tca/c log2N 
f = {11-2ffi 

Stca/c log2M ,-- --, 
The overhead decreases slowly to O if M is increased with N fixed. Figure 11-2 is a plot of f vs. logzN flog# , as measured on the Caltech Mark II hypercube. The points form a straight line corresponding to a value of tcommftca/c = 3.3 and show excellent agreement with the model. 

2.0 

.6 

1.2 

Overh~ad 

0.8 

0.4 

Fast Fourier Transform 

xV 

0 

00' 

v ,6 

O N~2 D N~4 " N=8 v N=16 x N =32 

<90 

I, , , ., 0.0 ~ ~ 

0.0 0.2 0.4 0.6 0.8 1.0 logN I logNx 
Figure 11.2 The overhead! = 11£-1 defined in Eq. (11-12) and discussed in Sec. 11-6. The points come fr<Xn measurements on the Caltech/JPL Malk n hypercube [Noerdlinger 86]. 

Even with tcomm ::: O the overhead does not vanish. This is a because both the pro- cessor and its complement compute the complex product in Code 11-5. Thus, we have introduced a small amount of algorithmic overhead in addition to the communication overhead. This overhead cannot be removed simply by computing the product in the odd processor, and transmitting only the result. In that case, although the even processor does no work, it is idle while the odd processor does the computation. A solution to this prob- lem does exist for the iterative algorithm (which appears in the Software Supplement). The difficulty in the recursive algorithm is that there is really nothing else for the even processor to do while the odd processor computes the product. In the iterative algorithm, 

196 

The Fast Fourier Transfonn Chap. II 

the two processors can agree each to compute half of the products at any given stage of the computation. 

11- 7 Permutation of the Indices by the FFT Algorithm 
In Sec. 11-1, we postulated a function permute which mapped the offset from the beginning of the partially transformed array, x. into the index of the element stored at that location. In this section, we derive the details of permute. 
We begin by giving a recursive description of permute U, M), where M is the number of items in the permutation. The simplest case occurs when M =1. Since there is only one item in the permutation, there is only one possible result: 
permuteU,I)=O (11-21) 
We can write permute U , 2M ) in terms of permute U, M) because when combining the even and odd M -point transformed arrays, we store the results in the 2M -point array as follows: 

(.ii'.ii+M)~(ei,ifi) (11-22) 
Let the pair of elements ( e i' if i) be at offset j in the M -element array. Then, by the .. 
definition of permute , 

io = permute (j, M) (11-23) 
After combining e i and if i, the values x i and x i +M are at offsets 2j and 2j + 1, respec- o o o o 
tively, in the 2M -element array. Thus, 
permute (2j, 2M) = io = permute (j, M) (11-24) permute(2j+1,2M)=io+M =permute(j,M)+M (11-25) Eq. (11-24) and (11-25) can be recast as a recursion formula for computing 
permute (j, M). 

permute U, M) = 

if j is even 

if 

is odd 

(11-26) 

Finally, we can translate the operations of Eq. (11-26) into words: 
1) Set the high-order bit (the one corresponding to the value of M /2) in the result to be the same as the low-order bitinj. 
2) Compute the rest of the bits by shifting both j and M one bit to the right and repeating the procedure. 
This procedure selects bits from j, starting at the right, and inserts them into the result, starting at the left. It simply reverses the bits ofj. 

Sec. 11-8 Inverting the Permutation in Parallel 

197 

11-8 Inverting the Permutation in Parallel 
One consequence of the permutation of indices is that elements of the final result are also permuted. Thus, ..f k is not found at an offset of k from the beginning of the array, x, after the completion of the FFf program. Instead, ..f k is found at an offset given by permute (k, M ), which is equal to the bit-reversal of the binary representation of k, as described in the previous section. The fact that ..f k is not where one might expect is usu- ally not a problem because it is simple enough to determine its true location. On the other hand, sometimes it is useful to invert the permutation and rearrange the elements of ..f k in order of increasing k. A sequential program which does this is trivial because per- mute is its own inverse. For each value of k , one simply exchanges the two elements at offsets k and permute (k , M ). 
The task is somewhat more complicated in the parallel case. Recall that the array, x, is distributed over all the nodes of the parallel processor. In Sec. 11-5, we found that the necessary communication topology for computation of the FFf is that of a hypercube. Thus, we restrict ourselves to a hypercube in the following discussion. The goal is to take the array, ..f , after the FFf has been performed, and to rearrange it so that ..f k is found at offset k mod n, in processor kIn, where n is the grain size, n == M IN. We divide the task into two parts: 1) Get the correct elements into the correct processors, regardless of order; and 2) Rearrange them locally without need for further communica- tion. 

To place elements in the correct processors, we need a way to identify the home processor of a given element. Recall that we want the number of the processor in which Xk is stored to correspond to the high bits of k , and that after the completion of the FFT , the high bits of k can be found, reversed, in the low bits of the offset of Xk. Thus, ele- ments destined for processor p are those with offsets from the beginning of their home processor whose low bits correspond to the bit-reversal of p .The bit-reversal is denoted by ii .These elements are evenly spaced in memory, separated by N elements. They are 
always at offsets ii, ii+N,ii+2N, Slightly trickier is determining the route that a given set of elements will follow to reach their destination. Remember that we are restricted to using the hypercube com- munication channels. As can be readily verified, destinations are not necessarily nearest neighbors in the hypercube topology. We must find a way to route the messages through several nodes, and we must be sure that the intermediate nodes are prepared to forward the data properly. This problem can be solved by insisting that all elements being transferred at any given time are taking the same route, where a route is defined as a sequence of traversed hypercube channels. A route can be uniquely identified by an integer, r, in the range O ~ r ~ (N -I ). Each bit in the binary representation of the integer represents a channel that is part of the route. Since all processors make use of the same route simultaneously, and hence the same channels, the polled communication scheme is perfectly adequate. Furthermore, since all processors are simultaneously working on the same route, there is no load imbalance or idle time. At any given time all the processors are active forwarding data to different final destinations. Organizing a long-range communication problem by enumeration of routes is a powerful idea that will 

198 

The Fast Fourier Transfonn Chap. 11 

be considered further in Chap. 21, where it is called the method of ordered cube geo- desics. The index algorithm of Chap. 21-3.2 is very similar to the algorithm of Code 11-6. In fact, Code 11-6 could be replaced by a local permutation followed by a call to index . 

Code 11.6 The Routine gloOOI-re\lerse 

With the communication phase complete. all that remains is to invert the permuta- tion locally. While straightforward in principle. the local permutation requires consider- able bit-oriented arithmetic to reverse bits. identify low-order and high-order bits. etc. Both C and FORTRAN implementations of this set of operations appear in the Software Supplement 
The cost of inverting the permutation is quite high. Almost every element of the array. x. must be transmitted over a long distance through the hypercube. The average distance traveled is half the diameter of the hypercube. or (log~)/2. We can estimate the time taken by Code 11-6 by counting the number of data transmissions. We find that each node executes (NlogzN)/2 exchanges. each one of M /N2 complex numbers. so that the total time taken is given approximately by: 
MlogzN 
Tinvert = (11-27) 

tcomm 

N 

The overhead due to inversion, f invert , which is a fonD of algorithmic overhead, is given by 

(11-28) T l(M) 10tcalc logzM 
Equation 11-28 has the same dependence on N and M as Eq. 11-20. Thus, the overhead due to permutation of indices in a parallel FFf has the same effect as slightly degraded communication speed. It does not affect the performance scaling with N or M except to modify the effective communication speed. 

finllert = 

Till1/ert (M) 

~ logzN 

= 

Sec. 11-9 Multidimensional Transfonns 

1QQ 

11-9 Multidimensional Transforms 
Multidimensional Fourier transforms are extremely useful in a wide variety of dif - ferential equation and image-processing problems. In a multidimensional transform (we will consider two dimensions for illustrative purposes), the input vector and the output transform are objects carrying two indices. The two-dimensional transform, with M I points in direction 1 and M 2 points in direction 2 is defined by: 
MI-l M2-1 
ijk = r. r. x/m exp(21ti(~+~» (11-29) 1=0 m=0 I 2 
We can factor the exponential and use the standard rules for removing a constant from a summation to obtain 

MI-l. M2-1 
ijk = r. eXp(21ti(1-» r. Xlm exp(21ti(-t;-» (11-30) 1=0 I m=O 2 
The inner summation in Eq.(11-30) requires the computation of MI one- dimensional Ffs, with one for each value of 1, each involving M 2 points. Similarly, the outer summation requires the computation of M 2 one-dimensional Ffs, one for each value of k, each over M I points. The total time predicted for these operations has the dependence: 

T oc M 1 (M 2log2M ~ + M 2(M llogzM J = M logzM 

In this expression, 

(11-31) 

M =MIM2 (11-32) is the total number of points in the two-dimensional mesh. The total time required for the multidimensional transform ofEq. (11-29) is exactly the same as for the one-dimensional case. The algorithm simply applies the one-dimensional algorithm many times in sequence, first to evaluate the inner summation and then to evaluate the outer summation. 
To compute multidimensional transforms in parallel, we still need a hypercube communication structure. We may consider a hypercube of dimension d to be a product, in a topological sense, of smaller hypercubes. The sum of the dimensions of the smaller cubes, d 1 and d2, in a two dimensional problem, must be equal to the dimension of the whole hypercube, d. When computing the inner summation, we use the d 1 communica- tion channels associated with the first subcube. Then, to compute the outer summation, we use the d2 communication channels associated with the second subcube. A detailed analysis shows that the results of Sec. 11-6 are unchanged for multidimensional Fourier transforms. In the formulas of Sec. 11-6, M should be considered as the total number of elements in the transform, regardless of the number of dimensions. 
An alternative decomposition for multi-dimensional transforms is the row decom- position. Again, for simplicity, consider a two-dimensional transform. We begin with one or more entire rows of the array stored within a single processor. The transforms of the rows can be computed using any sequential algorithm available. There is no 

200 

The Fast Fourier Transfonn Chap. 11 

communication within the transfoml phase. Then rows and columns of the partially transfomled array are transposed, so that complete columns are now stored within a sin- gle processor. Finally, the columns are transfomled, again without communication. This approach has some important advantages. 
Any sequential transfoml algorithm can be used. Non-binary transfomls, as well as algorithms tailored especially to known properties of the data like real-valuedness, may be used without difficulty. 
All communication is isolated into a single phase which can be optimized indepen- dently of the transfoml. Additionally, some advantage may be gained by transmit- ting large amounts of data at one time rather than in several small chunks. 
The algorithm to effect the transpose of the data is very similar to that used in Sec.II-8 to invert the bit-reversed transfoml. Each processor loops over routes through the hypercube, and sends data destined for the processor at the end of the route. As in Sec.II-8, the time needed for such communication is 
M logzN 
T (11-33) 

-IrunJpw" N "tcomm 
Since iliere is no oilier source of overhead, ilie total overhead in a row decomposed multi-dimensional FFr is given by 

- - 

(11-34) 10tcalc log2M 
This overhead is functionally the same as all others encountered in this chapter. How- ever, by performing all the communication at once, we have achieved a reduction in the coefficient of tco"."./tcalc that appears in the expression. This can be quite important on systems with a poor communication speed. Consolidating the communication can also be achieved for one-dimensional transforms, using the generalized collective communica- tion routine index described in Chap. 21-3.2. 

f= 

tcomm logzN -