Designing a SIMD Algorithm from Scratch
mcyoung.xyz
mcyoung.xyz
Although looking at the code, this article confirms my experience that Rust has rather poor ergonomics for simd and pointer-related work (and performance engineering in general).
But a "new" language like Rust would have the freedom to add special language syntax and operators for such things instead of going the C++ way of "just add it to the stdlib" (which IMHO is almost always the wrong place... it should either go into a regular cargo dependency, or into the language itself).
But yep portable SIMD in Rust is still not a good story, compared to C++. And getting down to raw byte regions, pointer & buffer manipulation, etc requires becoming comfortable with Pin, MaybeUninit, etc. Both portable simd and allocator_api have been sitting unstable for years. And the barrier to entry is a bit higher, and it's more awkward ... mostly on purpose
But there's nothing stopping one from building one's own abstractions (or using 3rd party crates etc) to make these things more ergonomic within one's own program?
I really wish compilers were better at auto vectorization. And some support added for annotations in the language to locally allow reordering some operations, ...
Floating point addition is not associative, so summing each value in order is not the same as summing every 8th element, then summing the remainder, which is how SIMD handles it. So even though this is an obvious optimization for compilers, they will prioritize the serial guarantee over the optimization unless you tell it to relax that particular guarantee.
It's a mess and I agree with janwas: use a library (and in particular: use Google Highway) or something like Intel's ISPC when your hot path needs this.
Problem is you need to enforce that requirement on all user compilations, and I don't know what for MSVC. In-language would be nice.
I tried it with:
double sum1(double arr[128]) {
double tot = 0.0;
for (int i=0; i<128; i++) {
tot += arr[i];
}
return tot;
}
__attribute__((optimize("-ffast-math")))
double sum2(double arr[128]) {
double tot = 0.0;
for (int i=0; i<128; i++) {
tot += arr[i];
}
return tot;
}
sum1(double*):
lea rax, [rdi+1024]
vxorpd xmm0, xmm0, xmm0
.L2:
vaddsd xmm0, xmm0, QWORD PTR [rdi]
add rdi, 32
vaddsd xmm0, xmm0, QWORD PTR [rdi-24]
vaddsd xmm0, xmm0, QWORD PTR [rdi-16]
vaddsd xmm0, xmm0, QWORD PTR [rdi-8]
cmp rax, rdi
jne .L2
ret
sum2(double*):
lea rax, [rdi+1024]
vxorpd xmm0, xmm0, xmm0
.L6:
vaddpd ymm0, ymm0, YMMWORD PTR [rdi]
add rdi, 32
cmp rax, rdi
jne .L6
vextractf64x2 xmm1, ymm0, 0x1
vaddpd xmm1, xmm1, xmm0
vunpckhpd xmm0, xmm1, xmm1
vaddpd xmm0, xmm0, xmm1
vzeroupper
ret
It can be very unreliable.
But that’s one of the points of a systems programming language (of which C++ is one) — it tries to be portably as efficient as possible but makes it easy to do target-specific programming when required.
> I really wish compilers were better at auto vectorization
FORTRAN compilers sure are, since aliasing is not allowed. C++ is kneecapped by following C’s memory model.
Here's a trivial algorithm that gcc borks over because 'bool' is used instead of 'int', even fixing that the codegen isn't great:
It also sucks that you have no way of knowing at build time whether autovectorisation succeeded. Even if you have it working initially it can silently break and you don't know where. So using clang, you're careful and you check the generated code and everything is great, then later it isn't and you have to go hunting to find out why.
The larger issue is that C and C++ have convoluted aliasing rules which makes reasoning and flagging non-aliased regions hard. The most innocent of cast can break strict aliasing rules and send you in horrible debugging scenarios. The sad side effect being people adding -fno-strict-aliasing to -O2/3.
Since with templating everything seems to be going header only anyways, I feel like we should give the compiler more power to reason about the memory layout, even through multiple levels of function invocation/inlining.
I wrote some library where basically all shapes can be infered at compile-time (through templates) but the compiler only rarely seems to take advantage of this.
Or ROCm (basically CUDA but for AMD).
I always was a fan of Microsoft's C++AMP though. I thought that was easiest to get into. Too bad it never stuck though.
ISPC is an interesting take on 'building a vector DSL which can easily be integrated in a C++ build-chain'.
Though calling GPUs SIMD nowadays is kind of reductive, since you also have to contain with a very constrained memory hierarchy, and you don't write or think your code as SIMD much but more. Also they don't give access to most bit-tricks, and the rigmarole of shuffles one would expect from a SIMD processor.
Oh rly?
* Popcnt is defined in ROCm and CUDA.
* Shifts, XORs, AND, OR, NOT are all defined in GPUs.
* You can actually do AES and a lot of other bit-level twiddling in GPU space very efficiently. (aka: see any mining software for highly-efficient bit-twiddling).
Yeah, not as many tricks (ex: pext / pdep) as a CPU. But popcnt is the "big" one. GPUs also have bit-reverse instructions (which reduces the need for LS1B instructions like CPUs because you can do LS1B tricks on "both ends" by just bit-reversing).
But honestly, GPUs are so parallel and Register-space is so plentiful that you can probably just build whatever you want out of AND/OR/NOT/XOR instructions alone.
And since GPUs have single-cycle "popcnt", any symmetric bit-twiddling function can be built off of popcnt within a few cycles.
> and the rigmarole of shuffles one would expect from a SIMD processor.
Abuse of __shared__ memory and the full crossbar of a GPU-core means that you can arbitrarily shuffle data through __shared__ memory in like 5 clock ticks (well... in certain ways... if you do it wrong it'd take many clock ticks).
GPUs are actually far superior to AVX512 in this front. NVidia's shfl.bfly instruction can implement FFTs, even-odd sorting, bitonic sorts, and more as high-speed communications.
AMD has similar instructions: https://gpuopen.com/learn/amd-gcn-assembly-cross-lane-operat...
Trust me on this case: GPU data-movement is far, far, far superior to CPU data-movement.
Intel does NOT have the bpermute or permute instructions like a GPU does. Meanwhile, both AMD and NVidia GPUs have bpermute / permute for both gather-and-scatter-like operation / distribution of data across lanes.
AMD/NVidia also have guarantees upon the speed of broadcasts across __shared__ memory. And butterfly-shuffles can theoretically implement any arbitrary shuffle you want in log2(SIMD-Width) steps, so we have incredible amounts of tools in GPU space that CPU-programmers have no idea about.
This sounds interesting, do you have a reference to this?
The gist is that if you have a symmetric function (ie: order doesn't matter) of say, 8-bits, then that means f(10101010) == f(11110000) == f(00001111), and all such combinations. Because this is the very definition of "order doesn't matter".
This necessarily implies that f(bits) can be rewritten into the form of f(bits) == g(popcnt(bits)).
Something to do with counting up all the possibilities of "order doesn't matter" and then pigeon-holing them into popcnt combinations or something really mathematical like that. I'm sorry I don't remember the proof, but hopefully this is enough to give you the gist of the idea.
------------------
For example: XOR is the simplest symmetric function. We can see that XOR can be rewritten as XOR(bits) == g(popcnt(bits)). Where G is "take the bottom-bit from popcnt".
A lot of the "power" of BDDs is their ability to uncannily decompose functions into symmetric and non-symmetric parts. Not consistently mind you, but if a BDD is well-behaved (low-memory space, high-speeds, etc. etc.), its likely because a large part of the calculations happened to be symmetric.
This can be beneficial with say addFourDWORDs(a, b, c, d), where a+b+c+d has a "partial" level of symmetry. (a_0 XOR b_0 XOR c_0 XOR d_0 determines output_0 bit... while the 1st bit is related to XOR (a1,b1,c1,d1), etc. etc.). So there's all kinds of 'hidden popcounts" that could pop up in practice for random functions (IE: "partially symmetric"), though its a puzzle on how to exactly decompose arbitrary bit-functions into this form.
And even if you do, its no guarantee that the popcnt form was actually faster either. But it does give you some "trick" to try when trying to optimize an arbitrary bitwise function into a high-speed routine.
-----------------
Other symmetric functions (assuming 8-bit numbers)
AND(bits) == (popcnt(bits) == 8).
OR(bits) == (popcnt(bits) > 0)
Etc. etc.
EDIT: Looks like Wikipedia has a link to this concept: https://en.wikipedia.org/wiki/Symmetric_Boolean_function
It seems I have lots of reading to do and lots of ways to improve my sorting networks / counting sort implementations.
Thanks again.
PTX is maintained over time. Its a high-level assembly so to speak, the full details of the machine remain abstracted so that code can be more portable.
SASS is not. SASS changes from architecture-to-architecture. SASS is the actual machine code of NVidia cards. There's an overall understanding of SASS in the GPU world but its not really documented and you "shouldn't" want to learn about it.
--------
I should note that Intel's "pshufb" instruction is very similar to the permute instruction in NVidia/AMD. So yeah, there's a high-speed generic shuffle that's key to Intel/AMD AVX512 code.
But having the backwards-direction (bpermute) available too, as well as __shared__ memory for all other cases is great.
That instruction is from AVX-512. There’s no AVX2 float to int cast. There’s in AVX1, the int is signed and the instruction is _mm256_cvtps_epi32.
I agree implementing dynamic dispatch is difficult, but Highway has taken care of that :)
I dip in and out of these waters (performance optimization, more systems bare metal engineering), but basically on personal projects. I wish I had jobs where it was more called for, but it's not what most industry work needs.
Not relevant to this example, which is `popcount`ing only a single number, but AVX512 did also introduce SIMD popcount instructions for getting many inputs at a time. Also true for other useful bit-funs like leading or trailing zeros. So if you're using zmm registers, you can do way better than that godbolt example.
pub fn popcnt(mut x: u32) -> u32 {
x.count_ones()
}
example::popcnt:
popcnt eax, edi
ret
32 bits is not large enough, and the code Rust produces is indeed comically bad.
[1]: https://github.com/llvm/llvm-project/blob/08a6968127f04a40d7... [2]: https://llvm.org/docs/LangRef.html#llvm-ctpop-intrinsic
It's dramatically different, C++ compiles down to 10 or so instructions even without the vectorization turned on.
Rust OTOH compiles down to 188 or 129 instructions with the vectorization. This is pretty bad.
Instruction counts are not everything. The rust assembly is completely branchless, and on my M2 MacBook I'm benching it at ~2ns per call vs ~10ns per call with the equivalent C++ code (running this: https://godbolt.org/z/nTanvfs76)
-----------------------------------------------------
Benchmark Time CPU Iterations
-----------------------------------------------------
BM_cpp 10.1 ns 10.1 ns 69503053
BM_rs 1.74 ns 1.72 ns 393022172
LLVM IR for the benchmark: https://godbolt.org/z/8PaoK4sYz
These results match a criterion benchmark: https://gist.github.com/orf/6dc9d6dbfdcd0df464c1c0a7274f6286
I could be doing something horribly wrong here - this was a nice intro into a few different things I haven't used before. But I ended up taking the rust function and outputting the LLVM IR for it, then compiling and linking that with the C++ benchmark. And it worked, after I realised how to get the optimizer to stop removing everything.
Similar results on Intel can be seen here: https://quick-bench.com/q/RxcWefXX5jo77NM-eK8u36W6ohk
.L4:
mov edx, edi
and edx, 1
cmp edx, 1
sbb eax, -1
shr edi
jne .L4
ret
Vectorized execution OTOH is not cheap and does not necessarily result in faster code so purely seeing it in assembly would not mean much without actually measuring it. Scalar versions of the same code may be faster than the compiler-generated auto-vectorization code.
Each SIMD instruction has a latency attached to it (see https://www.intel.com/content/www/us/en/docs/intrinsics-guid...) and the majority of those SIMD latencies are in between 1 and 7 CPU cycles so definitely not free as in a free beer.
That said, I am surprised by the results - in my mental model of how things work this does not add up. Rust assembly is ~50 instructions, roughly ~30 of them being vector (SIMD) instructions and the rest of ~20 instructions being the superset of (only) 6 instructions used in C++ assembly. I cannot see how those 6 instructions, even with the branch mispredict cost of ~15 cycles, can run ~2.5x slower than the pretty much convoluted version of the SIMD popcnt.
> I cannot see how those 6 instructions, even with the branch mispredict cost of ~15 cycles, can run ~2.5x slower than the pretty much convoluted version of the SIMD popcnt.
I fleshed out the benchmarks a bit, to compare `popcnt` itself and a no-op baseline: https://quick-bench.com/q/Bb36LyruJ1r8pzm8o_hjckU-xcI
The quick-bench tool gives some nice assembly-level timing information[1] but I couldn't get it to work with my janky inline-copy-paste-assembly, and I blew past my curiosity time budget for this. You seem pretty experienced and I'd love to know more about this, so maybe if you know how to get the "popcnt_rs" extern assembly function to be inlined in the quick-bench link below we could take a look in more detail?
1. Click "assembly" and select "BM_cpp" from here: https://quick-bench.com/q/Bb36LyruJ1r8pzm8o_hjckU-xcI
(I wrote this initially as an EDIT in my response above but since you replied at about the same time, I am copying it here)
This sort of thing is part of why you'd often like to target an intermediate vector representation and let the compiler figure out the details. E.g., on Haswell chips you had multiple floating point execution units per core, and the CPU could definitely execute more than one pipelined FP operation simultaneously, but only one of those could be an `add` instruction. If you had a bunch of additions to do which didn't rely on the previously computed results (avoiding stalls), you could double your addition throughput by also dispatching a fused-multiply-add instruction where the multiplicand was just one. That _could_ execute at the same time as a normal vector FP addition.
https://github.com/rust-lang/rust/issues/107617#issuecomment...
> Improvements to NumPy’s performance are important to many users. We have focused this effort on Universal SIMD (see NEP 38 — Using SIMD optimization instructions for performance) intrinsics which provide nice improvements across various hardware platforms via an abstraction layer. The infrastructure is in place, and we welcome follow-on PRs to add SIMD support across all relevant NumPy functions
"NEP 38 — Using SIMD optimization instructions for performance" (2019) https://numpy.org/neps/nep-0038-SIMD-optimizations.html#nep3...
NumPy docs > CPU/SIMD Optimizations: https://numpy.org/doc/stable/reference/simd/index.html
std::simd: https://doc.rust-lang.org/std/simd/index.html
"Show HN: SimSIMD vs SciPy: How AVX-512 and SVE make SIMD nicer and ML 10x faster" (2023-10) https://news.ycombinator.com/item?id=37808036
"Standard library support for SIMD" (2023-10) https://discuss.python.org/t/standard-library-support-for-si...
SIMD: Single instruction, multiple data: https://en.wikipedia.org/wiki/Single_instruction,_multiple_d...
Category:SIMD computing: https://en.wikipedia.org/wiki/Category:SIMD_computing
Vectorization: Introduction: https://news.ycombinator.com/item?id=36159017 :
> GPGPU > Vectorization, Stream Processing > Compute kernels: https://en.wikipedia.org/wiki/General-purpose_computing_on_g...
I recently helped my friend who needed a fast CCSDS Reed-Solomon decoder in Python. The existing Python library was really slow, like 40 blocks per second, while the same algorithm written in C# claimed 10,000+ blocks per second. It turned out that the Python version made use of Numpy but so badly that it was much faster to eschew Numpy. After some additional optimizations to avoid modulos, my version easily handled ~2,000 blocks per second without Numpy.