Vectorized and performance-portable Quicksort

What about performance-per-watt?

Some of us like to think of ourselves writing open source as serving the public interest. It's hard to do that if you're focusing on an ISA the public doesn't have. I haven't seen any consumer hardware that has AVX512.

He's just investigating how to reproduce the claimed result. Not sure where you got your take from.

But he’s right.

Thanks these are really great. That's a course I wish my CS program had.

These are really helpful too. Another reply to my question suggests a course that looks really good too.

This is helpful. Thank you

I agree with you about WFA. I was just very impressed when the striped SW first came out when I was still a fledgling PhD student.

Thank you

Thanks! Feel free to raise Github issues if you'd like to ask/discuss anything.

I'm also a huge fan of Godbolt/Compiler Explorer. Highway is integrated there, so you can just copy in your functions. Here's the last throwaway test that I was using: https://gcc.godbolt.org/z/5azbK95j9

> things might get better in the future but for now we have to implement it another way

There's several possible answers. 1) For the issue you pointed to, that we cannot have arrays of N vectors, it's a reasonable workaround to instead allocate an array of N vectors' worth of elements, and trust the compiler to elide unnecessary Load/Store. This often works using clang, and would avoid having to manually unroll here. I do prefer to minimize the amount of compiler magic required for good performance, though, so typically we're unrolling manually as shown in that code.

2) If there are compiler bugs, typically the workarounds have to stay in because in the open-source world people are still using that compiler years later.

Automatically detecting when things get better is an interesting idea but I am not yet aware of such infrastructure.

Thanks for sharing the pointers. I hadn't come across that, will read soon :)

:) If I'm reading 100Gb/sec correctly as Gigabits, a Skylake core can also reach such bandwidth usage. The issue is that the system can only sustain that for roughly 8 cores. It would be very interesting to see bench_parallel results for an M1 Max, if anyone would like to give that a try? I suspect it will be comparable to Skylake-X, because 8 (performance) cores are not quite enough to utilize all the M1 Max bandwidth in this app. M1 may be more power efficient, though. It is also great to see more bandwidth per core. A hypothetical M2 with say 256-bit SVE vectors and 16 cores could be very interesting.

L2 caches are typically partitioned and private to a core. If that also applies to the M1 Max, then each core would only access 3 MB, thus the working set is larger than cache as intended.

I do believe NEON is the limiting factor here. I haven't looked into how many IPC we should expect, but even if it is 4 (the number of 128-bit NEON units), Skylake is often reaching 2 (with 512-bit vectors), so the measurement that an M1 core is about half as fast as SKX for this use case seems plausible.

Thanks! Highway is included in Chrome and Firefox. Not sure what you mean by process? It's pretty easy for example using git submodules - you tell Git the Highway version you want to depend on, notify your CMake or Bazel build system of the hwy library dependency, and that's it.

One requirement is C++11 or later - some people have asked about supporting C but I believe that's infeasible.

Definitely, I have on my TODO to give turbostat a try.

SIMDe implements the exact semantics of a given ISA using either the native intrinsics, or other platform's intrinsics when possible, otherwise autovectorization. This is great when you've already written code using e.g. NEON intrinsics and want it to run on x86 (similar to neon2sse), or the reverse. It's an impressive undertaking to implement thousands of instructions for each platform :)

Highway instead aims for a path that each platform can implement efficiently. For example, ReorderWidenMulAccumulate bridges differences between NEON and x86 which user code need not care about if they just want a dot product, and CountTrue is efficient on both platforms without requiring NEON to go to the extra trouble of matching the exact x86 movmskb semantics. Also, Highway uses only intrinsics and does not rely on autovectorization (+), which makes for more predictable performance.

+ except in the EMU128 target, which is only used if SIMD/intrinsics are disabled

Good question. There are only a few functions implemented in .cc files, tagged with HWY_DLLEXPORT, notably memory allocation and detecting x86 CPU capabilities. If it were necessary, we could likely strip those out or do a header-only library. The ops/intrinsics called from user code are inlined in headers.

Yes, this is discussed in section 4.2 of our paper: https://arxiv.org/pdf/2205.05982.pdf

In short, it turns out not to help for single core with vectors, but a few initial passes of ips4o (with 64..256-way partitioning) is faster for parallel sorts.

Depends on your goals? If you'd like to contribute code under Apache2, I'd be happy to collaborate. Feel free to reach out to the email in the Highway README.

Some of the things on my short to medium-term todo list include:

- CompressStore for 8-bit lanes (before Icelake, that will likely require splitting up into runs of 8 bytes with one PSHUFB + store each)

- also ExpandLoad (analogous to CompressStore) and LeadingZeroCount

- the C++ STL isn't autovectorized much (http://0x80.pl/notesen/2021-01-18-autovectorization-gcc-clan...). Manual vectorization can help, we have some of the simpler algorithms already in hwy/contrib/algo, others such as unique, reverse, minmax_element, all_of are more interesting.

Thanks :) If I understand correctly, you have tuples of (number, T*) and want to sort by number? That can be done by storing these tuples interleaved, reinterpret_casting to hwy::uint128_t (making sure that the pointer bits are in the least-significant half), and sorting those.

Interesting, I did not know that, thanks for the pointer. There is also an effort underway to update LLVM std::sort to BlockQuicksort, which can use some SIMD instructions similarly to what you describe.

The blog post does seem to be quite explicit about the fact that their contribution is in the field of portability, not novelty.

nothing happened before google

https://dl.acm.org/doi/10.1145/113379.113380

:) FWIW we also use a constant-time sorting network as the base case of the Quicksort recursion, that's a widely used optimization.

We also compare with state of the art platform-specific code. The interesting thing here is that they turn out to be slower than our approach using portable intrinsics (github.com/google/highway) :)

Is that still true? I was under the impression that post-Haswell CPUs don’t have that particular issue.

The CPU has a certain upper bound on power, TDP. That limit is independent of whether it is reached by executing scalar code on lots of cores, or SIMD on a few. One little-known fact is that out of order CPUs burn lots of energy just on scheduling instructions, ~10x more than the operation itself. SIMD amortizes that per-instruction overhead over multiple data elements, and actually increases the amount of useful work done per joule. We're talking 5-10x increase in energy efficiency.

Right, though your use of the terms "different registers" and "properly pipelined" is quite nuanced for a non-specialist. It leaves a lot to the imagination or prior experience!

If it is a compute-heavy code that was spilling to a stack and now can fit in the SIMD registers, there could be a speed up with less memory traffic. This is analogous to a program's working set either fitting in RAM or spilling to swap while it runs, but at a different level of the memory hierarchy. In the extreme case, your working set could be a fixed set of variables in expensive loops for some sequential algorithm, and so the compiler register placement ability and number of registers available could be the difference between effectively O(1) or O(n) memory accesses with respect to the number of compute steps.

Of course, you can find such transitions between registers/L1 cache, L1/L2 cache, L2/L3 cache (if present), local/remote RAM (for modern NUMA machines), etc. This happens as you transition from a pure CPU-bound case to something with bigger data structures, where the program focus has to move to different parts of the data to make progress. Naively holding other things equal, an acceleration of a CPU-bound code will of course mean you churn through these different data faster, which means more cache spills as you pull in new data and have to displace something else. Going back to the earlier poster's question, this cache spilling can have a detrimental effect on other processes sharing the same cache or NUMA node, just like one extremely bloated program can cause a virtual memory swapping frenzy and hurt every other program on the machine.

One form of "pipelining" optimization is to recognize when your loop is repeatedly fetching a lot of the same input data in multiple iterations and changing that to use variables/buffers to retain state between iterations and avoid redundant access. This often happens in convolution algorithms on arrays of multi-dimensional data, e.g. image processing and neural nets and many cellular-automata style or grid-based simulations. Your algorithm slides a "sampling window" along an array to compute an output value as a linear combination of all the neighboring values within the window using some filtering kernel (a fixed pattern of multiplicative scalars/weights). Naive code keeps loading this whole window once for each iteration of the loop to address one output position. It is effectively stateless except for the loops to visit every output. Pipelined code would manage a window buffer and fetch and replace only the fringe of the buffer while holding and reusing inputs from the previous iterations. While this makes the loops more stateful, it can dramatically reduce the number of memory reads and shift a memory-bound code back into a CPU-bound regime.

Other major optimizations for these N-dimensional convolutions are done by the programmer/designer: some convolution kernels are "decomposable" and the expensive multi-dimensional operation can be strength-reduced to a sequence of 1-dimensional convolutions with appropriately chosen filter kernels to produce the same mathematical result. This optimization has nothing to do with SIMD but simplifies the work to embody the operation. Fewer input values to read (e.g. a 1D line of neighbors instead of 2D box) and the number of arithmetic operations to do all the multiply-and-add steps to produce the output. Imagine a 2D filter that operates on an 8x8 window. The 2D convolution has to sample and combine 64 input neighbors per output, while the decomposed filter does 8 neighbors in each axis and so one quarter as many steps total after one pass for the X axis and one pass for the Y axis.

Naively, decompsition is done as separate 1D passes over the data and so performs more memory writes and allocates an additional intermediate array between the original input and final output. This is often a big win in spite of the extra memory demands. It's a lot more coding, but this could also be "pipelined" to some degree in N-dimensions to avoid pushing the entire intermediate results arrays through main memory. Approaches differ, but you can make different tradeoffs for how many intermediate results to store or buffer versus redundant calculation in a less stateful manner.

Generally, as you make the above kinds of optimizations your code becomes more "block-based", loading larger chunks of data with fewer changes to new "random" memory locations. This is very much like how databases and filesystems optimized their access to disk storage in prior decades to avoid random seeks for individual bytes or blocks and to instead use clever structures to support faster loading of sets of blocks or pages. When this is successful, your program does most of its memory access sequentially and achieves higher throughput that is bound by total memory bus bandwidth rather than random access latency limitations.

The final form of "pipelining" matters once you really hit your stride with these optimized, block-based programs. All this access again can cause cache spilling as you sustain high bandwidth memory access. But if you can provide the right hints to the CPU, or if the CPU is clever enough to detect the pattern, it can shift into a different cache management regime. Knowing that you will only access each block of data once and then move on (because you are holding the reusable state in registers), the CPU can use a different prefetching and spilling policy to both optimize for your next memory access and to avoid churning the whole cache with these values you will only read once and then ignore. This reduces the impact of the code on other concurrent tasks which can still enjoy more reasonable caching behavior in spite of the busy neighbor.

Thanks for your post! Let me explain my problem, which fits right to it. It works as an example, too.

So, in freepascal there's some subpar pattern matching/string search function. It annoyed me so greatly, I've started researching into the "state of the art" and apparently they're all crap. My fixed size pattern matching solution fits into less than 256 bytes, not counting setup code, which is minimal.

And here start the problems for me, because I have no way of actually comparing my solution to "the state of the art", because "the state of the art" doesn't offer me downloadable, compiled benchmarks I can just run to compare stuff to my own and learning a new language (like, i don't speak C and certainly never will), installing a compiler, getting through all the walls that inevitably be in my way to getting it to run, etc ... totally breaks my timebudget and thus I drop the problem.

Basically, what I'm looking for is:

Give me a problem that's worth tackling, with a defined set of parameters and a compiled solution I can compare my own against and I will beat it eventually. I've tried doing this on my own, but that's just a big no-go on many levels, so I'm dependent on someone handing something to me.

Like ... the QuickSort. I have no way of comparing my potential own solution against "the state of the art", which is - as far as I know - pretty much all that's holding me back.

There are various Postgres extensions (such as Citus) that are columnar, so with them it should be possible.

How about AWS Redshift which is columnar (though that is also based on Postgres)...

The data in indexes is stored columnar by most RDBMSes, as far as I know.

The method exists. They could make it faster without breaking anything.

https://developer.mozilla.org/en-US/docs/Web/JavaScript/Refe...

My bad. Apologies.

I don’t see any mention in the paper of thermal clock throttling concerns, which can really neuter performance of tools that sustain use of AVX operations over a period of time. For the quick benchmarks presented in the paper, of course it will be faster. What if I’m continuously hammering my CPU with AVX operations? I expect it to severely downclock.

Do you expect these optimizations making their way into std::sort eventually?