In-place Parallel Super Scalar Samplesort (IPS$^4$o)

TL;DR

The paper introduces IPS$^4$o, an in-place, parallel, cache-efficient sorting algorithm that outperforms existing methods in speed and scalability on multi-core systems, combining theoretical innovations with practical performance.

Contribution

It presents novel in-place distribution-based sorting techniques that eliminate recursion stacks and improve parallel and cache-efficient performance.

Findings

IPS$^4$o scales well on multi-core machines.

Outperforms in-place competitors by up to 3x.

Faster than BlockQuicksort by up to 1.5x sequentially.

Abstract

We present a sorting algorithm that works in-place, executes in parallel, is cache-efficient, avoids branch-mispredictions, and performs work O(n log n) for arbitrary inputs with high probability. The main algorithmic contributions are new ways to make distribution-based algorithms in-place: On the practical side, by using coarse-grained block-based permutations, and on the theoretical side, we show how to eliminate the recursion stack. Extensive experiments show that our algorithm IPS$^4$o scales well on a variety of multi-core machines. We outperform our closest in-place competitor by a factor of up to 3. Even as a sequential algorithm, we are up to 1.5 times faster than the closest sequential competitor, BlockQuicksort.

Tables

Table 1. Table 1 . The first three rows show the speedups of IS 4 o relative to the fastest sequential in-place and non-in-place competitor on different input types executed on machine Intel2S , Intel4S , and AMD1S for n = 2 32 𝑛 superscript 2 32 n=2^{32} . The last rows show the speedups of IPS 4 o relative to the fastest parallel in-place and non-in-place competitor on different input types executed on different machine instances for n = 2 32 𝑛 superscript 2 32 n=2^{32} . Measurements in cells labeled with OOM ran out of memory.

Machine	Algo	Competitor	Uniform	Exponential	Almost	RootDup	TwoDup
Intel2S	IS4o	both	1.14	1.23	0.59	0.97	1.17
Intel4S	IS4o	both	1.21	1.54	0.77	1.65	1.44
AMD1S	IS4o	both	1.57	2.02	0.65	1.37	1.17
Intel2S	IPS4o	in-place	2.54	3.43	1.88	2.73	3.02
		non-in-place	2.13	1.79	1.29	1.19	1.86
Intel4S	IPS4o	in-place	3.52	4.35	3.62	3.19	2.89
		non-in-place	1.75	1.69	1.84	1.15	1.19
AMD1S	IPS4o	in-place	1.57	3.18	1.81	2.37	2.02
		non-in-place	OOM	OOM	OOM	OOM	OOM

Full text

1. Introduction

Sorting an array $A[1..n]$ of $n$ elements according to a total ordering of their keys is a fundamental subroutine used in many applications. Sorting is used for index construction, for bringing similar elements together, or for processing data in a “clever” order. Indeed, often sorting is the most expensive part of a program. Consequently, a huge amount of research on sorting has been done. In particular, algorithm engineering has studied how to make sorting practically fast in presence of complex features of modern hardware like multi-core (e.g., [30, 29, 4, 28]), instruction parallelism (e.g., [27]), branch prediction (e.g., [27, 18, 17, 9]), caches (e.g., [27, 6, 10, 4]), or virtual memory (e.g., [24, 16]). In contrast, the sorting algorithms used in the standard libraries of programming languages like Java or C++ still use variants of quicksort – an algorithm that is more than 50 years old. A reason seems to be that you have to outperform quicksort in every respect in order to replace it. This is less easy than it sounds since quicksort is a pretty good algorithm – it needs $\mathcal{O}\!\left(n\log n\right)$ expected work, it can be parallelized [30, 29], it can be implemented to avoid branch mispredictions [9], and it is reasonably cache-efficient. Perhaps most importantly, quicksort works (almost) in-place111In algorithm theory, an algorithm works in-place if it uses only constant space in addition to its input. We use the term strictly in-place for this case. In algorithm engineering, one is sometimes satisfied if the additional space is sublinear in the input size. We adopt this convention but use the term almost in-place when we want to make clear what we mean. Quicksort needs logarithmic additional space. which is of crucial importance for very large inputs. This feature rules out many contenders. Further algorithms are eliminated by the requirement to work for arbitrary data types and input distributions. This makes integer sorting algorithms like radix sort (e.g., [20]) or using specialized hardware (e.g., GPUs or SIMD instructions) less attractive, since these algorithms cannot be used in a reusable library where they have to work for arbitrary data types. Another portability issue is that the algorithm should use no code specific to the processor architecture or the operating system like non-temporal writes or overallocation of virtual memory (e.g. [26]). One aspect of making an algorithm in-place is that such “tricks” are not needed. Hence, this paper focuses on portable comparison-based algorithms and also considers how the algorithms can be made robust for arbitrary inputs, e.g., with a large number of repeated keys.

The main contribution of this paper is to propose a new algorithm – In-place Parallel Super Scalar Samplesort (IPS4o)222The Latin word “ipso” means “by itself”, referring to the in-place feature of IPS4o. – that combines enough advantages to become an attractive replacement of quicksort. Our starting point is super scalar samplesort (s3-sort) [27] which already provides a very good sequential non-in-place algorithm that is cache-efficient, allows considerable instruction parallelism, and avoids branch mispredictions. s3-sort is a variant of samplesort, which in turn is a generalization of quicksort to multiple pivots. The main operation is distributing elements of an input sequence to $k$ output buckets of about equal size. We parallelize this algorithm using $t$ threads and make it more robust by taking advantage of inputs with many identical keys. Our main innovation is to make the algorithm in-place. The first phase of IPS4o distributes the elements to $k$ buffer blocks. When a buffer becomes full, it is emptied into a block of the input array that has already been distributed. Subsequently, the memory blocks are permuted into the globally correct order. A cleanup step handles empty blocks and half-filled buffer blocks. The distribution phase is parallelized by assigning disjoint pieces of the input array to different threads. The block permutation phase is parallelized using atomic fetch-and-add operations for each block move. Once subproblems are small enough, they can be solved independently in parallel.

After discussing related work in section 2 and introducing basic tools in section 3, we describe our new algorithm IPS4o in section 4. section 5 makes an experimental evaluation. An overall discussion and possible future work is given in section 6. The appendix gives further experimental data and proofs.

2. Related Work

Variants of Hoare’s quicksort [14, 23] are generally considered some of the most efficient general purpose sorting algorithms. Quicksort works by selecting a pivot element and partitioning the array such that all elements smaller than the pivot are in the left part and all elements larger than the pivot are in the right part. The subproblems are solved recursively. A variant of quicksort (with a fallback to heapsort to avoid worst case scenarios) is currently used in the C++ standard library of GCC [23]. Some variants of quicksort use two or three pivots [31, 21] and achieve improvements of around 20% in running time over the single-pivot case. Dual-pivot quicksort [31] is the default sorting routine in Oracle Java 7 and 8. The basic principle of quicksort remains, but elements are partitioned into three or four subproblems instead of two. Increasing the number of subproblems (from now on called buckets) even further leads to samplesort [5, 4]. Unlike single- and dual-pivot quicksort, samplesort is usually not in-place, but it is well-suited for parallelization and more cache-efficient.

Super scalar samplesort [27] (s3-sort) improves on samplesort by avoiding inherently hard-to-predict conditional branches linked to element comparisons. Branch mispredictions are very expensive because they disrupt the pipelined and instruction-parallel operation of modern processors. Traditional quicksort variants suffer massively from branch mispredictions [18]. By replacing conditional branches with conditionally executed machine instructions, branch mispredictions can be largely avoided. This is done automatically by modern compilers if only a few instructions depend on a condition. As a result, s3-sort is up to two times faster than quicksort (std::sort), at the cost of $\mathcal{O}\!\left(n\right)$ additional space. BlockQuicksort [9] applies similar ideas to single-pivot quicksort, resulting in a very fast in-place sorting algorithm.

Super scalar samplesort has also been adapted for efficient parallel string sorting [3]. Our implementation is influenced by that work with respect to parallelization and handling equal keys. Moreover, we were also influenced by an implementation of s3-sort written by Lorenz Hübschle-Schneider. A prototypical implementation of sequential non-blocked in-place s3-sort in a student project by our student Florian Weber motivated us to develop IPS4o.

The best practical comparison-based multi-core sorting algorithms we have found are based on multi-way mergesort [29] and samplesort [28], respectively. The former algorithm is used in the parallel mode of the C++ standard library of GCC. Parallel in-place algorithms are based on quicksort so far. Intel’s Thread Building Blocks library [25] contains a variant that uses only sequential partitioning. The MCSTL library [29] contains two implementations of the more scalable parallel quicksort by Tsigas and Zhang [30].

There is a considerable amount of work by the theory community on (strictly) in-place sorting (e.g., [10, 11]). However, there are few – mostly negative – results on transferring these results into practice. Katajainen and Teuhola [19] report that in-place mergesort is slower than heapsort, which is quite slow for big inputs due to its cache-inefficiency. Chen [7] reports that in-place merging takes about six times longer than non-in-place merging. There is previous work on (almost) in-place multi-way merging or data distribution. However, few of these papers seem to address parallelism. There are also other problems. For example, the multi-way merger in [13] needs to allocate very large blocks to become efficient. In contrast, the block size of IPS4o does not depend on the input size. In-place data distribution, e.g., for radix sort [8], is often done element by element. Using this for samplesort would require doing the expensive element classification twice and would also make parallelization difficult.

3. Preliminaries

(Super Scalar) Samplesort.

Samplesort [12] can be viewed as a generalization of quicksort which uses multiple pivots to split the input into $k$ buckets of about equal size. A robust way for determining the pivots is to sort $\alpha k-1$ randomly sampled input elements. The pivots $s_{1}$,…$s_{k-1}$ are then picked equidistantly from the sorted sample. Element $e$ goes to bucket $b_{i}$ if $s_{i-1}\leq e<s_{i}$ (with $s_{0}=-\infty$ and $s_{k}=\infty$). The main contribution of s3-sort [27] is to eliminate branch mispredictions for element classification. Assuming $k$ is a power of two, the pivots are stored in an array $a$ representing a complete binary search tree: $a_{1}=s_{k/2}$, $a_{2}=s_{k/4}$, $a_{3}=s_{3k/4},\ldots$ More generally, the left successor of $a_{i}$ is $a_{2i}$ and its right successor is $a_{2i+1}$. Thus, navigating this tree is possible by performing a conditional instruction for incrementing an array index. We adopt (and refine) this approach to element classification but change the organization of buckets in order to make the algorithm in-place.

4. In-Place Parallel Super Scalar Samplesort (IPS4o)

IPS4o is based on the ideas of s3-sort. It is a recursive algorithm, where each step divides the input into $k$ buckets, such that each element of bucket $b_{i}$ is smaller than all elements of $b_{i+1}$. As long as problems with at least $\beta\frac{n}{t}$ elements exist, we partition those problems one after another with $t$ threads in parallel. Here, $\beta$ is a tuning parameter. Then we assign remaining problems in a balanced way to threads, which sort them sequentially.

The partitioning consists of four phases. Sampling determines the bucket boundaries. Local classification groups the input into blocks such that all elements in each block belong to the same bucket. Block permutation brings the blocks into the globally correct order. Finally, we perform some cleanup around the bucket boundaries. The following sections will explain each of these phases in more detail.

Sampling.

The sampling phase is similar to the sampling in s3-sort. The main difference is that we swap the sample to the front of the input array to keep the in-place property even if the oversampling factor $\alpha$ depends on $n$.

4.1. Local Classification

The input array $A$ is viewed as an array of blocks each containing $b$ elements (except possibly for the last one). For parallel processing, we divide the blocks of $A$ into $t$ stripes of equal size – one for each thread. Each thread works with a local array of $k$ buffer blocks – one for each bucket. A thread then scans its stripe. Using the search tree created in the previous phase, each element in the stripe is classified into one of the $k$ buckets, then moved into the corresponding local buffer block. If this buffer is already full, it is first written back into the local stripe, starting at the front. It is clear that there is enough space to write $b$ elements into the local stripe, since at least $b$ more elements have been scanned from the stripe than have been written back – otherwise no full buffer could exist.

In this way, each thread creates blocks of $b$ elements belonging to the same bucket. Figure 1 shows a typical situation during this phase. To achieve the in-place property, we do not track which bucket each block belongs to. However, we do keep count of how many elements are classified into each bucket, since we need this information in the following phases. This information can be obtained almost for free as a side effect of maintaining the buffer blocks. Figure 2 depicts the input array after local classification. Each stripe contains a number of full blocks, followed by a number of empty blocks. The remaining elements are still contained in the buffer blocks.

4.2. Block Permutation

In this phase, the blocks in the input array will be rearranged such that they appear in the correct order. From the previous phase we know, for each stripe, how many elements belong to each bucket. We perform a prefix sum operation to compute the exact boundaries of the buckets in the input array. In general, these will not coincide with the block boundaries. For the purposes of this phase, we will ignore this: We mark the beginning of each bucket $b_{i}$ with a delimiter pointer $d_{i}$, rounded up to the next block. We similarly mark the end of the last bucket $b_{k}$ with a delimiter pointer $d_{k+1}$. Adjusting the boundaries may cause a bucket to “lose” up to $b-1$ elements; this doesn’t affect us, since this phase only deals with full blocks, and any elements not constituting a full block remain in the buffers. Additionally, if the input size is not a multiple of $b$, some of the $d_{i}$s may end up outside the bounds of $A$. To avoid overflows, we allocate a single empty overflow block which the algorithm will use instead of writing to the final (partial) block.

For each $b_{i}$, a write pointer $w_{i}$ and a read pointer $r_{i}$ is introduced; these will be set such that all unprocessed blocks, i.e., blocks that still need to be moved into the correct bucket, are found between $w_{i}$ and $r_{i}$. During the block permutation, we maintain the following invariant for each bucket $b_{i}$, visualized in Figure 3:

• •

Blocks to the left of $w_{i}$ (exclusive) are correctly placed, i.e., contain only elements belonging to $b_{i}$.

• •

Blocks between $w_{i}$ and $r_{i}$ (inclusive) are unprocessed, i.e., may need to be moved.

• •

Blocks to the right of $\max(w_{i},r_{i}+1)$ (inclusive) are empty.

In other words, each bucket follows the pattern of correct blocks followed by unprocessed blocks followed by empty blocks, with $w_{i}$ and $r_{i}$ determining the boundaries. In the parallel case, we may need to establish this invariant by moving some empty blocks to the end of a bucket (see Appendix A for details); in the sequential algorithm, the result of the classification phase already has this pattern. The read pointers $r_{i}$ are then set to the first non-empty block in each bucket, or $d_{i}-1$ if there are none.

We are now ready to start the block permutation. Each thread maintains two local swap buffers. We define a primary bucket $b_{p}$ for each thread; whenever both its buffers are empty, a thread tries to read an unprocessed block from its primary bucket. To do so, it decrements the read pointer $r_{p}$ (atomically) and reads the block it pointed to into one of its swap buffers. If $b_{p}$ contains no more unprocessed blocks (i.e., $r_{p}<w_{p}$), it switches its primary bucket to the next bucket (cyclically). If it completes a whole cycle and arrives back at its initial primary bucket, there are no more unprocessed blocks and this phase ends. The starting points for the threads are distributed across that cycle to reduce contention.

Once it has a block, each thread classifies the first element of that block to find its destination bucket $b_{\text{dest}}$. There are now two possible cases, visualized in Figure 4:

• •

As long as $w_{\text{dest}}\leq r_{\text{dest}}$, write pointer $w_{\text{dest}}$ still points to an unprocessed block in bucket $b_{\text{dest}}$. In this case, the thread increases $w_{\text{dest}}$, reads the unprocessed block into its empty swap buffer, and writes the other one into its place.

• •

If $w_{\text{dest}}>r_{\text{dest}}$, no unprocessed block remains in bucket $b_{\text{dest}}$ but $w_{\text{dest}}$ now points to an empty block. In this case, the thread increases $w_{\text{dest}}$, writes its swap buffer to the empty block and then reads a new unprocessed block from its primary bucket.

We repeat these steps until all blocks are processed. We can skip unprocessed blocks which are already correctly placed: We simply classify blocks before reading them into a swap buffer, and skip as needed. We omitted this from the above description for the sake of clarity. In some cases, this reduces the number of block moves significantly.

It is possible that one thread wants to write to a block that another thread is currently reading from (when the reading thread has just decremented the read pointer, but has not yet finished reading the block into its swap buffer). To avoid data races, we keep track of how many threads are reading from each bucket. Threads are only allowed to write to empty blocks if no other threads are currently reading from the bucket in question, otherwise they wait. Note that this situation occurs at most once for each bucket, namely when $w_{\text{dest}}$ and $r_{\text{dest}}$ cross each other. In addition, we store each $w_{i}$ and $r_{i}$ in a single 128-bit word which we read and modify atomically. This ensures a consistent view of both pointers for all threads.

4.3. Cleanup

After the block permutation, some elements may still be in incorrect positions. This is due to the fact that we only moved blocks, which may span bucket boundaries. We call the partial block at the beginning of a bucket its head and the partial block at its end its tail.

We assign consecutive buckets evenly to threads; if $t>k$, some threads will not receive any buckets, but those that do only need to process a single bucket each. Each thread reads the head of the first bucket of the next thread into one of its swap buffers. Then, each thread processes its buckets from left to right, moving incorrectly placed elements into empty array entries. The incorrectly placed elements of bucket $b_{i}$ consist of the elements in the head of $b_{i+1}$ (or the swap buffer, for the last bucket), the partially filled buffers from the local classification phase (of all threads), and, for the corresponding bucket, the overflow buffer. Empty array entries consist of the head of $b_{i}$ and any (empty) blocks to the right of $w_{i}$ (inclusive). Although the concept is relatively straightforward, the implementation is somewhat involved, due to the many parts that have to be brought together. Figure 5 shows an example of the steps performed during this phase. Afterwards, all elements are back in the input array and correctly partitioned, ready for recursion.

4.4. The Case of Many Identical Keys

Having inputs with many identical keys can be a problem for samplesort, since this might move large fractions of the keys through many levels of recursion. We turn such inputs into easy instances by introducing separate buckets for elements identical to pivots (keys occurring more then $\frac{n}{k}$ times are likely to become pivots). Finding out whether an element has to go into an equality bucket (and which one) can be implemented using a single additional comparison [3] and, once more, without a conditional branch. Equality buckets can be skipped during recursion and thus are not a load balancing problem.

4.5. Analysis

Algorithm IPS4o inherits from s3-sort that it has virtually no branch mispredictions (this includes the comparisons for placing elements into equality buckets discussed in subsection 4.4). More interesting is the parallel complexity. Here, the main issue is the number of accesses to main memory. We analyze this aspect in the parallel external memory (PEM) model [1], where each of the $t$ threads has a private cache of size $M$ and access to main memory happens in blocks of size $B$. In Appendix B, we prove:

Theorem 4.1**.**

Assuming $b=\Theta\!\left(tB\right)$ (buffer block size), $M=\Omega(ktB)$, $n_{0}=\mathcal{O}\!\left(M\right)$ (base case size), $\alpha\in\Omega(\log t)\cap\mathcal{O}\!\left(t\right)$ (oversampling factor), and $n=\Omega\!\left(\max(k,t)t^{2}B\right)$, IPS4o has an I/O-complexity of $\mathcal{O}\!\left(\frac{n}{tB}\log_{k}\frac{n}{n_{0}}\right)$ block transfers with high probability.

Basically, Theorem 4.1 tells us that IPS4o is asymptotically I/O efficient if certain rather steep assumptions on cache size and input size hold. In particular, the blocks need to have size $b=\Theta\!\left(tB\right)$ in order to amortize contention on shared block pointers. Lifting those could be an interesting theoretical question and we would have to see how absence of branch mispredictions and the in-place property can be combined with previous techniques [1, 4]. However, it is likely that the constant factors involved are much larger than for our simple implementation. Thus, the constant factors will be the main issue in bringing theory and practice further together. To throw some light on this aspect, let us compare the constant factors in I/O-volume (i.e., data flow between cache and main memory) for the sequential algorithms IS4o (IPS4o with $t=1$) and s3-sort. To simplify the discussion, we assume a single level of recursion, $k=256$ and $8$-byte elements. In Appendix B, we show that IS4o needs about $48n$ bytes of I/O volume, whereas s3-sort needs (more than) $86n$ – almost twice that of IS4o. This is surprising since on first glance, the partitioning algorithm of IS4o writes the data twice, whereas s3-sort does this only once. However, this is more than offset by “hidden” overheads of s3-sort like memory management, allocation misses, and associativity misses.

Finally, we consider the memory overhead of IPS4o. In Appendix B, we show:

Theorem 4.2**.**

IPS4o requires additional space $\mathcal{O}\!\left(kbt+\log_{k}\frac{n}{n_{0}}\right)$.

In practice, the term $\mathcal{O}\!\left(kbt\right)$ (mostly for the distribution buffers) will dominate. However, for a strictly in-place algorithm in the sense of algorithm theory, we need to get rid of the $\mathcal{O}\!\left(\log n\right)$ term which depends on the input size. We discuss this separately in subsection 4.6.

4.6. From Almost In-Place to Strictly In-Place

We now explain how the space consumption of IPS4o can be made independent of $n$ in a rather simple way. We can restrict ourselves to the sequential case, since only $\mathcal{O}\!\left(\log_{k}t\right)$ levels of parallel recursion are needed to arrive at subproblems that are solved sequentially. We require the partitioning operation to mark the beginning of each bucket by storing the largest element of a bucket in its first entry. By searching the next larger element, we can then find the end of the bucket. Note that this is possible in time logarithmic in the bucket size using exponential/binary search. We assume that the corresponding function $\mathit{searchNextLargest}$ returns $n+1$ if no larger elements exists – this happens for the last bucket. The following pseudocode uses this approach to emulate recursion in constant space for sequential IS4o.

[TABLE]

4.7. Implementation Details

The strategy for handling identical keys described in subsection 4.4 is enabled conditionally: After the splitters have been selected from the initial sample, we check for and remove duplicates. Equality buckets are only used if there were duplicate splitters.

For buckets under a certain base case size $n_{0}$, we stop the recursion and fall back on insertion sort. Additionally, we use an adaptive number of buckets on the last two levels of the recursion, such that the expected size of the final buckets remains reasonable. For example, instead of performing two $256$-way partitioning steps to get $2^{16}$ buckets of 2 elements, we might perform two $64$-way partitioning steps to get $2^{12}$ buckets of about $32$ elements. Furthermore, on the last level, we perform the base case sorting immediately after the bucket has been completely filled in the cleanup phase, before processing the other buckets. This is more cache-friendly, as it eliminates the need for another pass over the data.

IPS4o has several parameters that can be used for tuning and adaptation. We performed our experiments using (up to) $k=256$ buckets, an oversampling factor of $\alpha=0.2\log n$, an overpartitioning factor of $\beta=1$, a base case size of $n_{0}=16$ elements, and a block size of about 2 KiB, or $b=\max\left(1,2^{\lfloor 11-\log_{2}s\rfloor}\right)$ elements, where $s$ is the size of an element in bytes. In the sequential case, we avoid the use of atomic operations on pointers. All algorithms are written in C++ and compiled with version 6.2.0 of the GNU compiler collection, using the optimization flags “-march=native -O3”. For parallelization, we employ OpenMP. Our implementation can be found at https://github.com/SaschaWitt/ips4o.

5. Experimental Results

We present the results of our in-place parallel sorting algorithm IPS4o. We compare the results of IPS4o with its in-place competitors, parallel sort from the Intel® TBB library [25] (TBB), parallel unbalanced quicksort from the GCC STL library (MCSTLubq), and parallel balanced quicksort from the GCC STL library (MCSTLbq). We also give results on the parallel non-in-place sorting algorithms, parallel samplesort from the problem based benchmark suite [28] (PBBS) and parallel multiway mergesort from the GCC STL library [29] (MCSTLmwm). We also ran sequential experiments and present the results of IS4o, the sequential implementation of IPS4o. We compare the results of IS4o with its sequential competitors, a recent implementation [15] of non-in-place Super Scalar Samplesort [27] (s3-sort) optimized for modern hardware, BlockQuicksort [9] (BlockQ), Dual-Pivot Quicksort [31] (DualPivot), and introsort from the GCC STL library (std-sort).

We ran benchmarks with nine input distributions: Uniformly distributed (Uniform), exponentially distributed (Exponential), and almost sorted (AlmostSorted), proposed by Shun et. al. [28]; RootDup, TwoDup, and EightDup from Edelkamp et. al. [9]; and Sorted (sorted Uniform input), ReverseSorted, and Ones (just ones). The input distribution RootDup sets $A[i]=i\mod\lfloor\sqrt{n}\rfloor$, TwoDup sets $A[i]=i^{2}+\frac{n}{2}\mod n$, and EightDup sets $A[i]=i^{8}+\frac{n}{2}\mod n$. We ran benchmarks with $64$-bit floating point elements and Pair, Quartet, and 100Bytes data types. Pair (Quartet) consists of one (three) $64$-bit floating point elements as key and one $64$-bit floating point element of associated information. 100Bytes consists of $10$ bytes as key and $90$ bytes of associated information. Quartet and 100Bytes are compared lexicographically. For $n<2^{30}$, we perform each measurement $15$ times and for $n\geq 2^{30}$, we perform each measurement twice. Unless stated otherwise, we report the average over all runs and use $64$-bit floating point elements.

We ran our experiments on machines with one AMD Ryzen $+1800$ 8-core processor (AMD1S), two Intel Xeon E5-2683 v4 16-core processors (Intel2S), and four Intel Xeon E5-4640 8-core processors (Intel4S). Intel2S and Intel4S are equipped with $512$ GiB of memory, AMD1S is equipped with $32$ GiB of memory. We use the taskset tool to set the CPU affinity for speedup benchmarks. We tested all parallel algorithms on Uniform input with and without hyper-threading. Hyper-threading did not slow down any algorithm. Thus, we give results of all algorithms with hyper-threading. Overall, we executed more than $12\,000$ combinations of different algorithms, input distributions and sizes, data types and machines. We now present a selection of our measurements and discuss our results. For the remaining (detailed) running time and hardware counter measurements, we refer to Appendix C.

Sequential Algorithms.

Figure 7 shows the running times of sequential algorithms on Uniform input executed on machine Intel2S. We see that IS4o is faster than its closest competitor, BlockQ, by a factor of $1.14$ for $n=2^{32}$. On machine Intel4S (AMD1S), IS4o outperforms BlockQ even by a factor of $1.22$ ($1.57$). DualPivot and std-sort, which do not avoid branch mispredictions, are at least a factor of $1.86$ slower than IS4o for $n=2^{32}$. The number of branch mispredictions of these algorithms for this input size is about $10$ times larger than that of IS4o. s3-sort is the slowest sequential sorting algorithm avoiding branch mispredictions and has fluctuations in running time for varying input sizes. Due to the initial overhead, IS4o is slower than BlockQ for $n\leq 2^{15}$.

As expected, the running times for inputs with a moderate number of different keys (TwoDup) are similar to the running times for Uniform. When the number of different keys decreases (Exponential, EightDup, and RootDup in decreasing order), IS4o becomes even faster by a factor of up to two on all machines. The running times of the competitors also decrease. However, only DualPivot on Intel2S with RootDup distributed input comes close for $n\geq 2^{28}$. Only input Ones and (almost) sorted input are hard for IS4o; for example, DualPivot outperforms IS4o on AlmostSorted input by a factor of $1.70$ for $n=2^{32}$ (Intel2S). For detailed measurements see also Figure 16-19 in Appendix C.

Parallel Algorithms.

Figure 8 (a-c) presents experiments of parallel algorithms on different machines for Uniform input. We see that IPS4o outperforms its closest competitors, e.g., for $n=2^{32}$ on Intel2S (AMD1S) by a factor of $2.13$ ($1.75$), and all but TBB and IPS4o fail to sort this input size on AMD1S due to memory limitations. For $n\geq 2^{26}$, IPS4o outperforms its closest non-in-place competitors on Intel2S (AMD1S) on average by a factor of $2.26$ ($1.69$) and its closest in-place competitors by a factor of $2.78$ ($1.98$). For the same input sizes, IPS4o outperforms its closest competitors on Intel4S in average just by a factor of $1.41$. We believe that the small difference in running time between IPS4o and its competitors on Intel4S is caused by two factors: The slower memory modules (DDR4 vs. DDR3), and the long load delays due to a ring interconnect between four sockets.

In Figure 8 (d-e), we present running times of parallel algorithms on input distributions with duplicates (TwoDup and RootDup) on machine Intel2S. For $n\geq 2^{26}$ and a moderate number of different keys (TwoDup), IPS4o still outperforms its in-place competitors on average by a factor of at least $2.88$ and its non-in-place competitors on average by a factor of at least $1.91$. Experiments have shown that the running times on EightDup and Exponential are similar to the running times on TwoDup. We also see that the non-in-place algorithms become almost as fast as IPS4o if we sort inputs which contain few different keys (RootDup). However, IPS4o still outperforms its in-place competitors by a factor of at least $3.43$ on this input for $n\geq 2^{20}$. Figure 8 (f) depicts the running times of parallel algorithms on AlmostSorted distributions on Intel2S. On AlmostSorted and ReverseSorted, the fastest non-in-place algorithm, PBBS, performs similarly to IPS4o for large input sizes. Only on Sorted and Ones, IPS4o is outperformed by TBB, an in-place competitor. This is because TBB detects these pre-sorted input distributions and terminates immediately. Further benchmarks on machines Intel4S and AMD1S show that IPS4o also outperforms its non-in-place competitors on any machine and that IPS4o is much faster than its in-place competitors except in the case of Sorted and Ones inputs. For detailed measurements see also Figure 9-11 in Appendix C.

In Figure 8 (g-h), we give running times of Pair and 100Bytes data types on machine Intel2S with uniformly distributed keys. We see that IPS4o outperforms its competitors, e.g., by a factor of $1.33$ (non-in-place competitor) and by a factor of $2.67$ (its in-place competitor) for $2^{29}$ 100Bytes elements. Further benchmarks on machines Intel4S and AMD1S show similar running times. For detailed measurements see also Figure 12-14 in Appendix C.

Figure 7 depicts the speedup of parallel algorithms executed on different numbers of cores relative to our sequential implementation IS4o on Intel2S, sorting Uniform input ($n=2^{30}$). We see that IPS4o outperforms its competitors on any number of cores. IPS4o outperforms IS4o on 32 cores by a factor of $28.71$, whereas its fastest non-in-place competitor, PBBS, outperforms IS4o just by a factor of $14.54$. The in-place algorithms, MCSTLubq and MCSTLbq, scale similarly to PBBS up to $16$ cores but begin lagging behind for larger numbers of cores. Further measurements show that IPS4o scales similarly on AMD1S. On Intel4S, IPS4o scales well on the first processor. However, as the input data is stored in the memory of the first processor, adding the second, third and fourth processors speeds up IPS4o by an additional factor of only 1.45; again caused by the slower memory modules (DDR4 vs. DDR3) and the long load delays due to a ring interconnect between four sockets. For detailed measurements see also Figure 15 in Appendix C.

6. Conclusion and Future Work

In-place super scalar samplesort (IPS4o) is among the fastest comparison-based sorting algorithms both sequentially and on multi-core machines. The algorithm can also be used for data distribution and local sorting in distributed memory parallel algorithms (e.g., [2]). Somewhat surprisingly, there is even an advantage over non-in-place algorithms because IPS4o saves on overhead for memory allocation, associativity misses and write allocate misses. Compared to previous parallel in-place algorithms, improvements by more than a factor of two are possible. The main case where IPS4o is slower than the best competitors (s3-sort and BlockQuicksort) is for sequentially sorting large objects (Quartet and 100Bytes, see Appendix C) because IPS4o moves elements twice in one distribution step. In this case, the overhead for the oracle information of s3-sort is small and we could try an almost-in-place variant of s3-sort with element-wise in-place permutation.

Several improvements of IPS4o can be considered. Besides careful adaptation of parameters like $k$, $b$, $\alpha$, and the choice of base case algorithm, one would like to avoid contention on the bucket pointers in the block permutation phase when $t$ is large. Perhaps the most important improvement would be to make IPS4o aware of non-uniform memory access costs (NUMA) depending on the memory module holding a particular piece of data. This can be done by preferably assigning pieces of the input array to “close-by” cores both for local classification and when switching to sequential sorting. In situations with little NUMA effects, we could ensure that our data blocks correspond to pages of the virtual memory. Then, one can replace block permutation with relabelling the virtual memory addresses of the corresponding pages.

Coming back to the original motivation for an alternative to quicksort variants in standard libraries, we see IPS4o as an interesting candidate. The main remaining issue is the code complexity. When code size matters (e.g., as indicated by a compiler flag like -Os), quicksort should still be used. Formal verification of the correctness of the implementation might help to increase trust in the remaining cases.

Acknowledgements.

We would like to thank the authors of [28, 9] for sharing their code for evaluation. Timo Bingmann and Lorenz Hübschle-Schneider [15] kindly provided code that was used as a starting point for our implementation.

Appendix A Details of the Algorithm

Empty block movement

An important observation is that within each stripe, all full blocks are at the beginning, followed by all empty blocks. This arrangement fulfils the invariant used during permutation, which is why there is no need to move empty blocks in the sequential algorithm. It also means that in the parallel algorithm, only the buckets crossing a stripe boundary need to be fixed.

To do so, each thread finds the bucket that starts before the end of its stripe, but ends after it. It then finds the stripe in which that bucket ends (which will be the following stripe in most cases) and moves the last full block in the bucket into the first empty block in the bucket. It continues to do this until either all empty blocks in its stripe are filled or all full blocks in the bucket have been moved.

In rare cases, very large buckets exist that cross multiple stripes. In this case, each thread will first count how many blocks in the preceding stripes need to be filled. It will then skip that many blocks at the end of the bucket before starting to fill its own empty blocks.

Appendix B Details of the Analysis

Proof B.1** (Proof of Theorem 4.1).**

It can be shown using Chernoff bounds that an oversampling ratio of $\alpha=\Omega\!\left(\log kt\right)$ is sufficient to produce (non-equality) buckets of size $\mathcal{O}\!\left(\frac{N}{k}\right)$ with high probability for subproblems of size $N$. Hence, $\mathcal{O}\!\left(\log_{k}\frac{n}{n_{0}}\right)$ levels of recursion suffice with high probability. On the other hand, for $\alpha=\mathcal{O}\!\left(t\right)$, even sequentially processing the sample does not become a bottleneck.

During the block distribution phase, each thread reads $\mathcal{O}\!\left(\frac{n}{tb}\right)$ logical data blocks, writes them to the buffers in its private cache, and eventually moves them back to main memory.

The same asymptotic cost occurs for moving blocks during block permutation. Each thread performs $\mathcal{O}\!\left(\frac{n}{tb}\right)$ successful acquisitions of the next block in a bucket. Charging $\mathcal{O}\!\left(t\right)$ I/Os for this accounts for possible contention with other threads. Overall, we get cost $\mathcal{O}\!\left(\frac{n}{b}\right)=\mathcal{O}\!\left(\frac{n}{tB}\right)$. Similarly, there are $k$ unsuccessful acquisitions before termination is determined, for which we charge an overall cost of $\mathcal{O}\!\left(kt\right)$ I/Os. Since we assume $n=\Omega\!\left(kt^{2}B\right)$, we have $k=\mathcal{O}\!\left(\frac{n}{t^{2}B}\right)$ and hence $\mathcal{O}\!\left(kt\right)=\mathcal{O}\!\left(\frac{n}{tB}\right)$.

In the cleanup phase, we consider a case distinction with respect to $k$ and $t$. If $k\leq t$, then each thread processes at most one bucket and it has to move elements from $t+2$ distribution buffers and bucket boundaries. This amounts to a cost of $\mathcal{O}\!\left(\frac{tb}{B}\right)=\mathcal{O}\!\left(t^{2}\right)$ I/Os. Since $n=\Omega\!\left(t^{3}B\right)$, we get $t^{2}=\mathcal{O}\!\left(\frac{n}{tB}\right)$. If $k>t$, then each thread processes $\mathcal{O}\!\left(\frac{k}{t}\right)$ buckets with a total cost of $\mathcal{O}\!\left(k/t\cdot t^{2}\right)=\mathcal{O}\!\left(kt\right)$. Since $n=\Omega\!\left(kt^{2}B\right)$, we have $kt=\mathcal{O}\!\left(\frac{n}{tB}\right)$.

Comparing the I/O volume of IS4o and s3-sort.

Both algorithms read and write the data once for the base case – $16n$ bytes (of I/O volume). IS4o reads and writes all data both during data distribution and block permutation phase – $32n$ bytes or $48n$ bytes overall. s3-sort reads the element twice but writes them only once in its distribution algorithm – $24n$ bytes. This sounds like a slight advantage. However, now we come to overheads unique to s3-sort. First, the algorithm reads and writes an oracle sequence that indicates the bucket for each element – $2n$ bytes. s3-sort has to copy the sorted result data back to the input array – $16n$ bytes. It also has to allocate the temporary arrays. For security reasons, that memory is zeroed by the operating system – $9n$ bytes.333In current versions of the Linux kernel this is done by a single thread and thus results in a huge scalability bottleneck. When writing to the temporary arrays or during copying back, there are so called allocate misses that happen when an element is written to a cache block that is currently not in memory – that block is read to the cache because the CPU does not know that none of the data in that block will ever be read. This amounts to an I/O volume of up to $17n$ bytes. Furthermore, s3-sort may suffer more associativity misses than IS4o – the relative positions of the buckets in the temporary array are not coordinated while IS4o essentially sweeps a window of size $\approx bk$ through the memory during the distribution phase. For an average case analysis refer to [22]. Even ignoring the latter overhead we get a total I/O volume of $86n$ byte – more than twice as much as IS4o. Much of this overhead can be reduced using measures that are non-portable (or hard to make portable). In particular, non-temporal writes eliminate the allocation misses and also help to eliminate the associativity misses. One could also use a base case sorter that does the copying back as as side-effect when the number of recursion levels is odd. When sorting multiple times within an application, one can keep the temporary arrays without having to reallocate them. However, this may require a different interface to the sorter. Overall, depending on many implementation details, IS4o may require slightly more I/O volume than s3-sort or significantly less.

Proof B.2** (Proof of Theorem 4.2).**

The main space overhead is for $k$ buffer blocks of size $b$ for each of $t$ threads. This bound also covers smaller amounts of memory for the search tree ($\mathcal{O}\!\left(k\right)$), swap buffers and overflow buffers ($\mathcal{O}\!\left(bt\right)$), read and write pointers ($\mathcal{O}\!\left(kB\right)$ if we avoid false sharing), end pointers, and bucket boundary pointers. All of these data structures can be used for all levels of recursion. The term $\mathcal{O}\!\left(\log_{k}\frac{n}{n_{0}}\right)$ stems from the space for the recursion stack itself.

Appendix C More Measurements