WT Framework and PRAM

DMM

T_par(P, n) = c_i * O(n * log n / P) + c_v * O(n * log P / P) + c_a * O(P).

The strong point of divide-and-conquer algorithms in general is that by their localizing nature the number of start-ups is low. For the above algorithm, the start-up costs are negligible if c_a * P << c_i * n * log n / P, that is when n * log n >> c_a / c_i * P^2. If we consider how large the typical problem is for which we will apply a parallel sorting algorithm, this condition will generally be satisfied.

The weakness of the parallel quicksort algorithm is that the routing volume is almost of the same order as the internal work: only when c_i * log n >> c_v * log P, the communication cost becomes negligible. This condition will rarely be satisfied, implying that the practical speed-up with parallel quicksort, though increasing almost linearly with P, will be small.

Convex Hull

The complexity of the planar convex-hull problem is closely related to that of sorting. For a set of integer values A = {a_i| 0 <= i < n}, we can determine in linear time the minimum and maximum. Then in linear time the set A can be turned into a set of points S: a_i is mapped to s_i = (a_i, f(a_i)), where f(a_i) is the unique positive value so that s_i lies on the circle through the minimum and maximum. The convex hull of S consists of all elements. Because they must be presented ordered, solving the convex-hull problem can be used to sort. So, T_sort(n) <= T_convex_hull(n) + O(n).

The planar convex-hull problem can also be solved by using sorting as a subroutine and is actually not much harder. The algorithm consists of log n rounds which each take O(log n) time. The total work is bounded by O(n * log n). The algorithm proceeds as follows:

1. Sort the points according to their x-coordinates.
2. Determine the points a and b with the smallest and largest x-coordinate, respectively.
3. Determine all points which lie above the line through a and b and those that lie under this line. Call these subsets S_upper and S_lower.
4. Determine the convex hull of S_upper and the hull of S_lower.
5. Construct the convex hull of S by enumerating the hull of S_upper followed by enumerating the hull of S_lower, enumerating a and b only once.

1. Select the point s so that half of the numbers have smaller and half of the numbers have larger x-coordinate.
2. Construct the sets S_left and S_right of points which lie to the left and right of s, respectively. Add s to both sets.
3. Recursively solve the problem for S_left and S_right.
4. Determine the upper common tangent of the convex hulls of S_left and S_right. For two finite sets of points, an upper common tangent is a line which intersects the sets in exactly one point and which has none of the points of the sets above it. Let s_left and s_right be the points in S_left and S_right, respectively, that lie on the upper common tangent.
5. Construct the convex hull of the whole set by enumerating the points in the hull of S_left until reaching s_left, then continue with the elements in S_right starting from s_right.

Later we will consider how to perform parallel sorting. It is rather simple however to perform it in O(log^2 n) time and O(n * log n) work on an EREW PRAM with n / log n PUs. All other steps of the algorithm except for determining the convex hulls of S_upper and S_lower can be performed in O(log n) time and O(n) work or less. Because the set of numbers has been sorted before, constructing S_left and S_right is trivial. The same is true for constructing the whole convex hull in the last step. For the time being we are happy with using the sequential algorithm for finding the upper common tangent. So, the time and work for the recursive part are given by:

T(n) = O(log n) + T(n / 2),
W(n) = O(n) + 2 * W(n / 2).

Partitioning

Merging is the process of turning two sorted arrays a[] and b[] of length n each into one array c[] of length 2 * n containing all values of a[] and b[] in sorted order. The rank of an element x in a set X, denoted rank(x, X), is defined to be the number of elements y in X with y < x. Here it is convenient to assume that all values are different, otherwise the index of the position in which a value is stored may be taken as a secondary criterion. The rank of an element x in A union B is given by rank(x, A) + rank(x, B). rank(a[i], A) = i and rank(b[i], B) = i, for all 0 <= i < n. Thus, it suffices to determine rank(a[i], B) and rank(b[i], A), for all 0 <= i < n, because once all ranks in A union B are known, it takes constant time and linear work to write the numbers of a[] and b[] in c[].

A trivial and fast way to determine all ranks is to perform binary search: for each i, rank(a[i], B) and rank(b[i], A) can be determined in O(log n) time and O(log n) work. Thus, merging can be performed in O(log n) time and O(n * log n) work. This time is satisfactory, but the superlinear work should be reduced.

The obvious improvement is to not determine the rank for all numbers but only a fraction. The largest fraction that can be used so that the work is still bounded by O(n) is n / log n. So, take all numbers a[j * log n] and b[j * log n], for 0 < j < n / log n, and rank them in the other set. This takes O(log n) time and O(n) work. This initial ranking partitions the problem of merging a[] and b[] into small subproblems: there are 2 * n / log n - 2 values for which we know the positions in both a[] and b[]. These subdivide each of the arrays in subsets of size at most log n which have to be merged pairwise. For each of these subproblems a sequential algorithm requires O(log n) time and work linear in the sum of the sizes of the subsets. So, the partitioning allows to solve the whole problem in O(log n) time and O(n) work.

Pipelining

As a further example we consider the problem of inserting k new numbers b_i, 0 <= i < k, into a 2-3 tree with n nodes with keys, a_i, 0 <= i < n, in sorted order. It is assume that k < n, otherwise it is more efficient to build a new 2-3 tree from scratch. Sequentially k numbers can be inserted in O(k * log n) time, and in general this cannot be done much faster, because on the one hand Omega(k * log k) is needed because as a result of the insertions the numbers get sorted, and if log k = o(log n), then it may happen that after the highest log k levels of the tree all paths are independent and must be followed to find the positions where to insert the elements. These remaining paths have length Omega(log n - log k) = Omega(log n) each. So, in our parallel algorithm we are happy with O(k * log n) work, but we would like to reduce the time as far as possible, that is, to O(log n).

First the smallest and largest number, which can easily be determined with linear work in O(log k) time, are inserted in the normal sequential fashion. This is done to eliminate some special cases: after this operation, we know that any element to insert fits between two consecutive leafs of the tree. The insertion takes O(log n) time and work. Then for all k - 2 other b_i the leaf a_j of the tree is determined after which b_i has to be inserted. A chain is the ordered subset B_j of the b_i which fits between a_j and a_{j + 1}. Let k_j = |B_j|. sum_j k_j = k - 2. First consider the special case k_j <= 1, for all j. In this case, the nodes of the tree may get degree larger than 3, but never more than 6. So, processing the tree bottom-up processing all insertions at the same level in parallel, the whole tree can be processed in O(log n) time with O(k * log n) work. Here it is essential that any node gets split only once and that therefore the property that between any two existing nodes at most one new node has to be inserted also holds at the higher levels of the tree. We also use that any node can be processed in constant time.

The above we call the basic algorithm which will be used as a subroutine for the general case. If k_j > 1, we do not insert all elements in one stroke, but rather first insert only the middle element of the chain. At this point the lists should have been sorted. If we are lucky all lists are short, but there might be one chain with k - 2 elements. In that case we need an optimal sorting algorithm running in logarithmic time. Inserting the middle element of a chain splits the chain in two. So, after at most log k of these operations all chains have length 1 or 0 and we are done. Thus, the general problem can be solved easily in O(log k * log n) time.

The final improvement comes by observing that there is no neeed to wait with the next application of the basic algorithm until the previous one has been completed. These operations can also be overlayed in such a way that the next insertion starts at the lowest level as soon as the previous insertion has finished there. All these waves run up to the root in O(log n) time each. So, the whole algorithm has running time O(log k + log n) = O(log n). The work is unchanged by this overlaying.

Accelerated Cascading

Consider the problem of computing the maximum of n numbers stored in an array a[]. On a Common CRCW PRAM with n * (n - 1) PUs this problem can be solved in O(1) time. The PUs are indexed PU_{i, j}, 0 <= i < n, 0 <= j < n - 1. In addition to the integer array a[] there is a boolean array b[]. The algorithm is simple:

  1. for (all i, 0 <= i < n, using  PU_{i, 0}) do
       b[i] = true;
  2. for (all i, j, 0 <= i < n, 0 <= i < n - 1, using PU_{i, j}) do
       if (a[i] < a[j] || (a[i] == a[j] && i < j))
         b[i] = false;
  3. for (all i, 0 <= i < n, using  PU_{i, 0}) do
       if (b[i])
         max_index = i;

The algorithm shows that the maximum can be computed very fast, but it is not optimal. However, we know optimal algorithms which can be combined with it. We summarize the algorithms:

• The sequential algorithm computes with 1 PU the maximum in O(n) time.
• The hypercube algorithm computes with n PUs the maximum in O(log n) time.
• The Common CRCW algorithm computes with n^2 PUs the maximum in O(1) time.

Suppose we have 2 * n PUs. Then for each subset of 2 elements there are 4 PUs, so the above algorithm can be applied. Hereafter there are still n / 2 numbers. For each subset of 4 PUs there are 16 PUs, so again the above algorithm can be applied. More generally, we have

Lemma: After t >= 0 rounds there are still n_t = n / 2^{2^t - 1} numbers from among which the maximum has to be determined.

Proof: The proof goes by induction. For t = 0, we have n_0 = n / 2^0 = n. So, assume the lemma holds for some t. The remaining n_t numbers are divided in subsets of size 2^{2^t}. For each of these subsets there are 2 * n * 2^{2^t} / (n / 2^{2^t - 1}) = (2^{2^t})^2 PUs available. So, the algorithm can be run again, and the number of candidates is reduced to 2^{2^t - 1} / 2^{2^t} = 2^{2^{t + 1} - 1}. End.

So, with n PUs the maximum can be computed in O(loglog n) time. This is much more efficient than before and still very fast, but not yet optimal. If we take n / loglog n PUs, then first each PU can compute the maximum of 2 * loglog n numbers in O(loglog n) time. Then the fast-but-not-optimal algorithm is used to complete the task.

Theorem: On a Common CRCW PRAM with n / loglog n PUs the maximum of n numbers can be computed in O(loglog n) time.

The given algorithm is not such a spectacular application of accelerated cascading. It can be made more interesting by performing several rounds of the hypercube algorithm before changing to the CRCW algorithm. This saves some PRAM steps, because the reduction with the CRCW algorithm is not particularly efficient when the subsets have size 2 or 4, but it has no impact on the asymptotic running time.

We want to point out that the given algorithm is not obtained by a simple application of the WT-scheduling principle: the work of the algorithm with n PUs is n * loglog n because all PUs are active in all rounds. So, it is not because of the scheduling principle that fewer PUs were taken. The decision rather was inspired as follows: considering the loglog algorithm, we wondered for what n' the work of this algorithm is still bounded by O(n). n' = n / loglog n is the largest such value. Then we wondered whether there is another algorithm which allows to reduce the problem size from n to n' with O(n) work in O(loglog n) time. In this case this other algorithm was the sequential algorithm applied to small subsets.

Symmetry Breaking

As an example we consider the problem of coloring the nodes of a directed cycle of length n. A node coloring of a graph with c colors, is an assignment of numbers c_i 0 <= c_i < c, to the nodes of the graph so that for any two adjacent nodes i and j, c_i != c_j. Coloring a list sequentially is simple: we can start at any node, coloring the nodes with colors 0 and 1 until we reach the last node, which has to be colored with color 2 if the number of nodes is odd.

When considering how to solve the problem in parallel, there is the problem that it does not appear to make sense to start coloring everywhere. In this section we show how to efficiently construct a three-coloring in a deterministic way. Of course randomization allows to solve the problem more easily, but it is good to also know how such problems can be solved deterministically.

The main idea of the parallel algorithm is to start with a correct coloring with many colors and to gradually reduce the number of used colors. A trivial initial coloring is given by c_i = i for all i, 0 <= i < n. The algorithm proceeds in rounds. In any round, all nodes are recolored in parallel as follows: for any node i with color c_i and successor j, the new color is given by c'_i = 2 * dif(c_i, c_j) + bit(dif(c_i, c_j), c_i). Here dif(x, y) gives the index of the least-significant bit which is not the same in x and y and bit(b, x) gives bit b of x, where bit 0 is the least significant one. The computation is simple. dif(,) is the only operation which is not a standard hardware operation, but it might have been so. Much harder operations, such as multiplication, can be performed in one clock cycle, so it is reasonable to assume that even dif(,) can be performed in constant time.

We must check that the new coloring is correct. The proof of this goes by contradiction. So, assume that c'_i = c'_j. Let node k be the successor of j, then we must have dif(c_i, c_j) = dif(c_j, c_k) and bit(dif(c_i, c_j), c_i) = bit(dif(c_j, c_k), c_j). Substituting the first in the second gives bit(dif(c_i, c_j), c_i) = bit(dif(c_i, c_j), c_j), a contradiction with the definition of dif(,), which in particular should give a bit position for which the two numbers are different. Because the previous coloring was correct, such a bit must exist.

If before the recoloring the largest color is c, requiring round_down(log c) + 1 bits, numbered 0, ..., round_down(log c), then afterwards the largest color is at most 2 * round_down(log c) + 1. From this it follows that in O(log* n) rounds the number of colors is reduced from n to 5. To further reduce the number of colors, we apply

  for (c = 5; c > 2; c--)
    for each node i with color c_i == c, recolor i with the 
    smallest color c'_i which does not conflict with the color
    of its predecessor or its successor.

Theorem: A three-coloring of a cycle can be constructed in O(log*^ n) time and O(log^* n * n) work.

A (purely theoretical) weakness of the above algorithm is its non-linear work. The algorithm can be made work-optimal at the expense of increasing the time consumption. The above algorithm is performed only once, reducing the number of colors to O(log n). For the special case that we have a set of n numbers with maximum value bounded by O(log n), the sorting problem can be solved in O(log n) time with O(n) work. Sorting means in this case creating buckets containing the indices of all nodes with a given color. After this sorting, the colors can be handled one-by-one as was done before to reduce the number of colors from 5 to 3. Because there are only O(log n) colors and each node has to be processed only once due to the sorting, this recoloring takes O(log n) time and O(n) work. If the super-fast algorithm is applied twice, the number of colors is reduced to O(loglog n). This has a positive impact on the cost of the final recoloring, but the sorting still takes O(log n) time and therefore this does not lead to a smaller asymptotic running time of the entire algorithm.

Exercises

1. Due to its simple structure parallel quicksort is very suitable for implementation on interconnection networks.
   • Specify missing details and analyze the algorithm for sorting n^2 numbers on an n x n mesh.
   • Specify missing details and analyze the algorithm for sorting n = 2^k numbers on a k-dimensional hypercube.
   • Consider making the hypercube algorithm work-optimal. In other words, consider sorting n * log n numbers on a k-dimensional hypercube. Which communication model do you need: one-port or full-port?

2. For numbers x and n, we want to compute x^i for all 0 <= i < n. Describe how to do this in O(log n) time and O(n) work. This leads to an efficient and fast way to evaluate a polynomial of degree n - 1 given by an array of coefficients a[i]. Write down the complete algorithm.

3. There are arrays a[] and b[] of length. a[i] contains integers and b[i] booleans. If b[i] = true, the interpretation is that a new segment starts for index i. The task is to compute the segmented prefix sums, that is, for any i, we want to compute the value s[i] = sum_{j = i_0}^i a[j], where i_0 <= i is the largest value with b[i_0] = true. Give an algorithm that solves this problem in O(log n) time and O(n) work on an EREW PRAM.

4. Given an array a[] of length n we want to compute the prefix minima, the values m[i], 0 <= i < n, defined by m[i] = min_{j <= i} a[j]. Give an algorithm that solves this problem in O(log n) time and O(n) work on an EREW PRAM.

5. Describe in some detail how to compute the upper common tangent of two upper hulls S_1 and S_2, with all elements of S_1 lying in the plane to the left of those of S_2. The algorithm should run in O(log n) time, where n = |S_1| + |S_2|.

6. Give the smallest numbers n_t, 1 <= t <= 5, for which the recoloring algorithm needs t rounds before reaching a situation with 5 colors.

7. Describe how to efficiently build a 2-3 tree from an array of n elements which stand in the array in sorted order in parallel. Specify the required PRAM submodel and the time and work of the algorithm.

8. Show that the maximum of n numbers can be determined in O(1) time on a PRAM with n^{1 + eps} PUs for any eps > 0. Which PRAM model do you need?

9. We consider the problem of sorting n numbers a[i], with 0 <= i < n, with values in a small range: 0 <= a[i] < m, for all i.
   • Give an algorithm for the case m = log n which runs in O(log n) time and O(n) work on an EREW PRAM. Hint: let S_x, 0 <= x < log n, be the subset of indices i, 0 <= i < n, with a[i] = x. First determine for all x and for all i in S_x, the rank of i within S_i.
   • Generalize the algorithm for arbitrary m and specify the corresponding time consumption and work in terms of n and m. Hint: use radix sort.

Sorting

• Transposition Sort
  • Linear Array
  • DMM
  • Mesh
• Merge Sort
  • Searching
  • Merging through Ranking
  • Bitonic Merge Sort
  • Odd-Even Merge Sort
• Sample Sort
  • Deterministic Sampling
  • Sampling through Randomization
  • Sampling through Shuffling
• Constant-Time Sort
  • PRAM
  • Reconfigurable Array
• Exercises

Transposition Sort

• Linear Array
• DMM
• Mesh

Linear Array

  for (t = 0; t < n / 2; t++)
  {
    for (i = 0; i < n / 2; i++) // Even
      compare_swap(2 * i, 2 * i + 1);
    for (i = 1; i < n / 2; i++) // Odd
      compare_swap(2 * i - 1, 2 * i);
  }

The correctness of the algorithm can be proven by induction. The idea underlying the argument is that the largest number, once it starts moving never gets delayed. It reaches its destination after at most n - 1 steps. The second largest number may be swapped with the largest number if the first step, and then it has to wait one step. Hereafter it will never be delayed anymore, so it arrives no later then after n steps. The operations in each of the inner loops can be executed in parallel. On an interconnection network, the simplest is to implement a compare-swap operation by exchanging the numbers to be compared and to discard one of them in each of the two involved PUs. So, on a linear array, the algorithm can be implemented with n steps in which each PU sends at most one number to a neighbor. Odd-even transposition sort is simple but inefficient: as a sequential algorithm it has complexity Theta(n^2) for sorting n numbers.

DMM

T(P, n) = c_i * O(n / P * log n + n) + c_v * n + c_a * P.

Mesh

1. Determine for all numbers their rank;
2. Route the numbers (and associated information) so that the number with rank r, 0 <= r < n, comes to stand in PU_{r / k}, for k = n / P. (Sometimes the number with rank r should be routed to PU_{r mod P}).

In the DMM model, the structure of the interconnection network is ignored. On interconnection networks, there is a structure, and the indexing scheme has impact on the efficiency with which sorting can be performed. Generally speaking, it is convenient if PUs for which the indices do not differ much also lie close to each other in the network. Nevertheless, sorting with respect to any indexing scheme which is fixed in advance is at least as hard as performing a k-k routing: if the hardest permutation pi routes a packet from PU_i to PU_j, then when sorting the value of the corresponding number in PU_i is set to j and this number will have to be routed from PU_i to PU_j in order to perform the sorting. The indexing scheme may have importance, but it should be noticed that the difference between the best and the worst scheme is limited: first sorting according to the best indexing scheme and then routing the numbers so that they come to stand according to the desired scheme costs at most twice as much as just sorting.

The simplest algorithm for sorting on a mesh is obtained by embedding a linear array into it. This algorithm sorts n^2 numbers on an n x n mesh in n^2 steps. This poor performance is due to the fact that the two-dimensional structure of the mesh is ignored. Ideally, if a problem can be solved in O(n) time on a linear array with n PUs it can also be solved on an n x n mesh in O(n) time. A natural way to obtain mesh algorithms is by combining operations on rows and columns. For sorting, it appears natural to consider algorithms which consist of an alternation of one-dimensional sorting phases, using odd-even transposition sort.

The shearsort algorithm, consists of an alternation of operations in which all columns are sorted in increasing order from top to bottom, and operations in which the even rows, that is the rows with index 0, 2, ..., n - 2, are sorted from left to right while the odd rows, with index 1, 3, ..., n - 1, are sorted from left to right. Only in the final sorting of the rows, the sorting direction can be chosen freely. Performing this operation as all others, the numbers are sorted in snake-like row-major order. When sorting all rows increasingly from left to right we get row-major order. The first alternative gives the following extremely simple algorithm:

  for (t = 0; t <= t_max(n); t++)
  {
    sort all columns from top to bottom;
    sort even rows from left to right and odd rows from right to left;
  }

0-1 Sorting Lemma: If a comparison-based sorting algorithm correctly sorts any distribution consisting of 0-s and 1-s, it correctly sorts any distribution with arbitrary numbers.

Proof: The proof goes by contradiction. Assume there is a set of numbers X = {x_0, ..., x_{n - 1}} for which the algorithm is not correct. Let pi'() be the inverse of the permutation realized by the algorithm. That is, afterwards the numbers have been rearranged to x_{pi'(0)}, ..., x_{pi'(n - 1)}. Let pi() be the inverse of the permutation giving a correct sorting of the numbers. That is, pi() is so that x_{pi(0)} <= x_{pi(1)} <= ... <= x_{pi(n - 1)}. Let k be the smallest index so that x_{pi'(k)} != x_{pi(k)}. All numbers with values smaller than x_{pi(k)} belong to the set {x_{pi(0)}, ..., x_{pi(k - 1)}} = {x_{pi'(0)}, ..., x_{pi'(k - 1)}}, so we must have x_{pi'(k)} > x_{pi(k)}. Because pi() is a permutation, this implies that there also must be an index l > k, so that x_{pi'(l)} = x_{pi(k)}. Now define the set of 0-1 values X' = {x'_0, ..., x_{n - 1}} by x'_i = 0 if x_i <= x_{pi(k)}, x'_i = 1 if x_i > x_{pi(k)}. Because for all 0 <= i, j < n, x_i < x_j => x'_i < x_j, the algorithm works on X' as it did on X (except that possibly some elements which now have the same value do not get exchanged). So, on X', the algorithm will produce the following sequence of values as output: ..., x'_{pi'(k)}, ..., x'_{pi'(l)}, ... = ..., 1, ..., 0, ... . That is, at position k there occurs a 1 while at position l > k there occurs a 0: the sorting is not correct, a contradiction with the assumption that it correctly sorts all 0-1 distributions. End.

  sort all columns from top to bottom;
  for (t = 0; t < t_max(n); t++)
  {
    sort even rows from left to right and odd rows from right to left;
    sort all columns from top to bottom;
  }
  sort even rows from left to right and odd rows from right to left;

Lemma: For any 0-1 number distribution on an n x n mesh, after t, 0 <= t <= log n, passes through the loop in shearsort there are still at most n / 2^t dirty rows.

Proof: The purpose of the first sorting of the columns is that it establishes the property that for all i, 0 <= i < n - 1, the number of 0-s in row is at least as large as in row i + 1. This is true because if row i + 1 would hold more 0-s than row i, there would be a position in row i with a 1 where row i + 1 has a 0, contradicting the sorting of the columns.

For t >= 0 we use induction. For t = 0 the claim holds. So, assume that the claim holds for some t >= 0. Sorting the clean rows does not change anything. So, we may focus on the n_t dirty rows. Let z_i, 0 <= i < n_t, denote the number of 0-s in the i-th of these rows. z_i is monotonously decreasing. Let m_0 = |{0 <= i < n_t| z_i >= n / 2}| and let m_1 = n_t - m_0. We claim that after pass t + 1, at least round_down(m_0 / 2) of these rows only contain 0-s, and at least round_down(m_1 / 2) only 1-s. This is a direct consequence of the the monotonicity of z_i and the fact that even rows are sorted differently from odd rows. Thus, any two consecutive rows with z_i >= n / 2 lead to one row with only 0-s. There are several cases to distinguish. If n_t is odd, then either m_0 is even or m_1 is even and thus round_down(m_0 / 2) + round_down(m_1 / 2) = n_t / 2. If n_t is even, then if m_0 and m_1 are even, round_down(m_0 / 2) + round_down(m_1 / 2) = (m_0 + m_1) / 2 = n_t / 2. Only the case that both m_0 and m_1 are odd deserves special consideration. So consider row m_0, which is sorted from right to left and row m_0 + 1 which is sorted from left to right. If z_{m_0} + z_{m_0 + 1} >= n, then this results in one more row with only 0-s, otherwise in one more row with only 1-s. So, in this case at least round_down(m_0 / 2) + round_down(m_1 / 2) + 1 = n_t / 2 new clean rows are generated. End.

Theorem: Shearsort sorts any distribution of n^2 numbers on an n x n mesh in (2 * log n + 2) * n steps.

The time consumption can be improved by almost a factor two, by realizing that the sorting in the columns does not need to be performed for n steps in the later passes of the loop: if there are only n_t rows which might contain 0-s and 1-s, it suffices to perform n_t steps. Thus, all sorting in the columns can be reduced to n + n + n / 2 + ... + 2 = 3 * n - 2. So, the total sorting time can be reduced to (log n + 4) * n - 2 steps. Even this is not such a good result, but it might serve as basis case of a recursive algorithm. The question is whether it is better than the result obtained by trivially embedding a linear array. The times that should be compared are given in the following table, the results are not very impressive:

Merge Sort

• Searching
• Merging through Ranking
• Bitonic Merge Sort
  • Algorithm
  • Work-Time Framework
  • Interconnection Networks
  • DMM
• Odd-Even Merge Sort
  • Algorithm
  • Work-Time Framework
  • Interconnection Networks
  • DMM

Searching

The algorithm is very simple. Initially x may lie in the whole interval ranging from 0 to n. This interval is reduced in several rounds. In each round the remaining interval is divided in P intervals I_i, 0 <= i < P, of approximately equal size. PU_i, 0 <= i < P, is in charge of I_i. Each of the PUs tests whether x lies in its interval. Assuming that the values are different (otherwise the smallest should be taken) there is a unique interval in which x lies. The PU in charge of this interval writes its index to a global variable c. All PUs read this variable and can start the next round. On an CREW PRAM, each round can be performed in O(1) time. So, the whole searching algorithm runs in O(log_P n).

Theorem: On a CREW PRAM with P PUs, an element x can be ranked in a sorted set of size n in O(log_P n) time.

Merging through Ranking

For two sorted arrays, a[] of length m and b[] of length n, let rank(a, b) denote the array (r_0, ..., r_{m - 1}) of length m with r_i = rank(a[i], b), the rank of a[i] in b[]. Assume that we have a CREW PRAM with n PUs, then if m = O(n^{1 - eps}), for some eps > 0, rank(a, b) can be determined in constant time, because for each of the m elements of a[] we can use n^eps PUs to perform the ranking.

Again we apply the partitioning strategy. For merging two sorted arrays a[] and b[] each of length n with n PUs, we select the subarrays a'[] and b'[] of length sqrt(n) with a'[i] = a[sqrt(n) * i] and b'[i] = b[sqrt(n) * i] (it is convenient to add one more element to each subarray with value infinity). With the available PUs rank(a', b) and rank(b', a) can be computed in constant time. The resulting partitioning reduces the merging problem to merging subarrays of a[] and b[] of length at most sqrt(n). Before we reduced the problem in one stroke to subproblems which could be solved sequentially, now we continue recursively until the size has become constant. Each phase of the recursion takes O(1), so the whole merging takes O(loglog n) time.

Lemma: On a CREW PRAM two sorted arrays of length n each can be merged in O(loglog n) time with linear work.

Proof: Above we have shown how to merge two sorted arrays a[] and b[] of length n on a CREW PRAM with n PUs in O(loglog n) time. The work of this algorithm is O(n * loglog n). This algorithm can be made work-optimal by some kind of partitioning. So, we first determine subarrays a'[] and b'[] each of length n / loglog n by taking a'[i] = a[i * loglog n] and b'[i] = b[i * loglog n]. a'[] and b'[] are merged with n / loglog n PUs in O(loglog n) time. Just as in all previous partitionings, as a result it remains to merge sections of a[] and b[] each of length at most loglog n. With one PU for each subproblem, each of these can be solved with the sequential algorithm in O(loglog n) time. End.

Using the above fast and efficient merging as a subroutine in a merge sort algorithm gives

Theorem: On a CREW PRAM with n / loglog n PUs n numbers can be sorted in O(log n * loglog n) time.

O(log n * loglog n) is not the best possible result: sorting can be solved in O(log n) time and O(n * log n) work on an EREW PRAM. This optimal algorithm, which is known in the literature as Cole's merge-sort algorithm, is considerably more complicated and practically a gap of O(loglog n) is hardly significant. A more serious drawback of the presented algorithm is that it is based on the parallel-search algorithm which really needs the CREW model. The parallel quick-sort algorithm can also be refined to run in O(log n) time with high probability.

Bitonic Merge Sort

• Algorithm
• Work-Time Framework
• Interconnection Networks
• DMM

Algorithm

  void bitonicMergeUp(int[] a, int n) {
    for (int s = n >> 1; s > 0; s = s >> 1)
      for (int i = 0; i < n; i += 2 * s)
        for (int j = i; j < i + s; j++)
          compareSwapUp(a, j, j + s); }

The correctness is not proven here. The invariant that must be proven is that after round r, k >= r >= 0, all numbers in a subset of size 2^r with bit r equal to 0 are smaller than all numbers in the corresponding subset with bit r equal to 1. Because later comparisons are performed only within the subsets, the earlier established properties are preserved. Because all properties together imply a sorted sequence, proving the invariant implies the correctness.

Work-Time Framework

T_merge(2^k) = k / 2 * 2^k, for all k >= 1.

For T_sort(2^k), the time to sort 2^k numbers based on this merge operation, we get

T_sort(2^1) = 1,
T_sort(2^k) = 2 * T_sort(2^{k - 1}) + T_merge(2^k), for all k > 1.

T_sort(2^k) = k * (k + 1) / 4 * 2^k, for all k >= 1.

Theorem: Using bitonic merge sort for sorting n numbers, the work is Theta(n * log^2 n). The time is Theta(log^2 n).

Interconnection Networks

1. PU p sends a copy of a to p', PU p' sends a copy of a' to p.
2. PU p and p' compare a and a' and decide which number to keep and which number to discard.

1. PU p sends a copy of a to p'.
2. PU p' compares a and a'. If it should keep a it sends a copy of a' to p, otherwise it simply discards a and keeps a'.

Bitonic merge sort is very suitable for implementation on a hypercube: in all steps number a[i] is compared with numbers a[j], where j = i +- 2^r, for some r, 0 <= r < k. This means that, when mapping the numbers of a[] in a straightforward way to the 2^k PUs of a k-dimensional hypercube, all communication takes place between neighbors, and that thus all steps can be performed in constant time.

Theorem: Using bitonic merge sort, 2^k numbers can be sorted in Theta(k^2) time on a k-dimensional hypercube.

Surprisingly bitonic merge sort is asymptotically optimal on linear arays and meshes. The analysis is most easy for linear arrays. Merging two subsequences standing each in m consecutive PUs takes m + m / 2 + + ... + 1 = 2 * m - 1 steps. So, for a linear array with n = 2^k PUs all merging rounds together take sum_{0 <= s < k} (2 * 2^s - 1) = 2 * sum_{0 <= s < k} 2^s - k = 2 * n - log n - 2 < 2 * n steps. This is not bad at all! Of course, for linear arrays odd-even transposition sort performs fewer routing steps, but at the expense of far more comparisons.

On a two-dimensional n x n mesh n^2 numbers can be sorted in 2 * n + o(n) steps, which is optimal, but the algorithm is complicated and the lower-order term makes this algorithm unsuitable for moderate values of n. How about bitonic merge sorting? It is essential that the PUs are indexed with the shuffled row-major indexing. In that case the distance between PU i and j with j = i +- 2^r, for some r, 0 <= r < k, is exactly 2^{round_down{r / 2}}. Because such PU i and j lie in the same row or column, the routing can be performed without delay. Thus, merging two subsequences of length 2^s each takes 1 + 1 + 2 + 2 + ... + 2^round_down{(s - 1) / 2} + 2^round_down{s / 2} < 4 * 2^{s / 2}. For the number of steps T(n) for the whole sorting this gives T(n) < 4 * sum_{0 <= s < k} 2^{s / 2} < 8 * sum_{0 <= s < k / 2} 2^s < 8 * n. It is worthwhile to compute T(n) exactly. T(2) = 1 + (1 + 1) = 3. T(4) = T(2) + (2 + 1 + 1) + (2 + 2 + 1 + 1) = 13. T(8) = T(4) + (4 + 2 + 2 + 1 + 1) + (4 + 4 + 2 + 2 + 1 + 1) = 37. More generally, T(n) = T(n / 2) + 3 1/2 * n - 4. The solution is T(n) = 7 * n - 4 * log n - 7. As noticed above, asymptotically this is not competitive, but for small n it might be the best achievable. It is better than transposition sort and shearsort for all n.

DMM

Odd-Even Merge Sort

• Algorithm
• Work-Time Framework
• Interconnection Networks
• DMM

Algorithm

  oddEvenMerge(int a[], int m) {
    if (m > 2) {
      int[] ae = new int[m / 2];
      int[] ao = new int[m / 2];
      for (int i = 0; i < m / 2; i++) {
        ae[i] = a[2 * i];
        ao[i] = a[2 * i + 1]; }
      oddEvenMerge(ae, m / 2);
      oddEvenMerge(ao, m / 2);
      a[0] = ae[0];
      for (int i = 1; i < m / 2; i++)
        if (ae[i] <= ao[i - 1]) {
          a[2 * i - 1] = ae[i];
          a[2 * i] = ao[i - 1]; }
        else {
          a[2 * i - 1] = ao[i - 1];
          a[2 * i] = ae[i]; } 
      a[m - 1] = ao[m / 2 - 1]; }
    else // m == 2
      if (a[0] > a[1]) {
        int x = a[0]; a[0] = a[1]; a[1] = x; } }

We consider the same example as before:

  a[]  initially: 10 12 14 16 20 25 32 34  8 13 24 26 28 36 38 40
  ae[] initially: 10 14 20 32  8 24 28 38
  ao[] initially: 12 16 25 34 13 26 36 40
  ae[] merged:     8 10 14 20 24 28 32 38
  ao[] merged:    12 13 16 25 26 34 36 40
  a[]  finally:    8 10 12 13 14 16 20 24 25 26 28 32 34 36 38 40

The correctness is not proven here, but using the 0-1 lemma this is not hard. The key observation is that the estimate of the position of a number in the global array can be predicted with an accuracy of 1 from the position in ae[] and ao[].

Work-Time Framework

T_merge(2^1) = 1,
T_merge(2^k) = 2 * T_merge(2^{k - 1}) + 2^{k - 1} - 1, for all k > 1.

T_merge(2^k) = (k - 1) / 2 * 2^k + 1, for all k >= 1.

T_sort(2^k) = ((k - 1) * k / 4 + 1) * 2^k - 1, for all k >= 1.

The following table gives the number of comparisons performed by the bitonic and the odd-even merge sort algorithms. First the difference between the algorithms increases, then it gradually decreases again. Due to their O(log^2 n * n) complexity, neither of the algorithms should be used for large values of n. On the other hand, for small n, odd-even merge sort is close to optimal.

Interconnection Networks

DMM

Sample Sort

• Deterministic Sampling
• Sampling through Randomization
  • DMM
  • Mesh
• Sampling through Shuffling
  • DMM
  • Mesh

Deterministic Sampling

1. Each PU sorts its k = n / P numbers.
2. Each PU selects out off its k numbers the x - 1 numbers with ranks i * k / x, for all i, 0 < i < x, as presplitters.
3. Gossip the presplitters, so that all PUs receive all P * (x - 1). presplitters.
4. Each PU sorts the received presplitters using P-way merge sort.
5. Each PU picks from the set of presplitters the elements with ranks i * (x - 1), 0 < i < P, as splitters.
6. Each PU partitions its numbers in P buckets using the splitters. Denote the set of numbers in bucket i of PU_j by B_{i, j}.
7. Perform an all-to-all routing so that PU_j, for all j, 0 <= j < P, receives all numbers in the buckets B_{i, j}, for all i, 0 <= i < P.
8. Each PU sorts all received numbers using P-way merge sort.

We analyze the algorithm for a DMM. The internal time consumption is dominated by the time for the several sorting and bucketing operations. It is given by c_i * O(x * log x + P * x * log P + k * log n). If we take care to select x so that P * x < k, then the time is dominated by the time for sorting k numbers. There are only two communication operations. Good performance can only be achieved if the all-to-all routing is performed in a direct way. Thus, the number of start-ups is log P + P - 1. The routing volume is (1 - 1 / P) * (P * x + k). If P * x < k, then this is dominated by the second term. So, under condition P * x = o(n / P), we get

T(P, n) = (1 + o(1)) * (T_sort(1, n / P) + c_v * n / P + c_a * P).

In the above analysis we silently assumed that the all-to-all routing is balanced and that finally all PUs receive approximately the same number of packets. This is not automatically true. There are two limitations:

• It may be the case that the numbers with values close together stand concentrated in a single PU. For example, if the numbers are standing in reversed order so that all numbers from PU_i have to be routed to PU_{P - 1 - i} for all i, 0 <= i < P.
• It may be the case that that the splitters do not divide the set of n numbers in P subsets each with close to n / P numbers. This would imply an unbalanced routing and unevenly loaded PUs in the final sorting.

The second problem is more specific and requires an analysis. We compute an upper bound on the number of numbers that finally may end up in a single PU. So, the question is how many numbers may lie between any two consecutive splitters. Consider two splitters s = s_{i - 1} s' = s_i. Let R_j and R'_j denote the rank of s and s', respectively, in the set of numbers initially stored in PU_j. Analogously, let r_j and r'_j denote the rank of s and s', respectively, in the set of presplitters selected by PU_j. By definition of the splitters, sum_{j = 0}^{P - 1} (r'_j - r_j) = x. The number S_j of numbers routed to PU_j is given by S_j = sum_{j = 0}^{P - 1} (R'_j - R_j). Using that in each PU exactly k / x numbers lie between any two presplitters, we get the following estimates:

S_j <= k / x * sum_{j = 0}^{P - 1} (r'_j - r_j + 1) = n / P + n / x.
S_j >= k / x * sum_{j = 0}^{P - 1} (r'_j - r_j - 1) = n / P - n / x.

So, there is an "unbalance" of O(n / x). Thus, in order to obtain the above claimed results we must take x = omega(P). Because at the same time we had the condition P * x = o(n / P), this gives the condition P^3 = o(n). This is quite a tough condition which is nevertheless satisfied for most practical values of P and n. By refining the procedure by which the splitters are selected from the presplitters, performing iterated subsampling the condition can be relaxed almost to P^2 = o(n), the typical condition which must be satisfied to obtain good performance on a DMM.

For optimal performance, x should be chosen so that the sum of all additional costs together is minimized. On the one hand there are the costs for handling the presplitters and splitters, on the other hand the costs due to an unbalanced routing and sorting. The first contribution increases with x in the same way the second decreases with x. In such cases an approximation of the optimal choice can be found by equating both cost contributions and solving x. This gives

(c_v + c_i * log P) * P * x = (c_v + c_i * log P) * n / x.

The given algorithm does not guarantee that finally each PU holds the same number of numbers. Furthermore, even though the numbers stand sorted, their ranks are not known. This can easily be overcome: counting the number of numbers in each PU and computing the prefix sum of these numbers, the ranks of the numbers can be determined cheaply. From this it follows where they have to be stored. In the normal case that n / x < n / P, this can be achieved by a routing in which a small fraction of the numbers in PU_i has to be routed to PU_{i - 1} and/or PU_{i + 1}. Of course this correction can be performed along with the final sorting operations.

Sampling through Randomization

• DMM
• Mesh

DMM

1. Each PU creates P buckets and uniformly and independently allocates each of its k = n / P numbers to one of the buckets at random. Denote bucket i of PU_j by B_{i, j}.
2. Perform an all-to-all routing so that PU_i, for all i, 0 <= i < P, receives all numbers in the buckets B_{i, j}, for all j, 0 <= j < P. Let k_i denote the number of numbers received by PU_i.
3. Each PU sorts all its numbers.
4. Each PU creates P buckets. PU_i allocates the number with rank r among its k_i numbers to bucket P * r / k_i. Denote bucket i of PU_j by B_{i, j}.
5. Perform an all-to-all routing so that PU_i, for all i, 0 <= i < P, receives all numbers in the buckets B_{i, j}, for all j, 0 <= j < P. Let k_i denote the number of numbers received by PU_i.
6. Each PU sorts all received numbers using P-way merge sort.
7. Each PU exchanges some packets with the PUs with index one smaller and one larger to correct minor sorting faults and to establish the desired final situation in which each PU holds the same number of numbers.

T(P, n) = (1 + o(1)) * (T_sort(1, n / P) + c_v * 2 * n / P + c_a * 2 * P).

We now analyze the minimal value of n as a function of P to assure that the above result holds with high probability. The following points must be satisfied:

• The all-to-all routings are balanced.
• Each PU has to process about the same number of numbers.
• All but o(n) numbers stand in the correct PU after the second all-to-all routing.

There are several ways to perform the final correction. Consider PU_i and PU_{i + 1}. Both have their numbers sorted, but possibly the values of the largest numbers in PU_i are larger than those of the smallest numbers in PU_{i + 1}. The above estimate has shown that we may assume that the number of these is at most x = c * P * sqrt(k * log n). Here c is a constant which can be determined by working out the details of the estimate. The simplest solution is to let both PUs exchange the x concerned numbers. Then each of them can merge the two subarrays of length x. PU_i keeps the smallest x and PU_{i + 1} the largest x. First exchanging samples, there is no need to make copy and the data exchange can be minimized.

Mesh

For k = 1, a good choice is m = n^{2/3} * log^{1/6} n. Doing that, n' = n / m = n^{1/3} / log^{1/6} n. So, there are P' = n'^2 = n^{2/3} / log^{1/3} n submeshes. k' = 1 * m^2 = n^{4/3} * log^{1/3} n. Any computation within the submeshes, corresponding to steps 1, 3, 4 and 6 of the algorithm above, takes O(m) = O(n^{2/3} * log^{1/6} n) time. The choice of m is made so that P' * (k' * log n)^{1/2} = n^{2/3} / log^{1/3} n * (n^{4/3} * log^{4/3} n)^{1/2} = n^{4/3} * log^{1/3} n = k'. So, in step of the algorithm the rearrangement takes place within subsets of constantly many submeshes. This takes O(m) time as well. Step 2 and 5 amount to routing a randomization and the inverse thereof. It can easily be shown that congestion is very limited. Each routing phase can be performed in n + O(sqrt(n * log n)) time. So, in total the sorting takes 4 * n + O(n^{2/3} * log^{1/6} n) time. This is twice as much as given by the distance lower bound.

For larger k optimality can be established. Now it is better to choose m smaller in dependency of k: the best choice is m = n^{2/3} * log^{1/6} n / k^{1/6}. This means that for very large k we will take m = 1, which sounds natural. The operations in the submeshes take O(k * m) = O(k^{5/6} * n^{2/3} * log^{1/6} n) time. Again the choice is made so that P' * (k' * log n)^{1/2} = n^{2/3} * k^{1/3} / log^{1/3} n * (n^{4/3} * k^{2/3} * log^{4/3} n)^{1/2} = n^{4/3} * k^{2/3} * log^{1/3} n = k'. In this case in step 2 and 5 half of the packets are routed XY and the other half YX. This makes that in each of the routing phases only k * n / 8 packets have to cross the bisection. Thus, together these steps take k * n / 2 + o(k * n) time.

Theorem: On a two-dimensional n x n mesh, k-k sorting can be performed in max{4 * n, k * n / 2} + o(k * n) steps with high probability.

Sampling through Shuffling

• DMM
• Mesh

DMM

Based on shuffles and unshuffles the following is a simple, versatile and efficient sorting algorithm which is first formulated for a DMM with P PUs and k = n / P numbers per PU:

1. Sort all numbers in each PU.
2. Shuffle the numbers.
3. Use P-way merge sort to sort the numbers in each PU.
4. Unshuffle the numbers.
5. Use P-way merge sort to sort the numbers in each PU.
6. Perform data exchange between pairs of PUs with consecutive indices to establish the desired global sorting with k numbers per PU.

The analysis of the algorithm goes analogously to those above. For a number x which after the sorting in step 3 has rank r in PU_i, the PU in which it currently resides, the global rank is estimated as R' = P * r. Let R denote the actual global rank of x. It is easy to put upper bounds on |R' - R|. For this we need some additional notation. Let r_j be the rank of x in PU_i among the numbers originating from PU_j. Among the numbers which originally resided in PU_j, x falls between the number with rank j + (r_j - 1) * P and the one with rank j + r_j * P. Summing over all PUs gives

P * (P - 1) / 2 + P * sum_{j = 0}^{P - 1} (r_j - 1) < R < P * (P - 1) / 2 + P * sum_{j = 0}^{P - 1} r_j

Mesh

The mesh is not square but has height n^2 and width n. Each PU holds one number. The numbers will be sorted in column-major order. The algorithm consists of a few simple one-dimensional operations:

1. Sort all columns from top to bottom.
2. Shuffle all numbers within the rows. That is, the number in PU_{i, j}, 0 <= i < n^2, 0 <= j < n, is routed to PU_{i, (i + j) mod n}.
3. Sort all columns from top to bottom.
4. Rotate the numbers in the columns. That is, the number in PU_{i, j}, 0 <= i < n^2, 0 <= j < n, is routed to PU_{(i + n * j) mod n^2, j}.
5. Unshuffle all numbers within the rows, that is, the number in PU_{i, j}, 0 <= i < n^2, 0 <= j < n, is routed to PU_{i, (i / n - 2 * j) mod n}.
6. Sort all even columns from top to bottom and all odd columns from bottom to top.
7. Sort pairs of adjacent row positions. That is, compare and swap the values in PU_{i, 2 * j} and PU_{i, 2 * j + 1}, 0 <= i < n^2, 0 <= j < n / 2.
8. Sort pairs of adjacent row positions. That is, compare and swap the values in PU_{i, 2 * j - 1} and PU_{i, 2 * j}, 0 < i < n^2, 0 <= j < n / 2.
9. Sort all columns from top to bottom.

Lemma: Starting with a zero-one distribution, there are at most n dirty rows after step 3.

Proof: Let z_j, 0 <= j < n, be the number of zeroes in column j. The shuffling is so that any column either receives round_down{z_j / n} or round_up{z_j / n} of these zeroes. So, a column receives at least l = sum_j round_down{z_i / n} and at most h = sum_j round_up{z_i / n} zeroes. So, after step 3, in row 0, ..., l - 1 there are only zeroes, while in row h, ..., n^2 - 1 there are only ones. Only in row l, ..., h - 1 there may be both zeroes and ones. Because h - l <= n, the claimed result follows. End.

Lemma: Starting with a zero-one distribution, there are at most 2 dirty columns after step 5.

Proof: Step 4 and 5 together are so that the contents of the columns is permuted among the rows so that the smallest n numbers in each column go to column 0, and more generally that all numbers standing in row j * n, ..., (j + 1) * n - 1 are ending up in column j. Let l and h be as in the proof of the previous lemma. There are two cases to distinguish. There is a j so that j * n <= l and h < (j + 1) * n. In that case only column j receives zeroes and ones. So, column j will be the only dirty column. Otherwise, because h - l <= n, there is a j so that j * n <= l and h < (j + 2) * n. In that case column 0, ..., j - 1 receives only zeroes, column j + 2, ..., n - 1 receives only ones, and column j and j + 1 may receive both zeroes and ones. End.

Theorem: The given column-sort algorithm correctly sorts any distribution of n^3 numbers on an n^2 x n mesh.

Proof: Given the situation as described in the proof of the above lemma, with zeroes and ones only in column j and j + 1 for some j, 0 <= j < n - 1, sorting the columns up and down alternatingly followed by sorting pairs of adjacent row positions reduces the number of dirty columns to one. This is a special case of the argument for shearsort that the number of dirty columns is reduced by a factor two in each round. The final sorting rearranges the numbers in this dirty column into sorted order. End.

Column-sort can also be formulated processor independent. In this generalized form it becomes immediately clear how to carry it over to interconnection networks. The numbers are divided over n^{1/3} buckets of size n^{2/3}. The role of the PUs is taken over by the buckets.

1. Sort the numbers in each bucket.
2. Shuffle the numbers.
3. Sort the numbers in each bucket.
4. Unshuffle the numbers.
5. Sort the numbers in each bucket.
6. Perform data exchange between pairs of buckets with consecutive indices to establish the desired global sorting with n^{2/3} numbers per bucket.

Sorting n numbers can be performed by a combination of a constant number of operations in which n^{1/3} buckets with n^{2/3} numbers each are sorted and some regular rearrangements of the elements between these buckets.

Theorem: On a 2-dimensional n x n mesh with k >= 8 numbers per PU, sorting can be performed in k * n / 2 + o(n) steps.

Proof: The buckets can be chosen to consist of the k^{2/3} * n^{4/3} numbers residing in each n' x n' submesh for n' = n^{2/3} / k^{1/6}. Sorting these submeshes can be performed with any sorting algorithm which achieves the optimal time order in O(k * n') = O(k^{5/6} * n^{2/3}). The final correction also has this time order. The shuffle and unshuffle send only half of the packets over the bisection and when exploiting the all-port facility they can be performed in k * n / 4 steps each. Thus, the whole sorting can be performed in k * n / 2 + O(k^{5/6} * n^{2/3}) steps. End.

Constant-Time Sort

• PRAM
• Reconfigurable Array

PRAM

The n! permutations can be enumerated systematically, and permutation p, 0 <= p < n!, is denoted pi_p. The PUs are indexed PU_{p, i}, 0 <= p < n!, 0 <= i < n. The algorithm is simple:

  1. for (all p, 0 <= p < n!, using  PU_{p, 0}) do
       b[p] = true;
  2. for (all p, i, 0 <= p < n!, 0 <= i < n - 1, using PU_{p, i}) do
       if ( a[p[i]]  > a[p[i + 1]] || 
           (a[p[i]] == a[p[i + 1]] && p[i] > p[i + 1]))
         b[i] = false;
  3. for (all p, i, 0 <= p < n!, 0 <= i < n, using PU_{p, i}) do
       if (b[p]) {
         c[p[i]] = a[i];
         a[i] = c[i]; }

Reconfigurable Array

Using O(n^3) hardware to sort n numbers is not particularly efficient. By how much can this be reduced while still bounding the sorting time to O(1)? We will see that O(n^2) hardware is enough. This is optimal because performing a permutation in O(1) requires a bisection bandwidth of n, while the grid also has a basis of length n.

The main conclusion from column-sort is that sorting n numbers can be performed by several times sorting O(n^{2/3}) numbers and several rearrangements. Of course each of these sortings themselves can also be performed using column-sort. So, the idea for sorting n numbers on a reconfigurable n x n mesh is to recursively apply column sort. The recursion is terminated as soon as the size n' of the subproblem to solve satisfies n'^2 < n, because then it can be solved in constant time using the presented algorithm for sorting on reconfigurable arrays. Let n_i, i >= 0, denote the size of the subproblems at level i of the recursion. n_0 = n, n_1 = n^{2/3} and n_2 = n^{4/9}. So, the recursion tree has constant depth. Because the number of subproblems increases only by a constant factor between two levels of the recursion, there are in total only constantly many subproblems of size n^{4/9} to solve. We conclude

Theorem: On a reconfigurable processor array with n^2 PUs, n numbers can be sorted in O(1) time.

Exercises

1. Give a transposition-sort algorithm for a linear array with an odd number n of PUs. How many steps do you need? Prove the correctness formally. Hint: use the 0-1 sorting lemma and induction.

2. Consider sorting n^2 numbers on an n x n mesh. Give an algorithm running in O(n) time. The final arrangement can be freely chosen but must be fixed.

3. Consider sorting n numbers on a k-dimensional full-port hypercube with P = 2^k PUs. Present an algorithm running in O(n / P) time. What is the minimum n for which your algorithm achieves this time order?

4. Generalize the shearsort algorithm for three-dimensional meshes. How long does your algorithm take to sort n^3 numbers on an n x n x n mesh?

5. Consider performing shearsort on the k-dimensional hypercube H_k. Distinguish two interpretations of such a hypercube: it can either be viewed as two interconnected copies of H_{k - 1} or as the direct product H_{k / 2} x H_{k / 2}. Give for both interpretations a rough estimate of the time consumption.

6. Consider sorting n^2 numbers on a DMM with P PUs using a variant of shearsort. Assume that n >= P. The mesh can be mapped on to the DMM in several ways:
   • Each PU taking care of a bundle consisting of n / P rows;
   • Each PU taking care of a bundle consisting of n / P columns;
   • Each PU taking care of an n / sqrt(P) x n / sqrt(P) submesh.
   For each of these views formulate the costs in the DMM model.

7. In the shearsort algorithm, the only purpose of the sorting in two-directions in the rows is to spread the values out. This can much more effectively be achieved by randomization. Analyze a variant of shearsort in which the sorting in the rows is replaced by performing a random permutation. A random permutation can efficiently be generated by choosing random numbers and sorting on them. The analysis should be based on the 0-1 sorting lemma. You may use that when n times picking white and black balls which are white uniformly with probability 1 / 2, that then the probability that the number of selected white balls exceeds n / 2 + h is bounded by e^{-h^2 / (c * n)}, for some small constant c.

8. The sorting algorithm based on first randomizing all numbers among the PUs, which was analyzed above for a mesh with k >= 8 numbers per PU, can also be used for the case that each PU holds exactly one number. Specify the time consumption of the algorithm for an n x n mesh.

9. In the text above the performance of bitonic merge sort was analyzed for linear arrays and two-dimensional meshes. Perform an analogous analysis for 3-dimensional n x n x n meshes for some n = 2^k. Most interesting is the constant of the term which is linear in n. How does this constant develop when the algorithm is applied for d-dimensional n x ... x n meshes?

10. Rewrite the method oddEvenMerge() in the given implementation of odd-even merge sort so that it becomes insitu and iterative. Hint: first study the above picture showing how 16 numbers are sorted. Compare the new version with the original one and with the implementation of bitonic merge sort, for a suitable range of problem sizes.

11. The presented three "sample-sort" algorithms (one with explicit sample, one randomizing the numbers and one shuffling them) can also be used for sorting on a hypercube with n PUs and one number per PU.
    • Formulate recurrency relations coarsely estimating the sorting time of each algorithm. It may be assumed that permutation routing can be performed in O(n) time.
    • Solve the recurrence relations and conclude that the differences, which before were hardly significant because they only affected the constant factor, now are more important.

12. Consider sorting n numbers on a reconfigurable array with n^2 PUs. How many steps does it take exactly when performing column-sort recursively? Assume that n is so large that c * n^{8/9} <= n for any constant c you may need and you may ignore problems of divisibility.

13. Assume that reconfigurable switches are relatively expensive, much more expensive than providing the same bisection bandwidth by conventional means. Then it may make sense to build a hybrid system: consisting of a reconfigurable grid which is used for computational purposes and some other network which is used for massive communication. We want to use this network for sorting n numbers on a system with n PUs. The time for the rearrangements is ignored. How big, in terms of n, must the reconfigurable array be in order to assure O(1) sorting time?

List Ranking

• Introduction
• Pointer Jumping
• Independent-Set Removal
• Sparse Ruling Sets
• One-by-One Cleaning
  • Two-Processor Algorithm
  • DMM
  • Circular Array
• Exercises

Two-Processor Algorithm

Introduction

We consider a set of lists with n nodes in total. This is internally represented by an array succ[] of length n. For a node u, succ[u] gives the successor of u. So, the lists are directed towards their final nodes. If u is the last node of its list, we have succ[u] = u. The task is to compute the values of the arrays last[] and dist[]. For a node u, last[u] indicates the last node of the list to which u belongs, and dist[u] gives the number of links that must be traversed in order to reach last[u] starting in u.

Sequentially list ranking can be solved by first figuring out the initial elements of the lists and then simply traversing the list while numbering the nodes.

  void sequential(int succ[], int[] frst, int dist[], int n)
  {
    boolean[] first = new boolean[n];
    for(int i = 0; i < n; i++)
      first[i] = true;
    for(int i = 0; i < n; i++)
      first[succ[i]] = false;
    for(int i = 0; i < n; i++)
      if (first[i]) // i is first node
      {
        int j = last[i] = i;
        int d = dist[i] = 0;
        while (succ[j] != j)
        {
          j = succ[j];
          d++;
          last[j] = i;
          dist[j] = d;
        }
      }
  }

There is no obvious parallel algorithm. Parallel list ranking is much harder than parallel sorting: when sorting any two elements can be compared and this adds some information, but while list ranking, the successor of a node is typically stored in another processor and just going there is expensive. Theoretically the parallel list-ranking problem has been solved: for ranking a list of length n, there are PRAM algorithms using n / log n PUs running in O(log n) time. But, practically it is hard to achieve reasonable speed-ups. Even for very large problems and a very good interconnection network it is not easy to achieve a speed-up of more than 30 on a parallel computer with 100 PUs. For sorting it would have been easy to achieve a speed-up of 80 or 90. This is mainly due to the hard, very non-local nature of list ranking, but also to the fact that its sequential complexity is only linear, while sequential sorting has running time O(n * log n).

The ideas in this section go back on many authors. Just to name a few: Wylie, Anderson, Cole, Miller, Vishkin, Reid-Miller, Ranade and Sibeyn.

Pointer Jumping

  void pointerJumping(int[] succ, int[] last, int[] dist, int n) 
  {
    // Initialization
    for (int i = 0; i < n; i++)
      if (succ[i] == i)
      {
        last[i] = succ[i] - n;
        dist[i] = 0;
      }
      else
      {
        last[i] = succ[i];
        dist[i] = 1;
      }

    // Performing rounds until all reached final node
    boolean flag = true;
    while (flag)
    {
      flag = false;
      for (int i = 0; i < n; i++)
        if (last[i] >= 0) // i has not yet reached last node
        {
          flag = true;
          dist[i] += dist[last[i]];
          last[i]  = last[last[i]];
        }
    }

    // Correction of last[] values
    for (int i = 0; i < n; i++)
      last[i] += n;
  }

Operations of the kind last[i] = last[last[i]], for all i, are called pointer-jumping steps. In a parallel algorithm, these are performed synchronizedly. So, as long as last[u] does not yet point to the last node of a list, the value dist[u] doubles in every pass. Therefore the while loop will be executed at most log n times. The work in each pass is O(n), so the total work is O(n * log n) and the time consumption is O(log n). Pointer jumping is simple and fast, but inefficient due to its super-linear work.

On an interconnection network, each pass requires two all-to-all routing operations, in which each question is resulting in two answers. So, in the DMM model this algorithm costs

T(P, n) = c_i * O(n / P * log n) + c_v * 3 * log n * n / P + c_a * 2 * log n * P.

Independent-Set Removal

The independent-set removal algorithm consists of a reduction and an expansion phase. In the reduction phase, repeatedly the nodes in an independent set are excluded. This continues until the problem has been reduced sufficiently to be solved in a more basic way. Then all excluded elements are reinserted in reversed order. For efficiently performing these operations it is practical to have a doubly linked list, that is, for every node i there is also data item pred[i] giving its predecessor next to succ[i] giving its successor.

During the reduction phase, the following steps are iterated: somehow an independent set S is constructed, then each node i in S, is excluded from the list by performing

  succ[pred[i]]  = succ[i]; 
  dist[pred[i]] += dist[i];
  pred[succ[i]]  = pred[i];

  last[i]  = last[succ[i]];
  dist[i] += dist[succ[i]];

In the work-time framework, each of the exclusion and inclusion operations can be performed in O(1) time and work that is proportional to the number of participating nodes. On an interconnection network each of these operations require a constant number of all-to-all routings, with a routing volume that is proportional to the number of participating nodes.

The main point that remains to be specified is the selection of the independent set. This can be done both randomly and deterministically. Investing some more effort in order to get a larger independent set might be a good idea, because this will result in a faster reduction of the problem size. However, the simplest idea is to perform random coin tossing: for each node 0 or 1 is chosen independently with probability 1/2. S is the set of all nodes with value 1 whose successor has value 0. For a list of length n, the expected number of nodes selected equals (n - 1) / 4. Using the Azuma inequalities it can even be shown that at least n / 4 - o(n) nodes are selected with high probability. This means that applying this selection method, the problem size is approximately multiplied by a factor 3 / 4 in each reduction round. In the work-time framework, S can be selected in constant time, on an interconnection network, it requires one all-to-all routing to inform all predecessors.

So, over all selection, reduction and expansion rounds, the sum of the the performed work is given by c * (1 + 3/4 + 9/16 + ...) * n <= 4 * c * n for some suitable constant c. So, it is linear in n. The time is proportional to the number of performed reduction rounds because for each reduction there is one corresponding selection and expansion. There is no need to continue reducing until the problem size is constant, though this would not be wrong. If the selection of the independent set is performed randomly, it is better to switch to another algorithm earlier. A reasonable choice is to reduce the problem size by a factor log n using independent-set removal in O(loglog n) time, and then to rank the remaining nodes using pointer jumping in O(log n) time. Both algorithms perform O(n) work.

Theorem: On an EREW PRAM with n / log n PUs the list-ranking problem for a set of lists with n nodes can be solved in O(log n) time.

Proof: The main problem is the allocation of work to the PUs. If initially each PU takes care of log n nodes, it is not guaranteed that this number decreases in a smooth way. However, in groups of log n PU, which initially take care of log^2 n nodes, such a smooth decrease is assured with high probability. These subsets of the PUs can be perceived as constituting (loglog n)-dimensional hypercubes. After reduction round r, each PU performs a load balancing with the PU at the other end of the (i mod loglog n)-axis of this conceptual hypercube. Doing this, it becomes highly unlikely that in round r a PU has to deal with substantially more than (3/4)^r * log n nodes. Neither independent-set removal nor pointer jumping require concurrent-memory access. End.

On an interconnection network with P PUs, we can proceed in the same way although in practice it is better to first reduce the problem size slightly more in order to reduce the cost of the pointer jumping. Alternatively, the problem size can first be reduced by a factor P and then all remaining nodes can be concentrated in PU 0, where the subproblem can be solved internally in O(n / P) time. This is not so elegant, but minimizes the number of start-ups. In terms of the DMM the second approach gives

T(P, n) = c_i * O(n / P) + c_v * O(n / P) + c_a * O(log P * P).

As observed in the section on "Routing" in the chapter on "Interconnection Networks", the number of start-ups for an all-to-all routing can be strongly reduced by simulating the routing of a d-dimensional meshes or hypercubes. In general this observation is useful only for problems that are too small for the used number of PUs, because such an indirect routing multiplies the routing volume and causes overhead for rearranging packets as well. So, this may only give an improvement if the time consumption is dominated by the time for the start-ups. However, in our case, the problem size is gradually reduced, and even when the overall time-consumption is dominated by the term (c_i + c_v) * O(n / P), it may happen that in some of the all-to-all routings the start-up time dominates. So, the best performance is obtained, when at any stage the best all-to-all routing subroutine is chosen as a function of P and the size of the packets. Working out this idea more carefully, shows that the factor log P form the term c_a * O(log P * P) can be omited without increasing the leading constants hidden in the other two terms.

The independent set can also be found deterministically. One possibility to obtain a large independent set is to deterministically construct a three-coloring as discussed in the section on "Symmetry Breaking" of the chapter on "Basic Techniques". The nodes with any of the three colors constitute an independent set. The largest of these sets has size at least n / 3 for a set with n nodes, The coloring and the computation of the sizes of the subsets can be performed in O(log n) time and O(n) work. Using this independent-set selection procedure until there is only one node left, the whole algorithm takes O(log^2 n) time and O(n) work. An immediate improvement is obtained by applying accelerated cascading: as soon as the problem size has been reduced to n / log n, that is, after at most loglog n / log (3/2) = O(loglog n) reduction rounds, we switch to pointer jumping. This idea reduces the total time consumption to O(loglog n * log n). A more refined algorithm runs in O(log n) time, but practically none of the deterministic algorithms can compete with the randomized one.

Sparse Ruling Sets

A ruling set in a graph is a subset S of the nodes so that each node is either in S or adjacent to a node in S. A k-ruling set is a subset S, so that each node has distance at most k from a node in S. The algorithm of this section is called the sparse-ruling-set algorithm, because it constructs a small subset of the nodes, which are called rulers, that are more or less regularly interspaced in the list. Once we have walked through the list from one ruler to the next, a new reduced list is constructed which has a node for each ruler and links with weight equal to the distance between the rulers. Once the ranks of the rulers are known, each non-ruler can ask the ruler from which it was reached for its rank and add an offset.

The following formulation of the algorithm is given for the DMM model. With minor modifications it can be turned into a PRAM or sequential algorithm.

1. Create an empty set of weighted links L. An entry (p, p', d) indicates a link from p to p' of length d.
2. For each node p, if p is a final node, mark each final node p as ruler and add a link (p, p, 0) to L. Otherwise mark p as unvisited.
3. Select s nodes randomly and uniformly from among the non-final nodes as rulers. Each of these rulers p initiates a wave: p prepares a packet (p, 0) to be sent to succ[p].

4.   while (there are still active waves) 
     {
       Send all prepared packets.
       for (each node p' that receives a packet (p, d)) 
         if (p' is a ruler) // wave of p ends in p'
         {
           Add a link (p, p', d + 1) to L.
           if (there are still unvisited non-ruler nodes)
           {
             Randomly select such a node p''. p'' initiates a new wave:
             p'' prepares a packet (p'', 0) to be sent to succ[p''].
           }
         }
         else // wave of p continues from p'
         {
           mark p' as visited with p as ruler and d + 1 as distance 
           thereof. Prepare a packet (p, d + 1) to be sent to succ[p'].
         }
     } 

5. Solve the weighted list-ranking problem defined by L.
6. For each ruler p extract the value of last[p] and dist[p] from the solution of the problem defined by L.
7. For each non-ruler node p', ask the ruler p of p' the values of last[p] and dist[p]. Set last[p'] = last[p] and dist[p'] = dist[p] - dist[p'], where dist[p'] is the previously marked distance.

In each round of the main while-loop exactly s previously unvisited nodes are reached, so all n nodes are visited in n / s rounds. Considering the number of unvisited nodes, it follows that the expected number s_t of waves that hits a ruler in round t, 0 <= t < n / s, is given by s_t = s * s / (n - t * s). Summing over t, gives the total number of newly created rulers:

  sum_t s_t = sum_{t = 0}^{n / s - 1} s * s / (n - t * s)
            = s * sum_{i = 1}^{n / s} 1 / i 
            = s * H_{n / s}.

In a DMM implementation, the algorithm performs n / s + O(1) all-to-all routings. During the main while-loop each node sends a packet containing the index of the ruler, the distance thereof and the destination, so the routing volume is 3 * n / p. In the last step, each non-ruler asks a question generating two answers, so the routing volume of the complete algorithm is 6 * n / p. So, the cost for a reduction of the problem size by a factor f = n / n' = n / s * log_e (n / s) is approximately

T(P, n) = c_i * O(n / P) + c_v * 6 * n / P + c_a * f * log_e f.

Taking s = o(n / (log n * loglog n)) gives n' = o(n / log n), which allows to solve the subproblem with pointer jumping at negligible costs. If s is chosen much smaller, s = o(n / (P * log P)), n' = o(n / P). In that case the data of the reduced list can be gossiped at negligible costs and the remaining problem can be solved internally. In that case the total routing volume is only 3 * n / P, but this goes at the expense of a much larger number of all-to-all routings. For small P this may be a very good idea though. If no new rulers are chosen in the course of the algorithm, the distance does not need to be sent along because it can be determined by counting the number of performed rounds. This has the disadvantage that we have to somehow deal with the initial nodes but allows to reduce the routing volume to just 2 * n / P, which appears to be the best achievable.

In practical applications, if P is somewhat larger, the focus should lie on reducing the number of routing operations. If s is chosen larger, for example, s = n / 50, the problem size is reduced by about a factor 10 with just over 50 all-to-all routings. Applying the same once more gives a reduction by about a factor 100, with 100 all-to-all routings and a routing volume of about 6.6 * n / P. In practice the choice of s should be made as a function of P, n, c_v and c_a in combination with estimates of the cost to solve the subproblem.

We remind that there are two main reasons why running a parallel program on a parallel computer with P PUs mostly gives much less than P speed-up:

• Overhead due to the communication.
• The parallel algorithm performs more work.

The weakness of sending too many packets can be alleviated by applying the following idea. Instead of performing all-to-all routings, at the end of each phase (which is the basic idea of the BSP approach!), it is better to intermix the computation and the communication: Round t is decomposed in P subrounds, denoted t_p, for 0 <= p < P. In subround t_p the algorithm running on PU i sends all it has to send to PU (i + p) mod P. Then it works on all the data it received and prepares the packet for PU (i + p + 1) mod P, and so on. Doing this, the routed packets become larger: a wave now has to wait only P / 2 subrounds before it is progressing again, going twice as fast. This idea does not only give improvements for list ranking, but for any problem in which a non-local structure has to be traversed.

The sparse-ruling-set algorithm may be considered to be a variant of independent-set removal. The difference is the speed at which the reduction is performed: independent-set removal gives a fixed small reduction factor. With the sparse-ruling-set algorithm the reduction factor can be chosen freely. A more interesting observation is that the structure of the sparse-ruling-set algorithm is very similar to that of sample sort: for a small subset of the elements the ranks are determined exactly, then the other elements are ranked between them.

Click here to see a sequential version of the sequential algorithm, the pointer-jumping algorithm and the sparse-ruling-set algorithm. It is important to observe that in the implemented variant of pointer jumping the distance may double by more than a factor two in one round because it may happen that node succ[i] gets processed before node i.

One-by-One Cleaning

• Two-Processor Algorithm
• DMM
• Circular Array

Two-Processor Algorithm

As before every node has fields succ[], last[] and dist[], where ultimately last[i] will point to the last node of the list of i and dist[i] at all times gives the distance from node i to last[i]. A final node i has succ[i] = i. There are two basic operations. The first is called autoclean, the other altroclean. If last[i] = j, we say that node i points to node j.

Autoclean(p) means that for all nodes i stored in some PU_p for which last[i] points to a node also stored in PU_p, we perform last[i] = last[last[i]] and dist[i] += dist[last[i]] until i is pointing to a node stored in another PU. This operation does not require any communication.

Altroclean(p, p') means that for all nodes i stored in some PU_p for which last[i] points to a node stored in PU_{p'}, we perform last[i] = last[last[i]] and dist[i] += dist[last[i]], by asking PU_{p'} for this information. An autoclean operation means that both involved PUs must send one packet. The first packet contains the index of the node the question is destined for, the second contains its last- and dist-value. In total three integers are exchanged.

Building on these two basic operations, we get the following algorithm for PU_p, 0 <= p < 2:

1.     for (each node i in PU_p)
         if (succ[i] == i) // i is a final node
         {
           last[i] = i - n;
           dist[i] = 0;
         }
         else
         {
           last[i] = succ[i];
           dist[i] = 1;
         } 

2.     autoclean(p); 

3.     altroclean(p, (p + 1) mod 2); 

4.     autoclean(p); 

5.     for (each node i in PU_p)
         last[i] += n; 

The algorithm is correct, because after the first autoclean no node points to a node in its own PU. Thus, after the altroclean, each node which does not point to a final node points to a node in the own PU. The second autoclean therefore finishes the problem, because if no node in PU_p points to a node in PU_{(p + 1) mod 2}, at the beginning of an autoclean, then it does not do so at the end of this operation. But by then no node points to a node in the own PU either, so all nodes must point to a final node.

The autocleans perform a linear number of internal operations. Essentially an autoclean is nothing but an internal list-ranking problem in which all nodes pointing to a node in another PU are considered to be final nodes. So, typically there will many very short lists.

T(2, n) = c_i * O(n) + c_v * 3 * n / 2 + c_a * 2.

DMM

The kernel of the algorithm consists of a loop in which all other PUs are addressed in a circular way:

1.     for (each node i in PU_p)
         if (succ[i] == i) // i is a final node
         {
           last_left[i] = last_rght[i] = i - n;
           dist[i] = 0;
         }
         else
         {
           last_left[i] = last_rght[i] = succ[i];
           dist[i] = 1;
         } 

2.     autoclean (p, last_left);
       autoclean (p, last_rght); 

3.     for (t = 1; t < P; t++)
       {
         altroclean(p, last_left, (p - t) mod P, last_rght);
         altroclean(p, last_rght, (p + t) mod P, last_left);
         autoclean (p, last_left);
         autoclean (p, last_rght);
       } 

4.     for (each node i in PU_p)
         last_rght[i] += n; 

We claim that at the end of the algorithm last_rght[i] gives the last node of the list of i for all i, 0 <= i < n. This can most easily be proven with invariants. We claim that after round t, 0 <= t < P, last_left[i] for all nodes i stored in PU_p is clean for p, p - 1, ..., p - t. By this we mean that last_left[i] does not point to any node stored in a PU with indices from {p, p - 1, ..., p - t}. In addition we claim that after round t, 0 <= t < P, last_rght[i] for all nodes i stored in PU_p is clean for p, p + 1, ..., p + t. For t = 0, we have performed autocleans, so in PU p, last_left and last_rght are clean for 0. More generally, assume the claims hold for some given t. So, at the beginning of round t + 1 last_left[i] is clean for p, p - 1, ..., p - t and last_rght[i] for p, p + 1, ..., p + t. Thus, if last_left[i] points to a node in PU_{p - t - 1}, putting last_left[i] = last_rght[last_left[i]] will return the index of a node which does not lie in a PU with indices in {p - t - 1, p - t, ..., p - 1}, but possibly it may lie in PU_p. Performing the autoclean establishes the claim for t + 1. In the above discussion, all PU indices are meant to be computed modulo P.

How about the time complexity? There are 2 * P - 2 altrocleans. Each operation consists of a question followed by an answer. Sending the second question in each pass of the loop together with the first answer, the number of start-ups is reduced to 3 * P - 3, this is only as much as in three all-to-all routings, without any indirect routing. On the other hand, if we have bad luck, there may be n / P questions in each altroclean. However, if the nodes indices are randomized, then we get much better performance: we may assume that in round t, the last values are the indices of nodes which are uniformly distributed over the remaining P - t PUs. If this assumption holds, the expected number of questions for an altroclean in round t equals (n / P) / (P - t). In total this gives sum_{i = 1}^{P - 1} (n / P) / (P - t) = n / P * H_P ~= n / P * log_e P. For every left-question we get back one answer, only a new last-value, for every right-question we get back two answers. So, for the presented version of the algorithm, the total routing volume is about 5 * (n / P) * log_e P. Altogether this gives

T(P, n) = c_i * O(n * log n) + c_v * 5 * n * log_e P + c_a * 3 * (P - 1).

Comparing with pointer jumping, we have slightly less work and routing volume because 5 * log_e P < 3 * log n, and much fewer start-ups. So, in practice pointer jumping will be worse. Comparing with the sparse-ruling-set algorithm we have lost almost a factor log_e P in work and routing volume, but the number or start-ups is far smaller. So, one-by-one cleaning will be better for small values of n / P^2, because then the start-up costs tend to dominate. There are other algorithms which are better for intermediate values, but it is correct to say that sparse ruling sets is best for large values of n / P^2 and one-by-one cleaning for small values of n / P^2. A combination of these two algorithms always gives reasonable results.

Circular-Array

The optimal PRAM algorithms can quite easily be turned into algorithms performing O(P) routing steps. This can be achieved on a linear array by geometrically reducing the problem size and the number of employed PUs until the reduced problem fits in one PU. Then it can be solved internally. Reversing the reduction steps solves the problem. The constant hidden in O(P) depends on how the details are worked out. It lies between 10 and 20. It is not clear how such algorithms that contract the problem in an ever smaller subset of PUs should profit from the wrap-around connection after the first reduction rounds, so on circular arrays, the number of routing steps is only slightly smaller.

The quite amazing thing is that on circular arrays list ranking can be solved much faster than this (and because circular arrays can be embedded on linear arrays with dilation 2, by extension even on linear arrays): P - 1 routing steps are sufficient! This can be achieved by an adaptation of the one-by-one cleaning algorithm. The idea is simple: the values of last_rght[] stay in the PU where they started, while the values of last_left[] travel around the network moving leftwards, that is, from PU_p to PU_{(p - 1) mod P}. In this way, the communication that is performed in the DMM algorithm can be resolved locally.

The packets sent in every step have size O(n / P) and the amount of computation is of the same order. Now one may complain that using P - 1 steps with O(n / P) work in each of them gives O(n) in total. This will not give much speed-up! This is true in general, but the important thing is the idea: this is an algorithm in which all communication is local. It is a kind of systolic list-ranking algorithm: two streams of information which are only interacting where meeting. The algorithm might even have practical importance. Maybe the "PUs" actually give highly connected clusters of PUs in a structure in which the clusters are more loosely connected (NUMA platforms). Or, we may have that n, the problem size, is greater than M, the size of the memory of a PU, but that n < P * M.

Exercises

1. Give a sequential list-ranking algorithm that performs at most 2 * n unstructured array accesses for a list of length n.

2. By a theoretical consideration we want to determine which algorithm in practice will be the best. The interesting point is that there is not a single algorithm dominating all others. There are several serious competitors: all of the presented algorithms and the sequential algorithm. Which algorithm is the best depends on the cost model and the values of P and n. Some algorithms which were not presented are even slightly better for certain ranges. We assume the DMM cost model with the following parameters: c_i = 10^{-9}, c_v = 10^{-8}, c_a = 10^{-5}. For a list with n nodes, estimate the number of internal operations performed by the sequential algorithm on 400 * n. The number of internal operations of the parallel algorithms can be estimated on 100 * routing_volume.
   • Which of the presented algorithms is best for very small problems? Draw a graph giving the maximum number of PUs that should be used for solving a problem of size n. Along the x-axis log_2 n increases from 14 to 34, along the y-axis log P increases from 0 to 10.
   • Compare pointer-jumping and one-by-one cleaning. Is one of them always better than the other, or do they both have their range of optimality? In the first case, indicate which algorithm is better, and prove this; in the second case, make a two-dimensional image like the one above and draw an approximate boundary.
   • Which of the presented algorithms is best for very large problems? Draw a graph of the maximum achievable speed-up as a function of n. Use again logarithmic scales along both axes. The cost of the sparse-ruling-set algorithm may be approximated by the cost for achieving a reduction by a factor P in one reduction round. Also indicate the number of PUs used for obtaining this maximum speed-up.

3. One of the tasks in the above exercise was to determine the maximum achievable speed-up in a practical setting. Now we take a theoretical stand point. Assume that we have arbitrarily many PUs in a DMM. Assume that c_i = c_v = c_a. The problem has size n. When answering the following questions you may use asymptotic notation, neglecting constant factors.
   • What is the maximum achievable speed-up? How many PUs do you need for this? What is the efficiency? Which algorithm do you use?
   • What is the maximum achievable speed-up if we insist that the efficiency should not converge to 0 in the limit for n goes to infinity? How many PUs do you need for this? Which algorithm do you use?

4. The list-ranking algorithm for a circular array with n PUs runs in n - 1 steps. This immediately gives an algorithm for a linear array running in 2 * n - 2 steps.
   • Show that even on a linear array n - 1 steps are sufficient.
   • Extend this idea to derive an improved algorithm for circular arrays running in n / 2 steps.

Trees

• Euler-Tour Technique
  • Construction
  • Applications
• Coloring Regular Bipartite Graphs
  • Sequential Algorithm
  • Parallel Algorithm
  • Unbalanced All-to-All Routing
• Tree Contraction
  • Problem Context
  • Raking
  • Expression Evaluation
  • Applications
• Lowest Common Ancestor
  • Definition and Applications
  • Special Cases
  • Reduction to Range-Minima Problem
  • Range-Minima Problem
• Exercises

Euler-Tour Technique

• Construction
• Applications
  • Tree Rooting
  • Preorder and Postorder
  • Depth
  • Size of Subtree
  • Inorder

Construction

suc(v_i, u) = (v, u_{(i + 1) mod deg_u}).

If the edges are numbered with different numbers, which can be performed by computing the prefix sums of the numbers of edges each PU is managing, then the defined successor relation simply gives a linked circular list. This can be turned into a non-circular list, by cutting it at a single place. For example, node 0 can cut the list between its first and last neighbor by setting suc(v_{deg_0 - 1}, 0) = null.

For a tree with n nodes and n - 1 edges, on a PRAM this construction takes O(log n) time with O(n) work for the prefix computation and O(1) time with O(n) work for the linking. On a DMM with P PUs, the prefix takes O(n / P) * c_i + log P * (c_v + c_a). If we assume that the information concerning edge (v_i, u) is stored in the PU responsible for node v_i, then the linking requires that the PU which is responsible for node u sends a packet to this other PU. Thus on a DMM the cost is dominated by that of an all-to-all routing in which each PU sends and receives one information for each edge of a node it is managing. If there are no extreme node degrees and the nodes are distributed evenly, this amounts to a k-k routing with k = 2 * n / P. More generally, the edges, after numbering them, should be distributed evenly over the PUs.

Basic Applications

• Tree Rooting
• Preorder and Postorder
• Depth
• Size of Subtree
• Inorder

Tree Rooting

Preorder and Postorder

Assume we have a rooted tree as input. Then to compute the preorder numbers, all edges pointing away from the root are given weight +1 the other edges are given weight 0. After running a weighted list-ranking algorithm, a node u can obtain its number from the rank of the edge coming from its parent. Switching the weights, that is giving weight +1 to the edges leading towards the root and 0 to the others computes the postorder numbers.

Depth

Size of Subtree

Inorder

• If there is a pair of consecutive edge in the list of the form (parent[u], u) followed by (u, parent[u]), this uniquely characterizes u as a leaf, then the second gets weight 1.
• If there is a pair of consecutive edges in the list of the form (v, u) followed by (u, w), with neither v nor w equal to parent[u], this uniquely characterizes u as an internal node, then the second gets weight 1.

Performing a weighted list ranking on the thus weighted list, a node u can find its inorder number from the rank of the first edge of the pair of edges that changed direction in u.

Coloring Regular Bipartite Graphs

• Sequential Algorithm
• Parallel Algorithm
• Unbalanced All-to-All Routing

Sequential Algorithm

Consider a bipartite graph with node sets V_1 and V_2, each with n nodes. Assume that the graph is regular of degree g. A first idea for constructing a coloring is to start allocating the colors to the first node of V_1, then the second and so on. When assigning the colors of node i, we should respect the conditions imposed by earlier assigned colors. If one is lucky this may work, but in general this approach will get stuck when we must assign a color c to an edge (i, j) while node j has an edge (i', j) for some i' < i, which was already assigned color c while coloring node i'. Another "greedy" strategy may also work, but not always. The idea is that one tries to assign the colors one by one. The algorithm gets stuck when, while assigning color c, one comes to a node i for which all the uncolored edges lead to nodes which already have an edge that was given color c before.

A working and efficient strategy is based on Euler splittings. For a g-regular graph G with g even, the graph can easily be split in two edge-disjoint subgraphs G_0 and G_1 each of degree g / 2. This is done by constructing an Euler tour of the graph, numbering the edges on the tour alternatingly 0 and 1. All edges which have been numbered 0 belong to G_0, all other edges to G_1. Because sequentially an Euler tour of a graph with n nodes and m edges can be constructed in O(n + m) time, this splitting costs time O(g * n).

In general, for a graph of degree g = 2^k, the algorithm consists of k rounds. In round i, 0 <= i < k, the algorithm is working on 2^i subgraphs in which all nodes have degree 2^{k - i}. Each of the 2^i operations in round i takes O(2^{k - i} * n) time, so the total amount of work in any of the rounds is O(2^k * n) = O(g * n) = O(m). Thus, in total over all rounds the algorithm takes O(k * m) = O(log g * m) time. For g which are not powers of 2, the construction of such colorings is considerably harder. After much research it has been established (by Cole, Ost and Schirra in Combinatorica, Vol. 21, 2001) that even the general case can be solved in O(log g * m) time.

Parallel Algorithm

The basic idea is to turn the problem of finding an Euler tour into a list-ranking problem. For every edge there is a node. These nodes have no predecessor and successor links, but for every edge e, there is a data element left[e] and rght[e]. left[e] gives the edge e' to which e is connected in its left endpoint i in V_1. This decision is made locally, without global coordination in node i. If left[e] = e', then left[e'] = e. rght[e] gives the edge e'' to which e is linked in its right endpoint j in V_2. A node of the list (corresponding to an edge of the graph) which is entered by a left[] link is left through a rght[] link and vice versa.

The list-ranking algorithms can, with minor modifications also be applied to figure out the minimum on a circular list and the distance thereof. This is now done for both lists (not knowing the global structure, it must be done for both). Of course we do not have a conventional list structure: each node has two "successors": one labeled "left" and one labeled "right". So, what is the appropriate successor? The answer is simple: the path starting over a "left" link must at even distances from the original node traverse "left" links and at odd distances "right" links. Further, notice that, because any cycle on a bipartite graph has even length, for any node e, if the distance to the node with the smallest index is even along one direction then it is also even along the other direction. This observation makes that we now can split the set of nodes (corresponding to the edges of the graph) in two: all those at even distances go into one subset, all those at odd distances into another.

Summarizing: first each node of the graph arbitrarily links up the edges connected to it in connected pairs; then a set of bidirectionally linked circular lists of a slightly modified nature with g * n nodes is constructed; then two list ranking problems are constructed. Clearly this whole problem can be solved in O(T_listrank(g * n)). Hereafter we have to solve two instances with each half as many edges. In a straightforward implementation this gives

T_color(P, n, g) = O(T_listrank(P, g * n) + 2 * T_color(P, n, g/2).

T_color(P, n, g) = O(log g * g * n / P * (c_i + c_v) + c_a * g * P).

This is a good result except for the number of startups. Organizing things slightly differently the number of startups can be minimized: after each splitting of the problem, each of the two subproblems should be concentrated in half of the available PUs. This gives a different recurrence relation:

T_color(P, n, g) = O(T_listrank(P, g * n) + T_color(P/2, n, g/2).

T_color(P, n, g) = O(log g * g * n / P * (c_i + c_v) + c_a * P).

If there are two independent problems which must be solved on a DMM or an interconnection network with P Pus, then it is mostly better to solve these problems in parallel each on P / 2 PUs than to solve them after each other each on P PUs.

On a two-dimensional m x m mesh the situation is even much better. Assuming that g * n > m^2, T_listrank(m^2, g * n) = O(T_route(m^2, g * n)) = O(g * n / m). So, we get

T_color(m^2, n, g) = O(g * n / m) + T_color(m^2 / 2, n, g / 2).

T_color(m^2, n, g) = O(g * n / m).

Of course the speed-up of this mesh algorithm is only m = sqrt(P), but this is the normal situation on meshes. There are very few problems like matrix multiplication and all-pairs shortest paths for which full speed-up can be achieved. All but very trivial problems require that at least all data are communicated once over the bisection of the network and therefore full speed-up can be achieved only if the sequential complexity of the problem is sufficiently high. A result like the obtained one is optimal up to constant factors in a two-dimensional world. The basic assumption of the DMM model that locality plays no role is definitely not valid without limitation on P.

Unbalanced All-to-All Routing

We have seen how to solve the unbalanced routing problem using randomization: sending all packets first to a random intermediate destination (or by shuffling the packets in a deterministic way), any routing pattern may be replaced by two almost balanced routing patterns. The main disadvantage of this approach is that it doubles the routing volume. A coloring approach offers an alternative solution which are mostly satisfied in practical contexts.

We consider a DMM with P PUs. Each PU has to send and receive in total k packets. Let a_{i, j} denote the number of packets PU_i has to send to PU_j. The conditions are

0 <= a_{i, j} <= k,
sum_{i = 0}^{P - 1} a_{i, j} = sum_{j = 0}^{P - 1} a_{i, j} = k.

Some form of sparsification is the key to the solution. Possibly we might group PUs together, but this inevitably leads to routing packets several times, and then we could just as well apply the simpler randomized routing. So, bundling the packets into superpackets is the only cost reduction we can try to apply. Let s be the size of the superpackets. That is, in PU_i the a_{i, j} packets which have to be routed to PU_j are packet into b_{i, j} = round_up(a_{i, j} / s) superpackets of size s, creating one partially filled superpacket. In total PU_i creates sum_{j = 0}^{P - 1} round_up(a_{i, j} / s) <= sum_{j = 0}^{P - 1} (a_{i, j} / s + 1) = k / s + P. So, this gives a bipartite graph with maximum degree k / s + P. The nodes do not need to have the same degree, but after computing prefix sums dummy packets can be added to make the graph regular. Because within certain limits s can be chosen freely, we may assume that k / s + P is a power of two.

On this reduced (k / s + P)-regular bipartite graph with P nodes on either side the presented Euler-tour-based coloring algorithm can be applied. This takes O(log(k / s + P) * (k / s + P) * (c_i + c_v) + c_a * P) time. This is a loss term. The second loss term comes from the fact that we have more than k / s superpackets: the subsequent routing takes (k / s + P) * (c_v * s + c_a). We want to determine a good choice for s and a condition on k and P so that all loss terms can be made negligible in comparison with the useful time c_v * k. Clearly P should be negligible in comparison with k / s. This allows to simplify the time for the coloring. Equating the loss terms gives the following equation:

log k * k / s = P * s.

Theorem: On a DMM with P PUs any k-k routing can be performed in (1 + o(1)) * k * c_v time if k = omega(P * log P).

Proof: Taking s = sqrt(k * log k / P), the loss terms are bounded by O(sqrt(k * log k * P) * c_v). If k = x * P * log P, sqrt(k * log k * P) = k * sqrt((1 + log x / log P + loglog P / log P) / x). Clearly for x = omega(1), this is o(k). End.

Tree Contraction

• Problem Context
• Raking
• Expression Evaluation
• Applications
  • Size of Subtree
  • Height

Problem Context

Evaluating arithmetic expressions in parallel does not appear to make sense because it is hard to imagine in which context sufficiently large expressions occur, but the expression-evaluation problem models other tree problems which are important and for which even have to be solved for large trees. An example of such a problem is that of computing the height of each node u in a tree, that is, the distance to the deepest lying leaf in the subtree of u. For this problem the Euler-tour technique does not appear to be suitable. As an expression-evaluation problem it can be solved by putting 0 in all leafs and using as operator max + 1, taking the maximum value of the children and adding 1. Another example is that of computing for each node u of a tree the minimum key value in the subtree of u. This problem can be solved as an expression-evaluation problem in which each node computes the minimum of its own value and that of its children.

For expression trees of moderate depth the evaluation problem can be solved efficiently in parallel by allocating a processor to each node connected to two leafs. This initial processor allocation can be performed using the Euler-tour technique. For each such a node its value can be computed in constant time. When as a result a new node connected to two leafs is generated, the allocated PU is taking over this job, otherwise it becomes passive. In each step the depth of the tree is reduced by one. Thus, for a tree of depth d with n nodes, the time is O(d + log n) and the work is O(n). If the tree is deep and skinny, for example a path, this way of evaluating gives very limited parallelization.

Raking

Many rakes can be performed in parallel, but it should be avoided to rake two leafs with the same grandparent. This is not hard to realize. The algorithm starts by consecutively numbering the leafs using the Euler-tour technique. Then the rake operation is first performed on all leafs with numbers 4 * k + 0, for some k, and then on all leafs with numbers 4 * k + 2, for some k. Hereafter each leaf with an odd numbers i is renumbered as (i - 1) / 2, and the process is repeated until there are only two leafs left. Because the raked leafs are always four numbers apart, they can never have the same grandparent. At this point we essentially use that the tree is binary. So, provided we know how to perform a rake operation, this correctly contracts the tree.

The above tree contraction can be performed very efficiently in parallel: the Euler-tour technique was analyzed before. On an EREW PRAM it can be performed in O(log n) time with O(n) work. Assuming that raking a leaf takes constant time, each contraction phase, consisting of two parallel rakings and a trivial renumbering, takes O(1) time. In each contraction phase the number of leafs is halved, so the whole contraction can be performed in O(log n) time. For some constant c, the work during the contraction is given by sum_t c * n / 2^t = O(n). In the time estimate it is assumed that there are no internal nodes of degree one, otherwise for a tree which has very few leafs only a few nodes disappear in every contraction phase until there is one leaf left. The easiest is to add a dummy leaf to any internal node with degree one. This leaf should get a neutral value so that the evaluated value remains the same for all real tree nodes.

Expression Evaluation

Invariant: At all times for any leaf u, a_u = 0 and b_u = val(u). For any internal node u, with children v and w, val(u) = (a_v * val(v) + b_v) op_u (a_w * val(w) + b_w).

• If op_v = +, then val(u) = a_w * val(w) + b_w + b_v. So, in this case we set (a_w, b_w) = (a_w, b_w + b_v).
• If op_v = *, then val(u) = a_w * a_v * val(w) + b_w * b_v. So, in this case we set (a_w, b_w) = (a_w * a_v, b_w * b_v).

The way the algorithm is specified, the correct value will be computed in the root of the tree, but not for all other nodes. In the following examples it will be discussed how to correct this. The idea is analogous to that for computing a prefix sum: the whole algorithm consists of two phases, one running upwards, one running downwards. The difference is now that operations may take place at different levels of the tree in the same parallel operation.

Applications

• Size of Subtree
• Height

Size of Subtree

This is an easy case. All nodes are given initial weight 1. In addition there are dummy leafs added to all internal nodes of degree one, these get weight 0. The only operation is the addition. It adds the value of both children to the weight. More precisely, when raking u with parent v and grandparent w which have weights a_u, a_v and a_w, respectively, then u and v are removed and a_w is replaced by a_w = a_w + a_v + a_u.

Notice that the above construction does not compute the size of the subtree for the node v which is skipped over. This problem can be solved easily. After the tree contraction there is a tree expansion, reversing the operations. So, all nodes which are excluded in phase i, i = 1, 2, ..., t, are reinserted in phase i = t, t - 1, ..., 1. After reinsertion phase i we may assume that all nodes which have been reinserted already know their correct values. So, by the time v is reinserted, u and its sibling u' at the time of the exclusion of v have been reinserted and obtained their final values a_u and a_{u'}. So, v can compute its value by setting a_v = a_v + a_u + a_{u'}.

Height

Lowest Common Ancestor

• Definition and Applications
• Special Cases
• Reduction to Range-Minima Problem
• Range-Minima Problem
  • Suboptimal Algorithm
  • Optimal Algorithm
  • Computing Most-Significant Different Bit

Definition and Applications

If the tree represents some kind of evolutionistic development, then the LCA of u and v gives the historical event where two "species" went their own independent way. If the tree T is a spanning tree of some undirected graph and we want to make it more like a DFS tree, then when discovering a cross edge (u, v) it may be clever to replace one of the edges from LCA(u, v) leading to u or v by the edge (u, v).

The distance from a node u to a node v in a tree is given by the sum of the distance from u to LCA(u, v) and that from v to LCA(u, v). The distance of u to LCA(u, v) is given by the difference of their depths in the tree. So, after computing the depths of all nodes, which can be performed in linear time, the distances between any pair of nodes can be computed in constant time plus the time for for computing the LCA. Replacing the depths by the sum of all weights on the path from the root to a node, this idea can even be used for computing the distances in weighted trees.

An analogous problem is when for a fixed weighted tree we frequently have to determine the heaviest edge on the path from a node u to another node v in the tree. So, instead of using the sum operator we use the maximum operator here. The difference is that the maximum operator has no inverse, so we cannot simply subtract the contribution from the path between the root and LCA(u, v). In this case we first have to construct an associated tree T'. We describe how this can be done sequentially. We need a union-find data structure in which all nodes are inserted. If T has n nodes, then T' has n leafs in one-one correspondence with the nodes. The leafs all get key value 0. The edges of T are sorted according to their weight and processed in order of increasing weight. For an edge (u, v) of weight weight(u, v), a new node w of T' with two children is created. w has key value weight(u, v). w is connected to the node corresponding to the component of u and the component of v. Hereafter these components are unified and w becomes the node representing this component. Hereafter, for any pair of node u and v the weight of the heaviest edge on the path from u to v in T is given by the key of the LCA of u and v in T'.

The above examples show that for several applications it is desirable to be able to answer the question "what is the LCA of u and v" for two nodes u and v of a tree. For any single such query the problem can be solved easily by walking towards the root from u and v marking the nodes on the path as visited. The first node that is reached for the second time is the LCA. If this walking is done by alternatingly making a step on the path from u and the path from v, the total number of visited nodes is proportional to the length of the path from u to v. The marks can be set back in the same time. This is clearly optimal and if there are several such queries, then these can be performed in parallel without problem (though concurrent reads may occur).

However, if there are very many such queries, it is better to preprocess the tree so that the queries can be faster than by walking along the tree. The problem of preprocessing the tree so that subsequent queries can be performed in constant time is called the LCA problem. So, the topic of this section is not how to parallelize the queries, which is trivial, but to show how this preprocessing, which sequentially can be performed in O(n) time for a tree with n nodes, can be performed in parallel.

Special Cases

More important is the case that the tree T is a perfect binary tree with n = 2^k - 1 nodes for some k > 0. The key idea is to perform an infix numbering of this tree, giving number 1 to the leftmost leaf. So, the node numbers are 1, ..., 2^k - 1. Let w = LCA(u, v) be the LCA of two nodes u and v. Let (u_{k - 1}, ..., u_0) and (v_{k - 1}, ..., v_0), be the binary expansions of the numbers of u and v respectively. Let x(u, v) be the index of the most-significant bit of u and v which are different. Using this notation, the binary expansion of the number of w can be expressed easily: w_i = u_i for all i > x(u, v); w_{x(u, v)} = 1; w_i = 0 for all i < x(u, v). Of course it is not necessary that the number of nodes is exactly a power of two minus one.

Theorem: For a perfect binary tree T with n = 2^k - 1 nodes a numbering of the nodes can be computed in O(log n) parallel time and O(n) work so that based on this numbering for any pair of nodes u and v of T a query LCA(u, v) can be solved in O(1) sequential time.

The idea for the perfect binary trees can be extended to all trees which are sufficiently like perfect binary trees. The sequential linear time algorithm is based on such an extension. We do not further consider this construction but assume the following as a fact:

Theorem: For any tree T with n nodes a tree T' of size O(n) can be constructed in O(n) time so that based on the information contained in T' for any pair of nodes u and v of T a query LCA(u, v) can be solved in O(1) time.

Reduction to Range-Minima Problem

The construction is again based on the Euler-tour technique. First an Euler tour is constructed running along all edges in the tree, starting and ending in the root node r. This tour is used to construct arrays index[] and depth[]. For a tree with n nodes each has length 2 * n - 1. index[i] gives the index of the node that is visited after traversing i edges of the tour. index[0] = index[2 * n - 2] = r. depth[i] gives the depth of node index[i]. It can be computed by weighting all list edges leading away from the root with +1 and all other edges with -1. Running a list-ranking algorithm, each of these arrays can be computed in O(log n) time and O(n) work.

From index[] and depth[] we compute the arrays first[] and last[]. These arrays have length n. For all u, 0 <= u < n, first[u] = min_{0 <= i < 2 * n - 1} {index[i] = u}, that is, it gives the first occurrence of u on the tour. Analogously, last[u] gives the last occurrence of j on the tour. Because the tree is connected all these values are well defined. For a leaf u, first[u] = last[u]. For an internal node u first[u] < last[u]. On an EREW PRAM these arrays can be constructed in O(1) time and O(n) work: For all i, 0 <= i < 2 * n - 1, it is checked whether depth[i - 1] < depth[i] or not. If this is the case, we know that the tour came to node u = index[i] from above. So, this must be the first visit to node u. In that case the testing PU sets first[u] = i. If depth[i + 1] > depth[i], this is the last visit to node u = index[i] and the testing PU sets last[u] = i.

For any pair of nodes u and v, there are five mutually exclusive possibilities:

1. u = v;
2. first[u] < first[v], last[v] < last[u]: u is an ancestor of v, LCA(u, v) = u;
3. first[v] < first[u], last[u] < last[v]: v is an ancestor of u, LCA(u, v) = v;
4. last[u] < first[v], LCA(u, v) is neither u nor v.
5. last[v] < first[u], LCA(u, v) is neither u nor v.

Lemma: If last[u] < first[v], LCA(u, v) = index[j], where j = min_{last[u] < i < first[v]}{depth[i] = d}, where d = min_{last[u] < i < first[v]}{depth[i]}.

Proof:

Let d and j be as defined and let w = index[j]. All nodes on the tour between last[u] and first[v] have depth larger than w. So, the tour has not left the subtree rooted at w. Thus, all these nodes, particularly u and v are descendants of w, that is, w is a common ancestor of u and v. Now consider another node w' which also lies on the subtour lying between last[u] and first[v]. Suppose w' is an ancestor of u. The tour is constructed so that once it leaves a node to a higher lying node it does not return. This holds in particular for w'. So, only after having visited all nodes in the subtree of w' the tour reaches w. These nodes in the subtree of w' will not be visited again. We know that hereafter the tour reaches v at least one more time because first[v] > j. Thus, v does not lie in the subtree of w' and thus w' is not a common ancestor. We conclude that on this subtour w is the only common ancestor of u and v. End.

Range-Minima Problem

• Suboptimal Algorithm
• Optimal Algorithm
• Computing Most-Significant Different Bit

Suboptimal Algorithm

p[i] = p_1[i], for all 0 <= i < n / 2;
p[i] = min{p_1[n / 2 - 1], p_2[i - n / 2]}, for all n / 2 <= i < n.

Assume that n = 2^k for some k >= 0. For other n, we can if necessary round up the length of the array to the next power of two. Let p_{l, m}[], 0 <= l <= k, 0 <= m < n / 2^l, be the array of prefix minima for the subarray of a[] of length 2^l starting at position m * 2^l. Define s_{l, m}[] to be the analogous array of suffix minima. By the above recursive construction all arrays p_{l, m} and s_{l, m} can be constructed in O(log n) time and O(n * log n) work:

T(n) = T(n / 2) + O(1);
W(n) = 2 * W(n / 2) + O(n).

We start with an example. Suppose we want to determine min_{393 <= i <= 438}{a[i]}. Then we first write these numbers u and v in binary: u = 393 = 00110001001, v = 434 = 00110110010. We have added a few leading zeroes to show the generality of the idea. Then we determine the most-significant different bit. In our example it occurs at position 5. The number w which is obtained by taking v and setting the 5 least-significant bits to zero is w = 0011010000 = 416 = 13 * 32. The number we are looking for is min{s_{5, 12}[9], p_{5, 13}[18]}. Here 23 = 393 - 12 * 32 and 18 = 434 - 13 * 32.

More generally, for numbers u and v, u <= v, the most-significant different bit x = x(u, v) is determined. Then the number w is constructed by taking v and setting the least-significant x bits to zero. Let d = 2^x and l = w / d. The definition of w and d is so that (l - 1) * d = w - d <= u and (l + 1) * d = w + d > v. The number we are looking for is min{s_{x, l - 1}[u - (l - 1) * d], p_{x, l}[v - l * d]}. This is so, because s_{x, l - 1}[u - (l - 1) * d] = min_{u <= j < w}{a[j]} and p_{x, l}[v - l * d] = min_{w <= j <= v}{a[i]}.

It is interesting to notice that the number w is just the index of the LCA of nodes u and v in a perfect binary tree with inorder indexing with smallest index 1, the special case of the LCA problem treated above. This shows the close relation between the LCA problem and the range-minima problem. In an implementation it is not necessary to make this tree structure explicit, the set of prefix and suffix arrays can be managed by a double indexed array of arrays of integers. The given computation shows which two values in this construct must be accessed to compute the range minimum from. On a RAM this takes constant time.

Lemma: For an array a[] of length n, in O(log n) time and O(n * log n) work a data structure requiring O(n * log n) storage can be created which allows it to solve range-minima queries on a[] in constant time.

Optimal Algorithm

To reduce the work and memory consumption, the array a[] of length n is divided in n / log n subarrays of length log n each. For each of these subarrays an optimal sequential algorithm is applied. For each of the subarrays this takes O(log n) time and work. So, the parallel time is O(log n) and the work is O(n). The created structure allows to perform range-minima queries inside any of the subarrays in constant time. Along with this construction we can also construct an array b[] of length n / log n consisting of the minima in each of the subarrays. The above given suboptimal algorithm is now applied to b[] rather than a[]. This takes O(log n) time and O((n / log n) * log(n / log n)) = O(n) work. So, the whole construction takes O(log n) time and O(n) work.

It remains to specify how this two-stage structure can be used to perform range-minima queries in constant time. For numbers u and v, u <= v, belonging to subarray i_u and i_v, respectively, we distinguish two cases (with a subcase in the second):

• i_v = i_u. The range-minimum query can be answered using the structure which was constructed for subarray i_u. subarray i_v.
• i_v > i_u. In general the range minimum is now computed as the minimum of three contributions: the minimum of the suffix minimum of u in subarray i_u; the prefix minimum of v in subarray i_v; and the minimum of b[] over the range from i_u + 1 up to i_v - 1. Notice that if i_v = i_u + 1, the last range is empty. In that case only the first two contributions should be taken.

Now the construction is perfect except that we have not told how to sequentially preprocess an array in linear time so that range-minima queries can be performed in constant time.

Theorem: For an array a[] of length n, in O(log n) time and O(n) work a data structure requiring O(n) storage can be created which allows it to solve range-minima queries on a[] in constant time.

Corollary: For a tree with n nodes, in O(log n) time and O(n) work a data structure requiring O(n) storage can be created which allows it to solve LCA queries for any pair of nodes in constant time.

Computing Most-Significant Different Bit

For k-bit numbers, any single computation of intlog() can trivially be performed in O(k) time. Using a binary-search strategy, based on comparisons and shifts, this can easily be reduced to O(log k) time. Using bit parallelism further improvements are possible, but we do not obtain a constant-time operation this way. If we have to perform many computations of intlog(), it may be better to first compute an array intlog[] of length n = 2^k containing all values. A later call intlog(i) can then be answered in constant time by returning intlog[i].

Sequentially all n values of intlog[] can easily be computed in O(n) time:

  void computeIntlog(int[] intlog, int n)
  {
    intlog[0] = -1; // Artificial value
    intlog[1] =  0; 
    for (int i = 2; i < n; i++)
      if ((i & (1 << (intlog[i - 1] + 1) == 0)
        intlog[i] = intlog[i - 1];
      else
        intlog[i] = intlog[i - 1] + 1;
  }

The given subroutine cannot easily be parallelized. There are several alternative constructions. If optimal work would not be a major concern, then we could apply any of the sequential algorithms. For a sequential time T(n), this gives work T(n) * n. We know that T(n) = O(loglog n) or less. It is very simple to make the algorithm optimal: first compute intlog[i] for all values 0 <= i < sqrt(n). This takes T(n) time and O(sqrt(n) * T(n)) = O(n) work. Then for any number i, sqrt(n) <= i < n, set intlog[i] = intlog[i >> (k / 2)] + k / 2. This takes constant time and O(n) work.

Theorem: If T(n) = O(loglog n) denotes the time for sequentially computing the value of intlog(i) for any number 0 <= i < n, then the array intlog[] with intlog[i] = intlog(i) for all i, 0 <= i < n, can be computed in parallel in O(T(n)) time and O(n) work.

Exercises

1. For a tree T the task is to number all leafs consecutively, that is, the numbering should be so that for any node u the set of all numbers of the leafs equals {l, l + 1, ..., h} for some l and h. Describe how the Euler-tour technique can be used to compute such a numbering.

2. Sequentially, coloring a bipartite graph can also be solved by repeatedly determining a maximum matching. Matching is more expensive than coloring, but different from the Euler-tour-based splitting technique it can even be applied for graphs of odd degree.
   • Why should Euler-tour splitting not be applied for graphs of odd degree?
   • The task is to color a bipartite graph with n nodes in both subsets of nodes and degree g with the minimum number of colors. If possible Euler-tour splitting should be applied, otherwise matching. How many matchings may be needed at most?
   • Reduce the number of matchings as far as possible. Hint: two disjoint regular subgraphs of degree g_1 and g_2 respectively can be combined to one subgraph of degree g_1 + g_2.

3. In the raking procedure it was used that the tree is binary. Trees of other degrees can be first made binary, but if the degree of any node is bounded by some constant d this is not necessary. Describe how to generalize the algorithm. Express the time consumption for a tree with n nodes in n and d.

4. Suppose we have a tree in which all nodes have degree at most two. Each node holds a key and for each node u we would like to compute the minimum value of the keys in the subtree rooted at u, u included. In the text it was suggested to use the tree-contraction method for this. Work out the details, particularly describe how to handle the internal nodes of degree one. For a tree with n nodes the algorithm should run in O(log n) time and O(n) work on an EREW PRAM independently of the tree structure.

5. Assume we have a black-box subroutine which allows to compute the minimum value m in a specified interval [l, h] of an array a[] of length n. Describe how this routine can be used to compute an index j, l <= j <= h so that a[j] = m. We assume that the range-minima subroutine only works on integer arrays, so that we cannot simply add the index as secondary key to the values of a[]. Distinguish the case that it is given that all values of a[] are different and the case that values may be the same.

6. This question deals with computing the least-significant different bit. All considered numbers have at most k bits for some even k > 0. n = 2^k. It is supposed that there is an array intlog[] of length n containing for all i, 0 <= i < n, the value of round_down(log_2 i).
   • Describe how to construct in O(n) sequential time an array reverse[] of length n containing for all i, 0 <= i < n, the value of the reversal of i, reverse(i). If i_j denotes bit j of i, then reverse(i)_j = i_{k - j - 1}. Hint: Start with i = 0 and reuse earlier computed values.
   • Use the values in intlog[] and reverse[] to define a constant-time operation lsdb(,), computing the least-significant different bit for any pair of numbers.
   • Show how to compute reverse[] in parallel with at most O(n) work. Try to obtain a running time of O(loglog n).

Graphs

Strings

Computational Geometry

Numerical Algorithms