18 

18-1 Introduction 

This chapter discusses sorting: the arrangement of data, which may consist of any of a variety of symbolic or numerical objects, into some set sequential order. With this generalized concept of order, we are not restricted to sorting only numerical values. In general, any data that can be keyed to values through some function are appropriate. This chapter is based upon work by Felten, Karlin, and Otto [Felten 85a]. The ideas presented here are appropriate for MIMD machines and are somewhat specific to hyper- cubes. Similar work appears in [Wagar 86]. Parallel sorting is discussed in other con- texts, in [Akl 85, Batcher 68, Baudet 78, Johnson 84, Li 86, Thompson 77]. 
Sorting is an algorithm of fundamental importance because it occupies a central part of many computational applications. To give just one example, we might consider a hypothetical symbolic manipulation program running on a concurrent machine. Having arrived at a complex set of algebraic expressions the program may be faced with the task of simplifying them to obtain a reduced result. This is done by expanding out all func- tions, combining expressions over common denominators, and so on. The process expands the original expression into many terms (perhaps hundreds of thousands), rapidly using up processors and memory throughout the parallel machine until, finally, only atomic expressions remain (for example, polynomials). The simplification is now a relatively straightforward process: like terms need only to be brought together, and their coefficients combined. This is where the parallel sorting algorithm enters. Sorting is also useful for load balancing parallel algorithms. In [Gottschalk 87] a parallel sort is employed to automatically decompose a missile tracking problem. Section 22-3 illus- trates a recursive bisection decomposer which is closely related to the parallel quicksort discussed in this chapter. 
It is not immediately clear how to define the canonical form for "sorted" data when dealing with a distributed-memory parallel computer. Should the resulting sorted list wind up in one of the processing nodes of the machine, or should the list remain distri- buted? The answer to this question depends somewhat upon the specific purpose for which the sort is intended. For the purposes of this chapter, we will be motivated by examples such as the above {load balancing, symbolic manipulation) and as a result, we will require that the sorted list remain distributed throughout the machine. In this way, it can easily form the input for further concurrent processing. Further pursuing the exam- ple of symbolic manipulation, consider the following: if the user, having simplified the expression down to some collection of relatively few terms, now desires the derivative of the expression, each processing node is able to compute the derivative of its part without preliminary manipulation. 

327 

328 

Sorting on the Hypercube Chap. 18 

In this discussion, we will suppose that the list to be sorted initially resides in a dis- tributed fashion in the parallel computer. The data items will be taken to be positive integers and the key will consist simply of their respective numeric values. We will require that at the end of the sorting, the data residing in each node are sorted internally and that these sublists are also sorted globally across the machine in some way. The natural way to specify this global order is to map the concurrent processor into a one- dimensional periodic auay or ring via the gridmap utility discussed in Chaps. 4 arid 14. The nodes along the ring are labeled by ring position, or ringpos. Global order is now defined by relative values of ringpos. Items in the node with ringpos=O precede items in the node with ringpos=l, which in turn precede the items in the node with ringpos=2, and so on. 
For some simple types of sorting strategies, such as bubblesort, a ring topology for the parallel computer would be sufficient. However, many faster sorting methods will turn out to require a much higher connectivity or dimension (actually, in the sense of Sec. 3-5, infinite system dimension) and we will assume that our parallel machine has a con- nectivity at least as rich as that of the hypercube. This is not to say that other types of connectivity (e.g., perfect shuffle) are not also attractive for sorting. Since we will see that the hypercube does a good job as compared with the best sequential algorithms, Shellsort and quicksort, no compelling argument for other connectivities is immediately obvious. 
The first algorithm we will present, Batcher's bitonic sort, does not leave the list in the desired one-dimensional ringpos order, but instead, the list ends up sorted according to the binary labeling scheme of the hypercube, (i.e., procnum ). Consecutive items may not even be near one another in the machine, but for many applications this presents no problem. Some applications of sorting simply require that one know where the items are, not that they be nearby. 
Three algorithms for sorting will be covered: In addition to Batcher's bitonic sorter [Batcher 68; Knuth 73a, pp. 232-33; Nassimi 79], we consider Shellsort and quicksort [Knuth 73a, pp. 379-81]. As discussed in Sec. 3-6, there are two ways of measuring the quality of a concurrent algorithm: "speed at any cost," or "speed per unit cost." It is interesting that the sorting algorithms to be discussed here illustrate these different cri- teria rather well. The bitonic sorter turns out to be the speed at any cost champion of the algorithms; if one fixes the list size and allows the hypercube to grow without bound (i.e., work on small lists) this algorithm is fastest On the other hand, bitonic is very inefficient For efficient algorithms, or equivalently, algorithms which run much faster on large lists, (holding the machine size fixed and allowing the list to grow) one needs to turn to Shellsort and quicksort. 
Before plunging into the first method, it should be pointed out that the subject of sorting is an extremely diverse one and we have chosen to address only a small fraction of the known methods. We have concentrated on the best general-purpose algorithms, so that no special properties of the sorted list are exploited. If the list to be sorted has spe- cial properties, such as a known distribution range (e.g., random numbers with a flat dis- tribution between O and 1) or high degeneracy (many redundant items, e.g., text files), 

Sec. 18-3 The Compare-Exchange Routine 

329 

then other strategies can be faster than Shellsort and quicksort. In the case of known data distribution, a bucketsort strategy (e.g., radix sort) is best, while the case of high degen- eracy is best handled by the distribution counting method [Knuth 73a, pp. 379-81]. We will comment on these methods at the end of the chapter and these form useful exten- sions for the interested reader to pursue. 

18-2 The Merge Strategy 
In the merging sttategy to be used in all our sorting algorithms, the first step is for each processor to sort its own sublist using some fast algorithm, which we will assume to be quicksort. Quicksort will be discussed later; for now let us just mention that an optimal choice is obtained through a combined quicksort/insertion sort and is described in detail as Algorithm Q by Knuth [Knuth 73a, pp. 118-19]. Once the local (processor) sort is complete, it must be decided how to merge all of the sorted lists in order to form one globally sorted list. This is done in a series of compare-exchange steps. In each step, two neighboring nodes exchange items so that each node ends up with a sorted list and all of the items in one node are greater than all of the items in the other. Thus, two sorted lists of m items each are merged into a sorted list of 2m items (stored collectively in the memory of the two nodes). 

18-3 The Compare-Exchange Routine 
The concept of the compare-exchange process between two processors is fairly simple, though we will see that what actually goes on between the processors becomes slightly tricky. The idea is illustrated in Fig. 18-1. The start is shown in Fig. 18-1(a). We have two processors which have each sorted their internal lists of five numbers. We wish to merge so that the processor on the right ends up with the high half of the collec- tive list, the one on the left ends up with the low half. The high processor will begin by sending its lowest number over to the empty buffer in the low processor (the empty box), offering to give this number over. The low processor will do the reverse, resulting in the situation shown in Fig. 18-1(b). 
The low processor now has what it needs to begin constructing the new list-it knows its own lowest member, and it has the other processor's lowest member also. The 7 and 10 are compared and it is decided that the new list should begin with the 7. The new list is constructed in a currently free section of memory and this is shown in Fig. 18-1(c). The high processor has done likewise, except that it works from the bottom up and has not yet used the offered number residing in the communications buffer. It will use it on the next cycle however, when the next pair of numbers are exchanged, leading to Fig. 18-1(d). The process continues, Figs. 18-1 (e)-(g), until the process terminates when the size of the new lists are the same as the old. 
Note that the new lists are constructed by actually copying the items from one memory location to another. If the items to be sorted were large data structures of many bits, we would, instead, deal with pointers to the objects and not the objects themselves. The objects would never move; lists of pointers would be manipulated. 

330 

(b) 

17 ~ : 14 10~ : 16 25 : 47 60 : r87 65 i@]92 

(e) 
m70f--@ 10 : 47 14 : 65 25 : 87 ~92 

Sorting on the Hypercube Chap. 18 

(c) 

(f) 

7 
j6 10 
14 25 60 

47 
£65 87 60 92 

Agure 18-1 Compare-exchange step between nodes. This figure illustrates a simple canpare-exchange betWeal two coonected nodes. Each processor contains a sorted list of five items and the goal of the canpare-exchange process is to put the largest five items in the right-hand processor and the smallest five it~s on the left. In (a), each processor sends an "unwanted" item to a register in the other processor. The right-hand processor, which is collecting the large items, sends its smallest it~, and the left-hand processor sends its largest. In (b), each processor inserts its new it~ into the list. One item is pushed" off the end" of each list, but that item is already residing safely in the other pro- cessor. This two-step process is repeated twice more (c, d, e, and f). After (f), the items are all in the correct final position, OOt the processors don't realize it yeL When the next exchange of items takes place, in step (g), neither processor will want the received item, so the canpare-exchange process will be terminated. 
In our description of the compare-exchange process, we were a bit sloppy in speci- fying when the items were actually transferred. The point is that, in a fast-polled com- munication environment, communications do not occur until both processors agree that more items need to be b"ansferred. This, in turn, means that sometimes one processor will wait while the other catches up, giving rise to load imbalance. In the above example, this occurs at the step of Fig. l8-l(c): the low processor has emptied its buffer and would like the next item to be b"ansferred, while the high processor is still busy. To avoid this, the actual implementation should b"ansfer many items at once, reducing the frequency of such occurrences. This is also a good approach because in almost any hardware imple- mentation, communications are cheaper (per bit) if many bits are sent in one large burst, than if sent as separate smaller messages. 

Sec. 18-3 The Compare-Exchange Routine 

331 

We next examine the implementation. The reader more interested in concepts than implementation details may wish to skip to the next section. Here, we will present the compare-exchange code and show how it works for the above example. This description is instructive in that it shows that relatively simple code can exhibit rather complex behavior and that one must think carefully about what is going on in the parallel machine in order to obtain a fast, efficient implementation. 
The routine cmp -ex-hi is called by a node which performs a compare-exchange step with a neighboring node, keeping the high part of the merged list. Similarly, cmp -ex-/ow is what the node on the other side runs so that it keeps the low half of the merged list Note that, in keeping with polled communication, if a given node calls cmp -ex-hi , its corresponding neighbor must call cmp -ex-/ow concurrently. 

proc-begin cmp -ex-hi ( sort items in 2 processors so this processor contains highest items ) declare-int nitems; ( number of items in this processor's list ) declare-int nbuj; ( number of items in communication buffer ) declare-buf old-Iist; ( room for nitems items, initial list ) 
declare-buf new-Iist; ( room for nitems items, sorted list ) 
declare-buf comm-but; ( room for nbutitems, communications buffer) 
[divide old-Ii;r1 into blocks of nbufitems, with smallest items in first block ] [exchange first block of old-Iist with block from other processor, store incoming block in comm-buj] [set i = (nitems- 1)] [setj = (nbuf-l)] 
for-begin ( k = (nitems -1) to O ) 
if begin ( ilh item in old list ?;jlh item in comm but) then -[klh item in new-Ii~ +-ilh item in old-Iist] [i +- (i-l)] 

else 

[klh item in new-list +- jlh item in comm-bu.f] 
if-begin ( j = O , i.e. all items in comm-bufhave been processed ) then [exchange next block of old-list with next block from other processor, 
store incomin~ block in comm-bu.f] [setj = (nbuf-l)] 

else 

u+- (j -I)] if end 
if end for end proc-end 

Code 18.1 The Routine cmp -ex-hi 
These routines are outlined in Code 18-1 and 18-2 respectively, while their opera- tion, for the example of Fig. 18-1, is illustrated in Fig. 18-2. The illustration is for the case of a buffer holding two items. The indices i and j are employed in this illustration to show in detail how load imbalance occurs for this particular choice of buffer size. By 

332 

Sorting on the Hypercube Chap. 18 

the time step (d) has been reached, the low processor has emptied its buffer. The proces- sor will have to wait at the exchange step within the k loop of cmp -ex-/ow. This waiting occurs because of the nature of polled communications: the low processor tries to exchange items, but the high processor has not yet reached the corresponding exchange, so the low processor must stand idle. The high processor continues processing through step (e), when it finally empties its buffer. It now reaches the exchange step and the items are communicated. The entire process stops when each processor has gone through the k loop nitems times; this means that nitems items have been assigned to the new list. 

proc begin cmp ex low ( sort items in 2 processors so this processor contains lowest items ) declire-lot nitems; ( number of items in this processor's list ) declare lot nbuf; ( number of items in communication buffer ) declare-bur old list,' ( room for nitems items, initial list ) 
declare-bur neW list,' ( room for nitems items, sorted list ) 
declare=bur coWn-buj;' ( room for nbufitems, communications buffer ) 
[divide old list into blocks of nbufitems, with largest items in first block] [exchange first block of old list with block from other processor, store incoming block in comm-bu.f] 

[set i = 0] [setj = 0] 
for begin ( k = oto (nitems -I) ) 
-If-begln ( il/o item in old-Iist ~jl/o item in comm-buf) then [kl/o item in new list +- il/o item in old list ] [i +- (i + I)] -- 

else 

[kIlt item in new list f- jIlt item in comm bu.f] 
if begIn ( j = ( ;ibuf- 1 ), i.e., all items in -comm bufhave been processed ) then -[exchange next block of old-Iist with next block from other processor, 
store incoming block in comm bu.f] [setj = 0] - 

else 

u ~ u + 1) ] If end 
If end - 
for end- proc-end 

Code 18-2 The Routine cmp -ex-/ow 
The number of times such processor idling takes place will be lessened if we increase the buffer size. In the exb'eme limit, if the buffer size is set equal to nitems, only one exchange need take place per compare-exchange and no waiting occurs. This approach, however, is memory-intensive for sorting large lists and it means that many items (on average, half of them) are transferred unnecessarily. 

Sec. 18-3 The Compare-Exchange Routine 

333 

Details of Compare-Exchange 'I..>w" Processor: "High" Proc~. 

Communicate . 

j~ l11J 

:1651 ~:~J 
, 

newlist: 

i --')0 10 25 (b) 60 65 87 
10 i --')025 (c) 60 65 
87 

10 j~25 (d) 60 65 
87 

~ 

7 [71 10 j -;o!1.4J 

10 14 ~ 

:;:J 

16"51 
~j 92 

~j 87 !§11 92 

67 ["6"5lE- j 87 !§1J 92 

(f) 10 i~25 60 
65 87 

(g) 10 i~25 60 
65 87 

7 10 14 

10 14 16 

(h) 10 25 i~ 65 
87 

7 10 14 16 
~25 

["i6\ ~ 

60 65 67 
~j 87 If!Q\ 92~ 

oldlisl: 7 14 16 67 92~ 

67~i 92 
7 14 16 67~i 92 
7 14 l~i 67 92 

7 14 I~ 

67 92 

16 67 92 

92 " 

67 92 

l'lgure 18.2 The time steps in the implemmtatioo of the compare exchange algorithm for two lists of 5 numbers to be sorted so that the low 5 md up in the left and the high 5 end up in the right processor. The numbers start in the arrayold and end in new. 

334 

Sorting on the Hypercube Chap. 18 

18-4 The Bitonic Algorithm 
Many algorithms for sorting on concurrent machines are based on Batcher's bitonic sorting algorithm [Batcher 68; Knuth 73a, pp. 232-33]. Bitonic sort is not a particularly good algorithm for sequential machines, but it is good for SIMD machines because it allows a high degree of concurrency without requiring independent flow control in the nodes. As mentioned earlier, the method sorts the list into ascending procnum order. 
The first step in the merge strategy is for each node to internally sort via quicksort. One is then left with the problem of constructing a series of compare-exchange steps which will correctly merge N = 2d sorted sublists. This problem is completely iso- morphic with the sequential problem of sorting a list of 2d items by pairwise comparisons between items. Each one of our sublist compare-exchange operations is equivalent to a simple compare-exchange between two individual items. We can construct a concurrent algorithm based on any sequential algorithm which uses the comparison strategy. In the algorithm specification of Code 18-3, cmp -ex-hi U) is called by a node which is to undergo a compare-exchange step with the neighboring node on channel j, keeping the high part of the merged list Similarly, cmp-ex-lowU) applies to a compare-exchange step with the node on channel j, with the calling processor keeping the low part This pattern of compare-exchanges is shown in Fig. 18-3 for the d =3 case. It is not intuitively obvious why the bitonic algorithm works at all. We merely give the d = 3 example and appeal to the references to give a proof that it does, indeed, sort. 

Code 18-3 Routine for Perfonning Bitooic Merge 

Sec. 18-4 The Bitonic Algorithm 

(a) 

(c) 

(e) 

335 

(b) 

(d) 

(t) 

l'lgure 18.3 Bitonic scheme for d=3. This figure illustrates the six C<Xnpare-exchange steps of the bitonic algorithm for d=3. Each diagram illustrates four cOOlpare-exchange processes which happen simuhaneously. A boldface arrow represents a compare- exchange between two processors. The largest items go the the processor at the point of the arrow. and the smallest items to the ooe at the base of the arrow. 

336 

Sorting on the Hypercube Chap.18 

18-5 Comments on the Efficiency of the Bitonic Algorithm 
It is fairly easy to analyze the expected perfonnance of this algorithm. If we are sorting M items, each node sorts M 12d items in the initial phase, taking time 
M M ts(M, d) ='t d-In(d) (18-1) 2 2 

where 't is an implementation-dependent time constant. In the merging phase, each node performs d (d + 1 )/2 cmp -ex-/ow or cmp -ex-hi operations. The time required by each step is a linear function of the number of items per node. However, we cannot assume that each of the steps takes the same amount of time; as the merge phase proceeds, the list begins to assume global order, so the steps should require less time. We write for the total merging time: 

tm(M, d) = f(d)t d 2 
The function f (d) is difficult to derive from a theoretical basis. We expect it to be pro- portional to da with 1 < a < 2. The lower bound is due to the minimum communication required to get each item to the proper node. Suppose that the algorithm is such that each item follows the optimum (shortest) path in going from the node in which it initially resided, to its ultimate destination. Since this distance grows linearly with d on a hyper- cube, a lower bound on a is unity .On the other hand, the bitonic algorithm involves a number of cycles quadratic in d, so that if all cmp -ex-hi and cmp -ex-low steps took the same amount of time, a would be 2. Since these comparisons actually become faster as the list becomes ordered, this value represents an upper bound on a. 
We can now derive an expression for the efficiency. By definition, we have 
1 ts(M,O) "0 ~, 

M 

(18-2) 

E= 
2d ts(M, d) + tm(M, d) 
Using the above equations, we obtain the result 

e= f(d)-d \..u~J 1+ 
lnM 
If we fix d and let M become arbitrarily large, the efficiency approaches 1. However, the approach to this limit is very slow due to the logarithmic dependence. If we increase d with M fixed, the efficiency goes to zero. We do not know f (d) a priori, so actual timings must be employed to evaluate it. 
Before presenting the performance measurements of the bitonic algorithm, it will be useful to explain the derivation of the efficiency measurements. The most robust method is to directly measure the speedup, S, by comparing the hypercube timing with the same sorting procedure running on a single processor of the hypercube, and then using the formalism described in Sec. 3-5 to find e. This was followed in the current examples, except for the larger sorts that would not fit within the memory (256 kbytes on 

Sec. 18-5 Comments on the Efficiency of the Bitonic Algorithm 

337 

the Caltech/JPL Mark II) of a single node. To circumvent this problem, single processor times were extrapolated via the theoretical functional form for sequential quicksort, M 1nM .This extrapolation is shown in Fig. 18-4. As seen from this figure, the fit to the theoretical form is rather good, so we can, with confidence, use this fit to extrapolate and give single-node times for large lists. 
Quicksort 

Time (see) 

20 

16 

8 

0 

log2M 
F1gure 18.4 The quicksort timing data for a single processor. Also shown is the M In M fit. 
Table 18-1 shows the actual times and efficiencies for our implementation of the bitonic algorithm. Results are shown for sorting lists of sizes (M) lk to 2048k items on hypercubes with dimensions, d, ranging from 1 (two nodes) to 7 (128 nodes). The same information is also shown, graphically, in Fig. 18-5. The efficiency data are fit fairly well by the functional form suggested above, with (X::: 1.9 . 
Clearly, the efficiencies fall off mpidly with increasing d. From the standpoint of cost effectiveness, this algorithm is a failure. On the other hand, a glance at Table 18-1 shows that for fixed size lists and increasing machine size, the execution times continue to decrease. Even for the very small 1024-item list, the 7-cube still outperforms the 6- cube: the classic rollover point (times actually increasing as machine size is increased) is not reached until there are just a few items per processing node. So, from the speed at any cost point of view, the algorithm is a success. We attribute the inefficiency of the bitonic algorithm partly to communication overhead and some load imbalance during the compare-exchanges, but mostly to nonoptimality of the algorithm itself. Remember, in our definition of efficiency we are comparing the parallel bitonic algorithm to the best available sequential algorithm-in this case, quicksort In bitonic, the number of cycles grows quadratically with d. This suggests that the efficiency can be improved greatly by using a parallel algorithm that sorts in fewer operations without sacrificing concurrency. 

338 

(a) 

Sorting on the Hypercube Chap. 18 

Bitonic vs Sequential Quicksort 

IOg2M 
Figure 18-5(3) The efficiency of the bitonic algorithm versus list size for various size cubes. 

Figure 18-S(b) Same as Fig. 18-5, but with the fit to Eq. (18-4) shown. 

Sec.18-6 Shellsort/Diminishing Increment Algorithm 

339 

Execution times (sec) 

18-6 Shellsort/Diminishing Increment Algorithm 
This algorithm again follows the merge strategy and is motivated by the fact that d compare-exchanges in d different directions of the hypercube result in a list with an almost sorted order. Global order is defined via ringpos, that is, the list will end up sorted on the embedded ring in the hypercube. After the d compare-exchanges, the algo- rithm switches to a simple mopping-up stage which is specially designed for almost sorted lists. This stage is optimized for moving relatively few items quickly through the machine and amounts to a parallel bubblesort or bucket brigade algorithm. Each node passes its largest items upward along the line and its smallest items downward. After each stage of this process, an efficient hypercube broadcast is done which decides whether any items remain to be propagated along the ring. 
Note that the mop-up algorithm takes advantage of the MIMD nature of the machine and that this characteristic is central to its speed-only the few items which need to be moved are examined and processed. The bitonic algorithm, on the other hand, is natural for a SIMD machine. It involves much extra work in order to be able to handle the worst case, which rarely occurs. 
We will refer to this algorithm as Shellsort [Shell 59; Knuth 73a, pp. 84-5, 102-5] or a diminishing increment sorting algorithm. This is not because it is a strict concurrent 

Sorting on the Hypercube Chap. 18 

implementation of the sequential namesake, but because the algorithms are similar in spirit. The important feature of Shellsort is that in early stages of the sorting process, items take very large jumps through the list, reaching their final destinations in few steps. As shown in Fig. 18-6, this is exactly what occurs in the concurrent algorithm, and, as we shall see below, is responsible for its speed. 

Shellsort 

(a) 

(b) 

On ring 

1~1 00000000 01326754 

J,~J,J,~J, 00000000 01326754 

/ (d) Othenp uh n o 
01326754 
Figure 18.6 The parallel Shellsort on a d=3 hypercube. The left side shows what the al- gorithm looks like on the cube, the right shows the same when the cube is regarded as a ring. 
The complete algorithm is outlined in Code 18-4. The for loop comprises the actual Shellsort process; the while is the mopping-up operation. 

~ 

~ 

0 

~ 

(c) 

~~~~ 00000000 01326754 

Sec. 18-7 Perfonnance of Shellsort 

341 

proc-begin shell-merge ( merge sorted sublists using Shellsort ) declare-int d; ( the dimension of the hypercube) 
declare-int ringyos; ( position in ring of processors ) for begin ( i = (d-l) to O ) 
-If-begln (ilh bit of ringyos is I) then proc-caU cmp -ex-hi (i ) 

else 

proc-caU cmp -ex-low (i ) if end for end 
[let done be FALSE] 
while begin ( done is FALSE) then 
ir-begin ( QIl' bit of ringyos is 1 ) then 
proc-caU cmp -ex-hi ( do compare-exchange with the previous processor in ring. keeping highest items in this processor ) 
proc-caU cmp -ex-low ( do compare-exchange with the next processor in ring. keeping lowest items in this processor ) 

else 

proc-call cmp -ex-/ow ( do compare-exchange with the next processor in ring, keeping lowest items in this processor ) 
proc-call cmp -ex-hi ( do compare-exchange with the previous processor in ring, keeping highest items in this processor ) 

if end 
if-begin ( received items need to be passed further ) then [broadcast FALSE to all processors ] 

else 

[broadcast TRUE to all processors] if end 
if-begin ( all processors broadcast TRUE) then [set done = TRUE] 
if end whIle end proc-end 

Code 18.4 Routine for Petfonning Shellsort Merge 

18.7 Performance or Shellsort 
We have not tried to predict the Shellsort efficiency theoretically, though this could certainly be done for both the average and worst cases. The algorithm was implemented and timed with the same starting problem used for the bitonic algorithm. The results appear in Table 18-2 and are also plotted in Fig. 18-7. 
This algorithm is much more efficient than the bitonic algorithm, and seems to offer the prospect of reasonable efficiency at large d .The inefficiency is the result of both communication overhead and algorithmic nonoptimality relative to quicksort. 

342 

Sorting on the Hypercube Chap. 18 

Executioo times (sec) 

Total Lisl Size 

, .~ .54 1.20 2.52 5.51 11.73 
.15 .31 .68 1.43 3.13 6.40 13.44 
.07 .16 .36 .76 1.62 3.30 7.03 14.96 
.04 .09 .19 .40 .85 1.82 3.70 7.91 16.66 
.05 .07 .11 .23 .44 .91 1.93 4.01 8.47 17.51 
.07 .08 .09 .16 .27 .52 1.05 2.12 4.36 9.15 18.81 
.21 .16 .16 .17 .23 .37 .67 1.30 2.44 5.23 10.22 21.03 
Efficiency £ 
.92 .94 .92 .95 .94 .94 
.77 .82 .81 .84 .82 .86 .88 
.82 .79 .77 .79 .80 .84 .84 .84 
.72 .70 .73 .75 .76 .76 .80 .79 .80 
.29 .45 .63 .65 .73 .76 .76 .78 .78 .80 
.10 .20 .38 .47 .60 .66 .70 .74 .76 .77 .78 
.02 .05 .11 .22 .35 .47 .55 .60 .68 .67 .72 .74 
Table 18-2: Shellsort 
The worst case mop-up phase may be summarized as follows. Suppose we seek to sort M items on a d -dimensional hypercube, so that there are n =M /2d items per node. Now, if more than n of the largest (or smallest) 2n items happen to occupy a particular pair of nodes in the first compare-exchange step, then some of the items will be put in the wrong plane of the hypercube in this first step. The subsequent d -1 stages will just rear - range the items in each (d-1)-dimensional subcube, but these wrong way items will stay in the wrong half of the hypercube. Consequently, in the mop-up phase, these items must travel a long distance to reach the top (or bottom) of the list, requiring 2d-1-1 steps. The probability that such a collision of items will occur, for random initial lists, is given by: 
p -(1)"-1 (2n)! 
(185)collisio" -- 
2d-1 II - n.n. 
This probability rapidly becomes negligibly small as n increases. In the other hypercube applications presented in this book, it is a general rule that algorithms become more efficient as the grain size, n , grows. That is also true here, but a qualitatively different effect has entered the picture. The probability of hitting a "bad" case decreases with increasing grain size. For the random lists which were run in the examples, the mop-up step was separately timed. The mop-up time, not surprisingly, increases with increasing d, but even for d =7 is only 10% to 20% of the total time for the larger lists. For the small lists on the d=7 cube, the mop-up time dominates. 

d 

lk 2 

4 

8 

16 

1 2 3 4 5 6 7 

32 

64 

128 

256 

512 

1024 2048k 

Sec. 18-8 Quicksort/Samplesort Algorithm 

343 

0.8 ~ 

0.6 

"Shellsort" vs Sequential Quicksort 

000000 
O O AOOOOAAA 
OA AA 'V 'V 'V X A'V'V~XXX 
'V 'V 'V X +++ 
+ 
+ o o 

0 

0 

log2M 

FIgure 18.7 Same as Fig. 18-5, but for the Shellsort algorithm. 18-8 QuicksortlSamplesort Algorithm 
The classic quicksort technique is a divide-and-conquer sorting method [Hoare 62; 
Knuth 73a, pp. 118-23]. As such it would seem to be amenable to a concurrent imple- mentation, and with a slight modification (actually, an improvement of the standard algo- rithm) this turns out to be the case. 
The standard algorithm begins by picking some item from the list and using this as the splitting key. A loop is entered which takes the splitting key and finds the point in the list where this item will ultimately end up once the sort is completed. This is the first splitting point. While this is being done, all items in the list which are less than the split- ting key are placed on the low side of the splitting point, all higher items are placed to the high side of the splitting point This completes the first divide. The list has now been broken into two independent lists, each of which still needs to be sorted. Now the algo- rithm is repeated recursively, splitting each of the two sublists, and so on, until the sort is complete. 
The essential idea of the concurrent (hypercube) quicksort is the same. The first splitting key is chosen (a global step to be described below) and then the entire list is split, in parallel, between two subcubes of the hypercube. All items higher than the split- ting key are sent in one direction in the hypercube, and all items less than the key are sent the other way. The procedure is then called recursively, splitting each of the subcubes , lists further. As in Shellsort, the ring-based labeling of the hypercube is used to define globalorder. Once d splits occur, there remain no further interprocessor splits to do, and the algorithm continues by switching to the internal quicksort mentioned earlier. An 

344 

Sorting on the Hypercube Chap. 18 

illustration of this algorithm is given in Fig. 18-8 

IL' 

I split 1 

Quicksort 

Cube space: 

List space: 

(a) -7 - 

-~- 

(b) 

-~- 

(c) 

~-- ~-- ~-- ~-- 

(d) 

then, internal quicksort. 

Figure 18.8 An illustratioo of the parallel quicksort. 
So far, we have described the standard quicksort. For quicksort to perfonn well, even on sequential machines, it is essential that the splitting keys be chosen so that the splitting points land somewhere near the median of the list. If this isn't true, quicksort may behave rather poorly, the usual example being the quadratic time that standard quicksort takes on almost-sorted lists. To counteract this, it is a good idea to choose the splitting keys with some care, so as to make evenhanded splits of the list. 
This point becomes much more important on the concurrent computer. In this case, if the splits are done haphazardly, not only will an excessive number of operations be 

Sec. 18-8 Quicksort/Samplesort Algorithm 

345 

necessary , but large load imbalances will also occur. If half of the hypercube has a sub- list three times as long as the other half, much idle time will be incurred. Therefore, in the concurrent algorithm, the splitting keys are chosen with some care. The present ver- sion of the concurrent algorithm randomly samples a subset of the entire list (giving an estimate of the true distribution of the list) and then picks splitting keys based upon this sample. To save time, all 24-1 splitting keys are found at once. This modified quicksort should perhaps be called a parallel samplesort and contains the following steps: 
(1) each processor picks sample of I items at random. (2) sort the sample of I 24 items using the Shellsort algorithm. (3) choose splitting keys as if this was entire list. 
(4) broadcast splitting keys so that all nodes know all splitting keys. (5) perform the splits in the d directions of the hypercube. 
(6) each node quicksorts its sublist Correct choice of the sample size, I, is important A large value of I increases the time spent sorting the sample, but tends to reduce the load imbalance in the splitting phase. Experimentation with various values of I has shown that the best value depends on the number of items per processor. Times and efficiencies for the parallel quicksort are shown in Table 18-3. The efficiencies are also plotted in Fig. 18-9. The times shown are for fairly good choices of the sample size, I. The first part of Table 18-4 shows the choices for I (in terms of number of items per node sampled) which correspond to the runs of Table 18-3. 

Executioo times (sec) 

Total List Size 

.27 11.37 
.18 6.10 12.71 
.12 3.32 6.69 13.25 
.09 1.83 3.80 7.46 15.64 
.06 1.02 2.02 3.94 8.26 19.45 
.08 .61 1.11 2.12 4.26 10.11 25.87 
.20 .45 .76 1.29 2.38 5.59 13.66 37.91 
Efficiency £ .85 .84 .94 .95 .98 .97 
.64 .75 .78 .85 .88 .91 .93 
.48 .63 .64 .72 .81 .83 .88 .95 
.32 .49 .58 .68 .67 .76 .78 .84 .85 
.24 .35 .46 .53 .63 .68 .73 .80 .80 .72 
.09 .14 .27 .37 .47 .57 .66 .74 .78 .69 .57 
.02 .05 .10 .19 .34 .38 .49 .61 .70 .63 .54 .41 
Table 18-3: Quicksort 
There are two main causes of inefficiency in this algorithm. The first is a result of the time wasted sorting the sample. The second is due to load imbalance in the 

lk 

2 ~ 

2 

4 

6 '7 

.ro .34 .20 .13 .00 .11 .16 

8 

1.18 .71 .43 .24 .15 .13 .18 

16 

2 3 4 5 6 7 

2.51 1.41 
.83 .44 .28 .20 .20 

32 

5.24 2.93 1.59 .96 .51 .34 .24 

64 

128 

256 

512 

1024 

346 

Sorting on the Hypercube Chap. 18 

Parallel Quicksort vs Sequential Quicksort 

<> <> 

+ o 0 
O . 10 I: 

20 

22 

og2M 

F1gure 18.9 Efficiency data for the parallel quickson described in the texL 
splitting/sorting phase, as a result of uneven distribution of items among the processors. B y varying the sample size 1 , we achieve a trade-off between these two sources of inefficiency. 
The second part of Table 18-4 conveys some information about the amount of load imbalance for each run. As an algorithm becomes load imbalanced, the efficiency is degraded by the ratio of average processor sort time to the maximum of all processor times. This relation holds as long as the sum of all processor times remains a constant with varying configurations. In other words, it holds as long as the computation time of a node depends linearly on the number of items in the node. This isn't precisely correct in our case (the dependence includes logarithmic factors) but the balance factor, average/maximum, is still generally indicative of the amount of load imbalance. This fac- tor is determined in the examples by employing timing routines independently in each node, yielding a timing profile of the concurrent machine. 
The striking feature of the quicksort results, in comparison to those of Shellsort, is that while the parallel quicksort does surpass the latter for a large range of list and hyper- cube sizes, for the largest lists on the largest cubes, it exhibits a sharp rolloff in perfor- mance. This behavior is a direct consequence of the statistical nature of the sampling procedure for estimating the medians of the list. Because the method is statistical, the estimated median converges to the true value only as 1/..![ .This is reasonably effective for the smaller lists, and does the job well enough so that the performance surpasses that of Shellsort. For the largest lists on the largest cubes, however, precise load balancing becomes critical and the algorithm rapidly succumbs to imbalance. It is clear that a more 

Sec. 18-8 Quicksort/Samplesort Algorithm 

Sample size per node 

Sorting on the Hypercube Chap. 18 

18-9 Concluding Remarks 
There seem to be several reasonable algorithms for sorting on a hypercube. Batcher's bitonic method is the leader of the low-efficiency, speed at any cost algorithms and it is the best method for small lists (measured relative to the machine size). The Shellsort- and quicksort-inspired algorithms, on the other hand, can give close to linear speedups for large lists. An important enhancement of the quicksort approach is the incorporation of a deterministic algorithm for finding the splitting points. 
It is interesting to note that the two most efficient algorithms utilize the MIMD nature of the hypercube. As in sequential sorting, good algorithms make dynamical adjustments during execution. On a SIMD machine, one would be stuck with something like the bitonic algorithm, sacrificing a great deal of possible speed. 
Quicksort is a very natural algorithm for the hypercube. The first d divide-and- conquer splits correspond nicely with the d dimensions of the hypercube. In addition, each datum arrives at its ultimate destination along a cube geodesic, that is, along a shor- test path. The communication pattern of quicksort on a hypercube is the same as that of the crystal-router, an important communication utility to be discussed in Chap. 22. 
This chapter has concentrated on sorting lists of unknown data distribution and low degeneracy (few identical items). Other cases require specialized algorithms. When the input data distribution is known, bucketsort (with small buckets) is most efficient on a sequential machine. A simple generalization of the quicksort (remember, it can be thought of as a bucketsort) to many buckets per processor would handle this case. One would simply replace the sequential quicksort phase with a sequential bucketsort in each processor. 
The high degeneracy case is handled best by an algorithm which simply counts the number of items of each type. In this case, we could perform the same counting within each node, and then use an efficient global sum algorithm to add up the total result. 

19 

Scalar Products: 

19-1 Introduction 
Since the topological structures of many interesting algorithms can be represented by various types of b"ees, we will investigate the mapping of b"ees onto a hypercube in this chapter. Among the applications that use b"ee algorithms are computer chess, sym- bolic manipulation, and more general artificial intelligence applications, such as expert systems. Although a single b"ee, such as the binary b"ee shown in Fig. 19-1, maps natur- ally onto a hypercube, significant load imbalance occurs in this decomposition. Rather surprisingly, we will see that a collection of many b"ees maps onto a hypercube both naturally and efficiently. Many b"ees arise in a concurrent b"ee-based algorithm in which the number of "leaves" is greater than the number of nodes in the concurrent processor. If a b"ee algorithm with a topology as shown in Fig. 19-1 is to be run concurrently on eight nodes, it could be decomposed as il1usb"ated by the dashed line. Each of the eight nodes could be placed at the top of one of the eight b"ees below the dashed line. Since relatively little of the b"ee is above the dashed line, most of the computation occurs below the dashed line and can be run concurrently. 

;\ 

A 

A Binary Tree Topology 

..-::. 

)\ 

/ \ 

" 

i. 

A 

1\I\l\Al~I~I~I~I~I~ 1~1~1~1~1~1~ 
l'lgure 19.1 A binary tree with 32 temlinal nodes discussed in Sec. 19-1. The artificial intelligence applications mentioned above are quite complicated, so that in general they would require a sophisticated operating system capable of maintain- ing load balance dynamically. Since such a system is beyond the scope of this book, we will examine a much simpler application that uses a fixed tree. As will be explained in 
349 

A A A 

~ 

~ 

A 

i ~ 

350 

Scalar Products: Mapping the Hypercube into a Tree Chap. 19 

Chap. 21, many matrix algorithms involve the calculation of several global scalar pro- ducts as a subcomponent Examples include matrix-vector multiplication and the deter- mination of the eigenvalues of full matrices. In this chapter we examine a method of cal- culating multiple scalar products concurrently on a hypercube. In Chap. 21 we will gen- eralize the addvec algorithm described here in the routines fold and expand, and will argue that these represent a class of distinct and interesting communication algorithms. We will show in Chap. 22 how they generalize to the crystal-accumulator and are very important in some circuit simulations [Funnanski 87b]. 
The scalar product of two vectors, A and B , each of which contains Q elements is given by: 

Q R = r. Ai Bi (19-1) i=l 
where R is the resulting scalar product We see that calculating a scalar product involves multiplying the corresponding 'elements of the two vectors A and B and adding the pro- ducts to obtain a single value. On a sequential computer, several scalar products are cal- culated by simply computing each one in succession. However, on a concurrent proces- sor we would like an algorithm that distributes the computation among the nodes so that many scalar products can be genemted in parallel. We will assume that each of the vec- tors being multiplied is distributed evenly among the nodes with a contiguous subset of both vectors' elements in each node. Such a distribution of the vectors is natural and recalls the treatment of Chap. 8 as well as the example matrix algorithms to be found in Chap.21. 
We first examine the calculation of a single scalar product on a concurrent proces- sor consisting of N nodes. We rewrite the vectors A and B as matrices in which column j contains the vector elements assigned to node j .For simplicity we will assume that each node has q = Q IN elements. We then rewrite Eq. (19-1) as 
N q 
R = r. r. Aij Bij (19-2) j=l i=l 
From Eq. (19-2) we see that each node can multiply its pairs of matrix elements Aij and Bij in parallel with all the other nodes. The multiplication step involves no commun- ication since each node already contains the pairs of matrix elements to be multiplied. In Eq. (19-2) the resulting products are added in two steps in contrast to Eq. (19-1), which has only one addition step. The summation over the index i adds the products within each node and runs completely concurrently since no communication is required. The summation over the index j adds the results from each node to obtain the final result. Since the second addition step combines the partial results from each node, communica- tion is required and for a single scalar product has inevitable load imbalance. 
We now discuss the calculation of multiple scalar products on a concurrent proces- sor. We start by generalizing Eq. (19-2) to the case of many scalar products by adding another index, k, which labels each of the p scalar products: 
N q 
Rk = r. r. Aijk Bijk (19-3) j=l i=l 

Sec. 19-2 Description of Algorithm 

351 

Calculating p scalar products is very similar to calculating a single scalar product. As before, both the multiplication step and the summation over the index i for each of the scalar products requires no communication and so runs completely concurrently. Label- ing the partial results from each node as rjk and rewriting Eq. (19-3) gives: 
N 
Rk = L rjk (19-4) j=l 
which is the formula for adding the vectors formed by the rows of the matrix r jk , which has N rows and p columns. Thus, we can think of the p scalar products as forming a vector, Rk' which is the sum of the vectors rk containing the p partial results from each of the N nodes. Accordingly, the algorithm that we will study in this chapter adds a vec- tor from each of the nodes and then distributes the resulting vector to all of the nodes. 
In Sec. 4-4, we introduced the routine vm-combine that takes a block of data items from each node and combines them into a single block using a combining function specified as an argument. By using a combining function that returns the sum of its two arguments, we can form the sum of each node's vectors. However, the routine presented in Sec. 14-5 that implements combine on a hypercube uses a tree algorithm that calcu- lates the same result in each node. When a large number of data items are to be com- bined, a more efficient algorithm would calculate separate results in each node and then distribute the results to all of the nodes. By distributing the calculation among the nodes, many separate results can be calculated concurrently instead of performing the same cal- culation in each node. When the calculation is finished, the separate results must be col- lected so that each node contains all of the results. 
The routine we describe in this chapter is more specific than combine because it only calculates the sum of a vector from each of the nodes. The vector is an array of floating-point values and is the same length in each of the nodes. After discussing the implementation of the routine addvec, which adds the vectors; we will analyze its efficiency and compare it with combine. Although we only examine the specific routine addvec, the underlying algorithm can be easily generalized to provide an alternate imple- mentation of combine that would be faster for other tasks that involve combining many data items. 

19-2 Description of Algorithm 
As we mentioned in the previous section, the algorithm we will use to implement addvec is similar to that used by combine but without its redundancy of computation. As with the algorithm used to implement combine in Sec. 14-5, we will map the hypercube onto a collection of binary trees with each node at the top of a different tree, but in the addvec algorithm each tree calculates separate results. The method we will describe was developed by Furmanski following a remark by Fox [Fox 86b, Cisneros 85] . 
The essential feature of the algorithm is that the communication of partial results is correctly interleaved so that the computation and communication loads are balanced optimally. Since we do not want to impose any conditions on the number of elements in the vectors being added, the algorithm should disperse the computation as evenly as 

352 

Scalar Products: Mapping dte Hypercube into a Tree Chap. 19 

possible even when the number of elements in the vectors is not a multiple of the number of nodes. In addition, the results must also be distributed to all nodes efficiently at the end of the calculation so that each node contains the complete result. 
We begin the description of the algorithm by considering the case of a 1-cube (two nodes) in which both nodes contain a vector consisting of two elements that we label a and b to avoid confusion. Instead of exchanging both elements, as was dulle in combine , node 0 sends element b to node 1, while node 1 sends element a to node 0. After both nodes add their pair of corresponding elements, each contains one of the two results. The two trees that comprise the algorithm are illustrated in Fig. 19-2. It is evident that each node is at the top of one of the trees. The individual results are then distributed by going back down the trees, with node 0 sending the result for element a to node 1, and node 1 sending the result for element b to node 0. When the bottoms of both trees are reached after the results are exchanged, both nodes contain the sum of the two original vectors. 

Addvec Algorithm on a l-cube 

Element a 

Element b 

/ "' / ~ node O node 1 node O node 1 
F1gure 19.2 The algorithm for addvec described in Sec. 19-2 for the case of two trees (scalar products) on two nodes. 
The algorithm is easily generalized to vectors containing p elements by first group- ing the elements into blocks containing two elements. If p is odd, the last block will have a missing element which will not cause any problem. For each block of two ele- ments the previous procedure is repeated, except that all the necessary exchanges for all the blOCks are performed together. Each node then calculates the result for one of the elements in each block and exchanges its results with the other node. If one of the blocks contains a missing element, the steps involving the missing element are skipped. Thus, both nodes contain the complete result of adding the two vectors of p elements, but each did only half the number of additions that would have been required by the implementa- tion of combine in Sec. 14-5. We recall that combine exchanges all of the elements in the vector once, rather than half of the elements twice as in the current algorithm, so both algorithms require the same total number of elements to be exchanged. 
We can easily extend this algorithm to implement addvec on hypercubes with more than two nodes. From the description of the algorithm on a 1-cube, we see that applying 

node O /"" 

node 1 /"\ 

Sec. 19-2 Description of Algorithm 

353 

the algorithm to a block of N elements, where N is the number of nodes, fully defines the algorithm since the generalization to several such blocks is straightforward. If the vector does not contain enough elements to completely fill the last block, the steps involving the missing elements are merely ignored. Thus, we will outline the extension of the algo- rithm to a 2-cube (where N = 4) by only considering vectors that contain four elements. 
As in the previous case, the nodes will be mapped onto a collection of binary trees with one tree per element of the vectors being added, and a different node at the top of every tree. The four resulting trees are shown in Fig. 19-3, in which the nodes are labelled with their binary processor numbers. As before, the algorithm traverses each of the trees twice, once from bottom to top to calculate the result for each element, and again from top to bottom to distribute each of the results to all of the nodes. The first step in the algorithm is similar to the previous case: node 0 sends elements b and d to node I, which in turn sends elements a and c to node 0. After adding their pairs of corresponding elements, node 0 contains partial results for elements.nr 99 10 a and c while node 1 contains partial results for elements b and d. Nodes 2 and 3 perform the same first step as if node 2 were node 0 and node 3 were node 1. In fact, if we use only the lowest bit of the processor numbers to label the nodes, the numbering is identical to that in a 1-cube in which the processor numbers only have one relevant bit. We also note that in the first step all communication occurs across channel 0. 
The second step consists of node 0 sending its partial result for element c to node 2, while node 2 sends its partial result for element a to node 0. Nodes 1 and 3 per- form a similar exchange involving elements b and d. Each node again sums its pair of corresponding elements, so that each of the nodes contains the result for one of the four elements. In the second step we note that both exchanges use channell and involve half as many elements being exchanged and half as many additions as in the first step. After the second step, the algorithm is at the top of each of the trees in Fig. 19-3. 
We now travel back down each of the trees in order to distribute the result for each of the elements to all the nodes. In the first step, which involves exchanging elements across channel 1, node 0 sends its result for element a to node 2, while node 2 sends its result for element c to node 0. Nodes 1 and 3 perform a similar exchange of elements b and d. Next, the results are exchanged across channel 0, as in the algorithm for a 1-cube, with node 0 sending elements a and c to node 1, while node 1 sends elements b and d to node 0. Nodes 2 and 3 perform a similar exchange. The algorithm is again at the bottom of the trees in Fig. 19-3, where it started, and all four nodes contain the result for each of the four elements. 
We now generalize the algorithm used in addvec to a v-cube with N = 2v. As in the previous two examples, the basic algorithm operates on blocks of N elements and maps the hypercube onto N binary trees, placing each node at the top of one of the trees. The trees are traversed from bottom to top to calculate the result for each element and then from top to bottom to distribute each result to every node. The number of steps in both directions is v and the number of the channel on which the exchanges occur con- veniently labels each step, where channel 0 corresponds to the bottom of the trees and channel v-I to the top. 

354 

Scalar Products: Mapping the Hypercube into a Tree Chap. 19 

/ 00 

Addvec Trees on a 2.cube 

Element a 

Element b 

/ 00 
/\ 

/ 

, 01 

00 /""" 

... 
10 
/\ 10 11 

/ 00 

, 01 1\ 

/ 

\ 

\ 01 

01 /' 

"" 

, 

, 
11 
1\ 

I 10 

I 

Element c 

Element d 

\ 11 

/ "' / " 00 10 01 11 
/\ /\ /\ /\ 00 01 10 11 00 01 10 11 FIgure 19.3 The algorithm for addvec described in Sec. 19-2 for the case of four trees 
(scalar products) on four nodes. 
We begin a typical step C by grouping the N elements into blocks containing 2c+l elements. Every node exchanges a single element from each group with the node across channel c. The element from each group that is received is labeled by the rightmost c + 1 bits of the node's processor number, while the element that is transmitted is labeled by the same c+1 bits of the destination node's processor number. The elements in each group are numbered from O to (2c+l-1). Figure 19-3 illustrates the pattern of exchanges by conveniently labeling the processor numbers in binary. In the first part of the algo- rithm, when the trees are being traversed from bottom to top, the received element for each block is added to the corresponding element in the block. Since each step involves half as many blocks as the previous step, each step involves half as many additions and requires half as many elements to be exchanged as the previous step. In the second part of the algorithm, where the order of the steps is reversed, each step involves exchanging twice as many elements as the previous step. However, the steps in the second part of the algorithm require no additions since the received values are the results for each element, and so are simply stored in their appropriate locations. 

10 /"" 

/ 

11 /" 

Sec. 19-3 Implementation of Algorithm 

355 

19-3 Implementation of Algorithm 
We will now examine the actual implementation of the routine addvec, which adds up N contributions to p vectors and returns the p results to all N of the nodes. Even though the algorithm is rather complicated, its implementation is relatively simple because the sequence of exchanges and additions is so regular and fits the hypercube topology perfectly. As we noted in the previous section, the number of elements in a block is different for each step in the procedure, so the subset of elements being exchanged is different in every step. Since the elements to be exchanged cannot be con- veniently gathered into a contiguous block at the beginning of the routine, we have a good example of a case in which the communication routines that transmit and receive vectors are useful. As usual, the nodes will always be transmitting and receiving together, so we use vshift in the implementation of addvec. 

proc begin addvec ( add vector, V, of length p over hypercube ) 
declare bur T; ( same length as V, for storing data received from other nodes ] declare-lnt S; ( the number of elements sent ) 
declare-lnt R; ( the number of elements received ) 
rOr-begln ( each dimension, c, from O to (doc-l) ) 
[ divide Vand T into blocks of length b =2c+l ] [ r a = the lowest c+ 1 bits of procnum ] [ sa = r a with bit number c flipped ] 
[ R = number of blocks having an element at position r a ] [ S = number of blocks having an element at position s a ] 
proc-caU v shift ( send S elements on channel c, one from position s a of each block of V; receive R elements on channel c, and store in T at position r a in each block ) 
[ add to the element at position r a of each block of V the corresponding element in T ] 

for end 
for=begln ( each dimension c. from (doc-l) to O ) 
[ divide V and T into blocks of length b =2c+l ] [ So = the lowest c+l bits of procnum ] [ r o = s o with bit number c flipped ] 
[ R = number of blocks having an element at position r o ] [ S = number of blocks having an element at position s o ] 
proc-caU v shift ( send S elements on channel c. one from position s o of each block of V; receive R elements on channel c. and store in Vat position r o in each block ) 
for end proc-end 

Code 19.1 The addvec Algorithm in Pseudocode 
In order to use vshift we must specify the location of the first element to be sent, the size of the elements, the separation between elements that we want to send, and the total number of elements to be sent In addition, the same information must be provided for the elements to be received by vshift. The size of each of the elements is given by sizeoj(float), which returns the number of bytes in a variable defined to be a float. The offset between both the elements to be sent and those to be received is sizeoj(float) times 

356 

Scalar Products: Mapping the Hypercube into a Tree Chap. 19 

the number of elements in the block. This is just 2c+l, where c is the current channel number. The index of the first element to be received is obtained by extracting the right- most c + 1 bits of procnum, which we will recall is the node's processor number. The index of the first element to be transmitted is obtained from the exclusive-OR of the channel mask of channel c and the index of the first element to be received. The number of elements to be received is the number of blocks containing the element labeled by the receive index. Since the number of elements to be transmitted can differ from the number to be received by one if the last block is not full, the number of elements to be transmitted must be calculated separately. 
The arguments of addvec consist of but, a pointer to the array containing the vector to be added to the vectors of the other nodes, tmp, a pointer to a scratch buffer that is at least as large as the vector, and P, the number of elements in the vector. The scratch buffer is needed to hold the vector elements received from other nodes until they are combined with the partial results. The routine addvec makes the resultant vector avail- able by overwriting the contents of buf. In addition, addvec returns a value of -1 if vshitt returns an error. Otherwise, it returns a value of 0. The routine addvec is outlined in pseudocode in Code 19-1, and in the C language in Code 19-2. 
Although Code 19-2 contains quite a few lines of code, most are involved with housekeeping the various sizes and indices. The variable mask contains a mask for the bits of procnum that are used in each iteration of the loop. Thus, mask is one less than the number of elements in the current block, which is stored in block. The variable chan is the channel mask of the channel being used in the current iteration of the loop, which is also the highest bit in mask. The variable size contains the size in bytes of the floating-point variables and offset is size times the number of elements in the current block. After each iteration of the loop, chan is shifted to the next higher channel, one more bit is included in mask, and block and offset are multiplied by two. The bitwise left shift operator is used in the C implementation because it is normally faster than multiply- ing by two and it gives the same result. 
The index of the first vector element to be received, which is stored in recv, is obtained by using only the bits in procnum indicated by mask. The index of the first vec- tor element to be transmitted to the node across channel c, whose channel mask is chan, is obtained from the exclusive-OR of recv and chan, and is stored in send. The variables recvn and sendn contain the numbers of elements to be received and transmitted, respec- tively. The correctness of the formulas for obtaining recvn and sendn are easily verified by applying them to test cases. 
The actual exchange of vector elements is accomplished through the call to vshitt with the various arguments that have already been discussed. The incoming elements are placed in a scratch buffer because the original buffer may not have enough room for them. Even if we were to overwrite the outgoing elements with the incoming elements, we would not be guaranteed sufficient room, since the number of elements received may be one more than the number transmitted. Since we do not want to impose any restric- tions on P , the number of elements in the vectors, we simply require that the application program must furnish a suitable scratch buffer. 

Sec. 19-3 Implementation of Algoritlun 

357 

int addvec(buf,bnp,P) float bu£[],bnp[]; int P; 

int cJ.error,mask.chan,size. block.offset.recv .send,recvn.sendn; error = 0; 
mask = chan = 1; 
size = sizeof(float); block = 2; 
offset = size«l; 
for( c=O;c<doc; +-K: ) { 
recv = procnum&mask; send = recv.chan; 
recvn = (P+mask-recv»>(c+1); sendn = (P+mask-send»>(c+1); 
if( vshift(tmp+recv ,chan.size.offset,recvn. buf+send.chan,size,offset.sendn) == -1 ) error= -1; 
for(j=recv;j<P; j += block) buffj] += tmpfj]; chan = chan«l; 
mask = masklchan; block = block«l; offset = offset«l ; 

chan = 1«(doc-1); mask = nproc-1; 
offset = size«doc; block = l«doc; 
for( c=(doc-1);c>=0;--c ){ 
send = procnum&mask; recv = sendAchan; 
sendn = (P+mask-send»>(c+1); recvn = (P+mask-recv»>(c+1); 
if( vshift(buf +recv .chan.size.offset.recvn, buf+send.chan.size.offset.sendn) = -1) error= -1; 
mask = mask&(-chan); chan = chan»l; 
block = block»l; offset = offset» 1 ; 

} 

Code 19.2 The addvec Algorithm in C 
After exchanging the appropriate vector elements, the corresponding pairs of ele- ments are added and then the various parameters are adjusted for the next step in the algorithm. When the exchange of elements and addition of the corresponding pairs has been repeated for all the channels in the hypercube, the algorithm is at the top of each of 

358 

Scalar Products: Mapping the Hypercube into a Tree Chap. 19 

the N trees. As discussed in the previous section, each of the nodes now contains the result for the element whose index within each of the blocks of N elements is given by the node's processor number. 
The second part of the algorithm involves traversing the trees from top to bottom in order to distribute the results obtained in the first part. Thus, the order of the channels is reversed and the elements to be transmitted and received are interchanged. The various pmcuneters, such as chon, mask, and offset, are the same for each step as in the first part of the algorithm, but they are generated in reverse order. When the second series of calls to vshift is complete, but contains the result of summing the original vectors from all of the nodes. 
In Chap. 21 we will consider a further generalization of the first step. This calcula- tion of scalar products will be encountered in the routine fold. The second step, the dis- tribution of results, will be found in the routine expand in the library of generalized com- munication utilities of Chap. 21. As explained there, expand is closely related to concat introduced in Chap. 14. 

19-4 Efficiency Analysis 
Having described the algorithm for adding vectors from each of the nodes and its implementation, we are now ready to examine the concurrent efficiency of the algorithm. Regarding the routine addvec as only a part of the overall algorithm we have developed to calculate scalar products concurrently, we should properly calculate the efficiency of the total concurrent scalar product algorithm, rather than analyzing addvec in isolation. As usual, we assume that all floating-point operations require a characteristic time of tcalc and that all communication operations, such as exchanging a floating-point value with a neighboring node, involve a characteristic time of tcomm .The efficiency of the concurrent algorithm follows from Eq. (5-5): 

N TN(P,Q) (19-5) 
where T I (P ,Q) is dIe time taken by dIe sequential algorithm, T N (P ,Q ) is dIe time taken by dIe parallel algoridIm, and N is dIe number of nodes in dIe concurrent processor. 
In order to calculate T 1, we count dIe number of floating-point operations required to calculate dIe p scalar products and multiply dIat quantity by tcalc .From Eq. (19-1) we see dIat calculating each scalar product requires Q multiplications and Q -1 additions. Thus, dIe total number of floating-point operations required to calculate p scalar pro- ducts is 

£= 

T1(P,Q) 

ncalc = p (2Q -1) 
and total amount of time taken by the sequential computer is T 1 = ncalc tcalc = p (2Q -1) tcalc 

(19-6) 

(19-7) 

Sec. 19-4 Efficiency Analysis 

359 

In order to calculate T N , we break the concurrent scalar product algorithm into two com- ponents. The first component involves all of the computation that is performed without need for communication, while the second component involves the summing of partial results from each node via addvec to give the final result. We will denote the time taken by the first component as T a and the time taken by the second component as T b . 
Q is not restricted to be a multiple of N in this analysis, so the amount of computa- tion required by the first component of the algorithm is not necessarily the same in every node. While the average number of multiplications per node for each scalar product is Q IN, the maximum number is r Q IN 1 , where r x 1 is the integer part (ceiling) of x rounded to the next higher integer. Since the communication between nodes in addvec will force the nodes that perform less computation to wait for those that perform more, load imbalance occurs when Q is not an exact multiple of N. To account for the load imbalance, we must use the maximum number of multiplications per node when we cal- culate the efficiency. Accordingly, the maximum number of multiplications per node for each scalar product is r Q INl and the maximum number of additions is r Q INl -I, so the total time for the p products is: 

T a = p [2r Q 1 -1] tcalc (19-8) N 
We now examine the amount of time taken by addvec to sum the vector containing the p partial results from each node as given in Eq. (19-4). As discussed in Sec. 19-2 and 19-3, the first part of the algorithm consists of log2N steps, one for each dimension of the hypercube. In the first step half of the elements from each node's vector are exchanged with the node connected by channel O and the resulting pairs are summed by each node. For each of the channels in the hypercube, addvec performs the same type of exchange and addition using half as many elements as in the previous step. When the first part of the algorithm is complete, the resulting scalar products are distributed to each of the nodes in the second part of the algorithm. The second part of the algorithm, which involves another log2N steps, uses the same amount of communication as the first part but involves no computation. 
The number of floating-point additions in addvec is given by: 
logN . 
ncalc = L r2-' pl i=l 

(19-9) 

where we have used the ceiling function r x to account for the effect of any load imbal- ance at each step. The load balance of ~ch step must be considered separately because each step involves communication calls which resynchronize the nodes. The number of elements that are exchanged is twice nca/c since the second part of addvec requires no floating-point additions. As a result, the total amount of time taken by addvec is given by: 

logN . 
T b = (tcalc + 2tcomm ) L r 2-1 pl 

(19-10) 

360 

Scalar Products: Mapping the Hypercube into a Tree Chap. 19 

Adding the two contributions from Eqs. (19-8) and (19-10) gives the total time for the calculation of the scalar products on the concurrent processor: 
Q logN . 
TN =p [2r N 1 -1] tca/c + (tca/c + 2tcomm) :r. r2-1 pl (19-11) i=1 
Using Eq. (19-5), the efficiency of the concurrent algorithm is predicted to be: p 
(2 
Q-l)t 
/cca (19-12) 

£= 

Although Eq. (19-12) is quite complicated and cannot easily be simplified due to the integer ceiling function r x 1 used to account for load imbalance, we can simplify it considerably if we assume that the loads are balanced. In particular, we consider the case in which Q and p are multiples of N .In this case, 
1 2(N -1) tcomm 
-= 1 + (19-1~) 

e -2Q -1 tcalc ,-- --, 
As expected, me efficiency given by Eq. (19-13) is independent of me number of scalar products. When bom Q and P are multiples of N, me problem is already completely balanced, so mat increasing P yields no further improvement in me efficiency. If tcomm ::: 2tcalc and N is reasonably large, me efficiency is about 30% when each node con- tains only one element from each vector (i.e., Q = N) and increases as me number of ele- ments per node increases. wim only ten elements from each vector per node, i.e., Q = ION , me efficiency rises to about 80%. Thus, me concurrent scalar product algo- rimm yields good efficiency, even when each node contains only a comparatively small number of vector elements. 
It is desirable, however, not to impose any restrictions on P , such as requiring it to be a multiple of N in order to obtain an efficient algorithm. Therefore, we next examine me efficiency when P is not a multiple of N. Using tcomm = 2tca/c and N = 32, we have plotted me efficiency given by Eq. (19-12) for several values of Q as a function of P in Fig. 19-4. Since N = 32, me maximum efficiency is obtained at P = 32 and P = 64. The efficiency of me scalar product algorimm rapidly approaches me upper bound, and mis limiting value increases wim Q .The jumps in me plots illustrate me sensitivity of me efficiency to me load balance of me algorithm, but also show mat me sensitivity declines as P increases. Thus, even in me arbitrary case given by Eq. (19-12) and illustrated in Fig. 19-4, me scalar product algorimm implemented wim addvec yields good perfor- mance. Figure 19-5 illustrates me insensitivity of mese results to me value of P . TN(P,Q) is plotted as a function of P for me fixed values N=Q=16. These measured timings were obtained on me Caltech Mark II hypercube using an optimized assembly code version of addvec for which tcomm/tcalc = 1.5 in 64-bit arimmetic. The results fit an essentially straight line wim only minor dips at me special (optimal) points where P =16, 32, 48 (multiples of N). This graph also breaks down me timings by component, as defined at me end of Sec. 19-3. 

Sec. 19-4 Efficiency Analysis 

361 

Performance of ADDVEC 

100% 90% 80% t 70% 60% 
Efficiency 
50% 
40% 30% 20% 10% 0% ' .I. , .I, , , I, , , I, , , I, ., I, ., I. , , I 
0 8 16 24 32 40 48 56 64 Number of Scalar Products p 
Flgure 19.4 The efficiency formula Eq. (19-12) plotted as a function of P, the number of scalar products, for various values of Q, the number of terms in each scalar product. The number of nodes N is fixed at 32. 
We now calculate the efficiency that would have been obtained using the imple- mentation of combine in Sec. 14-5 in place of addvec. Each node separately calculates all of die p scalar products in combine, so that all nodes perform p log2N floating-point additions (die combining algoridim takes IOg2N steps to sum each scalar product). In addition, all p partial results from each node are exchanged across each of the log2N channels in die hypercube, so that the total number of elements exchanged by each node is also p IOg2N. Thus, die total amount of time taken by combine is 
Tb =p log2N (tca/c +tcomm) (19-14) 
For simplicity, we will assume diat Q is a multiple of N and use Eqs. (19-5), (19-7), and (19-8) to obtain the efficiency: 

~ 

N=32 

" 

5 ~ ~--- 

r.J 

I 

tcomm~ 2tcalc 

2 ~~ 

~ 

, 10 

-I 

QiN Number of vector components per node. 

(19-15) (2Q -N) tcalc + N log2N (tcalc + tcomm/ 
It is apparent from Eq. (19-15) that when combine is used in place of addvec the efficiency is independent of P , regardless of whether p is a multiple of N .For N > > 1, the simplified expression is obtained: 
1 log2N -1 log2N tcomm -=1+ 

, (19-16) £ tcalc 
where q = Q IN. Noteworthy is the appearance of a factor of log2N in both the calcula- tion and communication loads of combine, in contrast with Eq. (19-13). While the 

e= 

(2Q -I) tcalc 

+-c--- 

2q 

\ 

2q 

362 

Scalar Products: Mapping the Hypercube into a Tree Chap. 19 

20 

20-1 Introduction 
In Chap. 10 we showed how concurrent computers, and the hypercube in particu- lar, can be used effectively for matrix multiplication. Although multiplication is the sim- plest, it is not the most important matrix algorithm. Furthermore, it is this simplicity that makes the decomposition easier and the concurrent performance better than may be the case for other matrix problems. However, as stated in Chap. 10, we believe that con- current computers like the hypercube can perform reasonably well on essentially all matrix problems. In the present chapter, we will illustrate this point by discussing in detail the solution of linear systems of equations, represented in matrix notation by: 
A,!. = .Q. (20-1) 
We will examine the LU decomposition method for obtaining the solution :!. This approach does not produce the inverse of the matrix A, which is a more costly computa- tion. In Chap. 21 we will generalize the results of Chaps. 10 and 20 and discuss other key matrix algorithms more quantitatively. 
The numerical solution of Eq. (20-1) on any computer depends crucially on the structure of A. We will choose the case in which A is banded, as illustrated in Fig. 20-1, with the nonzero elements of A clustered about the diagonal. Such matrices are charac- terized by the values M and w, the matrix order and bandwidth, respectively. Here we will consider the case: 

M » w » 1 

This is an important parameter range, as the solution of many partial differential equations as discussed in Chaps. 7 and 8 can be formulated as the solution of Eq. (20-1 ) with parameters satisfying (20-2). For instance, if we consider the finite difference matrix for a N x x N y grid in two dimensions, then: 
M =Nx Ny ; w -2 min(Nx,Ny) (20-3) 
As both N x and N y are large in problems of practical importance, the inequalities (20-2) are satisfied. 
We have chosen this example because it illustrates the handling of both linear equation solving and banded matrices on concurrent processors. In the next section, we discuss conventional sequential algorithms by way of introduction to Sec. 20-3, in which we discuss the concurrent decomposition. This discussion is extended in Sec. 20-4 to a 
363- 

364 

Matrix Algorithms II: Banded Matrix LU Decomposition Chap. 20 

Structure of the Matrix A 

(a) Full Matrix 

~ 

E:m ~ (b) Banded Matrix 

~ 

EEm ~ 

Shaded area contains all nonzero matrix elements. 
Unshaded area corresponds to 
zero matrix elements. 

Figure 20-1 Comparison of zero structure of (a) Full and (b) Banded Matrix A in Eq. (20-1). 
concuuent algorithm for LU decomposition along with performance analyses and meas- urements from the implementation on Caltech hypercubes. In Sec. 20-5 we extend the algorithm to include pivoting. In Sec. 20-6 we complete the treatment of equation solv- ing by presenting the forward reduction and backsubstitution steps. These are also accompanied by experimental performance measurements. In the final sections con- cerned with systems of equations (Sec. 20-7 and 20-8), we discuss various "end effects" and then summarize the lessons from the previous sections. 
The discussion of the concuuent LU decomposition algorithm is based on work by T. Aldcroft. A. Cisneros, W. Furmanski, D. Walker, and G. Fox [Aldcroft 87]. Similar issues are discussed in [Chan 86, Chu 87, Geist 86, Geist 87, O'Neil 87, Utku 86] with particular attention to the BBN Butterfly and INTEL iPSC implementations. 

Sec. 20-2 The Sequential Algorithm 

20-2 The Sequential Algorithm 

365 

20-2.1 Overview 
In general, we will seek to solve Eq. (20-1) for a number nb of distinct right hand sides !2. .The nb solutions:! may be regarded as the columns of a rectangular solution matrix X. We can write the generalized problem in the form: 
AX =B (20-4) 
where X and B are each Mx nb matrices. Our approach [Martin 67] involves three stages which are described below: 
Stage I: LU Decomposition: 
Factorize the coefficient matrix according to 
A =LU (20-5a) 
where L is a lower triangular matrix and U is upper triangular. As we will see, L and U share a band structure similar to that of A. Equation (20-5a) is a unique factorization if we restrict L to have unit diagonal elements. Other choices are possible, as in Cholesky's method in which L and U have equal valued diagonal matrix elements. The factorization (20-5a) is useful because it is essentially trivial to invert L and U . 
Stage II: Forward Reduction Stepping through the succession of right-hand sides in B, we next make the replacement of B by BFR where: 
BFR =L -lB (20-5b) 

Stage III: Backsubstitution 
Finally, the solution vectors comprising X are obtained through: 
X = U-1BFR (20-5c) 
Following this cursory outline of the basic method, we will now proceed to detail the standard sequential algorithm, which is actually a modification of the Gaussian elimi- nation method for inversion of general matrices. We begin by noting that it is convenient and natural to store L and U in exactly the same space used by A. L overlaps the lower triangular part of A , and U occupies the remaining space. The diagonal elements of U replace those of A , and there is no need to explicitly store the unit diagonal of L .Simi- larly, B , B FR ' and X can make use of the same memory locations. Thus the algorithm amounts to a modification of the original A into L and U, while B is successively transformed to obtain B FR and finally the solution X . 

366 

Matrix Algorithms ll: Banded Matrix LU Decomposition Chap. 20 

20-2.2 LU Decomposition The sequential algorithm for LU decomposition involves M steps where step k, (k=O, ..., M -1), is illustrated in Fig. 20-2. This step involves an m x m square "win- dow" of the matrix A, where the bandwidth is given by w = 2m-l. Elements of A above the window have already undergone factorization. At the kth step, we subtract a multiple of the first row of the window, which is the k th row of the original matrix A , from the m-l following rows. In Gaussian elimination, this corresponds to the familiar strategy of zeroing the coefficients of Xk in all but the first row. Notationally this pro- cedure is expressed as follows: 
Find the kth column of L with elements given by: 
Lkk = 1 
Lk+i.k =A.tklAk+i.k 1 ~ i < m (20-6) Find the kth row of U with elements given by: 
Uk.k+j =Ak.k+j O~j <m (20-7) Modify elements of A inside the window but not in the first row or column by: 
Ak+i.k+j =Ak+i.k+j -Lk+i.k Uk,k+j O<i <m, O<j <m (20-8) 

.1 

A 

-E--nb-7 

x 

= 

= 

~nb~ 

B 

Step k of Sequential Algorithm Dark Shaded Areas are Active 
Agure 20.2 Step k of sequential algorithm showing current window. 

Sec. 20-3 Concurrent Decomposition 

367 

20-2.3 Forward Reduction 
We may choose to obtain the modification B f- BFR simultaneously with the for- mation of L and U .However, for clarity we choose hereto complete the factorization first, and then apply the following algorithm. Forward reductiori consists of M steps labeled by first row k as above. At step k, we modify elements of B according to: 
Bk+i,i = Bk+i,i -Bk,i Lk+i,k 1 ~ i < m, 0 ~ j < nb (20-9) 

20-2.4 Backsubstitution 
The backsubstitution step works from the last row to the first. We label the M steps by k, an index which now runs from M -I to 0. The algorithm now proceeds with the steps: 

Xk,i 

=Bk,j U.tk1 O~j <nb. l~i <m BL ..=BL ..-XL .UL .L 
~-J.J ~-J,J ~,J ~-1.~ 

(20-10) 

20-3 Concurrent Decomposition 
In Chap. 10 we showed that a simple square subblock decomposition was optimal. As shown in Fig. 20-3(a), this presents two serious difficulties for the problem at hand. Firstly, many of the processors are wasted, by operating on zeros outside the band. Secondly, each of the steps of the algorithms in Eqs. (20-6) to (20-10) involve only a lim- ited number of the processors. Clearly we need to ensure that each of the windows involve all processors in as load-balanced a fashion as possible. This is achieved with the scattered square decomposition illustrated in Fig. 20-3(b) [Fox 84d]. Adjacent matrix elements are assigned to different processors. We shall assume for the present discussion that the number of processors N is a perfect square. For the hypercube topology this res- tricts us to consideration of cubes of even dimension. We impose a periodic template over the matrix as shown in Fig. 20-3(c). In this, and all subsequent figures in this chapter, the value associated with a given cell of the template indicates the processor in which the corresponding matrix element is to be stored. The labeling of processors corresponds to that appropriate for a row-ordered, two-dimensional mesh topology, and is not the same as the procnum labeling scheme for hypercubes described in Sec. 14- 7. 
We note that the scattered square decomposition would in fact have been as good as the /oca/ square decomposition in the application of Chap. 10. The locality property of contiguous decompositions is not particularly important in matrix multiplication prob- lems. We previously employed the scattered decomposition in the discussion of Chap. 17 and the utility of this approach has been further demonstrated for irregular geometry finite element problems [Morison 86]. We will assume here that the topology of our con- current processor includes the periodic two-dimensional mesh connectivity. A straight- forward modification of the scattered square decomposition is required for other topolo- gies. In the case of a one-dimensional machine topology, a scattered ring decomposition would be a viable approach. 

368 

Matrix Algoritluns ll: Banded Malrix LU Decomposition Chap. 20 

Two Decompositions ofa16 by 16 Matrix Onto a 16 Node Concurrent Computer 

(a) Local Square 

(b) Scattered Square 

Each cell is a matrix element. 
The Number in each cell is the node in which this element is stored. 
The nodes are numbered 0...15. 
Flgure 20.3 C<Xnparison of (a) local square and (b) scattered square decompositioos of a 16 x 16 matrix onto a 16-node coocurrent processor. Each cell represents a matrix ele- ment and the number in the cell indicates the node in which this element is stored. 
We note that the existence of the matrix band structure does not alter the applica- bility of the scattered decomposition. The same decomposition template is used what- ever the band size. Similarly. we may define the scattered decompositions for the matrices X and B. The resultant decompositic~ for a sample problem is shown in Fig. 20-4. 

20-4 Concurrent LU Algorithm 
We now consider the concurrent implementation of the basic LU decomposition of Eqs. (20-6.7.8). We consider the kth step illusttated in Fig. 20-5. This figure shows the current window for the parameter value m=12. We now introduce the quantity: 
m= m/{f[ (20-11) 

Sec. 20-4 Concurrent LU Algorithm 

369 

(c) Scattered Decomposition Templates 

[QIiJ ~ 

4 Nodes 

16 Nodes 

Figure 20-3(c) Shows the basic template whose periodic placement leads to the scattered decomposition. 
Each processor holds n = m2 matrix elements. These elements do not comprise a con- tiguous m by m submatrix, but are rather scattered uniformly throughout the window. This illustration brings out a useful feature of the scattered decomposition. By compar- ing the windows for steps k and k+1, we see that row k and column k are active at step k but not afterward. The processors assigned to these now inactive elements obtain a compensating amount of new work (from row k+m, column k+m) in the form of the new matrix elements introduced at the next step. This illustrates the importance of the periodic property of this decomposition scheme. 
The algorithm of Eqs. (20-6), (20-7) and (20-8) consists of five stages in its con- current implementation. These are illustrated in Fig. 20-6 and are described here: 
a) Invert coefficient Au in the processor at the top left of the current window. In the examples of Figs. 20-5 and 20-6 this is taken to be processor 0. 
b) The elements Ak.k+j (0 ~j < m) are now the appropriate elements of the kth row of U with the diagonal element inverted as needed below and in the backsubstitu- tion algorithm. Broadcast this row downward with the pipe utility that was described in Sec. 10-3.2. 
c) We can now form the multipliers of Eq. (20-6). This only involves the first "column of processors." These are processors 0, 4, 8, and 12 in the example of Fig.20-6(c). The other processors are idle at this time. 
d) Transmit these multipliers from the first column of processors to their row neigh- bars. After this stage, all processors in the same row of a template (such as that of Fig. 20-3(d» contain identical subsets of the multipliers Lk+i k (0 < i < m ). 
e) Modify the matrix elements of A in the window according t~ Eq. (20-8). All infor- mation necessary to do this for a given matrix element is contained within the pro- cessor holding that element. 
Before examining this algorithm in further detail, let us mention some interesting points. We first note that in steps (a) and (c) serious load imbalance occurs. In fact, in stage (a) only one processor is doing any work. As there is but one floating-point calcu- lation (Au ~ A.tk1), this is a somewhat trivial example of a sequential bottleneck. We can attribute the load imbalance to the lack of available work. This is an example of 
~ 

370 

Matrix Algorithms ll: Banded Matrix LU Decomposition Chap.20 

1 2 3 
5 6 7 4 
91011 8 9 
131415121314 1 2 3 
lU012 3j~ Typical window 6 7 4 5 6 7 4 4 x 4 submatrix 
11 8 9 1011 8 9 
12131415121314 
230123 
6745674 
11 8 9 1011 8 9 12131415121314 
1230123 6745674 
11 8 9 1011 8 9 12131415121314 
1 2 3 0 1 2 674567 
11 8 9 10. 121314 

Bandwidth 7 

A 

! n,,= 8 

Q 

nb= 8 

16 Processors 
Figure 20.4 A banded matrix equation AX = B corresponding to the matrix order M=20 , bandwidth w=7, submatrix size m=4 and nx = nb = 8 equations. We show the Scattered Square decomposition onto a 16 node concurrent processor. Each cell represents a matrix element and the numbers in the cell indicates the node in which this element is stored. 

Amdahl's law introduced in Sec. 3-6 which states that the concurrency attainable in an algorithm is limited by the sequential parts of the problem [Amdahl 67]. Fortunately, Amdahl's law does not pose a serious problem in this case, and it is worth understanding why that is the case both here and quite generally for coarse-grained algorithms. This gives a specific example of the general analysis in Sec. 3-6. 
In Table 20-1, we analyze the three stages ofLU decomposition. Stage (a) is very inefficient but this becomes irrelevant when we compare the net concurrent time of the dominant stage (e). Stage(c) is more efficient than (a) but in fact contributes rather more to the net inefficiency given in the last column of Table 20-1. The lesson in this discus- sion is that so-called sequential bottlenecks correspond to the stages of the calculation involving very few degrees of freedom. These particular steps are implemented poorly from a concurrency standpoint but are unimportant because the concurrent stages require 

Sec. 20-4 Concurrent LU Algorithm 

371 

~ 7 !?iagonal 

Begin band 

Row k+m-l ~~ 

End band 

-step k v Next window -step k+l 

Processor number 

9 
'3 '4 
, 2 3 
5 6 7 4 
9 '0" 8 9 
'3 '4 '5 '2 '3 '4 '230'23 56745674 
9 '0" 8 9 '0" 8 9 
'3 '4 '5 '2 '3 '4 '5 '2 '3 '4 '230'230'23 

Matrix elements to be used 

FIgure 20.5 Scattered decomposition of a banded matrix of full bandwidth 23 onto 16 processors. We show the window at step k of the LU decanposition. 

times of order n , where n = ni2 is the grain size of the dominant part of the calculation. Sequential components tend to contribute terms proportional to l/n to the inefficiency, and these terms are usually small. 
A second general remark may be made with regard to comparison of the communi- cation in this algorithm with that in the algorithm of Chap. 10. In the present case, it is necessary only to broadcast rows and columns, in strong resemblance to the row broad- cast strategy of Chap. 10. On the other hand, we find here nothing analogous to the sim- ple rolling shift of the matrix B discussed in Sec. 10-3.1. The broadcast algorithm is cen- tra1 to the discussion in this chapter because we are only transmitting single rows or 

372 

Mabix Algorithms ll: Banded Matrix LU Decomposition Chap. 20 

5 Stages of LU Decomposition 

Shaded Matrix Elements Unshaded Matrix Elements 

li:j Active D Inactive 

(a) Invert Coefficient Akk of Variable to be Eliminated 

Row k ~ 12 1 2 3 O 1 2 3 0 1 2 " 

/' 

k = Ukk -7 At 

/ 

456745674561 8 9 1011 8 9 1011 8 9 1011 
1121314151213141512131415 0 1 2 3 0 1 2 3 0 1 2 -, 
456745674567 18 9 1011 8 9 1011 8 910 11 12131415 121314151213 1415r 012301230123 456745674567 
18 9 1011 8 910 11 8 9 1011 121314151213 141512131415 

Flgure 20.6(a) The first of five stages involved in step k of LU deromposition. Active matrix elements are shaded. 
columns with a message length of m= ~ , as opposed to the message lengths of m2 = n in Chap. 10. As a result, the startup time given in Eq. (10-3) is proportionately more important in the Eqs. (20-12) and (20-13) below. In addition to this machine- independent performance impact, the employment of shorter messages contributes more to overhead on machines which incur significant fixed costs for message initiation. 
A third comparison between the present algorithm and that of Chap. 10 is that both involve approximately equal amounts of communication among the matrix rows and columns. Thus the square subblock (albeit scattered rather than local) decomposition remains preferable to rectangular decompositions. In particular, pure row or column decompositions are distinctly inferior. The analysis relevant to this point follows by analogy with that in Sec. 10-4 for multiplication and need not be repeated here. 
We next analyze the performance of the LU algorithm discussed above. Up to this point we have only considered the cases in which {F[ is an exact divisor of the window order m. However, as we have seen previously in Chap. 10, significant performance degradation results in the case that m is not exactly divisible by {F[ .The working matrices (the window in A and active parts of B and X) must be padded with zeros in order to provide each processor node with a square submatrix. The worst case occurs when m exceeds an exact multiple of {F[ by one. In this case, m assumes the value [m /{F[ ], where the square-bracketed quantity is truncated to the integer part. The 

Sec. 20-4 Concurrent LU Algorithm 

373 

(b) Transmit Top Rows Downwards 

Processors 0 4 8 12 now contain 1 5 9 13 same subset 2 6 10 14 of top row 3 7 1115 
F1gure 20-6(b) The se(:ond of five stages involved in step k of LU decomposition. Ac- tive matrix elements are shaded. 

( c ) Calculate Multipliers in First Column 

101 230 1 230 1 23 56745674567 1011 8 9 1011 8 91011 14151213141512131415 12301230123 
5 6 7 4 5 6 7 4 5 6 7i 1011 8 9 1011 8 9 1011i 14151213141512131415 12301230123 56745674567 
9 1011 8 910 11 8 9 1011 14 15 12 1314 15 12 13 14 151 

Multipliers 4+;, k = A;kAtk 

F1gure 20-6(c) The third of five stages involved in step k of LU decomposition. Active matrix elements are shaded. 

374 

Matrix Algorithms ti: Banded Matrix LU Decomposition Chap. 20 

(d) Transmit Multipliers From First Column 

Processors 

O 4 8 12 

1 5 9 13 

2 6 10 14 

3 7 11 15 

now contain same subset 

of multipliers (first column) 

FIgure 20-6( d) The fourth of five stages involved in step k of LU decomposition. Active matrix elements are shaded. 
padded decomposition now assigns mabices of size (ni+1) x (ni+1) to each node. This procedure is illustrated in Fig. 20-7 for a simple example. In Table 20-2 the sequential and concurrent theoretical execution times are listed for the five stages of concurrent LU decomposition. The results are tabulated separately for the best case (m = ni-.[f[) and worst case (m = niN + 1) decompositions described above. 
We may now proceed to derive the net overhead and concurrent efficiencies, fol- lowing the definitions established in Chap. 5. Actually, the true algorithm is somewhat more complicated than the simple model used to derive Table 20-2 and may not be fully represented. With this in mind, we begin by noting that the timings obtained for stages (b) and (d) apply to the Caltech/JPL Mark II hypercubes which incorporate no hardware assist for broadcast operations. The already quoted formulae apply to the communication algorithm given in Code 10-3. This is the algorithm also used in the hypercube timings of the current chapter. It is possible however, to devise a better algorithm for the same hardware environment, which we will refer to as the split pipe. The initiating node in the scheme sends information in both directions (e.g., "up" and "down" in stage (b), with the "up" direction wrapped periodically to the bottom of the window). This enhancement approximately halves the startup time previously given by (-.[f[ -2)tcomm .Example pseu- docode illustrating this method is given in Code 20-1. 

Sec. 20-4 Concurrent LU Algorithm 

375 

(e) Correct Remaining Matrix Elements 

l' New column ofL 

At+;, k+j -;;. At+;, k+j -4+;, k Uk, k+j i, j =I O 

New row of U 

Figure 20-6(e) The fifth of five stages involved in step k of LU decomposition. Active matrix elements are shaded. 

a) InvertAkk 0(1) 0(1) O(I/N) O(I/m-) c) Fonn Multipliers m~ -1 m 0(1/~) O(I/m) IKe) Update Matrix (m~ -1)£1 m2 --0(I/m~) O(I/m~) 

Table 20-1 Utilization of Hypercube for Three Stages of LU Decomposition 

In Table 20-2, we have assumed that the times required for multiplications and additions are both equal to tcalc. Floating-point scalar inversion is assumed to take a time tiny, although this term is irrelevant to the asymptotic form of the concurrent overhead. 
Returning to the simple predictions of Table 20-2, we may derive the theoretical overhead relations: 

f(b)0.5 't LU est case = A -liT + --;;- + --;0 mVlV m m 

1.5 

fLU (worst case) = ~ + ~ 't m .[jil A + -;;- 
m m 

1 + Jj!!:;:?:l m 

1 + 1:1f!:;:Jl m 

(20-12) 

(20-13) 

376 

Matrix Algoritluns ll: Banded Matrix LU Decomposition Chap. 20 

Effect of Imbalanced Decomposition 

Full Bandwidth 25 16 Processors 

0 1 2 3 0 1 2 3 0 1 2 3 rdj ,;;:Cc 
456 74567456 7"~' 9101189101189 1011cc~ 
121314151213141512131415~~ 0 1 2 3 0 1 2 3 0 1 2 3 "~ 4 5 6 7 4 5 6 7 4 5 6c 7 "I 8 9 1011 8 9 1011 8 9 1011 S;; 121314151213141512131415~~ ~ ~ : ~ ~ ~ : ~ ~ ~ : ~ "i 8910118910118 91011"IJ 

13 x 13 window 

Processor O 16 elements 1,2,3,4,8,12 12 elements 
Remaining processors 9 elements Figure 20.7 An example of an imbalanced dec<Xnposition. We show the 13 x 13 win. 
dow corresponding to a matrix of bandwidth 25. We shade the matrix elements which give processors 0, 1, 2, 3 additional work. 
These formulae have ignored the nonasymptotic contribution, but more impor- tantly, they have ignored some issues related to concrete implementation. As an exam- ple, stages (c) and (d) for the case of the Caltech machines are almost completely over- lapped (concurrent). In addition, there enters significant software overhead of the type discussed in Sec. 10-5.1. This effect will alter the effective value of tcomm and call into question the assumption tstal1 = tcomm used in Eq. (10-3) to obtain Table 20-2. As a result of these departures from ideality, we can at best expect Eqs. (20-12) and (20-13) to pro- vide an indicative, qualitative description of the measurements. 
Figure 20-8 presents the measurements on the Caltech/JPL Mark II hypercube in comparison with the above th-eoretical predictions. From this plot we derive the ratio: 
t 
comm =4.7 (20-14) 

't= 

tca/c