The Block Distributed Memory Model.
AbstractWe introduce a computation model for developing and analyzing parallel algorithms on distributed memory machines. The model allows the design of algorithms using a single address space and does not assume any particular interconnection topology. We capture performance by incorporating a cost measure for interprocessor communication induced by remote memory accesses. The cost measure includes parameters reflecting memory latency, communication bandwidth, and spatial locality. Our model allows the initial placement of the input data and pipelined prefetching.We use our model to develop parallel algorithms for various data rearrangement problems, load balancing, sorting, FFT, and matrix multiplication. We show that most of these algorithms achieve optimal or near optimal communication complexity while simultaneously guaranteeing an optimal speed-up in computational complexity. Ongoing experimental work in testing and evaluating these algorithms has thus far shown very promising results.
Introduction
Parallel processing promises to offer a quantum leap in computational power that
is likely to have a substantial impact on various aspects of the computing field, and
that in particular can be exploited to investigate a wide range of what has been called
"grand challenge" problems in science and engineering. It is widely recognized [23]
that an important ingredient for the success of this technology is the emergence of
computational models that can be used for algorithms development and for accurately
predicting the performance of these algorithms on real machines. We take a similar
view as in [24] in that the computation model should be a "bridging model" that links
the two layers of hardware and software. Existing computation models tend to be
biased towards one or the other layer, except for very few exceptions. The Bulk Synchronous
Parallel (BSP) model advocated by Valiant [24] is one of the few exceptions.
In this paper, we introduce a computation model that specifically attempts to be a
bridging model between the shared memory (single address) programming model and
the distributed-memory message passing architectures. Distributed memory systems
configured as a single address space are usually referred to as (scalable) shared memory
multiprocessors. These machines achieve the scalability of distributed memory
architectures and the simple programming style provided by the single address space.
Our model can also be used for predicting performance of data parallel algorithms
running on distributed memory architectures.
Since a computation model should predict performance on real machines, we start
with a discussion on the basis of our measure of communication costs incurred by
accessing remote data. As indicated in [8], the hardware organizations of massively
parallel processors (MPPs) seem to be converging towards a collection of powerful
processors connected by a communication network that can be modeled as a complete
graph on which communication is subject to the restrictions imposed by the
latency and the bandwidth properties of the network. According to this common or-
ganization, the communication between the different processors is handled by point-
to-point messages whose routing times are controlled by parameters related to the
network latency, processor communication bandwidth, overhead in preparing a mes-
sage, and network capacity. Such a model avoids a description of the exact structure
of the network since algorithms exploiting specific features of the network are less
likely to be robust enough to work well on a variety of architectures and to adapt
easily to possible future technological changes. However programming the machine
at the message-passing level imposes a heavy burden on the programmer and makes
algorithms development and evaluation quite complicated. On the other hand, the
data-parallel and the shared-memory programming models are appealing in terms of
their ease of use and in terms of their close relationship to sequential programming.
Both models assume a single address space.
The Block Distributed Memory (BDM) model introduced in the next section captures
the performance of shared memory (single address space) algorithms by incorporating
a cost measure for interprocessor communication caused by remote memory
accesses. The cost is modeled using the latency and the communication bandwidth of
each processor. Since a remote memory access involves the transmission of a packet
that typically contains a number of consecutive words, our model encourages the use
of spatial locality by incorporating a parameter m that represents a cost associated
with accessing up to m consecutive words; this cost will be incurred even if a single
word is needed. Our model allows the initial placement of input data and includes
the memory latency hiding technique of pipelined prefetching. Since we measure the
amount of local computation and the amount of communication separately, we are
able to normalize the communication cost and drop one parameter so as to make
the analysis of the corresponding algorithms simpler. We use our model to develop
parallel algorithms for various data rearrangement problems, load balancing, sort-
ing, the Fast Fourier Transform (FFT) computation, and matrix multiplication. We
show that most of these algorithms achieve optimal or near optimal communication
complexity while simultaneously guaranteeing an optimal speed-up in computational
complexity.
In the next section, we provide the details of our model while Section 3 describes
a collection of algorithms for handling data rearrangements that occur frequently
in shared memory algorithms. The load balancing problem is addressed in Section 4
where a communication efficient algorithm is presented, and Section 5 is devoted to the
presentation of efficient algorithms for sorting, FFT, and matrix multiplication. Most
of the resulting algorithms seem to share a common structure with high-performance
algorithms that have been tested on real machines.
2 The Block Distributed Memory (BDM) Model
Our computation model, the Block Distributed Memory (BDM), will be defined in
terms of four parameters p, - , oe, and m. As we will see later, the parameter oe can
be dropped without loss of generality. The parameter p refers to the number of pro-
cessors; each such processor is viewed as a unit cost random access machine (RAM).
In addition, each processor has an interface unit to the interconnection network that
handles communication among the different processors. Data are communicated between
processors via point-to-point messages; each message consists of a packet that
holds m words from consecutive locations of a local processor memory. Since we are
assuming the shared memory programming model, each request to a remote location
involves the preparation of a request packet, the injection of the packet into
the network, the reception of the packet at the destination processor, and finally the
sending of a packet containing the contents of m consecutive locations, including the
requested value, back to the requesting processor. We will model the cost of handling
the request to a remote location ( read or write) by the formula -
is the maximum latency time it takes for a requesting processor to receive the appropriate
packet, and oe is the rate at which a processor can inject (receive) a word
into (from) the network. Moreover, no processor can send or receive more than one
packet at a time. As a result we note the following two observations. First, if - is
any permutation on p elements, then a remote memory request issued by processor
destined for processor P -(i) can be completed in - +moe time for all processors
simultaneously. Second, k remote access requests issued by k
distinct processors and destined to the same processor will require k(- + moe) time to
be completed; in addition, we do not make any assumption on the relative order in
which these requests will be completed.
Most current interconnection networks for multiprocessors use several hardware
and software techniques for hiding memory latency. In our model, we allow pipelined
prefetching for hiding memory latency. In particular, k prefetch read operations issued
by a processor can be completed in -
The underlying communication model for BDM is consistent with the LogP and
the postal models [8, 13, 5] but with the addition of the parameter m that incorporates
spatial locality. However, our model does not allow low-level handling of message
passing primitives except implicitly through data accesses. In particular, an algorithm
written in our model can specify the initial data placement among the local memories
of the p processors, can use the processor id to refer to specific data items, and can use
synchronization barriers to synchronize the activities of various processors whenever
necessary. Remote data accesses are charged according to the communication model
specified above. As for synchronization barriers, we make the assumption that, on
the BDM model, they are provided as primitive operations. There are two main
reasons for making this assumption. The first is that barriers can be implemented in
hardware efficiently at a relatively small cost. The second is that we can make the
latency parameter - large enough to account for synchronization costs. The resulting
communication costs will be on the conservative side but that should not affect the
overall structure of the resulting algorithms.
The complexity of a parallel algorithm on the BDM model will be evaluated in
terms of two measures: the computation time T comp , and the communication time
Tcomm . The measure T comp refers to the maximum of the local computation performed
on any processor as measured on the standard sequential RAM model. The
communication time Tcomm refers to the total amount of communication time spent
by the overall algorithm in accessing remote data. Our main goal is the design of parallel
algorithms that achieve optimal or near-optimal computational speedups, that
is,
is the sequential complexity of the problem under
consideration, in such a way that the total communication time Tcomm is minimized.
Since Tcomm is treated separately from T comp , we can normalize this measure by
dividing it by oe. The underlying communication model for BDM can now be viewed
as the postal model [5] but with the added parameter m reflecting spatial locality.
Hence an access operation to a remote location takes - +m time, and k prefetch read
operations can be executed in - time. Note that the parameter - should now
be viewed as an upper bound on the capacity of the interconnection network, i.e.,
an upper bound on the maximum number of words in transit from or to a processor.
In our estimates of the bounds on the communication time, we make the simplifying
(and reasonable) assumption that - is an integral multiple of m.
We believe that locality is an important factor that has to be taken into consideration
when designing parallel algorithms for large scale multiprocessors. We have
incorporated the parameter m into our model to emphasize the importance of spatial
locality. The notion of processor locality also seems to be important in current multiprocessor
architectures; these architectures tend to be hierarchical, and hence the
latency - is much higher for accessing processors that are further up in the hierarchy
than those that are "close by". This feature can be incorporated into our model by
modifying the value of - to reflect the cost associated with the level of hierarchy that
needs to be used for a remote memory access. This can be done in a similar fashion
as in the memory hierarchy model studied in [2] for sequential processors. However
in this paper we have opted for simplicity and decided not to include the processor
locality into consideration.
Several models that have been discussed in the literature, other than the LogP
and the postal models referred to earlier, are related to our BDM model. However
there are significant differences between our model and each of these models. For
example, both the Asynchronous PRAM [9] and the Block PRAM [1] assume the
presence of a shared memory where intermediate results can be held; in particular,
they both assume that the data is initially stored in this shared memory. This makes
data movement operations considerably simpler than in our model. Another example
is the Direct Connection Machine (DCM) with latency [14] that uses message passing
primitives; in particular, this model does not allow pipelined prefetching as we do in
the BDM model.
3 Basic Algorithms for Data Movements
The design of communication efficient parallel algorithms depends on the existence of
efficient schemes for handling frequently occurring transformations on data layouts.
In this section, we consider data layouts that can be specified by a two-dimensional
array A, say of size q \Theta p, where column i of A contains a subarray stored in the local
memory of processor transformation \Pi on the layout A
will map the elements of A into the layout \Pi(A) not necessarily of the same size. We
present optimal or near optimal algorithms to handle several such transformations
including broadcasting operations, matrix transposition, and data permutation. All
the algorithms described are deterministic except for the algorithm to perform a
general permutation.
We start by addressing several broadcasting operations. The simplest case is to
broadcast a single item to a number of remote locations. Hence the layout A can
be described as a one-dimensional array and we assume that the element A[0] has
to be copied into the remaining entries of A. This can be viewed as a concurrent
read operation from location A[0] executed by processors . The next
lemma provides a simple algorithm to solve this problem; we later use this algorithm
to derive an optimal broadcasting algorithm.
Lemma 3.1 Given a p-processor BDM and an array A[0 resides
in processor P j , the element A[0] can be copied into the remaining entries of A in
time.
Proof: A simple algorithm consists of rounds that can be pipelined. During
the rth round, each processor P j reads A[(j
only A[0] is copied into A[j]. Since these rounds can be realized with
prefetch read operations, the resulting communication complexity is -
We are now ready for the following theorem that essentially establishes the fact
that a k-ary balanced tree broadcasting algorithm is the best possible for
(recall that we earlier made the assumption that - is an integral multiple of m).
Theorem 3.1 Given a p-processor BDM, an item in a processor can be broadcast to
the remaining processors in - 2-d log p
log( -
e communication time. On the other hand,
any broadcasting algorithm that only uses read, write, and synchronization barrier
instructions requires at least (- log p
log( -
+m log p) communication complexity. 4
Proof: We start by describing the algorithm. Let k be an integer to be determined
later. The algorithm can be viewed as a k-ary tree rooted at location A[0];
there are dlog k pe rounds. During the first round, A[0] is broadcast to locations
using the algorithm described in Lemma 3.1, followed by a
synchronization barrier. Then during the second round, each element in locations
is broadcast to a distinct set of locations, and so on. The
communication cost incurred during each round is given by
3.1). Therefore the total communication cost is Tcomm
we set
m)d log p
log( -
log( -
e.
We next establish the lower bound stated in the theorem. Any broadcasting
algorithm using only read, write, and synchronization barrier instructions can be
viewed as operating in phases, where each phase ends with a synchronization barrier
(whenever there are more than a single phase). Suppose there are s phases. The
amount of communication to execute phase i is at least - is the
maximum number of copies read from any processor during phase i. Hence the total
amount of communication required is at least
m). Note that by the
end of phase i, the desired item has reached at most
remote locations. It follows that, if by the end of phase s, the desired item has
reached all the processors, we must have 1). The
communication time
m) is minimized when k
log(k+1)
and the communication time is at least
log(k+1)
. We complete the proof of this theorem by proving the
following claim.
log(k+1)
log( -
+m, for any k - 1.
Proof of the Claim: Let
log(k+1)
log(k+1)
, and
log(k+1)
. Then,
(Case decreasing and f 2 (k) is increasing in this
range, the claims follows easily by noting that f 1
log( -
and
4 All logarithms are to the base 2 unless otherwise stated.
(Case We show that f(k) is increasing when k ? r+1 by showing that
1. Note that since k -
is at least as large as
which is positive for all nonzero integer values of k. Hence f(k) - f( -
and the
claim follows. 2
The sum of p elements on a p-processor BDM can be computed in at most 2-d log p
log( -
e
communication time by using a similar strategy. Based on this observation, it is easy
to show the following theorem.
Theorem 3.2 Given n numbers distributed equally among the p processors of a
BDM, we can compute their sum in O( n
log -
time and at most
2-d log p
log( -
e communication time. The computation time reduces to O( n
Another simple broadcasting operation is when each processor has to broadcast
an item to all the remaining processors. This operation can be executed in minf-
time as shown in the next lemma.
Lemma 3.2 Given a p-processor BDM and an array A[0 distributed one element
per processor, the problem of broadcasting each element of A to all the processors
can be done in minf- communication time.
Proof: The bound of - follows from the simple algorithm described in
Lemma 3.1. If p is significantly larger than m, then we can use the following strategy.
We use the previous algorithm until each processor has m elements. Next, each block
of m elements is broadcast in a circular fashion to the appropriate (d p
processors.
One can verify that the resulting communication complexity is 2-
Our next data movement operation is the matrix transposition that can be defined
as follows. Let q - p and let p divide q evenly without loss of generality. The data
layout described by A is supposed to be rearranged into the layout A 0 so that the first
column of A 0 contains the first q=p consecutive rows of A laid out in row major order
form, the second column of A 0 contains the second set of q=p consecutive rows of A,
and so on. Clearly, if this corresponds to the usual notion of matrix transpose.
An efficient algorithm to perform matrix transposition on the BDM model is
similar to the algorithm reported in [8]. There are rounds that can be fully
pipelined by using prefetch read operations. During the first round, the appropriate
block of q=p elements in the ith column of A is read by processor P (i+1)modp into the
appropriate locations, for During the second round, the appropriate
block of data in column i is read by processor P (i+2)modp , and so on. The resulting total
communication time is given by
pm
the amount of local computation is O(q). Clearly this algorithm is optimal whenever
pm divides q. Hence we have the following lemma.
Lemma 3.3 A q \Theta p matrix transposition can be performed on a p-processor BDM
pm e; this bound reduces to -
We next discuss the broadcasting operation of a block of n elements residing on
a single processor to p processors. We describe two algorithms, the first is suitable
when the number n of elements is relatively small, and the second is more suitable
for large values of n. Both algorithms are based on circular data movement as used
in the matrix transposition algorithm. The details are given in the proof of the next
theorem.
Theorem 3.3 The problem of broadcasting n elements from a processor to p processors
can be completed in at most 2[2-d log p
log( -
using
a k-ary balanced tree algorithm. On the other hand, this problem can be solved in
at most 2[-
pm e] communication time by using the matrix transposition
algorithm.
Proof: For the first algorithm, we use a k-ary tree as in the single item broadcasting
algorithm described in Theorem 3.1, where
+1. Using the matrix transposition
strategy, distribute the n elements to be broadcast among k processors, where each
processor receives a contiguous block of size n
k . We now view the p processors as
partitioned into k groups, where each group includes exactly one of the processors
that contains a block of the items to be broadcast. The procedure is repeated within
each group and so on. A similar reverse process can gradually read all the n items
into each processor. Each forward or backward phase is carried out by using the
cyclic data movement of the matrix transposition algorithm. One can check that the
communication time can be bounded as follows.
km
log p
log( -
pm, we can broadcast the n elements in 2[-
pm
communication
time using the matrix transposition algorithm of Lemma 3.3 twice, once to distribute
the n elements among the p processors where each processor receives a block of size
, and the second time to circulate these blocks to all the processors. 2
The problem of distributing n elements from a single processor can be solved by
using the first half of either of the above two broadcasting algorithms. Hence we have
the following corollary.
Corollary 3.1 The problem of distributing n elements from one processor to
processors such that each processor receives n
p elements can be completed in at most
minf(2-d log p
log( -
pm e)g communication time. 2
We finally address the following general routing problem. Let A be an n
array
of n elements initially stored one column per processor in a p-processor BDM machine.
Each element of A consists of a pair (data,i), where i is the index of the processor to
which the data has to be relocated. We assume that at most ff n
elements have to be
routed to any single processor for some constant ff - 1. We describe in what follows
a randomized algorithm that completes the routing in 2(-
and O( n
computation time, where c is any constant larger than maxf1+ 1
ffg.
The complexity bounds are guaranteed to hold with high probability, that is, with
probability positive constant ffl, as long as
6
the logarithm to the base e.
The overall idea of the algorithm has been used in various randomized routing
algorithms on the mesh. Here we follow more closely the scheme described in [20] for
randomized routing on the mesh with bounded queue size.
Before describing our algorithm, we introduce some terminology. We use an auxiliary
array A 0 of size cn
\Theta p for manipulating the data during the intermediate stages
and for holding the final output, where c ?
ffg. Each column of
A 0 will be held in a processor. The array A 0 can be divided into p equal size slices,
each slice consisting of cn
consecutive rows of A 0 . Hence a slice contains a set of cn
consecutive elements from each column and such a set is referred to as a slot. We are
ready to describe our algorithm.
Algorithm Randomized Routing
Input: An input array A[0 : n
such that each element of A consists of
a pair (data,i), where i is the processor index to which the data has to be routed. No
processor is the destination of more than ff n
p elements for some constant ff.
Output: An output array A
holding the routed data, where c is
any constant larger than maxf1
ffg.
begin
Each processor P j distributes randomly its n
elements into the p slots of the
jth column of A 0 .
so that the jth slice will be stored in the jth processor, for
Each processor P j distributes locally its - cn
elements such that every element
of the form (*,i) resides in slot i, for
Perform a matrix transposition on A 0 (hence the jth slice of the layout generated
at the end of Step 3 now resides in P j ).
The next two facts will allow us to derive the complexity bounds for our randomized
routing algorithm. For the analysis, we assume that
Lemma 3.4 At the completion of Step 1, the number of elements in each slot is no
more than cn
with high probability, for any c ?
pProof: The procedure performed by each processor is similar to the experiment of
throwing n
bins. Hence the probability that exactly cn
balls are placed
in any particular bin is given by the binomial distribution
, and
. Using the following Chernoff bound for estimating
the tail of the binomial distribution
we obtain that the probability that a particular bin has more than cn
balls is upper
bounded by
Therefore the probability that any of the bins has more than cn
balls is bounded by
and the lemma follows. 2
Lemma 3.5 At the completion of Step 3, the number of elements in any processor
which are destined to the same processor is at most cn
with high probability, for any
ffProof: The probability that an element is assigned to the jth slice by the end of
Step 1 is 1
. Hence the probability that cn
elements destined for a single processor
fall in the jth slice is bounded by b( cn
no processor is the destination of
more than ffn
elements. Since there are p slices, the probability that more than cn
elements in any processor are destined for the same processor is bounded by
pe \Gamma( c
and hence the lemma follows. 2
?From the previous two lemmas, it is easy to show the following theorem.
Theorem 3.4 The routing of n elements stored initially in an n
array A of a
p-processor BDM such that at most ff n
p elements are destined to the same processor
can be completed with high probability in 2(-
computation time, where c is any constant larger than maxf1
ffg, and
Remark: Since we are assuming that
6 the effect of the parameter m is
dominated by the bound c n
(as n
Balancing
Balancing load among processors is very important since poor balance of load generally
causes poor processor utilization [19]. The load balancing problem is also important
in developing fast solutions for basic combinatorial problems such as sorting,
selection, list ranking, and graph problems [12, 21].
This problem can be defined as follows. The load in each processor P i is given
by an array A represents the number of useful elements in P i
such that max
n. We are supposed to redistribute the
elements over the p processors such that n
p elements are stored in each processor,
where we have assumed without loss of generality that p divides n.
In this section, we develop a simple and efficient load balancing algorithm for the
BDM model. The corresponding communication time is given by Tcomm - 5- +M
. The overall strategy is described next.
Then, the overall strategy
of the load balancing algorithm can be described as follows: First, the load balancing
problem of the (n 0;1 +n stored in the p arrays A
considered and hence an output array A 0
The array A 0
may have km or
pm (steps 2 and
3). Next, processors with elements from appropriate
processors (Step 4). The details are given in the next algorithm. We assume for
simplicity that all of n, are powers of two.
Algorithm Load Balancing
Each processor P i contains an input array A i [0
The elements are redistributed in such a way that n
p elements are stored
in the output array A 0
begin
Each processor P i reads the held in
the remaining processors. This step can be performed in at most 2-
communication time by Lemma 3.2.
performs the following local computations:
2.1 for
2.2 Compute the prefix sums s
2.3 if (i ! t) then
else
Remark: The index t is chosen in such a way that, for
read n
elements, and for i - t, P i will read n
elements. The indices l i and
r i will be used in the next step to determine the locations of the n
p or n
elements that will be moved to P i . Notice that this step takes O(p) computation
time.
reads appropriate numbers of
elements from A l i
and A r i
respectively.
Remark: This step needs a special attention since there are cases when a set of
consecutive processors read their elements from one processor, say P
that h processors, P have to read some appropriate elements
from P i . Notice that h - d M
n=p\Gammam
1. Then this step can be divided into two
substeps as follows: In the first substep,
read their elements from each such this substep can be done in -
communication time by applying Corollary 3.1. In the second substep, the rest
of the routing is performed by using a sequence of read prefetch operations since
the remaining elements in each processor are accessed only by a single processor.
Hence the total communication time required by this step is (- +M)+(-+ n
.
reads the remaining m elements from the appropriate
processors; the corresponding indices l 0
i can be computed locally as in
Step 2.
Remark: This step can be completed in (-
each processor reads its m elements from at most m processors, and these reads
can be prefetched.
Thus, one can show the following theorem.
Theorem 4.1 The load balancing of n elements over p processors, such that at most
elements reside in any processor, can be realized in 5-
communication
time . 2
5 Sorting, FFT, and Matrix Multiplication
In this section, we consider the three basic computational problems of sorting, FFT,
and matrix multiplication, and present communication efficient algorithms to solve
these problems on the BDM model. The basic strategies used are well-known but the
implementations on our model require a careful attention to several technical details.
5.1 Sorting
We first consider the sorting problem on the BDM model. Three strategies seem to
perform best on our model; these are (1) column sort [15], (2) sample sort (see e.g.
[6] and related references), and (3) rotate sort [18]. It turns out that the column sort
algorithm is best when n - and that the sample sort and the rotate sort
are better when
The column sort algorithm is particularly useful if n - can be implemented
in at most (4-
time. When , the column sort algorithm is not practical since its constant
term grows exponentially as n decreases. The sample sort algorithm is provably
efficient when
6 can be implemented in (3- +(p \Gamma 1)d 5
time and O( n log n
probability. The rotate sort
algorithm can be implemented in 8(-
computation time, whenever n - 6p 2 .
We begin our description with the column sort algorithm.
Column Sort
The column sort algorithm is a generalization of odd-even mergesort and can be
described as a series of elementary matrix operations. Let A be an q \Theta p matrix of
elements Initially, each entry of the
matrix is one of the n elements to be sorted. After the completion of the algorithm,
A will be sorted in column major order form. The column sort algorithm has eight
steps. In steps 1, 3, 5, and 7, the elements within each column are sorted. In steps
2 and 4, the entries of the matrix are permuted. Each of the permutations is similar
to matrix transposition of Lemma 3.3. Since
in this case, these two steps
can be done in 2(-
communication time. Each
of steps 6 and 8 consists of a q-shift operation which can be clearly done in -
communication time. Hence the column sort algorithm can be implemented on our
model within (4-
Thus, we have the following theorem.
Theorem 5.1 Given n - elements such that n
p elements are stored in each
of the local memories of a set of p processors, the n elements can be sorted in column
major order form in at most (4-
computation time. 2
The second sorting algorithm that we consider in this section is the sample sort
algorithm which is a randomized algorithm whose running time does not depend on
the input distribution of keys but only depends on the output of a random number
generator. We describe a version of the sample sort algorithm that sorts on the BDM
model in at most (3-
em
computation time whenever
. The complexity bounds are guaranteed with
high probability if we use the randomized routing algorithm described in Section 3.
The overall idea of the algorithm has been used in various sample sort algorithms.
Our algorithm described below follows more closely the scheme described in [6] for
sorting on the connection machine CM-2; however the first three steps are different.
Algorithm Sample Sort
Input: n elements distributed evenly over a p-processor BDM such that
Output: The n elements sorted in column major order.
begin
Each processor P i randomly picks a list of 5 ln n elements from its local memory.
Each processor P i reads all the samples from all the other processors; hence
each processor will have 5p ln n samples after the execution of this step.
Each processor P i sorts the list of 5p ln n samples and pick (5
samples as the
Each processor P i partitions its n
such that the elements in set S i;j belong to the interval between jth pivot and
(j 1)st pivot, where 0th pivot is \Gamma1, pth pivot is +1, and
Each processor P i reads all the elements in the p sets, S 0;i
using Algorithm Randomized Routing.
Each processor P i sorts the elements in its local memory.
The following lemma can be immediately deduced from the results of [6].
Lemma 5.1 For any ff ? 0, the probability that any processor contains more than
elements after Step 5 is at most
ne
Next, we show the following theorem.
Theorem 5.2 Algorithm Sample Sort can be implemented on the p-processor BDM
in O( n log n
and in at most (3-
em
time with high probability, if
6
Proof: Step 2 can be done in -
em communication time by using a
technique similar to that used to prove Lemma 3.2. By Lemma 5.1, the total number
of elements that each processor reads at Step 5 is at most 2 n
elements with high
probability. Hence, Step 5 can be implemented in
communication time with
high probability using Theorem 3.4. The computation time for all the steps is clearly
O( n log n
6 and the theorem follows. 2
Rotate Sort
The rotate sort algorithm [18] sorts elements on a mesh by alternately applying transformations
on the rows and columns. The algorithm runs in a constant number of
row-transformation and column-transformation phases (16 phases). We assume here
that n - 6p 2 .
Naive implementation of the original algorithm on our model requires 14 simple
permutations similar to matrix transpositions, and 14 local sortings within each pro-
cessor. We slightly modify the algorithm so that the algorithm can be implemented
on our model with 8 simple permutations and at most 14 local sortings within each
processor. Since each such simple permutation can be performed on our model in
communication time, this algorithm can be implemented in 8(-
communication time and O( n log n
computation time on the BDM model.
For simplicity, we assume that n
2t. The results can
be generalized to other values of n and p. A slice is a subarray of size p p \Theta p, consisting
of all rows i such that l
block is a subarray
of size p p \Theta p p, consisting of all positions (i; j) such that l
and r
We now describe the algorithm briefly; all the details appear in [18]. We begin
by specifying three procedures, which serve as building blocks for the main algo-
rithm. Each procedure consists of a sequence of phases that accomplish a specific
transformation on the array.
Procedure BALANCE: Input array is of size v \Theta w.
(a) Sort all the columns downwards.
(b) Rotate each row i (i mod w) positions to the right.
(c) Sort all the columns downwards.
Procedure UNBLOCK:
(a) Rotate each row i ((i mod p p) p p) positions to the right.
(b) Sort all the columns downwards.
Procedure SHEAR:
(a) Sort all even-numbered columns downwards and all the odd-numbered columns
upwards.
(b) Sort all the rows to the right.
The overall sorting algorithm is the following.
Algorithm ROTATESORT
1. BALANCE the input array of size n
p \Theta p.
2. Sort all the rows to the right.
3. UNBLOCK the array.
4. BALANCE each slice as if it were a p \Theta p p array lying on its side.
5. UNBLOCK the array.
6. "Transpose" the array.
7. SHEAR the array.
8. Sort all the columns downwards.
For the complete correctness proof of the algorithm, see [18]. We can easily prove
that this algorithm can be performed in at most 14 local sorting steps within each
processor. We can also prove that each of the simple permutations can be done in
time in a similar way as in Lemma 3.3. Steps 1, 3, 5,
and 6 can each be done with one simple permutation. Steps 2 and 4 also can each be
done with one simple permutation by overlapping their second permutations with the
first permutations of steps 3 and 5 respectively. Originally, Step 7 is "Repeat SHEAR
three times" which is designed for removing six "dirty rows" that are left after Step
5; hence, this step requires 6 simple permutations on our model. Since we assumed
- 6p, the length of each column is larger than that of each row, and we can reduce
the number of the applications of SHEAR procedure in Step 7 by transposing the
matrix in Step 6. Thus, since the assumption n
- 6p implies that there are at most
two dirty columns after the execution of Step 6, one application of procedure SHEAR
is enough in Step 7 for removing the two dirty columns and we have the following
theorem.
Theorem 5.3 Given n - 6p 2 elements, n
elements in each of the p processors of the
BDM model, the n elements can be sorted in column major order form in 8(-
communication time and O( n log n
Notice that if we need to repeat SHEAR dlog(1
)e times at
Step 7 for removing the dirty columns, and the communication time for Algorithm
ROTATESORT is at most k(-
)e.
Other Sorting Algorithms
When the given elements are integers between 0 and p O(1) , the local sorting needed in
each of the previous algorithms can be done in O( n
applying
radix sort. Hence we have the following corollary.
Corollary 5.1 On a p-processor BDM machine, n integers, each between 0 and p O(1) ,
can be sorted in column major order form in O( n
computation time and in (1) (4-
communication
time
n )e, or (3) (3-
em
communication time with high probability, if
Two other sorting algorithms are worth considering: (1) radix sort (see e.g. [6]
and related references) and (2) approximate median sort [2, 22]. The radix sort can
be performed on our model in O(( b
r
r )Tcomm (n; p) communication
time, where b is the number of bits in the representation of the keys, r is such
that the algorithm examines the keys to be sorted r-bits at a time, and Tcomm (n; p) is
the communication time for routing a general permutation on n elements (and hence
the bounds in the above corollary apply). The approximate median sort, which is
similar to sample sort with no randomization used but with
from each processor after sorting the elements in each processor, can be done on
our model in O( n log n
and in at most 3-
communication time, if n -
5.2 Fast Fourier Transform
We next consider the Fast Fourier Transform (FFT) computation. This algorithm
computes the discrete Fourier transform (DFT) of an n-dimensional complex vector
x, defined by y
cos 2-
\Gamma1, in O(n log n) arithmetic operations.
Our implementation on the BDM model is based on the following well-known fact
(e.g. see [17]). Let the n-dimensional vector x be stored in the n
p \Theta p matrix X in
row-major order form, where p is an arbitrary integer that divides n. Then the DFT
of the vector x is given by
where W n is the submatrix of W n consisting of the first n
rows and the first p columns
(twiddle-factor scaling), and   is elementwise multiplication. Notice that the resulting
output is a p \Theta n
matrix holding the vector y = W n x in column major order form.
Equation (1) can be interpreted as computing DFT( n
p ) on each column of X, followed
by a twiddle-factor scaling, and finally computing DFT(p) on each row of the resulting
matrix.
Let the BDM machine have p processors such that p divides n and n - p 2 . The
initial data layout corresponds to the row major order form of the data , i.e., the
local memory of processor P i will hold x 1. Then the
algorithm suggested by (1) can be performed by the following three stages. The first
stage involves a local computation of a DFT of size n
p in each processor, followed by
the twiddle-factor scaling (elementwise multiplication by W n ). The second stage is a
communication step that involves a matrix transposition as outlined in Lemma 3.3.
Finally, n
local FFTs each of size p are sufficient to complete the overall FFT computation
on n points. Therefore we have the following theorem.
Theorem 5.4 Computing an n-point FFT can be done in O( n log n
divides n, the
communication time reduces to -
Remark: The algorithm described in [8] can also be implemented on our model
within the same complexity bounds. However our algorithm is somewhat simpler to
implement. 2
5.3 Matrix Multiplication
We finally consider the problem of multiplying two n\Thetan matrices A and B. We assume
that
log n
. We partition each of the matrices A and B into p 2
3 submatrices, say
each of A i;j and B i;j is of size
\Theta n
assuming without loss of generality that p 1
3 is an integer that divides n. For
simplicity we view the processors indices are arranged as a cube of size p 1
3 \Theta p 1
3 \Theta p 1
that is, they are given by P i;j;k , where 0 -
We show that product
can be computed in O( n 3
computation time and 6(2-d log p
communication time. The overall strategy is similar to the one used in [1, 10], where
some related experimental results supporting the efficiency of the algorithm appear
in [10]. Before presenting the algorithm, we establish the following lemma.
Lemma 5.2 Given p matrices each of size n \Theta n distributed one matrix per proces-
sor, their sum can be computed in O(n 2 ) computation time and 2(2-d log p
log( -
communication time.
Proof: We partition the p processors into p
k groups such that each group contains
processors, where
1. Using the matrix transposition algorithm, we put
the first set of n=k rows of each matrix in a group into the first processor of that
group, the second set of n=k rows into the second processor, and so on. Then for each
processor, we add the k submatrices locally. At this point, each of the processors in a
group holds an n
k \Theta n portion of the sum matrix corresponding to the initial matrices
stored in these processors. We continue with the same strategy by adding each set of
k submatrices within a group of k processors. However this time the submatrices are
partitioned along the columns resulting in submatrices of size n
\Theta n
. We repeat the
procedure d log p
log( -
e times at which time each processor has an n
portion of
the overall sum matrix. We then collect all the submatrices into a single processor.
The complexity bounds follow as in the proof of Theorem 3.3. 2
Algorithm Matrix Multiplication
Input: Two n \Theta n matrices A and B such that p - n 3
log n . Initially, submatrices A i;j
and B i;j are stored in processor P i;j;0 , for each 0 -
Output: Processor P i;j;0 holds the submatrix C i;j , where
and each C i;j is of size n
\Theta n
begin
blocks A i;j and B j;k that are initially stored in processors
respectively. Each such block will be read concurrently by p 1=3
processors. This step can be performed in 4(2-d log p
communication
time by using the first algorithm described in Theorem 3.3.
Each processor multiplies the two submatrices stored in its local memory. This
step can be done in O( n 3
computation time.
For each 0 -
the sum of the product submatrices in the p 1processors P i;j;k , for
computed and stored in processor P i;j;0 .
This step can be done in O( n 2
computation time and in 2(2-d log p
communication time using the algorithm described in Lemma 5.2.
Therefore we have the following theorem.
Theorem 5.5 Multiplying two n \Theta n matrices can be completed in O( n 3
time and 6(2-d log p
communication time on the p-processor BDM model,
log n . 2
Remark: We could have used the second algorithm described in Theorem 3.3 to
execute Steps 1 and 3 of the matrix multiplication algorithm. The resulting communication
bound would be at most 6[-



--R

On Communication Latency in PRAM Computations
Hierarchical Memory with Block Transfer
APRIL: A Processor Architecture for Multiprocessing
Designing Broadcasting Algorithms in the Postal Model for Message-Passing Systems
Multiple Message Broadcasting in the Postal Model
A Comparison of Sorting Algorithms for the Connection Machine CM-2
Overview of the KSRI Computer System
LogP: Toward a Realistic Model of Parallel Computation
Asynchronous PRAM Algorithms
Scalability of Parallel Algorithms for the Matrix Mul- tiplication
The Cache-Coherence Protocol of the Data Diffusion Machine
Load Balancing and Routing on the Hypercube and Related Networks
Optimal Broadcast and Summation in the LogP Model
A Complexity Theory of Efficient Parallel Algorithms
Tight bounds on the Complexity of Parallel Sorting
The Stanford Dash Multiprocessor
Computational Frameworks for the Fast Fourier Transform
Sorting in Constant Number of Row and Column Phases on a Mesh
A Probabilistic Analysis of a Locality Maintaining Load Balancing Algorithm
Optimal Routing Algorithms for Mesh-connected Processor Arrays
Efficient Algorithms for List Ranking and for Solving Graph Problems on the Hypercube
Parallel Sorting by Regular Sampling
Report of the Purdue Workshop in Grand Challenges in Computer Architecture for the Support of High Performance Computing
A Bridging Model for Parallel Computation
--TR

--CTR
David R. Helman , David A. Bader , Joseph JJ, Parallel algorithms for personalized communication and sorting with an experimental study (extended abstract), Proceedings of the eighth annual ACM symposium on Parallel algorithms and architectures, p.211-222, June 24-26, 1996, Padua, Italy
Assefaw Hadish Gebremedhin , Mohamed Essadi , Isabelle Gurin Lassous , Jens Gustedt , Jan Arne Telle, PRO: a model for the design and analysis of efficient and scalable parallel algorithms, Nordic Journal of Computing, v.13 n.4, p.215-239, December 2006
parallel algorithms;parallel model;personalized communication;matrix multiplication;load balancing;broadcasting;sorting;Fast Fourier Transform