Optimal Broadcasting in Mesh-Connected Architectures ^{∗}

Michael Barnett Department of Computer Sciences The University of Texas at Austin Austin, Texas 78712–1188 mbarnett@cs.utexas.edu

David G. Payne Supercomputer Systems Division Intel Corporation 15201 N.W. Greenbrier Pkwy. Beaverton, Oregon 97006 payne@ssd.intel.com

Robert van de Geijn Department of Computer Sciences The University of Texas at Austin Austin, Texas 78712–1188 rvdg@cs.utexas.edu

December 4, 1991

Abstract

In this paper, we disprove the common assumption that the time for broadcasting in a mesh is at best proportional to the square root of the number of processors, at least in the presence of worm-hole routing. We present an optimal algorithm for broadcasting in mesh-connected distributed-memory architectures with worm-hole routing. By organizing the processing nodes in a logical spanning tree, the algorithm executes in time proportional to the logarithm of the number of nodes without inducing contention in the communication network. We restrict the number of nodes in each dimension of the processor mesh to be a power of two. Our method provides insight into how to avoid and/or reduce network contention on meshes for other communication operations. Experimental results on the Intel Touchstone Delta system are included.

Keywords: distributed-memory, mesh-connected, broadcast, parallel processing, worm-hole routing

1 Introduction

We investigate broadcast algorithms for mesh-connected distributed-memory architectures with worm-hole routing. While broadcast algorithms on hypercubes have been studied for quite some time [6, 7, 10], such algorithms can cause network contention on mesh-connected architectures. We will show that, under reasonable assumptions, there is a hypercube-optimal algorithm that is also optimal for meshes. Furthermore, this algorithm creates no network conflicts and has the same logarithmic time complexity on meshes and hypercubes.

We conducted our experiments on the Intel Touchstone Delta, a distributed-memory, messagepassing multicomputer of the Multiple Instruction Multiple Data (MIMD) class developed jointly by the Defense Advanced Research Projects Agency (DARPA) and the Intel Corporation [8]. It consists of 520 i860-based nodes, interconnected via a communications network having the topology of a two-dimensional rectangular grid. (Scaling is not restricted to a power-of-two increment typical of hypercube topologies.) It has a peak double-precision floating-point system performance of ≈ 32 Gigaflops and an aggregate system memory of ≈ 8 Gigabytes. The interconnection network employs a Mesh Routing Chip (MRC), developed at the California Institute of Technology, at each system node. Each MRC provides five channels, one for its associated node and four for its adjacent neighbors in the two-dimensional mesh. The channels consist of two unidirectional buses: one for data flow into the MRC, one for data flow out of the MRC. The peak interprocessor communications bandwidth is ≈ 30 MBytes/s in each direction. The system supports explicit message-passing, with a latency of ≈ 75 microseconds via worm-hole routing using a packet-based protocol. Interconnect blocking is minimized by interleaving packets associated with distinct messages which need to traverse the same interconnect path.

∗This research was supported in part by Intel Corporation.

The Delta system is owned by the Concurrent Supercomputing Consortium (an association of 13 research institutions, government laboratories and government agencies) and resides at the California Institute of Technology.

This research was motivated by the need for efficient broadcast algorithms for codes being developed as part of an effort to port a subset of algorithms included in LAPACK to distributedmemory MIMD architectures [2, 1]. In particular, the broadcast discussed in this paper has been incorporated in an implementation on the Delta of a proposed communication library, the Basic Linear Algebra Communication Subprograms (BLACS) [3], and is used by parallel implementations of the LU factorization [12] and matrix reduction algorithms [5].

Section 2 discusses linear arrays and the extension to meshes is in Section 3. Each section includes timing results on the Intel Touchstone Delta system. Conclusions and future work are discussed in Section 4.

2 Broadcasting on Linear Arrays

We start by concentrating on the one dimensional case. For the busy reader, we briefly explain the principles in Section 2.1, after which sufficient background will have been gained to skip the theory developed in Section 2.2-2.4. Notation is explained in Section 2.2. The basic algorithm is presented in Section 2.3. We present its theoretical properties in Section 2.4. Results from experiments on the Intel Touchstone Delta are given in Section 2.5.

2.1 In a Nut Shell

On hypercube architectures, when the send is a blocking operation (i.e. the sending process does not return until the message has been received by the destination node), optimal implementations can be derived by embedding a minimum spanning tree in the hypercube. The broadcast progresses as illustrated in Figure 1 for a network of eight nodes. Nodes are numbered 0 through 7 (labeled in binary). The minimum spanning tree is derived by toggling successive bits. Note that the order in which bits are toggled, in Figure 1 least significant bit to most significant bit, is not important for hypercubes.

Imagine the same broadcast on a linear array of p nodes with worm-hole routing, as illustrated in Figure 2. Notice that when the bits are toggled least-significant bit first, network conflicts become more and more pronounced as the broadcast proceeds. During step k, when messages of n bytes are sent, there is a link carrying 2^{k−1} messages of length n. If the bandwidth between nodes equals τ, the time complexity for step k is at least 2^{k−1}n/τ. Hence, the total broadcast of log p steps requires time at least (p − 1)n/τ, if step k + 1 cannot start until step k has completely finished. Even if steps are allowed to pipeline and interleave, during the last step there exists a link that must accomodate 2^{d−1} = p/2 messages of length n. This yields a lower bound for the time of at least (p/2)n/τ. To summarize then, it can be shown that when the message-length n is fixed and the minimum spanning tree is embedded by toggling the least-significant bit first, the broadcast is a O(p) operation. Moreover, this result is independent of the root of the spanning tree.

Next, consider the broadcast on a linear array when the order in which the bits are toggled is reversed, as is illustrated in Figure 3. In this case, all network conflicts are resolved, and for fixed n, the broadcast becomes a O(log p) operation. We prove this formally in the next few sections.

2.2 Notation and Assumptions

There are p processing nodes numbered from zero to p − 1. For simplicity, we assume p = 2^{d} for some d, 0 ≤ d. A network is a set of nodes and the associated communication network for passing messages between the nodes. A sub-network is a contiguous subset of the network: [i . . j] denotes the set of nodes {k | k ∈ N ∧ i ≤ k ≤ j}.

We make the following assumptions about the architecture. They are realistic for the class of architectures with which we are concerned.

1. Consecutively numbered nodes are physical neighbors.

2. The communication network is not ring connected; a message between processor i and processor j requires the (sub)-network [min(i, j) . . max(i, j)].

3. The communication network uses cut-through (worm-hole) routing, i.e., intermediate nodes along a message path do not have to handle the message. The result is that multi-hop messages do not take noticeably longer than single-hop messages.

4. All links have the same fixed bandwidth.

Each processor’s identity can be represented with d bits. bit(i, x) is the value of bit i in the binary representation of the number x. The least significant bit is bit number zero. The value of x ⊕ i is the same as that of x, but with bit i “toggled”:

(A j : 0 ≤ j < d ∧ j ≠ i : bit(j, x ⊕ i) = bit(j, x)) ∧ bit(i, x ⊕ i) ≠ bit(i, x)  (1)

The procedure send(m, y), executed at a processor x, sends the message m to processor y. The send does not return until the message has been received (i.e. communication is blocking). x is the source of the message, y the target. Messages that contain more than one item are indicated by tuples using angled brackets, with the elements separated by commas. The execution of recv(m), by a processor x, blocks until a message matching the structure of m is sent to x. Messages not matching the structure of m cause an error.

2.3 The Algorithm

The basic idea is that each message creates two partitions of the network, after which each partition continues the broadcast independently. No message ever crosses a partition boundary, therefore no contention is ever created in the network. The number of partitions doubles with each step of the algorithm, thus the time complexity is logarithmic in the number of nodes. The binary representation of the nodes’s identity is used to calculate the identity of the recipient of each

[figure: Step 1; 000; 001; 010; 110; 100; 111; 101; 011; Step 2; 000; 001; 010; 110; 100; 111; 101; 011; Step 3; 000; 001; 010; 110; 100; 111; 101; 011]

Figure 1: Broadcasting a message from node 0 on a hypercube (least significant bit to most significant bit).

[figure: Step 1; 000; 001; 010; 011; 100; 101; 110; 111]

[figure: Step 2]

[figure: Step 3]

Figure 2: Broadcasting a message from node 0 on a linear array (least significant bit to most significant bit).

[figure: Step 1; 000; 001; 010; 011; 100; 101; 110; 111]

[figure: Step 2]

[figure: Step 3]

Figure 3: Broadcasting a message from node 0 on a linear array (most significant bit to least significant bit).

bcast (nat source, msg∗ M ) :  (1)
nat i, level  (2)
msg M  (3)
begin  (4)
if source ≠ x  (5)
then recv (< M, level >)  (6)
else level := 0  (7)
for i from d − 1 − level downto 0 by −1 do  (8)
send (< M, d − i >, (x ⊕ i))  (9)
end bcast  (10)

Figure 4: Basic Algorithm

message; it is by accessing the bits from the most significant to the least significant that the algorithm yields no contention.

The procedure bcast in Figure 4 is defined at each node. (Line numbers are shown on the right.) The parameter source identifies the node which initiates the broadcast: at that node M is a pointer to the message to be broadcast, while at all other nodes, it points to the location at which to store the message. The variable x is the identity of the processor executing the procedure. All nodes execute the procedure concurrently; level indicates the level of node x in the minimum spanning tree rooted at the source of the broadcast. To initiate a broadcast of a message M at processor j, the procedure is executed by all nodes as bcast(j, M).

2.4 Properties

This section presents the theoretical properties of the algorithm. First, in Section 2.4.1, we prove that the algorithm does indeed implement a broadcast. In Section 2.4.2, a schedule is shown for the broadcast that is logarithmic in the size of the network, thus establishing its time optimality. Finally, Section 2.4.3 proves that the broadcast produces no contention in the communication network.

2.4.1 Correctness

To show that the algorithm implements a global broadcast, we provide a proof by induction on d, the number of bits in the binary representation of the identities of the nodes. The case when d equals zero is trivially satisfied.

Theorem 1 The algorithm in Figure 4 implements a global broadcast.

Proof:

d = 1: When d is equal to one, there are two nodes. Whichever one is the source of the broadcast will send one message since the left and right bounds of the loop in line 8 are the same. The message will be sent to the other processor in line 9, effecting the broadcast. d > 1: Let s be the processor that initiates the broadcast. Then on the first iteration, when

i = d − 1, s will send a message to a processor y whose most significant bit is the opposite of s’s most significant bit. From Equation 1, all messages sent from y will be to nodes with a most significant bit the same as y’s. Similarly, all messages sent from those nodes will also be to nodes with the same most significant bit. By the same reasoning, all further messages that s sends will be to nodes that have the same most significant bit as s, and those nodes will also only send messages to nodes with the same most significant bit as s. (d > 1 implies s will send at least two messages.) Therefore, the first message partitions the network into two halves, one comprising nodes with the same most significant bit as s, the other half’s most significant bit having the other bit value. Each partition has 2^{d−1} nodes: each partition can be regarded as a distinct network, each processor’s identity consisting of the low-order d − 1 bits from its “real” identity. There are now two networks initiating a broadcast, one from s, the other from y, since the number of bits in their identities is d−1, by the inductive hypothesis, each process in the two partitions will receive the message.

□

2.4.2 Time Optimality

The algorithm’s time complexity, O(d), i.e., O(log p), is optimal. We assume that sending a message is a constant time operation (see Section 2.2); clearly every other operation in the algorithm can be performed in constant time.

Obviously, a send is not a constant time operation in the same sense as an add or multiply, as it depends on the distance between the source and target. But from assumption 3 (Section 2.2), multi-hop messages do not take noticably longer than single hop messages. So although the amount of time, especially on current architectures, is much more than operations within the CPU, it is still some constant amount of time. Even if one assumes that all messages travel the maximum distance possible in the entire network, our algorithm is still logarithmic in the number of nodes.

Theorem 2 The algorithm in Figure 4 is completed in O(d) steps

Proof: We show that there exists a schedule of O(d) steps in which the broadcast is completed, with each step containing only constant time operations. Let {s_{k} | 0 ≤ k < d} be the set of steps. Each step, s_{k}, consists of the sends that are executed by a node with the loop index, i, equal to (d − 1) − k. Each step, s_{k}, contains 2^{k} sends. The proof is again by induction. Step s_{0} contains the one send the source node begins the broadcast with. If step s_{k} contained 2^{k} sends, then in step s_{k+1}, each of the 2^{k} nodes that received these messages did so with the value of level equal to k + 1. Each begins to execute the loop in lines 8 and 9, thus sending a message, and the 2^{k} nodes that sent the 2^{k} messages are executing their next iteration. Thus, there are 2^{k+1} nodes sending messages in step s_{k+1}.

After d steps, there have been ∑_{k=0}^{d−1} 2^{k} sends, which is 2^{d} − 1 = p − 1 sends. The broadcast is complete with p − 1 sends, so d steps is actually a tight bound on the algorithm. □

2.4.3 Contention

We define contention as two messages that both need the same communication link between any two nodes. To show that the algorithm induces no contention, we use another proof by induction.

Theorem 3 The algorithm in Figure 4 induces no contention.

Proof: The proof is by induction on d, the number of bits in the identities of the nodes.

d = 1 With two nodes, there is only one message, which by definition cannot cause contention.

d > 1 As before, the first message partitions the network. No message is sent by another process until the first message is received. No message ever crosses between the partitions, so each partition can be considered separately. By the inductive hypothesis, the ensuing broadcasts in each partition induce no contention. Therefore, the entire broadcast induces no contention.

□

2.5 Experimental Results

In this section, we discuss experimental results from implementation of the algorithm in Figure 4 on the Intel Touchstone Delta system.

The actual code is slightly different than that given in Figure 4: since all nodes, a priori, know the index of the source, they pre-compute the value of level, making it unnecessary to include this information in the messages. Also, since d was at most the number of bits in a machine word, the variable i was implemented as a bit vector mask: the mask was all zeros, except for the bit number i. On each iteration, the destination of the message is the exclusive or of the mask and x. The mask is then shifted one position to the right for the next iteration. The broadcast routine was coded in C, using the blocking message-passing calls provided by the Intel communication environment.

The largest subnetwork of the Delta that forms a linear array consists of 32 nodes. We collected data both for our algorithm, which we will denote M2L for most-to-least and the alternative algorithm, L2M, that embeds the minimum spanning tree by toggling bits least-significant bit first. The collected data for various size linear arrays is reported in Figures 5 and 6. In these graphs, data from M2L and L2M are indicated by dashed and solid lines, respectively.

In the left panel of Figure 5, we report the elapsed time as a function of message length on a linear array of 32 nodes, averaging over all possible sources for the broadcast. Clearly, the algorithm of choice is the one which toggles bits from most-significant to least-significant, indicated by the dashed line.

Some insight can be gained by comparing the two curves in Figure 5 (left). The ratio of the slopes of the two curves in the left panel of Figure 5 indicates the extent to which message contention influences the time complexity of the broadcast for L2M. From the arguments presented in Section 2.1 we would expect that the solid line be steeper than the dashed line by a factor approximately equal to (p − 1)/log_{2}p = 31/5 = 6.2. Empirically we find the solid line to be steeper than the dashed line by a factor approximately equal to 2.75. Hence message contention does increase the time required to carry out a broadcast, but not as significantly as the theory suggests.

The difference between what is observed in practice and the theoretical result can be explained by noting that the time required to complete the broadcast for the two methods is a function of the message size, link bandwidth, bandwidth from node to link (the rate with which a node can send message), number of contending messages, and the way the network resolves contention. In particular, for the Delta, the bandwidth with which a node can send a message (using the standard communication primitives), the limiting factor for M2L, is approximately half the bandwidth of a link, the limiting factor for L2M. Furthermore, a limited amount of pipelining between steps of the broadcast can occur for L2M.

The right panel of Figure 5 further illustrates the effects of message contention on the broadcast. The dashed line appears logarithmic in p, as expected. We believe that the solid line would have become more linear if it had been possible to use a larger linear array. Again, the elapsed time is obtained by averaging over all possible sources of the broadcast.

Figure 5: Elapsed time of the minimum spanning tree broadcast on a subnetwork of the Delta system that forms a linear array. Left: execution time as a function of message size when the number of compute nodes p = 32. Right: execution time as a function of the number of compute nodes p when messages of n = 10000 bytes are sent.

The graph in Figure 6 illustrates how the broadcasts are affected by changing the source of the broadcast. While the elapsed time for M2L is, as expected, invariant with respect to the source, time for L2M is very much affected. We believe this can be attributed to the way the communication network of the Delta resolves conflicts.

3 Broadcasting on Mesh Arrays

The same technique can be used for broadcasting on meshes with more than one dimension. Order the dimensions of the mesh arbitrarily: the originator of the broadcast performs first a broadcast along a line of nodes, defined by those nodes whose identities differ from the identity of the originator only in the first dimension. This partitions the network along the first dimension; each partition is a distinct network with one less dimension. The same strategy continues to recursively partition the network; when the sub-networks are one dimensional, the algorithm reduces to that for a linear array.

Extensions to the previous notation are presented in Section 3.1. Section 3.2 presents the algorithm for meshes of any dimension. That the extended algorithm enjoys the same properties as the linear case is shown in Section 3.3 The results of an implementation of the algorithm for the two-dimensional Delta are in Section 3.5.

Under certain conditions, the original algorithm can be used without change. As long the mapping of logical nodes to physical locations is done properly, the broadcast is still optimal and free of contention. Section 3.4 discusses these conditions.

3.1 Notation and Assumptions

The mesh has n dimensions; in each dimension the number of nodes is a power of two. Nodes are identified by a n-tuple with n components. For a n-tuple x, x_{i} is component i. For each component, i, 0 ≤ x_{i} < p_{i}. Here, p_{i} is the number of nodes in dimension i; we restrict p_{i} to be a power of two, Figure 6: Execution time as a function of the source of the broadcast when p = 32 and n = 10000 bytes are sent.

i.e., (A i : 0 ≤ i < n : p_{i} = 2^{d_{i}}) for some d_{i} ≥ 0. Let x, y, and z be n-tuples. A sub-network (a contiguous subset of the mesh) [x . . y] now denotes the set of nodes:

{z | (A i : 0 ≤ i < n : x_{i} ≤ z_{i} ≤ y_{i})}

With each n-tuple x, we can associate a unique natural number. Let d be the maximum of all the d_{i}. Then the number is equal to ∑_{i=0}^{n−1} x_{i} ∗ 2^{i∗d}. Note that the binary representation of this number is the concatenation of those for each component of x, when they are regarded as d bit numbers.

Assumptions made in Section 2.2 are modified as follows:

1. Nodes x and y are physical neighbours exactly when they differ in exactly one component, and the difference is exactly one. (That is, if x and y are looked at as vectors, their difference is a unit vector; the mesh is a rectangular lattice.)

2. The communication network is not torus connected; for n-tuples x and y, let min(x, y) be the component-wise minimum of x and y; similarly for max(x, y). A message between nodes x and y requires the sub-network [min(x, y) . . max(x, y)].

The operation of “toggling” a bit now also specifies the component in which the bit is located: x⊕ < i, j > is the result of toggling bit j in component i. More formally:

(A k : k ≠ i : (x ⊕ < i, j >)_{k} = x_{k}) ∧ (x ⊕ < i, j >)_{i} = (x_{i}) ⊕ j  (2)

3.2 The Algorithm

The algorithm for any dimension mesh is shown in Figure 7. The type id stands for n-tuples of n natural numbers. Considering the numbering of the nodes as row-major order, the mesh is first partitioned into rows, then each row performs the broadcast independently. Note that the comparison in line 3 is between two n-tuples. The destination of the message sent in line 10 is computed by flipping bit i in component j of the node’s identity, x.

bcast (id source, msg∗ M ) :  (1)
nat i, j, dim, level  (2)
msg M  (3)
begin  (4)
if source ≠ x  (5)
then recv (< M, dim, level >)  (6)
else dim := level := 0  (7)
for j from dim to n − 1 do  (8)
for i from d − 1 − level downto 0 by −1 do  (9)
send (< M, j, d − i >, (x ⊕ < j, i >))  (10)
end bcast  (11)

Figure 7: Basic Mesh Algorithm

3.3 Properties

The proofs for correctness and time optimality in Sections 2.4.1-2.4.2 are independent of the underlying topology of the network. They merely depend on the fact that all nodes can be assumed to be connected. Hence these proofs also apply to any dimension mesh.

That the algorithm generates no network contention is a straight forward result of Theorem 3:

Theorem 4 The algorithm in Figure 7 induces no contention.

Proof: Follows immediately from Theorem 3 since the communication in each stage is within independent linear arrays.

□

3.4 Extensions

The key to the success of the algorithm in Figure 4 is that as the algorithm proceeds, the mesh is automatically partitioned into partitions which are rectangular in shape. Further communications are restricted to the interior of each of these partitions. (Messages within rectangular sub-networks do not need to use links outside the sub-network.) Further, for certain numbering schemes which map multi-dimensional node identities to natural numbers, the original algorithm will still work.

This condition holds for the case where the binary representation of the logical node index is a concatenation of the binary representations of the row and column index, as we assumed before (Figure 8 (a)). Since the order of concatenation doesn’t matter, a column-major numbering would be acceptable as well (Figure 8 (b)). Indeed, even a serpentine-like mapping (Figure 8 (c)) works. Not all mappings conform to the condition, as is illustrated in Figure 8 (d). This final mapping embeds a ring so that all neighbors in the ring are physical neighbors in the mesh.

It is important to realize that if the nodes of the mesh are identified by a single index, execution of the algorithm in Figure 4 does not require the nodes to know the physical dimensions of the mesh.

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 (a) row 0 4 8 12 1 5 9 13 2 6 10 14 3 7 11 15 (b) column

0 1 2 3
7 6 5 4
8 9 10 11
15 14 13 12

(c) serpentine

0 1 2 3 15 6 5 4 14 7 8 9 13 12 11 10 (d) ring

Figure 8: Numbering Patterns

Figure 9: Left: Execution time as a function of the grid dimensions when the number of nodes p = 32. Right: Execution time as a function of the number of nodes p when packets of n = 10000 bytes are sent and the grid dimensions are taken to produce minimal execution time.

In a later paper, we will modify the algorithm to include the case where the number of processors in a dimension of the network is not restricted to be a power of two.

3.5 Experimental Results

For the experimental results, the algorithm in Figure 7 was used for a two dimensional mesh with row and column dimensions equal to powers of two. Since the nodes were logically numbed as illustrated in Figure 8(a), no change needed to be made to either M2L or L2M.

The introduction to Section 3 argues that M2L recursively performs linear array broadcasts in the different dimensions of the mesh. A similar argument holds for L2M. This, together with the argument in Section 2.1, would yield an expected execution time for L2M that is O(p_{r} × p_{c}), where p_{r} and p_{c} are the row and column dimension of the mesh, respectively. Since p_{c} = p/p_{r}, the expected optimal mesh dimensions are p_{r} = p_{c} = \sqrt{(}p) if p is an even power of two. If p is an odd power of two, the optimal grid sizes would be those where p_{r} and p_{c} differ by a factor two. Hence, the expected execution time of L2M is O(\sqrt{(}p)). This is illustrated in Figure 9 (left) for fixed p = 32 and various row and column dimensions. By contrast, performance of M2L is not noticeably changed.

In Figure 9 (right), execution time is reported for the optimal choice of grid size when the number of nodes is varied. The execution times appear to be O(\sqrt{(}p)) and O(log p) for L2M and M2L, respectively. Clearly, the method that toggles the most significant bit first is the method of choice.

4 Conclusions

Our analysis shows that under reasonable assumptions, a minimum spanning tree broadcast can be performed as efficiently on mesh-connected architectures with worm-hole routing as on hypercubes. However, the order in which the minimum spanning tree is generated can affect execution time considerably.

Our results can be trivially extended to fan-in operations, where a minimum spanning tree is used to combine results from all nodes to a single node. Since the communication associated with this operation is essentially a broadcast in reverse, in this case the least significant bit should be toggled first. Moreover, similar techniques can be applied to more complicated communications, e.g., total exchange [6, 10] and more complicated combine algorithms [6, 11, 13], in order to reduce network conflicts.

In a future paper, we will generalize our algorithm to allow the number of processors in a dimension of the mesh to be any number, not just a power of two.

Acknowledgments

We acknowledge useful conversations with Greg Regnier and Joe Brandenburg. This research was performed in part using the Intel Touchstone Delta System operated by the California Instutute of Technology on behalf of the Concurrent Supercomputing Consortium. Access to this facility was provided by the California Institute of Technology and Intel Supercomputer Systems Division.

References

[1] E. Anderson, Z. Bai, C. Bischof, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen. LAPACK: A portable linear algebra library for highperformance computers. In Proceedings Supercomputing ’90, pages 2–11. IEEE Computer Society Press, Los Alamitos, California, 1990.

[2] E. Anderson, A. Benzoni, J. Dongarra, S. Moulton, S. Ostrouchov, B. Tourancheau, and R. van de Geijn. LAPACK for distributed memory architectures: Progress report. In Parallel Processing for Scientific Computing, Fifth SIAM Conference. SIAM, 1991.

[3] E. Anderson, A. Benzoni, J. Dongarra, S. Moulton, S. Ostrouchov, B. Tourancheau, and R. van de Geijn. Basic Linear Algebra Communication Subprograms. In Sixth Distributed Memory Computing Conference Proceedings, pages 287–290. IEEE Computer Society Press, 1991.

[4] E. W. Dijkstra and C. S. Scholten. Predicate Calculus and Program Semantics. Texts and Monographs in Computer Science, Springer-Verlag, 1990.

[5] J.J. Dongarra and R.A. van de Geijn. Reduction to condensed form for the eigenvalue problem on distributed memory architectures. LAPACK Working Note 30, Technical Report CS-91-130, University of Tennessee, 1991.

[6] G. Fox, M. Johnson, G. Lyzenga, S. Otto, J. Salmon, and D. Walker. Solving Problems on Concurrent Processors, volume I. Prentice Hall, 1988.

[7] Ching-Tien Ho and S. Lennart Johnsson. Distributed routing algorithms for broadcasting and personalized communication in hypercubes. In Proceedings of the 1986 International Conference on Parallel Processing, pages 640–648. IEEE, 1986.

[8] S.L. Lillevik. The Touchstone 30 Gigaflop DELTA Prototype. In Sixth Distributed Memory Computing Conference Proceedings, pages 671–677. IEEE Computer Society Press, 1991.

[9] C.L. McCreary, M.E. McArdle, and J.D. McCreary. Broadcast communication delay metric for the iPSC/2 and iPSC/860 hypercubes. Technical Report CSE-91-28, Auburn University, 1991.

[10] Y. Saad and M.H. Schultz. Data communciation in parallel architectures. Research Report YALEU/DCS/RR-461, Yale University, 1986.

[11] R.A. van de Geijn. Efficient global combine operations. In Sixth Distributed Memory Computing Conference Proceedings, pages 291–294. IEEE Computer Society Press, 1991.

[12] R.A. van de Geijn. Massively parallel LINPACK benchmark on the Intel Touchstone Delta and iPSC/860 systems. Computer Science report TR-91-28, Univ. of Texas, 1991.

[13] R.A. van de Geijn. On global combine operations. LAPACK Working Note 29, Technical Report CS-91-129, University of Tennessee, 1991.