All those don't matter if they hardly make a 10-15% difference in everyday use, and require specialized code to take advantage of.
All those don't matter if they hardly make a 10-15% difference in everyday use, and require specialized code to take advantage of.
But my work (statistics) seems poorly represented in benchmarks and the community at large. Ryzen is in the later group, and gets rave reviews.
While most press I've seen talking about avx512 has been negative. Yet I get a lot more flops out of the 10-core i9 7900x than the similarly priced 16-core Threadripper 1950x in simulations.
A big part of that is probably doing most of my work in Julia, where because the code you run is compiled locally the first time a function/argument type combination is called (a) code will be generated for all your CPU's features, and (b) it is easy to encourage aggressive optimizations to vectorized code by the compiler, eg by parameterizing arrays by their bounds.
R benchmarks will not see benefits at all, as much of R's math libraries are terribly unoptimized -- the default BLAS doesn't even use any runtime dispatch! (Thankfully it's easy to change, but the benchmarks all reflect that reality -- and most people's personal experiences too).
Similar likely applies to many NumPy benchmarks, but depending on the distributor of those binaries, the situation likely isn't better. Either generic unoptimized binaries, or shipping with something like Intel MKL which, while the fastest linear algebra library for Intel processors, dispatches to slow code paths for AMD processors -- a very unfair choice for benchmarking. I'm less less familiar with NumPy though, so perhaps there's a booth fair and well optimized distribution.
I fear the future will move against my use cases.
But a CPU has latency benefits. Most noticeably, it takes only a few clock cycles for an Intel CPU to transfer data from its RAX register -> AVX vector registers, compute the solution, and transfer it back to the "scalar world".
This level of data-transfer is done within single-digit nanoseconds.
In contrast, any GPU memory transfer takes hundreds-of-nanoseconds to microseconds... an order of 100x to 10,000x slower than transferring between Scalar-code and AVX-code on a CPU.
----------
So GPUs are good for bulk compute (big matrix multiplies, Deep Learning, etc. etc.)... but I'm definitely interested in the low-latency uses of SIMD.
For example: Consider a SIMD-based Bloom Filters to augment a Hash Table (or any data-structure, really). Multi-hash systems like Cuckoo-hashes have also been successfully implemented with CPU-SIMD acceleration.
You can use CPU-accelerated SIMD to accelerate CPU-based algorithms. And CPUs certainly can benefit from fatter, and better designed, SIMD instruction sets.
unfortunately, it takes millions of cycles for the CPU to switch to the higher-power state necessary for AVX instructions, so your latency really isn't that much better.
the radfft guy did a good writeup https://gist.github.com/rygorous/32bc3ea8301dba09358fd2c64e0...
But I've found getting into GPU computing extremely difficult. Admittedly, part of that may be having bought an AMD Vega graphics card instead of NVidea (not a fan of vendor locking, and AMD is building an open source stack with ROCm/HIP).
My code does lots of different things, many small, but given a billion iterations, it can add up. If you're running a Monte Carlo simulation over millions of different data sets, fitting each iteration with Markov Chain Monte Carlo, with a model that has a few long for loops and requires (automatic) differentiation and some number of special functions... It's not hard to rack up slow execution times out of small pieces.
Part of the problem is I need to actually find a GPU project that's accessible for me, so I can gain some experience, and learn what's actually possible. How I could actually organize computation flow. It's easy with a CPU. Break up Monte Carlo simulations among processors. Depending on the model being fit, as well as the method, different computations can vary dramatically in width. But optimizing is as simple as breaking up all those computations into the appropriately sized AVX units (and libraries + compilers auto-vectorization will often handle most of that!), so wider units directly translate to faster performance.
Part of the problem is I don't really know how to think about a GPU. Can I think of the Vega GPU as a processor with 64 cores, a SIMD width of 64 (64 * 64 = advertised "4096), and more efficient gather loads/scatter store operations?
If there were a CPU like that, compiler + library calls and macros (including libraries you've written or wrapped yourself) go a long way so you can quickly write well optimized code.
I really need to dedicate the time to learn more. My 7900x processor is about 4x faster at gemm than my 1950x, but 10x slower than the Vega graphics card that cost less money. I see the potential.
To start, I need to find a GPU project that lets me really get experimenting and figure out how programs can even look.
Is there an Agner Fog of GPGPU?
From my limited experience, that's the wrong way of looking at things.
SIMD is good at looking at "converged" instruction streams, and bad at "divergent" instruction streams.
"Converged" is when you've got your bitmask (in AMD / OpenCL: the scalar bitmask register) executing all 64-threads as much as possible. This means all 64-threads take "if" statements together, they loop together, etc. etc.
"Diverged" is when you have a tree of if-statements, and different threads take different paths. The SIMD processor has to execute both the "Then" AND the "Else" statements. And if you have a tree of "if / else" statements, the threads have less-and-less in common.
-------
That's it. You try to group code together such that as much code converges as possible, and diverges as little as possible. It helps to know that the work-group size on AMD can be anywhere from 64 to 256 (so keep your "thread groups" together as large as 256 at a time).
The OpenCL compiler will automatically translate most code into conditional moves and do its best to avoid "real if/else branches". But as the programmer, you just have to realize that nested if-statements and nested-loops can cause thread-divergence.
---------
Matrix maths are the easiest for SIMD, because of how regular their structure is. When you get to "real" cases like... well...
> My code does lots of different things, many small, but given a billion iterations, it can add up. If you're running a Monte Carlo simulation over millions of different data sets, fitting each iteration with Markov Chain Monte Carlo, with a model that has a few long for loops and requires (automatic) differentiation and some number of special functions... It's not hard to rack up slow execution times out of small pieces.
Okay, well that's why its hard to program properly on a GPU. Because they're all diverging. Unless you figure out a way to "converge" these if statements, it won't run on a GPU correctly.
Chances are, if you have a million-wide Monte-carlo, a lot of those threads are going to be converging. Can you split up the steps and "regroup" tasks as appropriate?
Lets say your million threads run into a switch statement:
* 106,038 threads will take path A.
* 348,121 threads take path B.
* 764 threads take path C.
Etc. etc.
Can you "regroup" so that all your threads in pathA start to execute together again? To best take advantage of the SIMD-architecture?
Think about it: a gang of 64 may have thread #0 take PathA, thread#1 take PathB, thread #2 take PathA, etc. etc. So it all diverges and you lose all your SIMD.
But if you "rebatch" everything together... eventually you'll get hundreds of threads ready for PathA. At which point, you gang up your group of 64 threads and execute PathA all together again, taking full advantage of SIMD.
SIMD is an architecture that allows "similar" threads to run at the same time, in groups of 64 or more (up to 256 threads at once). And they truly execute simultaneously as long as they are all taking the same if/else branches and loops. The real hard part is designing your data-structures to maximize this "convergence" of threads.
Something like Chess would practically be impossible to SIMD on a GPU. Its just impossible to handle how divergent the typical chess analysis engine gets due to the precision of chess board positions. But something like Monte-Carlo surely will have hundreds, or thousands, of threads taking any "particular execution path". So I'd have hope for a SIMD-based execution on Monte-Carlo simulations.
You still need some understanding of gpus, but at least you can write your code in Julia rather than CUDA or OpenCL.
Concretely, I’ve used this for parameter sweeps of large systems of ODEs and SDEs.
Maybe. Sparsity of the matrix in question matters a lot. Matrixes with very little sparsity, like say an image, do well. Others may not. It’s more like GPUs do well on algorithms with predictable branching.
The future will probably include more avx512, but not yet.
Video Editing, Video Games, Graphics, 3d Modeling, Photoshop. Even Stockfish Chess uses new instructions (not SIMD: but the Bit-board popcnt and pext / pdep instructions) to accelerate chess computations.
At a bare minimum, AVX grossly accelerates memcpy and memset operations. (Setting 256-bits per assembly instruction instead of 64-bits per operation is a big improvement). And virtually every program can benefit from faster memcpy and faster memsets.
"Standard Software" isn't written to be very fast. But anything that's even close to CPU-bound is being upgraded to use more and more SIMD instructions.
How often are memcpy and memset CPU-bound, though?
L3-cache on Skylake architectures has 18-bytes-per-clock sustained, which is still greater than 128-bits per clock bandwidth to your L3 cache.
Soooo... I'd say roughly for any memset or memcpy smaller than 256kB or so... you're going to benefit to an AVX-based memset or memcpy.
At ~8MB or so, where you're hitting L3 cache, it probably is a benefit but not much of one. The important thing is Skylake only has one store unit, so you can only write once per clock cycle.
So... do you write one 64-bit value, or do you write one 256-bit (32-byte) value?
----------
To be fair: I'm pretty sure every compiler's default settings today outputs SSE-based memcpy and memsets (128-bit). So AVX doubles that to 256-bit.
Still, it seems like the documentation elsewhere says that AVX 128-bit doesn't cause any clocking issues. So AVX512 (applied to 128-bit) should still lead to good speeds without any clocking problems.
void memcopy(double* restrict a, double* restrict b, long N){
for (long n; n < N; n++){
a[n] = b[n];
}
}
You can check assembly with: gcc -march=skylake-avx512 -O2 -ftree-vectorize -mprefer-vector-width=128
-shared -fPIC -S memcopy.c -o memcopy128.s
gcc -march=skylake-avx512 -O2 -ftree-vectorize -mprefer-vector-width=512
-shared -fPIC -S memcopy.c -o memcopy512.s
to confirm they use `xmm` and `zmm` registers, respectively. (The 512 version also uses xmm to finish off a remainder. Ideally, it would use masked load/stores, but I haven't seen any auto-vectorizers actually do that.I'm using gcc 8.2.1. I believe the option `-mprefer-vector-width` was added recently.
Now I compiled both into shared libraries (drop the `-S` and choose an appropriate file name), and benchmarked from within Julia.
julia> using BenchmarkTools, Random
julia> memcopy128!(a, b) = ccall((:memcopy,
"/home/chriselrod/Documents/progwork/C/libmemcopy128.so"),
Cvoid, (Ptr{Cdouble},Ptr{Cdouble},Clong), pointer(a), pointer(b), length(a));
julia> memcopy512!(a, b) = ccall((:memcopy,
"/home/chriselrod/Documents/progwork/C/libmemcopy512.so"),
Cvoid, (Ptr{Cdouble},Ptr{Cdouble},Clong), pointer(a), pointer(b), length(a));
julia> b = randn(32); a = similar(b);
julia> b = randn(32); a = similar(b);
julia> all(a == b)
false
julia> memcopy128!(a, b); all(a == b)
true
julia> randn!(a); all(a == b)
false
julia> memcopy512!(a, b); all(a == b)
true
julia> @btime memcopy128!($a, $b)
7.517 ns (0 allocations: 0 bytes)
julia> @btime memcopy512!($a, $b)
5.277 ns (0 allocations: 0 bytes)
julia> b = randn(64); a = similar(b);
julia> @btime memcopy128!($a, $b)
10.980 ns (0 allocations: 0 bytes)
julia> @btime memcopy512!($a, $b)
7.300 ns (0 allocations: 0 bytes)
julia> b = randn(100); a = similar(b);
julia> @btime memcopy128!($a, $b)
15.060 ns (0 allocations: 0 bytes)
julia> @btime memcopy512!($a, $b)
14.031 ns (0 allocations: 0 bytes)
julia> b = randn(200); a = similar(b);
julia> @btime memcopy128!($a, $b)
33.867 ns (0 allocations: 0 bytes)
julia> @btime memcopy512!($a, $b)
19.546 ns (0 allocations: 0 bytes)
julia> b = randn(400); a = similar(b);
julia> @btime memcopy128!($a, $b)
53.923 ns (0 allocations: 0 bytes)
julia> @btime memcopy512!($a, $b)
28.376 ns (0 allocations: 0 bytes)
julia> b = randn(800); a = similar(b);
julia> @btime memcopy128!($a, $b)
98.277 ns (0 allocations: 0 bytes)
julia> @btime memcopy512!($a, $b)
57.960 ns (0 allocations: 0 bytes)
julia> b = randn(2_000); a = similar(b);
julia> @btime memcopy128!($a, $b)
232.138 ns (0 allocations: 0 bytes)
julia> @btime memcopy512!($a, $b)
133.229 ns (0 allocations: 0 bytes)
julia> b = randn(20_000); a = similar(b);
julia> @btime memcopy128!($a, $b)
3.688 μs (0 allocations: 0 bytes)
julia> @btime memcopy512!($a, $b)
2.759 μs (0 allocations: 0 bytes)
julia> b = randn(200_000); a = similar(b);
julia> @btime memcopy128!($a, $b)
112.598 μs (0 allocations: 0 bytes)
julia> @btime memcopy512!($a, $b)
117.437 μs (0 allocations: 0 bytes)
The advantage wasn't very impressive here, but it persisted for vectors of length 20,000. At 200,000 we saw the memory bottleneck hit in force, when it busts the per core L2 cache. The 7900x's per core L2 cache is 1048576, which translates to 131,072 doubles.Of course there's a performance advantage from using 512 bit registers for a memcpy - but a memcpy is rarely a major performance bottleneck by itself and is usually surrounded by other code. Unless that code is also AVX-512, you've just made it slower by optimizing the memcpy. My point was that a compiler can't usually decide whether it's worth making the optimization in light of the broader context.
The other point was whether using AVX-512 while sticking to xmm registers is faster than just using xmm SSE/AVX code. I don't have an AVX-512 capable machine at the moment, perhaps you'd like to check if your 128 bit version is any faster than just doing "gcc -march=skylake -O2 -mprefer-vector-width=128" (thereby retaining the microarchitecture optimizations, but sticking to AVX2)?
Not necessarily, for large operations and depending on processor generation a simple rep stos/movsb will be simpler (no alignment requirements) and saturate your memory bandwidth just as well as any AVX sequence will with less icache pressure.
Seems to be a feature in Ivy Bridge and later, which happens to be around the time AVX2 started.