A Novel Hybrid Quicksort Algorithm Vectorized using AVX-512 on Intel Skylake

TL;DR

This paper introduces a new hybrid Quicksort algorithm optimized with AVX-512 vectorization on Intel Skylake CPUs, significantly improving sorting performance for various data types compared to existing libraries.

Contribution

It presents a novel two-part hybrid Quicksort leveraging AVX-512, combining a new partitioning method with a branch-free Bitonic sort for small arrays.

Findings

Achieves 4x faster sorting than GNU C++ std::sort

Outperforms Intel IPP library by 1.4x

Effective adaptation of classical algorithms with AVX-512

Abstract

The modern CPU's design, which is composed of hierarchical memory and SIMD/vectorization capability, governs the potential for algorithms to be transformed into efficient implementations. The release of the AVX-512 changed things radically, and motivated us to search for an efficient sorting algorithm that can take advantage of it. In this paper, we describe the best strategy we have found, which is a novel two parts hybrid sort, based on the well-known Quicksort algorithm. The central partitioning operation is performed by a new algorithm, and small partitions/arrays are sorted using a branch-free Bitonic-based sort. This study is also an illustration of how classical algorithms can be adapted and enhanced by the AVX-512 extension. We evaluate the performance of our approach on a modern Intel Xeon Skylake and assess the different layers of our implementation by sorting/partitioning integers, double floating-point numbers, and key/value pairs of integers. Our results demonstrate that our approach is faster than two libraries of reference: the GNU \emph{C++} sort algorithm by a speedup factor of 4, and the Intel IPP library by a speedup factor of 1.4.

Figures

Full text

I Introduction

Sorting is a fundamental problem in computer science that always had the attention of the research community, because it is widely used to reduce the complexity of some algorithms. Moreover, sorting is a central operation in specific applications such as, but not limited to, database servers [1] and image rendering engines [2]. Therefore, having efficient sorting libraries on new architecture could potentially leverage the performance of a wide range of applications.

The vectorization — that is, the CPU’s capability to apply a single instruction on multiple data (SIMD) — improves continuously, one CPU generation after the other. While the difference between a scalar code and its vectorized equivalent was ”only” of a factor of 4 in the year 2000 (SSE), the difference is now up to a factor of 16 (AVX-512). Therefore, it is indispensable to vectorize a code to achieve high-performance on modern CPUs, by using dedicated instructions and registers. The conversion of a scalar code into a vectorized equivalent is straightforward for many classes of algorithms and computational kernels, and it can even be done with auto-vectorization for some of them. However, the opportunity of vectorization is tied to the memory/data access patterns, such that data-processing algorithms (like sorting) usually require an important effort to be transformed. In addition, creating a fully vectorized implementation, without any scalar sections, is only possible and efficient if the instruction set provides the needed operations. Consequently, new instruction sets, such as the AVX-512, allow for the use of approaches that were not feasible previously.

The Intel Xeon Skylake (SKL) processor is the second CPU that supports AVX-512, after the Intel Knight Landing. The SKL supports the AVX-512 instruction set [13]: it supports Intel AVX-512 foundational instructions (AVX-512F), Intel AVX-512 conflict detection instructions (AVX-512CD), Intel AVX-512 byte and word instructions (AVX-512BW), Intel AVX-512 doubleword and quadword instructions (AVX-512DQ), and Intel AVX-512 vector length extensions instructions (AVX-512VL). The AVX-512 not only allows work on SIMD-vectors of double the size, compared to the previous AVX(2) set, it also provides various new operations.

Therefore, in the current paper, we focus on the development of new sorting strategies and their efficient implementation for the Intel Skylake using AVX-512. The contributions of this study are the following:

• •

proposing a new partitioning algorithm using AVX-512,

• •

defining a new Bitonic-sort variant for small arrays using AVX-512,

• •

implementing a new Quicksort variant using AVX-512.

All in all, we show how we can obtain a fast and vectorized sorting algorithm 111 The functions described in the current study are available at https://gitlab.mpcdf.mpg.de/bbramas/avx-512-sort. This repository includes a clean header-only library (branch master) and a test file that generates the performance study of the current manuscript (branch paper). The code is under MIT license. .

The rest of the paper is organized as follows: Section II gives background information related to vectorization and sorting. We then describe our approach in Section III, introducing our strategy for sorting small arrays, and the vectorized partitioning function, which are combined in our Quicksort variant. Finally, we provide performance details in Section IV and the conclusion in Section V.

II Background

II-A Sorting Algorithms

II-A1 Quicksort (QS) Overview

QS was originally proposed in [3]. It uses a divide-and-conquer strategy, by recursively partitioning the input array, until it ends with partitions of one value. The partitioning puts values lower than a pivot at the beginning of the array, and greater values at the end, with a linear complexity. QS has a worst-case complexity of O(n2), but an average complexity of O(n log n) in practice. The complexity is tied to the choice of the partitioning pivot, which must be close to the median to ensure a low complexity. However, its simplicity in terms of implementation, and its speed in practice, has turned it into a very popular sorting algorithm. Figure 1 shows an example of a QS execution.

We provide in Appendix A-A the scalar QS algorithm. Here, the term scalar refers to a single value at the opposite of an SIMD vector. In this implementation, the choice of the pivot is naively made by selecting the value in the middle before partitioning, and this can result in very unbalanced partitions. This is why more advanced heuristics have been proposed in the past, like selecting the median from several values, for example.

II-A2 GNU std::sort Implementation (STL)

The worst case complexity of QS makes it no longer suitable to be used as a standard C++ sort. In fact, a complexity of O(n log n) in average was required until year 2003 [4], but it is now a worst case limit [5] that a pure QS implementation cannot guarantee. Consequently, the current implementation is a 3-part hybrid sorting algorithm i.e. it relies on 3 different algorithms 222See the libstdc++ documentation on the sorting algorithm available at https://gcc.gnu.org/onlinedocs/libstdc++/libstdc++-html-USERS-4.4/a01347.html#l05207. The algorithm uses an Introsort [6] to a maximum depth of 2 $\times$ log2 n to obtain small partitions that are then sorted using an insertion sort. Introsort is itself a 2-part hybrid of Quicksort and heap sort.

II-A3 Bitonic Sorting Network

In computer science, a sorting network is an abstract description of how to sort a fixed number of values i.e. how the values are compared and exchanged. This can be represented graphically, by having each input value as a horizontal line, and each compare and exchange unit as a vertical connection between those lines. There are various examples of sorting networks in the literature, but we concentrate our description on the Bitonic sort from [7]. This network is easy to implement and has an algorithm complexity of O(n log(n)2 ). It has demonstrated good performances on parallel computers [8] and GPUs [9]. Figure 2(a) shows a Bitonic sorting network to process 16 values. A sorting network can be seen as a time line, where input values are transferred from left to right, and exchanged if needed at each vertical bar. We illustrate an execution in Figure 2(b), where we print the intermediate steps while sorting an array of 8 values. The Bitonic sort is not stable because it does not maintain the original order of the values.

If the size of the array to sort is known, it is possible to implement a sorting network by hard-coding the connections between the lines. This can be seen as a direct mapping of the picture. However, when the array size is unknown, the implementation can be made more flexible by using a formula/rule to decide when to compare/exchange values.

II-B Vectorization

The term vectorization refers to a CPU’s feature to apply a single operation/instruction to a vector of values instead of only a single value [10]. It is common to refer to this concept by Flynn’s taxonomy term, SIMD, for single instruction on multiple data. By adding SIMD instructions/registers to CPUs, it has been possible to increase the peak performance of single cores, despite the stagnation of the clock frequency. The same strategy is used on new hardware, where the length of the SIMD registers has continued to increase. In the rest of the paper, we use the term vector for the data type managed by the CPU in this sense. It has no relation to an expandable vector data structure, such as std::vector. The size of the vectors is variable and depends on both the instruction set and the type of vector element, and corresponds to the size of the registers in the chip. Vector extensions to the x86 instruction set, for example, are SSE [11], AVX [12], and AVX512 [13], which support vectors of size 128, 256 and 512 bits respectively. This means that an SSE vector is able to store four single precision floating point numbers or two double precision values. Figure 3 illustrates the difference between a scalar summation and a vector summation for SSE or AVX, respectively. An AVX-512 SIMD-vector is able to store 8 double precision floating-point numbers or 16 integer values, for example. Throughout this document, we use intrinsic function extension instead of the assembly language to write vectorized code on top of the AVX-512 instruction set. Intrinsics are small functions that are intended to be replaced with a single assembly instruction by the compiler.

II-B1 AVX-512 Instruction Set

As previous x86 vectorization extensions, the AVX-512 has instructions to load a contiguous block of values from the main memory and to transform it into a SIMD-vector (load). It is also possible to fill a SIMD-vector with a given value (set), and move back a SIMD-vector into memory (store). A permutation instruction allows to re-order the values inside a SIMD-vector using a second integer array which contains the permutation indexes. This operation was possible in since AVX/AVX2 using permutevar8x32 (instruction vperm(d,ps)). The instructions vminpd/vpminsd return a SIMD-vector where each value correspond to the minimum of the values from the two input vectors at the same position. It is possible to obtain the maximum with instructions vpmaxsd/vmaxpd.

In AVX-512, the value returned by a test/comparison (vpcmpd/vcmppd) is a mask (integer) and not an SIMD-vector of integers, as it was in SSE/AVX. Therefore, it is easy to modify and work directly on the mask with arithmetic and binary operations for scalar integers. Among the mask-based instructions, the mask move (vmovdqa32/vmovapd) allows for the selection of values between two vectors, using a mask. Achieving the same result was possible in previous instruction sets using the blend instruction since SSE4, and using several operations with previous instruction sets.

The AVX-512 provides operations that do not have an equivalent in previous extensions of the x86 instruction sets, such as the store-some (vpcompressps/vcompresspd) and load-some (vmovups/vmovupd). The store-some operation allows to save only a part of a SIMD-vector into memory. Similarly, the load-some allows to load less values than the size of a SIMD-vector from the memory. The values are loaded/saved contiguously. This is a major improvement, because without this instruction, several operations are needed to obtain the same result. For example, to save some values from a SIMD-vector v at address p in memory, one possibility is to load the current values from p into a SIMD-vector v’, permute the values in v to move the values to store at the beginning, merge v and v’, and finally save the resulting vector.

II-C Related Work on Vectorized Sorting Algorithms

The literature on sorting and vectorized sorting implementations is extremely large. Therefore, we only cite some of the studies that we consider most related to our work.

The sorting technique from [14] tries to remove branches and improves the prediction of a scalar sort, and they show a speedup by a factor of 2 against the STL (the implementation of the STL was different at that time). This study illustrates the early strategy to adapt sorting algorithms to a given hardware, and also shows the need for low-level optimizations, due to the limited instructions available at that time.

In [15], the authors propose a parallel sorting on top of combosort vectorized with the VMX instruction set of IBM architecture. Unaligned memory access is avoided, and the L2 cache is efficiently managed by using an out-of-core/blocking scheme. The authors show a speedup by a factor of 3 against the GNU C++ STL.

In [16], the authors use a sorting-network for small-sized arrays, similar to our own approach. However, instead of dividing the main array into sorted partitions (partitions of increasing contents), and applying a small efficient sort on each of those partitions, the authors perform the opposite. They apply multiple small sorts on sub-parts of the array, and then they finish with a complicated merge scheme using extra memory to globally sort all the sub-parts. A very similar approach was later proposed in [17].

The recent work in [18] targets AVX2. The authors use a Quicksort variant with a vectorized partitioning function, and an insertion sort once the partitions are small enough (as the STL does). The partition method relies on look-up tables, with a mapping between the comparison’s result of an SIMD-vector against the pivot, and the move/permutation that must be applied to the vector. The authors demonstrate a speedup by a factor of 4 against the STL, but their approach is not always faster than the Intel IPP library. The proposed method is not suitable for AVX-512 because the lookup tables will occupy too much memory. This issue, as well as the use of extra memory, can be solved with the new instructions of the AVX-512. As a side remark, the authors do not compare their proposal to the standard C++ partition function, even so, it is the only part of their algorithm that is vectorized.

III Sorting with AVX-512

III-A Bitonic-Based Sort on AVX-512 SIMD-Vectors

In this section, we describe our method to sort small arrays that contain less than 16 times VEC_SIZE, where VEC_SIZE is the number of values in a SIMD-vector. This function is later used in our final QS implementation to sort small enough partitions.

III-A1 Sorting one SIMD-vector

To sort a single vector, we perform the same operations as the ones shown in Figure 2(a): we compare and exchange values following the indexes from the Bitonic sorting network. However, thanks to the vectorization, we are able to work on the entire vector without having to iterate on the values individually. We know the positions that we have to compare and exchange at the different stages of the algorithm. This is why, in our approach, we rely on static (hard-coded) permutation vectors, as shown in Algorithm 1. In this algorithm, the compare_and_exchange function performs all the compare and exchange that are applied at the same time in the Bitonic algorithm i.e. the operations that are at the same horizontal position in the figure. To have a fully vectorized function, we implement the compare_and_exchange in three steps. First, we permute the input vector v into v’ with the given permutation indexes p. Second, we obtain two vectors wmin and wmax that contain the minimum and maximum values between both v and v’. Finally, we selects the values from wmin and wmax with a mask-based move, where the mask indicates in which direction the exchanges have to be done. The C++ source code of a fully vectorized branch-free implementation is given in Appendix A-B (Code 1).

III-A2 Sorting more than one SIMD-vectors

The principle of using static permutation vectors to sort a single SIMD-vector can be applied to sort several SIMD-vectors. In addition, we can take advantage of the repetitive pattern of the Bitonic sorting network to re-use existing functions. More precisely, to sort $V$ vectors, we re-use the function to sort $V/2$ vectors and so on. We provide an example to sort two SIMD-vectors in Algorithm 2, where we start by sorting each SIMD-vector individually using the bitonic_simd_sort_1v function. Then, we compare and exchange values between both vectors (line 5), and finally applied the same operations on each vector individually (lines 6 to 11). In our sorting implementation, we provide the functions to sort up to 16 SIMD-vectors, which correspond to 256 integer values or 128 double floating-point values.

III-A3 Sorting Small Arrays

Each of our simd-Bitonic-sort functions are designed for a specific number of SIMD-vectors. However, we intend to sort arrays that do not have a size multiple of the SIMD-vector’s length, because they are obtained from the partitioning stage of the QS. Consequently, when we have to sort a small array, we first load it into SIMD-vectors, and then, we pad the last vector with the greatest possible value. This guarantee that the padding values have no impact on the sorting results by staying at the end of the last vector. The selection of appropriate simd-Bitonic-sort function, that matches the size of the array to sort, can be done efficiently with a switch statement. In the following, we refer to this interface as the simd_bitonic_sort_wrapper function.

III-B Partitioning with AVX-512

Algorithm 3 shows our strategy to develop a vectorized partitioning method. This algorithm is similar to a scalar partitioning function: there are iterators that start from both extremities of the array to keep track of where to load/store the values, and the process stops when some of these iterators meet. In its steady state, the algorithm loads an SIMD-vector using the left or right indexes (at lines 19 and 24), and partitions it using the partition_vec function (at line 27). The partition_vec function compares the input vector to the pivot vector (at line 47), and stores the values — lower or greater — directly in the array using a store-some instruction (at lines 51 and 55). The store-some is an AVX-512 instruction that we described in Section II-B1. The initialization of our algorithm starts by loading one vector from each array’s extremities to ensure that no values will be overwritten during the steady state (lines 12 and 16). This way, our implementation works in-place and only needs three SIMD-vectors. Algorithm 3 also includes, as side comments, possible optimizations in case the array is more likely to be already partitioned (A), or to reduce the data displacement of the values (B). The AVX-512 implementation of this algorithm is given in Appendix A-B (Code 2). One should note that we use a scalar partition function if there are less than 2 $\times$ VEC_SIZE values in the given array (line 3).

III-C Quicksort Variant

Our QS is given in Algorithm 4, where we partition the data using the simd_partition function from Section III-B, and then sort the small partitions using the simd_bitonic _sort_wrapper function from Section III-A. The obtained algorithm is very similar to the scalar QS given in Appendix A-A.

III-D Sorting Key/Value Pairs

The previous sorting methods are designed to sort an array of numbers. However, some applications need to sort key/value pairs. More precisely, the sort is applied on the keys, and the values contain extra information and could be pointers to arbitrary data structures, for example. Storing each key/value pair contiguously in memory is not adequate for vectorization because it requires transforming the data. Therefore, in our approach, we store the keys and the values in two distinct arrays. To extend the simd-Bitonic-sort and simd-partition functions, we must ensure that the same permutations/moves are applied to the keys and the values. For the partition function, this is trivial. The same mask is used in combination with the store-some instruction for both arrays. For the Bitonic-based sort, we manually apply the permutations that were done on the vector of keys to the vector of values. To do so, we first save the vector of keys k before it is permuted by a compare and exchange, using the Bitonic permutation vector of indexes p, into k’. We compare k and k’ to obtain a mask m that expresses what moves have been done. Then, we permute our vector of values v using p into v’, and we select the correct values between v and v’ using m. Consequently, we perform this operation at the end of the compare_and_exchange in all the Bitonic-based sorts.

IV Performance Study

IV-A Configuration

We asses our method on an Intel(R) Xeon(R) Platinum 8170 Skylake CPU at 2.10GHz, with caches of sizes 32K-Bytes, 1024K-Bytes and 36608K-Bytes, at levels L1, L2 and L3 respectively. The process and allocations are bound with numactl –physcpubind=0 –localalloc. We use the Intel compiler 17.0.2 (20170213) with aggressive optimization flag -O3.

We compare our implementation against two references. The first one is the GNU STL 3.4.21 from which we use the std::sort and std::partition functions. The second one is the Intel Integrated Performance Primitives (IPP) 2017 which is a library optimized for Intel processors. We use the IPP radix-based sort (function ippsSortRadixAscend_[type]_I). This function require additional space, but it is known as one of the fastest existing sorting implementation.

The test file used for the following benchmark is available online333 The test file that generates the performance study is available at https://gitlab.mpcdf.mpg.de/bbramas/avx-512-sort (branch paper) under MIT license. , it includes the different sorts presented in this study plus some additional strategies and tests. Our simd-QS uses a 3-values median pivot selection (similar to the STL sort function). The arrays to sort are populated with randomly generated values.

IV-B Performance to Sort Small Arrays

Figure 4 shows the execution times to sort arrays of size from 1 to 16 $\times$ VEC_SIZE, which corresponds to 128 double floating-point values, or 256 integer values. We also test arrays of size not multiple of the SIMD-vector’s length. The AVX-512-bitonic always delivers better performance than the Intel IPP for any size, and better performance than the STL when sorting more than 5 values. The speedup is significant, and is around 8 in average. The execution time per item increases every VEC_SIZE values because the cost of sorting is not tied to the number of values but to the number of SIMD-vectors to sort, as explained in Section III-A3. For example, in Figure 4(a), the execution time to sort 31 or 32 values is the same, because we sort one SIMD-vector of 32 values in both cases. Our method to sort key/value pairs seems efficient, see Figure 4(c), because the speedup is even better against the STL compared to the sorting of integers.

IV-C Partitioning Performance

Figure 5 shows the execution times to partition using our AVX-512-partition or the STL’s partition function. Our method provides again a speedup of an average factor of 4. For the three configurations, an overhead impacts our implementation and the STL when partitioning arrays larger than 107 items. Our AVX-512-partition remains faster, but its speedup decreases from 4 to 3. This phenomena is related to cache effects since 107 integers values occupy 40M-Bytes, which is more than the L3 cache size. In addition, we see that this effect starts from 105 when partitioning key/value pairs.

IV-D Performance to Sort Large Arrays

Figure 6 shows the execution times to sort arrays up to a size of 109 items. Our AVX-512-QS is always faster in all configurations. The difference between AVX-512-QS and the STL sort seems stable for any size with a speedup of more than 6 to our benefit. However, while the Intel IPP is not efficient for arrays with less than 104 elements, its performance is really close to the AVX-512-QS for large arrays. The same effect found when partitioning appears when sorting arrays larger than 107 items. All three sorting functions are impacted, but the IPP seems more slowdown than our method, because it is based on a different access pattern, such that the AVX-512-QS is almost twice as fast as IPP for a size of 109 items.

V Conclusions

In this paper, we introduced new Bitonic sort and a new partition algorithm that have been designed for the AVX-512 instruction set. These two functions are used in our Quicksort variant which makes it possible to have a fully vectorized implementation (at the exception of partitioning tiny arrays). Our approach shows superior performance on Intel SKL in all configurations against two reference libraries: the GNU C++ STL, and the Intel IPP. It provides a speedup of 8 to sort small arrays (less than 16 SIMD-vectors), and a speedup of 4 and 1.4 for large arrays, against the C++ STL and the Intel IPP, respectively. These results should also motivate the community to revisit common problems, because some algorithms may become competitive by being vectorizable, or improved, thanks to AVX-512’s novelties. Our source code is publicly available and ready to be used and compared. In the future, we intend to design a parallel implementation of our AVX-512-QS, and we expect the recursive partitioning to be naturally parallelized with a task-based scheme on top of OpenMP.

Appendix A APPENDIX

A-A Scalar Quicksort Algorithm

A-B Source Code Extracts