Physics Sorting

C HAPTER 9 Sorting 1 The correctness of this compare-split operation can be proved using induction on the number of compareexchange operations. Processor Pi has a sequence of k elements x = x 1 , x 2 , . . . , x k  sorted in increasing order and processor P j has k elements y = y1 , y2 , . . . , yk  sorted in decreasing order. During the nth step of the compare-split operation processors Pi and P j exchange their nth elements. Processor Pi keeps the max{x n , yn } and P j keeps min{x n , yn }. The inductive proof is as follows: 1. After the first step of the compare-split operation, Pi holds z = max{x 1 , y1 } and P j holds w = min{x1 , y1 }. For the operation to be correct z must be larger than at least k elements of the combined sequences x and y, and w must be smaller than at least k elements of the combined sequence x and y. Since sequence x is sorted in increasing order, we know that x 1 is smaller than the remaining k − 1 elements of sequence x. Since w = min{x 1 , y1 }, w is either x 1 or y1 . If w = x 1 , then x 1 < y1 ; thus, y1 and x 2 , x 3 , . . . , x k are the k elements that are larger than w. If w = y1 , then the elements of sequence x are the k elements that are larger than w. Therefore w is correctly assigned to processor P j . Similarly we can show that z is also correctly assigned to processor Pi . 2. Assume that up to step n −1 of the bitonic-split operation, elements were correctly assigned to processors Pi and P j . 3. We need to prove that w = min{x n , yn }, and z = max{x n , yn } are assigned to processors Pi and P j correctly. Again, consider the case where w = x n . w is smaller than the remaining k − n elements of sequence x. Also, w is smaller than yn . However, since sequence y is sorted in decreasing order, y is smaller than the first n − 1 elements of sequence y. Thus, there are k elements that are larger than w. If w = yn , then w is smaller than the n − k + 1 elements of x (that is, x n , x n+1 , . . . , x k ) and also is smaller than the previous n − 1 elements of y (that is, y1 , y2 , . . . , yn−1 ). Therefore, in either case, w is correctly assigned to P j . Similarly, it can be shown that z is assigned to Pi correctly. As mentioned in the problem, the compare-split operation can terminate as soon as x n ≥ yn . The reason is that xi ≥ yi , ∀i > n, as the x sequence is sorted in increasing order and the y sequence is sorted in decreasing order. Consequently, if the compare-exchange operations are performed, all elements x 1 ∀i > n go to processor Pi and all the yi elements ∀i > n go to processor P j . Because, the compare-exchanges can stop as soon as x n ≥ yn , the run time of this operation varies. If we have to sort the sequences after each compare-split operation, the run time of this formulation will be (k log k). However, at the end of the compare split operation, processor Pi has two subsequences, one from x and one from y. Furthermore, each of these subsequences is sorted. Therefore, they can be merged in linear time to obtained the final sorted sequence of size k . An advantage of this operation is that requires no extra memory. Fox et al. [fox] provide further details of this compare-split operation. Table 9.1 Communication requirements of row-major mapping. 1st iteration 2nd iteration 3rd iteration .. . (d/2)th iteration (d/2 + 1)st iteration (d/2 + 2)nd iteration (d/2 + 3)rd iteration .. . dth iteration 1 1 1 + + 2 2 + 4 1 1 1 1 + + + + 2 2 2 2 + + + + 4 4 4 4 + + + + ··· ··· ··· ··· + + + + 2 2 d/2−1 2 d/2−1 d/2−1 2 + + + 1 1 1 + + 2 2 + 4 1 + 2 + 4 + ··· + 2d/2−1 + 1 + 2 + 4 d/2−1 + ··· + 2d/2−1 2 A relation R on A is • reflexive if for every a ∈ A (a , a ) ∈ R. • antisymmetric if for every a, b ∈ A with (a, b) ∈ R and (b, a) ∈ R, then a = b. • transitive if for every a, b, c ∈ A whenever (a, b) ∈ R, and (b, c) ∈ R, then (a, c) ∈ R . From the above definitions the fact that ≤ is a partial ordering can be easily shown. 3 To prove (a) note that there is an element ai (0 < i < n/2) of sequence s such that for all j < i, min{a j , an/2+ j } belongs to {a0 , a1 , . . . , an/2−1 } and for all i < j < n/2 min{a j , an/2+ j } belongs to {an/2 , an/2+1 , . . . , an−1 }. Similarly, for all j < i, max{a j , an/2+ j } belongs to {an/2 , an/2+1 , . . . , an−1 } and for all n/2 > j > i, max{a j , an/2+ j } belongs to {a0 , a1 , . . . , an/2−1 }. Therefore, s1 = {a0 , a1 , . . . , ai , an/2+i +1 , . . . , an−1 } and s2 = {an/2 , an/2+1 , . . . , an/2+i , ai +1 , . . . , an/2−1 }. Note that both s1 and s2 are bitonic sequences. Parts (b) and (c) can be proved in a similar way. For more detailed proof of the correctness of the bitonic split operation refer to Knuth [Knu73], Jaja [Jaj92], Quinn [Qui87], or Batcher [Bat68]. 5 From Figure 9.10 (page 389) we see that bitonic merge operations of sequences of size 2k take place during the kth stage of the bitonic sort. The elements of each bitonic sequence merged during the kth stage correspond to wires whose binary representation have the same (log n − k) most significant bits (where n is the number of element to be sorted). Because each wire is mapped to a hypercube processor whose label is the same as that of the wire’s, each sequence is mapped onto a k-dimensional subcube. Since the most significant (log n − k) bits of each bitonic sequence of length 2 k are different, these sequences are mapped onto distinct subcubes. 7 Let n be the number of processors such that d = log n . (a) From Figure 9.11 (page 391) we can see that the communication requirements are shown in Table 9.1. Let A be the communication required during the first d/2 iterations. Then d/2−1 A=  i =0 i  j=0 2j = d/2−1   22i − 1  √ 12 =2 n−2− log n. i =0 Let B be the communication required during the first d/2 steps of the last d/2 iterations. Then B= d d/2−1 i 1 2 = log n(√n − 1) = 1 √n log n − 1 log n. 2  2 2 2 i =0 Table 9.2 Communication requirements of snakelike mapping. 1st iteration 2nd iteration .. . (d/2)th iteration (d/2 + 1)st iteration (d/2 + 2)nd iteration (d/2 + 3)rd iteration .. . dth iteration 1 1 + 2 1 1 1 1 + + + + 2 2 2 2 + + + + ··· ··· ··· ··· + + + + 2d/2−1 2 d/2−1 d/2−1 2 2 d/2−1 + + + √ 1+ n √ 1+ n √ 1+ n 1 + 2 + ··· + 2d/2−1 + 1+ √ n + + 2 2 + 4+ + 2 + 4+ √ √ n n + ··· + 2d/2−1 Table 9.3 Communication requirements of the row-major shuffled mapping. 1 1 1 1 1 1 1 1st iteration 2nd iteration 3rd iteration 4th iteration 5th iteration 6th iteration .. . (d/2 − 1)th iteration dth iteration 1 1 2 ... d −1 + 1 + 2 + 2 + 4 + 4 + ··· + 1 + 2 + 2 + 4 + 4 + ··· + 2 d/2−1 d/2−1 + 2 + + + + + 1 1 1 1 1 3 + + + + 4 2 2 2 2 + + + 5 2 2 2 + + 4 4 6 + d 4 + 2 d/2−1 The total communication performed by the row-major formulation of the bitonic sort is 1√ √ 3 n log n + 4 n − log n − 4. 2 2 From Figure 9.11 (page 391) we can see that the communication requirements are shown in Table 9.2. Note that the communication requirements for the first d/2 iterations are similar to the row-major mapping. However, during the last d/2 iterations, the communication requirements of snakelike mapping is higher than that of the row-major mapping. The additional communication is approximately B + 2A = (b) √ d/4 n2  i ≈ i =1 √ n log2 n . 16 Therefore, the total communication requirements is √ (c) √ n log2 n + 1 √ 3 n log n + 4 n − log n − 4. 16 2 2 The communication requirements of the row-major shuffled mapping is shown in Table 9.3. To derive the total communication required we first concentrate on the odd numbered columns. Note that the sum of the odd numbered columns is d/2−1 2  i =0 i √ i  2 = 4 n − log n − 4. j=0 Now note that in Table 9.3, column two differs from column one by 1, column four differs from column two by 2 and so on. Therefore, the sum of the even numbered columns is equal to that of the odd numbered columns minus d/2−1  √ 2i = n − 1. i =0 Therefore, the total communication requirements of row-major shuffled mapping is √ √ √ 2(4 n − log n − 4) − n + 1 = 7 n − 2 log n − 7. Therefore, row-major shuffled mapping requires less communication than the other two mappings. However, Nassimi and Sahni [NS79] have developed variants of bitonic sorting for the row-major and the snakelike mapping that have the same communication requirement as the row-major shuffled mapping. 8 The solution is presented by Knuth [Knu73] (merging network theorem). 9 In a ring-connected parallel computer, input wires are mapped onto processors in a way similar to the hypercube mapping. The wire with label a is mapped onto processor with label a . The total amount of communication performed by the processors is: j log n−1   j=0 2i = 2n − log n − 2 ≈ 2n i =0 The parallel run time, speedup, and efficiency when p processors are used to sort n elements are: TP S E 10 =  = =  p log n  p n log n   n + log n)log2 p + (n) (n p p n p     +  1pnlog2 p + (n)   2     log p log p +  log +  logp n 1 −  log n n For this formulation to be cost optimal p/ log n = O(1). Thus, this formulation can efficiently utilize up to p = (log n) processors. The isoefficiency function is ( p2 p ). This proof is similar to that presented in Solution 8. 11 Recall that in the hypercube formulation of shellsort, each processor performs compare-split operations with each of its neighbors. In a d-dimensional hypercube, each processors has d neighbors. One way of applying the spirit of shellsort to a mesh-connected parallel computer is for each processor to perform compare-split operations with each of its four neighboring processors. The power of shellsort lies in the fact that during the first few iterations, elements move significant closer to their final destination. However, this is not the case √ with the mesh formulation of shellsort. In this formulation elements move at most p processors. An alternate way of performing shellsort on a mesh-connected computer is shown in Figure 9.1. The mesh is divided in to equal halves along the y-dimension, and corresponding processors compare-exchange their elements. Then, the mesh is halved and the same is repeated for the smaller meshes. This is applied recursively until communication takes place between adjacent mesh rows. A similar sequence of compare-exchanges is performed along the x-dimension. Note that the processors paired-up for compare-split are similar to those paired-up in the hypercube formulation of shellsort. 12 The sequential shellsort algorithm consists of a number of iterations. For each iteration i there is an associated distance di . The n elements to be sorted are divided into groups, such that each group consists of elements whose indices are multiples of di . These groups are then sorted. In each iteration di decreases. During the last step, di = 1, and the algorithms sorts the entire list. Shellsort can be parallelized by performing each iteration in parallel. In order to take advantage of the fact that as the sort progresses sequences tend to become almost sorted, a sorting algorithm such as the odd-even transposition must be used. Step 1 Step 2 Step 3 Figure 9.1 An example of performing shellsort on a mesh with 16 processors. 13 [fox] The worst case can be summarized as follows. Let n be the number of elements to be sorted on a d dimensional hypercube. Each processors is assigned n/2d = m elements. Now if more than m of the largest 2m elements occupy a particular pair of processors in the first compare-split step, some of the items will be put in the wrong subcube in this first step. The subsequent (d − 1) stages will just rearrange the items in each (d − 1)-dimensional subcube, but these elements will stay in the wrong subcube. Consequently, in the final odd-even transposition sort, these items must travel a long distance to reach the top of the list, requiring 2d−1 − 1 steps. The probability that such a collection of items will occur, for random initial lists, is given by: m−1 (2m)! 1 m !m ! d−1 2 This probability becomes small rapidly as m increases.  14 The CREW formulation of quicksort that assigns each subproblem to a different processor can be extended for a p processor system by stopping the recursive subdivision as soon as p subproblems have been created. At this point, each of the p processors internally sorts one of the p subsequences. At the end of this step, the n elements are sorted. If we assume that the pivot selection leads to perfectly balanced splits, then, if n = 2d and p = 2k , the parallel run time is: n n k−1 n n n 1 ). log + = log + 2n(1 − i p p p p p  2 i =0 The efficiency of this formulation is: E = = n log n n log n − n log p + 2np − 2n 1 1 − log p/ log n + 2 p/ log n − 2/ log p From this equation, we see that for a cost efficient formulation p = O(log n). Therefore, this formulation cannot use more that O(log n) processors. The isoefficiency function is ( p2 p ). 15 The algorithm is described by Chlebus and Vrto [CV91]. 16 The expected height of the tree is less than or equal to 3 log n + 6. The analysis is discussed by Chlebus and Vrto [CV91]. 17 The solution is provided in the first edition of this book. 18 The proof that the hypercube formulation of quicksort correctly sorts n elements on p processors, is as follows. After the first split, p/2 processors will end up having elements greater than the pivot while the remaining p/2 processors have elements smaller than the pivot. In particular, the processors whose processor label has a zero at the most significant bit get the smaller, while the other processors get the larger elements. Thus, the first split, puts the smaller elements in the right half-cube. Assume that after k splits, the elements are sorted according to the k most significant bits of the processors’ labels. That is, if we treat each subcube of size 2 d−k as a single processor, and the elements assigned to this subcube, as a block of elements, then these blocks are sorted. Now, during the (k + 1)st split each block is partitioned into two smaller blocks. In each subcube of size 2d−k , each of these two smaller blocks is assigned to one of the two 2d−k −1 half-cubes contained into a 2d−k subcube. This assignment is such that the block with the smaller elements is assigned to the half-cube that has a zero at the (d − k − 1)st bit, while the block with larger elements is assigned to the half-cube that has a one. Therefore, the label of the half-cube determines the ordering. Furthermore, since this is a local rearrangement of the elements in each 2d−k subcube, the global ordering has not changed. 20 This pivot selection scheme was proposed by Fox et al., [fox]. The performance of the algorithm depends on the sample size l; the larger the sample the better the pivot selection. However, if the sample is large, then the algorithm spends a significant amount of time in selecting the pivots. Since l cannot be too large, a good choice seems to be l = (log n). Furthermore, if the distribution of elements on each processor is fairly uniform, this pivot selection scheme leads to an efficient algorithm. 21 This problem was motivated by the (n) worst case selection algorithm is described in [CLR89]. Following a proof similar to that in [CLR89] it can be shown that this scheme selects a pivot element such that there are (n) elements smaller than the pivot, and (n) elements larger than the pivot. Initially, the sequences stored in each processor are sorted. For each processor pair, computing the median of the combined sequences takes (log n) time. Sorting 2i −1 medians can be done using bitonic sort, which takes (i 2 ) time. Finally, the median of the median is broadcasted in (i) time. 22 The barrier synchronization is not necessary. Removing it does not affect neither the correctness nor the asymptotic performance. 23 The parallel run time of the hypercube formulation of quicksort if we include the tc , ts , and tw costs is approximately n n TP = log   p n tc + tw + ts log p + log2 p(ts + tw ) The efficiency is: p p E = n log ntc n log ntc − n log ptc + (ntw + pt s ) log p + p log2 p(ts + tw ) log p tw p log p ts p log2 p + (9.1) + −1  log n tc n log n tc ts + tw = 1− t c n log n If K = E/(1−E), then from Equation 9.1 we have that the isoefficiency function due to the (log p/ log n)(tw /tc ) term is  logpn log + log p tw tw 1 ⇒ log n = = K log p ⇒ n = 2(K t w /tc ) log p = p K t tc log n tc K w /tc ⇒ n log n = tw K p K tw /tc log p. tc The isoefficiency function due to the ( p log p)/(n log n)(ts /tc ) term is K (ts /tc ) p log p, due to the ( p log2 p)/(n log n)(ts /tc ) term is K (ts /tc ) p log2 p, and due to the ( p log2 p)/(n log n)(tw / tc ) term is K (tw / tc ) p log2 p. Therefore, the most significant isoefficiency functions are: tw ts tw K p K tw /tc log p, K p log2 p. K p log2 p, and tc tc tc Note that the asymptotic complexity of the first term depends on the value of K and the tw /tc ratio. (a) tc = 1, tw = 1, ts = 1 E = 0.5 K = E /(1 − E ) = 0.5/(1 − 0.5) = 1. The isoefficiency functions are: p log p, p log2 p, and p log2 p. Therefore, the overall isoefficiency function is p log2 p. E = 0.75 K = 0.75/(1 − 0.75) = 3. The isoefficiency functions are: 3 p 3 log p, 3 p log2 p, and 3 p log2 p. Therefore, the overall isoefficiency function is 3 p 3 log p . E = 0.95 K = 0.95/(1 − 0.95) = 19. The isoefficiency functions are: 19 p 19 log p, 19 p log2 p, and 19 p log2 p. Therefore, the overall isoefficiency function is 19 p 19 log p . (b) tc = 1, tw = 1, ts = 10 E = 0.5 K = E /(1 − E ) = 0.5/(1 − 0.5) = 1. The isoefficiency functions are: p log p, p log2 p, and 10 p log2 p. Therefore, the overall isoefficiency function is 10 p log2 p. E = 0.75 K = 0.75/(1 − 0.75) = 3. The isoefficiency functions are: 3 p 3 log p, 3 p log2 p, and 30 p log2 p. Therefore, the overall isoefficiency function is 3 p 3 log p . E = 0.95 K = 0.95/(1 − 0.95) = 19. The isoefficiency functions are: 19 p 19 log p, 19p log2 p , and 190 p log2 p . Therefore, the overall isoefficiency function is 19 p 19 log p . (c) tc = 1, tw = 10, ts = 100 E = 0.5 K = E /(1 − E ) = 0.5/(1 − 0.5) = 1. The isoefficiency functions are: 10 p 10 log p, 10p log2 p , and 100 p log2 p . Therefore, the overall isoefficiency function is 10 p 10 log p . E = 0.75 K = 0.75/(1 − 0.75) = 3. The isoefficiency functions are: 30 p 30 log p, 30p log2 p , and 300 p log2 p . Therefore, the overall isoefficiency function is 30 p 30 log p . E = 0.95 K = 0.95/(1 − 0.95) = 19. The isoefficiency functions are: 190 p 190 log p, 190 p log2 p, and 1900 p log2 p. Therefore, the overall isoefficiency function is 190 p 190 log p . Therefore, as the efficiency increases, the hypercube formulation becomes less scalable. The scalability also decreases, as the tw /tc ratio increases. 24 The solution is provided in the first edition of this book. 25 [SKAT91] Let n be the number of elements to be sorted on a p-processor mesh-connected computer. Each processor is assigned n/ p elements. In order to extend this algorithm to apply to the case of multiple elements per processor, the following changes have to be made: (a) (b) (c) One of the elements in the first processor is chosen as the pivot. Pivot distribution remains the same. This √ step takes ( p) time, since the worst case distance traveled in the vertical and horizontal directions is √ p hops. On receiving the pivot, each processor divides its elements into sets of elements smaller and larger than the pivot. Also, it maintains the number elements in each of these two sets. Information propagated from the leaves of the binary tree to the root takes into account the number of smaller and larger elements at each node. Since there are n/ p elements per processor, it takes (n/ p) time to divide the set. √ Propagating the information back takes ( p) time. The information propagated from the root to the leaves about the next free processor is modified to account for multiple elements per processors. In particular, the two partitions are separated at a processor boundary. Therefore, the number of elements per processor may differ somewhat in the two partitions. The next free processor is specified along with the number of elements that can be added to the proces√ sors. This takes ( p) time. (d) Elements are moved in large messages. Notice that a message from one processor may have to be split across two destination processors. Then the portion for the other destination processor can be sent in one √ hop. This takes (n/ p) time. Since we assume perfect pivot selection, the algorithm requires (log p) iterations, the parallel run time is: TP = 26 n n log p  √ ( p) +      n + √ (log p) p p p √ √ p log p p log p). The isoefficiency function of the algorithm is (2 The CRCW PRAM algorithm can be easily adapted to other architectures by emulation. + n CREW PRAM The main difference between a CREW and a CRCW PRAM architecture is the ability to perform a concurrent write. However, a concurrent write operation performed by n processors can be emulated in a CREW PRAM in (log n) time. Therefore, the parallel run time of the enumeration sort is (log n). EREW PRAM In this model, in addition to requiring (log n) time to emulate a concurrent write operation, it needs a (log n) time to perform a concurrent read operation. Therefore, the parallel run time of the enumeration sort is (log n). Hypercube The hypercube formulation is similar to the EREW PRAM formulation. Each read or write operation takes (log n) time. Also, the elements can be permuted to their final destination in (log n) time, since n 2 processors are used. Therefore, the parallel run time of the enumeration sort is (log n). √ Mesh In a mesh-connected computer we can emulate a concurrent write or read operation in ( n) time. √ Therefore, the parallel run time of the enumeration sort is ( n). When p processors are used, we assign n/ p elements to each processor. Each processor now performs computations for each of its elements, as if a single element is assigned to it. Thus, each physical processor emulates n/ p virtual processors. This results in a slowdown by a factor of n/ p . 27 Since the sequential run time is (n), the speedup and efficiency are: S = E = = (n) (n/ p ) + ( p log p ) (n) (n) + ( p 2 log p) 1 1 + (( p 2 log p)/n) Therefore, the isoefficiency function of bucket sort is ( p 2 log p ). Comparing the isoefficiency function of bucket sort with the other sorting algorithms presented in the chapter, we see that bucket sort has better isoefficiency function than all but the quicksort algorithm for the hypercube. Recall that the bucket sort algorithm assumes that data is uniformly distributed. Under the same assumption the hypercube formulation of quicksort has good performance, and its isoefficiency functions is ( p log2 p). Note also, that the isoefficiency function of bucket sort is similar to that of sample sort, under the assumption of uniform distribution. However, if the distribution is not uniform, then the parallel run time of bucket sort becomes 28   n TP =  + (n) + ( p log p). p sort is cost optimal only when p = (1). Since the processor-time product is (np), bucket The proof is presented by Shi and Schaeffer [SS90]. 29 When the sizes of the subblocks are roughly the same, the speedup and efficiency of sample sort are: S = E = (n log n) ((n / p) log(n / p )) + ( p2 log p) + ( p log(n/ p )) + (n/ p ) + ( p log p ) (n log n) (n log(n / p )) + ( p3 log p ) + ( p2 log(n / p)) + (n ) + ( p2 log p)       −1   p 3 log p p2 p 2 log p log p 1 + + + +  log n log n n log n n n log n The isoefficiency function is ( p 3 log p). When the size of the subblocks can vary by a factor of log p, then there is a processor that has a block of size  =  1+ n/ p + log p. The time spent in communication the is (n/ p + log p ) + O( p log p ). However, this new term will only add a ( p log p)/(n log n) term in the denominator of the efficiency function. Therefore, it does not affect the asymptotic scalability. 30 The algorithm can be easily modified to sort the sample of p( p − 1) elements in parallel using bitonic sort. This can be done by using the hypercube formulation of bitonic sort presented in Section 9.2.2 (page 392). Each processor has ( p − 1) elements; thus, sorting the p( p − 1) elements takes ( p log2 p) time. At the end of the bitonic sort, each processor stores ( p − 1) elements. To select the ( p − 1) splitter elements, all the processors but the last one (that is, processor P p−1 ), select their ( p − 1)st element. These selected elements are the ( p − 1) splitters. The splitters are send to all the processors by performing an all-to-all broadcast operation (Section 4.2 (page 157)). This operation takes ( p) time. Assuming a uniform distribution of data, the parallel run time is: n  p log    n n TP =  + ( p log2 p) +  p log + (n/ p ) + O( p log p ) p The isoefficiency function of thispformulation is ( p2 log p). 31 Recall that the parallel runtime is b 2r ((log n ) + (n )) (9.2) r The optimal value of r is such that it minimizes Equation 9.2. This value can be found if we compute the exact parallel run time of the algorithm, using tc , tw , and ts , and then differentiate the obtained parallel run time. By setting the derivative to zero and solving for r we obtain the value that minimizes the run time. TP = 32 Let n be the number of elements to be sorted using p processors. We can use the radix sort algorithm described in Algorithm 9.8 (page 416) to perform the task, by using virtual processors. This is done by emulating n/ p processors on each physical processor. The run time of this formulation is slower, at most, by a factor of n/ p. However, the procedures parallel sum and prefix sum can be optimized to use the fact that each processor stores n/ p elements. Similarly, the communication step can be optimized taking into account that n/ p processors are assigned to each physical processor. These optimizations are described by Blelloch et al. [BLM+ 91] that is based on ideas by Cole and Vishkin [CV86].