Samplesort

Pseudocode

The following listing shows the above-mentioned three step algorithm as pseudocode and shows how the algorithm works in principle.^[5] In the following, A is the unsorted data, k is the oversampling factor, discussed later, and p is the number of splitters.

function sampleSort(A[1..n], k, p)
    // if average bucket size is below a threshold switch to e.g. quicksort
    if n / k < threshold then smallSort(A) 
    /* Step 1 */
    select S = [S1, ..., S(p−1)k] randomly from // select samples
    sort S // sort sample
    [s0, s1, ..., sp−1, sp] <- [-∞, Sk, S2k, ..., S(p−1)k, ∞] // select splitters
    /* Step 2 */
    for each a in A
        find j such that sj−1 < a <= sj
        place a in bucket bj
    /* Step 3 and concatenation */
    return concatenate(sampleSort(b1), ..., sampleSort(bk))

The pseudo code is different from the original Frazer and McKellar algorithm.^[3] In the pseudo code, samplesort is called recursively. Frazer and McKellar called samplesort just once and used quicksort in all following iterations.

Complexity

The complexity, given in Big O notation, for a parallelized implementation with ${\displaystyle p}$ processors:

Find the splitters.

${\displaystyle O\left({\frac {n}{p}}+\log(p)\right)}$

Send to buckets.

${\displaystyle O(p)}$ for reading all nodes

${\displaystyle O(\log(p))}$ for broadcasting

${\displaystyle O\left({\frac {n}{p}}\log(p)\right)}$ for binary search for all keys

${\displaystyle O\left({\frac {n}{p}}\right)}$ to send keys to bucket

Sort buckets.

${\displaystyle O\left(c\left({\frac {n}{p}}\right)\right)}$ where ${\displaystyle c(n)}$ is the complexity of the underlying sequential sorting method.^[1] Often ${\displaystyle c(n)=n\log(n)}$.

The number of comparisons, performed by this algorithm, approaches the information theoretical optimum ${\displaystyle \log _{2}(n!)}$ for big input sequences. In experiments, conducted by Frazer and McKellar, the algorithm needed 15% fewer comparisons than quicksort.

Oversampling

The oversampling ratio determines how many times more data elements to pull as samples, before determining the splitters. The goal is to get a good representation of the distribution of the data. If the data values are widely distributed, in that there are not many duplicate values, then a small sampling ratio is sufficient. In other cases where there are many duplicates in the distribution, a larger oversampling ratio will be necessary. In the ideal case, after step 2, each bucket contains ${\displaystyle n/p}$ elements. In this case, no bucket takes longer to sort than the others, because all buckets are of equal size.

After pulling ${\displaystyle k}$ times more samples than necessary, the samples are sorted. Thereafter, the splitters used as bucket boundaries are the samples at position ${\displaystyle k,2k,3k,\dots ,(p-1)k}$ of the sample sequence (together with ${\displaystyle -\infty }$ and ${\displaystyle \infty }$ as left and right boundaries for the left most and right most buckets respectively). This provides a better heuristic for good splitters than just selecting ${\displaystyle p}$ splitters randomly.

Bucket size estimate

With the resulting sample size, the expected bucket size and especially the probability of a bucket exceeding a certain size can be estimated. The following will show that for an oversampling factor of ${\displaystyle S\in \Theta \left({\dfrac {\log n}{\epsilon ^{2}}}\right)}$ the probability that no bucket has more than ${\displaystyle (1+\epsilon )\cdot {\dfrac {n}{p}}}$ elements is larger than ${\displaystyle 1-{\dfrac {1}{n}}}$.

To show this let ${\displaystyle \langle e_{1},\dots ,e_{n}\rangle }$ be the input as a sorted sequence. For a processor to get more than ${\displaystyle (1+\epsilon )\cdot n/p}$ elements, there has to exist a subsequence of the input of length ${\displaystyle (1+\epsilon )\cdot n/p}$, of which a maximum of S samples are picked. These cases constitute the probability ${\displaystyle P_{\text{fail}}}$. This can be represented as the random variable: $${\displaystyle X_{i}:={\begin{cases}1,&{\text{if }}s_{i}\in \left\langle e_{j},\dots ,e_{j}+(1+\epsilon )\cdot {\dfrac {n}{p}}\right\rangle \\0,&{\text{otherwise}}\end{cases}},X:=\sum _{i=0}^{S\cdot p-1}X_{i}}$$

For the expected value of ${\displaystyle X_{i}}$ holds: $${\displaystyle E(X_{i})=P(X_{i}=1)={\dfrac {1+\epsilon }{p}}}$$

This will be used to estimate ${\displaystyle P_{\text{fail}}}$: $${\displaystyle P(X<S)\approx P(X<(1-\epsilon ^{2})S)=P(X<(1-\epsilon )E(X))}$$

Using the Chernoff bound now, it can be shown: $${\displaystyle P_{\text{fail}}=n\cdot P(X<S)\leq n\cdot \exp \left({\dfrac {-\epsilon ^{2}\cdot S}{2}}\right)\leq n\cdot {\dfrac {1}{n^{2}}}{\text{ for }}S\geq {\dfrac {4}{\epsilon ^{2}}}\ln n}$$

Determining buckets

In a comparison based sorting algorithm the comparison operation is the most performance critical part. In Samplesort this corresponds to determining the bucket for each element. This needs ${\displaystyle \log k}$ time for each element.

Super Scalar Sample Sort uses a balanced search tree which is implicitly stored in an array t. The root is stored at 0, the left successor of ${\displaystyle t_{i}}$ is stored at ${\displaystyle t_{2i}}$ and the right successor is stored at ${\displaystyle t_{2i+1}}$. Given the search tree t, the algorithm calculates the bucket number j of element ${\displaystyle a_{i}}$ as follows (assuming ${\displaystyle a_{i}>t_{j}}$ evaluates to 1 if it is true and 0 otherwise):

j := 1
repeat log2(p) times
    j := 2j + (a > tj)
j := j − p + 1

Since the number of buckets k is known at compile time, this loop can be unrolled by the compiler. The comparison operation is implemented with predicated instructions. Thus, there occur no branch mispredictions, which would slow down the comparison operation significantly.

Partitioning

For an efficient partitioning of the elements, the algorithm needs to know the sizes of the buckets in advance. To partition the elements of the sequence and put them into the array, we need to know the size of the buckets in advance. A naive algorithm could count the number of elements of each bucket. Then the elements could be inserted to the other array at the right place. Using this, one has to determine the bucket for each elements twice (one time for counting the number of elements in a bucket, and one time for inserting them).

To avoid this doubling of comparisons, Super Scalar Sample Sort uses an additional array ${\displaystyle o}$ (called oracle) which assigns each index of the elements to a bucket. First, the algorithm determines the contents of ${\displaystyle o}$ by determining the bucket for each element and the bucket sizes, and then placing the elements into the bucket determined by ${\displaystyle o}$. The array ${\displaystyle o}$ also incurs cost in storage space, but as it only needs to store ${\displaystyle n\cdot \log k}$ bits, these cost are small compared to the space of the input array.

Local classification

In a first step, the input array is split up into ${\displaystyle p}$ stripes of blocks of equal size, one for each processor. Each processor additionally allocates ${\displaystyle k}$ buffers that are of equal size to the blocks, one for each bucket. Thereafter, each processor scans its stripe and moves the elements into the buffer of the according bucket. If a buffer is full, the buffer is written into the processors stripe, beginning at the front. There is always at least one buffer size of empty memory, because for a buffer to be written (i.e. buffer is full), at least a whole buffer size of elements more than elements written back had to be scanned. Thus, every full block contains elements of the same bucket. While scanning, the size of each bucket is kept track of.

Block permutation

Firstly, a prefix sum operation is performed that calculates the boundaries of the buckets. However, since only full blocks are moved in this phase, the boundaries are rounded up to a multiple of the block size and a single overflow buffer is allocated. Before starting the block permutation, some empty blocks might have to be moved to the end of its bucket. Thereafter, a write pointer ${\displaystyle w_{i}}$ is set to the start of the bucket ${\displaystyle b_{i}}$ subarray for each bucket and a read pointer ${\displaystyle r_{i}}$ is set to the last non empty block in the bucket ${\displaystyle b_{i}}$ subarray for each bucket.

To limit work contention, each processor is assigned a different primary bucket ${\displaystyle b_{prim}}$ and two swap buffers that can each hold a block. In each step, if both swap buffers are empty, the processor decrements the read pointer ${\displaystyle r_{prim}}$ of its primary bucket and reads the block at ${\displaystyle r_{prim-1}}$ and places it in one of its swap buffers. After determining the destination bucket ${\displaystyle b_{dest}}$ of the block by classifying the first element of the block, it increases the write pointer ${\displaystyle w_{dest}}$, reads the block at ${\displaystyle w_{dest-1}}$ into the other swap buffer and writes the block into its destination bucket. If ${\displaystyle w_{dest}>r_{dest}}$, the swap buffers are empty again. Otherwise the block remaining in the swap buffers has to be inserted into its destination bucket.

If all blocks in the subarray of the primary bucket of a processor are in the correct bucket, the next bucket is chosen as the primary bucket. If a processor chose all buckets as primary bucket once, the processor is finished.

Cleanup

Since only whole blocks were moved in the block permutation phase, some elements might still be incorrectly placed around the bucket boundaries. Since there has to be enough space in the array for each element, those incorrectly placed elements can be moved to empty spaces from left to right, lastly considering the overflow buffer.