8 ms·

BLAS-level CPU Performance in 100 Lines of C

216 points by henning 5y ago

pstrike 5y ago
cool approach to implement matrix multiplication
nyanpasu64 5y ago
Is OpenBLAS faster for large matrices because of optimizations or asymptotically faster algorithms? https://stackoverflow.com/a/11421344 https://stackoverflow.com/a/11421344 claims BLAS doesn't use Strassen's algorithm.
- jleahy 5y ago
  Due to optimisations similar to those in this article. The asymptotically faster algorithms are not useful in practice.
- quanto 5y ago
  Strassen's algorithm is not numerically stable and thus is not normally used in floating point operations. Usually, the performance differences between numerical libraries, especially for something as routine as matrix computation, are due to engineering implementations (unless explicitly stated otherwise).
- sampo 5y ago
  I've read that Strassen's algorithm, and fast matrix multiplication algorithms in general, accumulate too much numerical errors to be useful with single precision (f32) floats. That's why they were not in use before year 2000. With f64 it's better, and in the past 20 years there's been research on how to augment the fast algorithms to decrease the numerical errors. But those extra steps also make the fast algorithms slower. I don't know what the current situation is.
olliej 5y ago
The core of their multiplication is in https://github.com/HazyResearch/blocking-tutorial/blob/master/matrix_kernel_vectorized.cpp#L64 https://github.com/HazyResearch/blocking-tutorial/blob/maste... Where they implement the matrix inner product. And their example implementations using that start at https://github.com/HazyResearch/blocking-tutorial/blob/master/compare_blas.cpp#L90 https://github.com/HazyResearch/blocking-tutorial/blob/maste... The inner multiply is all AVX2(w/FMA extension) so it's probably making some degree of LoC win from not needing to handle hardware capabilities or being portable, but still it is impressive. It is also worth noting that the inner product is c++ templates - it could probably be turned into a horrific macro to make it pure C, but I'd say it isn't worth it. But then I prefer C++ to C :D
raverbashing 5y ago
"It can be done in 100 lines of C" sure, but if I use the library that's 100 lines less that I have to think about and write While it is good to understand the techniques to make it faster, for most cases using BLAS is probably better
- cyber_kinetist 5y ago
  OpenBLAS is quite a heavy dependency and can be quite cumbersome to integrate into your build system. In that case that 100-liner can be really useful if you just need some large-matrix multiplication somewhere in your code.
- Sesse__ 5y ago
  On the other hand: If you have those 100 lines, it's one library less you have to think about and carry around.
trombonechamp 5y ago
Once I needed a single LAPACK function (tridiagonal matrix multiplication) in a C Python extension. I spent an eternity trying to link to LAPACK reliably across platforms until giving up and deciding to implement it myself. My implementation ended up being about 10 lines of code and ran slightly faster (!!!) than the LAPACK function. (I think this is because LAPACK provided numerical stability for special cases which didn't apply to my use case.)
- stephencanon 5y ago
  It’s because tridiagonal multiplication doesn’t do enough work to benefit from most of the BLAS techniques discussed here. You can’t do non-trivial vectorization or cache blocking or threading because there simply isn’t very much work to be done unless your matrix is enormous, so a lot of BLAS implementations will basically just use scalar code for this operation.
  - stephencanon 5y ago
    That said, it _is_ optimizable, it's just that: - the gains are the sort of typical 2-8x speed improvements from vectorization, not the multiple-orders-of-magnitude gains that you can get on dense GEMM. - the absolute number of flops performed is O(n^2) rather than O(n^3) for GEMM, so even if you could make tridiagonal operations infinitely fast, that optimization effort would be better spent on even small speedups to the O(n^3) work that probably comprises other parts of your algorithm.
miraculixx 5y ago
Sure. Many a thing that we use libraries for are 100 lines or less. It's just if you have 100 times 100 lines that's O(10000) more work. Ups.
- ant6n 5y ago
  Well, O(10000) = O(1), so it's not really a problem.
  - miraculixx 5y ago
    :-) yeah I was kind of expecting that comment
eachro 5y ago
There's "knowing C syntax" and then there's knowing C well enough to be able to read this code and understand what's going on. One of these days I'd like to get to the point where I can conceivably read this code and understand it well enough to build my own matrix library using similar optimizations for the rest of the matrix operation primitives. Some day ...
- saagarjha 5y ago
  This code is less about knowing C and more about knowing vector intrinsics which are not standard C. The actual code is mostly straightforward constructs.
  - qayxc 5y ago
    Interestingly knowing vector intrinsics is exactly what optimising compilers are good at. The article even shows that - look at the C++ code - no instrisics in sight, just pure basic C++ getting the optimisation and performance from the compiler alone.
    - saagarjha 5y ago
      You probably want to look at the kernels themselves: https://github.com/HazyResearch/blocking-tutorial/blob/master/matrix_kernel16_6.c https://github.com/HazyResearch/blocking-tutorial/blob/maste...
      - qayxc 5y ago
        I was referring to https://github.com/HazyResearch/blocking-tutorial/blob/master/compare_blas.cpp https://github.com/HazyResearch/blocking-tutorial/blob/maste... which upon closer inspection uses intrinsics in the inner loop, but to a much lesser degree and arguably specifying -mavx2 -mfma compiler options would yield comparable results...
        saagarjha 5y ago
GistNoesis 5y ago
There is something, I found quite liberating, was writing matrix operations in the naive way. Forget about performance for a moment, and just write those imbricated loops. They can run everywhere from javascript in the web, to your raspberry pi, to your arduino where you want to have a linear classifier, windows, linux, embedded platform. No need for complex to compile library. When you realize that most deep learning library, "einsum is all you need", and that einsum is just loops and sum, it kind of screams you are doing it the wrong way when using numerical libraries ranging from BLAS to cuda. Loops can also be automatically differentiated. Loops can also be reordered. Loops can be tiled. There are some techniques (polyhedral loop optimization) to do it automatically. Loops can be parallelized and vectorized. Most of the optimization are always the same and due to the memory hierarchy. The problem space is well defined. Candidates solutions can be benchmarked automatically. Some architectures like to have their memory access in a Structures of Arrays, other in Array of Structures, or even Arrays of structures of Arrays. Benchmark depend of sizes in various dimensions. Are you really always going to write the optimizations manually, and do the complex code transforms by hand. This feels like one area ripped for disruption. The other day there was this thread about not writing your own computer language, but here one domain specific language can have a huge impact. Take the naive form and automatically generate the forward, backward, (nth order autodiff) optimized version given an instruction set whether for cpu, gpu or tpu, and for various memory usage profile. Optimizing code, even though it feels quite fun and rewarding to do manually is a job for compilers. Write naive algorithms but smart compilers. Extra kudos if you handle the mathematical function : For example if you want an exponential : write it as an infinite taylor sum. If you want matrix exponential write it as a taylor sum also ( and if your compiler automatically derive the Pade approximation for improved numerical stability you know you are on the right track).
- MaxBarraclough 5y ago
  > Optimizing code, even though it feels quite fun and rewarding to do manually is a job for compilers. > Write naive algorithms but smart compilers. I'm not sure where you're coming from here. The article's charts show that in this case, naive code gets you less than 2% the performance of an optimised implementation.
d_tr 5y ago
I wrote a matrix multiplication routine on a Haswell a few years ago, for no particular reason. I chose an 8x8 minimum block size since it seemed like the most convenient one, and I had one or two more block sizes for cache utilization IIRC. It was a few tens of lines of code and the performance in GFLOPS was 80% of peak which was actually a pleasant surprise. I think it was because of fmaddps latency. My code would like four but on Haswell it is five cycles. So the other types of instructions (prefetch, broadcast, loops) basically did not have any impact on speed.
Traster 5y ago
I guess that's an thorough way of saying "The problem with building a library is supporting all the use cases of a library, not solving the simple underlying problem". Sure, you've got BLAS level performance on this one CPU, but you've literally hard coded stuff according to the exact cache size of each of your caches and the exact low level ISA extensions your CPU has. I think it would have been far more interesting to take this and then put it on a Threadripper to see "Given I've optimized the hell out of this for one CPU, how portable is it".
- d_tr 5y ago
  True, and this operation is close to as simple as it gets.
- paulsutter 5y ago
  The solution is code generation based on CPU parameters, not general portable code
- jcranberry 5y ago
  Cache sizes can be gotten programatically. I don't see that as an impediment. Would need to be benchmarked of course to see how that affects performance.
  - henrydark 5y ago
    Eigen does this
  - mkbosmans 5y ago
    Sure, but that doesn't change the fact that you cannot simply plug that into a simple formula to get the e.g. the optimal block size for cache-blocking. To be more precise: you can try to figure out how cache size would affect optimal block size in theory and in practice this gets you pretty decent performance. But to get the last bit of performance out of the hardware you need to benchmark several configurations on the hardware and select the best one based on run time parameters like u-arch and cache size.
throwamon 5y ago
How is Julia doing in that regard?
- adgjlsfhk1 5y ago
  Really well with LoopVectorization The following code is 4x faster than openblas for small matrices, and 10% faster at 1000x1000. It still loses to MKL at larger sizes, but if you look at Octavian.jl, it's a pure Julia library that beats or is equal to MKL at pretty much all sizes. using LoopVectorization function mul!(C,A,B) @tturbo for i in 1:size(C,1) for j in 1:size(C,2) C[i,j] = zero(eltype(C)) for k in 1:size(A,1) C[i,j] += A[i,k]*B[k,j] end end end return C end
kazinator 5y ago
If you're "using Intel intrinsics", it's not C. It's effectively a form of inline assembly. BLAS-level performance in C means you write ISO C, and the compiler spits out the required instructions.
- bigdict 5y ago
  Any pointers to docs showing how to get GCC/Clang to emit vectored code based on scalar C source?
  - deleted 5y ago
    [deleted]
  - gnufx 5y ago
    -ftree-vectorize, but you'll typically want something like -Ofast to allow lax numerics; and -mcpu=native, of course. You also want to ensure you use C "restrict" appropriately, though GCC will multi-version for possible aliasing. That's not the real problem with level 3 BLAS, though, it's taking account of the micro-architecture and matrix dimensions with appropriate blocking and packing. You can get >80% of hand-optimized performance with plain C [1]. Leaving pure C, you should be able to improve that with prefetching, and I should try to figure why that didn't work in my tests. 1. https://fx.srht.site/posts/2022-01-03-pure-c-gemm.html https://fx.srht.site/posts/2022-01-03-pure-c-gemm.html Edit: You probably want -funroll-loops too.
    - moonchild 5y ago
      > -Ofast to allow lax numerics Please don't do this if you want your results to be reproducible.
      - stephencanon 5y ago
        -Ofast doesn't license any transforms that can break GEMM. Everything that it licenses is already used universally by GEMM kernels, because those transforms have no effect on the standard backwards stability bounds for matrix multiplication.
twoodfin 5y ago
This reminds me of a well-read paper some years ago that was titled something like “The Missing 100X”, on the disparity between raw CPU computational throughput gains and measured software performance + an evaluation of techniques to close it. My Google-fu has failed me, however.
- deleted 5y ago
  [deleted]
nspattak 5y ago
I would be very happy if we lived in a world where a code using SIMD instructions and hand crafted cache-blocking and threaded code was a "naive implementation" but we are not. a "naive implementation" is a simple nested loop. The rest are complex and they become even more complex when one tries to implement them in a way which compile efficiently on multiple CPUs (from different caches to different architectures). Interesting article but the authors' conclusion is completely misguided.
mainland 5y ago
An interesting, older paper that explores how to bridge the "ninja gap" between naive code and hand-optimized code is Can Traditional Programming Bridge the Ninja Performance Gap for Parallel Computing Applications? <https://web.eecs.umich.edu/~msmelyan/papers/isca-2012-paper.pdf https://web.eecs.umich.edu/~msmelyan/papers/isca-2012-paper....>. They reduce the ninja gap for a range of benchmarks from 24x to 1.3x by relying on "smart" compilers and using only basic program transformations---no intrinsics!
gnufx 5y ago
I'm unconvinced by the need for intrinsics in many cases (other than possibly for prefetch). I took out those out of the final version of [1] and got code that actually ran faster on a contemporary processor, targeting the same old SIMD. 1. https://github.com/flame/how-to-optimize-gemm/tree/master/src https://github.com/flame/how-to-optimize-gemm/tree/master/sr...
eternityforest 5y ago
Pls don't actually do this in production though. It won't support other CPUs, nor newer or older opcodes on the same CPU, unless you add that too, but that's just adding fragmentation.

216 points by henning 5y ago

cool approach to implement matrix multiplication

Is OpenBLAS faster for large matrices because of optimizations or asymptotically faster algorithms? https://stackoverflow.com/a/11421344 https://stackoverflow.com/a/11421344 claims BLAS doesn't use Strassen's algorithm.

jleahy 5y ago
Due to optimisations similar to those in this article. The asymptotically faster algorithms are not useful in practice.
quanto 5y ago
Strassen's algorithm is not numerically stable and thus is not normally used in floating point operations. Usually, the performance differences between numerical libraries, especially for something as routine as matrix computation, are due to engineering implementations (unless explicitly stated otherwise).
sampo 5y ago
I've read that Strassen's algorithm, and fast matrix multiplication algorithms in general, accumulate too much numerical errors to be useful with single precision (f32) floats. That's why they were not in use before year 2000. With f64 it's better, and in the past 20 years there's been research on how to augment the fast algorithms to decrease the numerical errors. But those extra steps also make the fast algorithms slower. I don't know what the current situation is.

olliej 5y ago

The core of their multiplication is in https://github.com/HazyResearch/blocking-tutorial/blob/master/matrix_kernel_vectorized.cpp#L64 https://github.com/HazyResearch/blocking-tutorial/blob/maste... Where they implement the matrix inner product. And their example implementations using that start at https://github.com/HazyResearch/blocking-tutorial/blob/master/compare_blas.cpp#L90 https://github.com/HazyResearch/blocking-tutorial/blob/maste... The inner multiply is all AVX2(w/FMA extension) so it's probably making some degree of LoC win from not needing to handle hardware capabilities or being portable, but still it is impressive. It is also worth noting that the inner product is c++ templates - it could probably be turned into a horrific macro to make it pure C, but I'd say it isn't worth it. But then I prefer C++ to C :D

raverbashing 5y ago

"It can be done in 100 lines of C" sure, but if I use the library that's 100 lines less that I have to think about and write While it is good to understand the techniques to make it faster, for most cases using BLAS is probably better

cyber_kinetist 5y ago
OpenBLAS is quite a heavy dependency and can be quite cumbersome to integrate into your build system. In that case that 100-liner can be really useful if you just need some large-matrix multiplication somewhere in your code.
Sesse__ 5y ago
On the other hand: If you have those 100 lines, it's one library less you have to think about and carry around.

trombonechamp 5y ago

Once I needed a single LAPACK function (tridiagonal matrix multiplication) in a C Python extension. I spent an eternity trying to link to LAPACK reliably across platforms until giving up and deciding to implement it myself. My implementation ended up being about 10 lines of code and ran slightly faster (!!!) than the LAPACK function. (I think this is because LAPACK provided numerical stability for special cases which didn't apply to my use case.)

stephencanon 5y ago
It’s because tridiagonal multiplication doesn’t do enough work to benefit from most of the BLAS techniques discussed here. You can’t do non-trivial vectorization or cache blocking or threading because there simply isn’t very much work to be done unless your matrix is enormous, so a lot of BLAS implementations will basically just use scalar code for this operation.
- stephencanon 5y ago
  That said, it _is_ optimizable, it's just that: - the gains are the sort of typical 2-8x speed improvements from vectorization, not the multiple-orders-of-magnitude gains that you can get on dense GEMM. - the absolute number of flops performed is O(n^2) rather than O(n^3) for GEMM, so even if you could make tridiagonal operations infinitely fast, that optimization effort would be better spent on even small speedups to the O(n^3) work that probably comprises other parts of your algorithm.

miraculixx 5y ago

Sure. Many a thing that we use libraries for are 100 lines or less. It's just if you have 100 times 100 lines that's O(10000) more work. Ups.

ant6n 5y ago
Well, O(10000) = O(1), so it's not really a problem.
- miraculixx 5y ago
  :-) yeah I was kind of expecting that comment

eachro 5y ago

There's "knowing C syntax" and then there's knowing C well enough to be able to read this code and understand what's going on. One of these days I'd like to get to the point where I can conceivably read this code and understand it well enough to build my own matrix library using similar optimizations for the rest of the matrix operation primitives. Some day ...

saagarjha 5y ago
This code is less about knowing C and more about knowing vector intrinsics which are not standard C. The actual code is mostly straightforward constructs.
- qayxc 5y ago
  Interestingly knowing vector intrinsics is exactly what optimising compilers are good at. The article even shows that - look at the C++ code - no instrisics in sight, just pure basic C++ getting the optimisation and performance from the compiler alone.
  - saagarjha 5y ago
    You probably want to look at the kernels themselves: https://github.com/HazyResearch/blocking-tutorial/blob/master/matrix_kernel16_6.c https://github.com/HazyResearch/blocking-tutorial/blob/maste...
    - qayxc 5y ago
      I was referring to https://github.com/HazyResearch/blocking-tutorial/blob/master/compare_blas.cpp https://github.com/HazyResearch/blocking-tutorial/blob/maste... which upon closer inspection uses intrinsics in the inner loop, but to a much lesser degree and arguably specifying -mavx2 -mfma compiler options would yield comparable results...
      - saagarjha 5y ago

GistNoesis 5y ago

There is something, I found quite liberating, was writing matrix operations in the naive way. Forget about performance for a moment, and just write those imbricated loops. They can run everywhere from javascript in the web, to your raspberry pi, to your arduino where you want to have a linear classifier, windows, linux, embedded platform. No need for complex to compile library. When you realize that most deep learning library, "einsum is all you need", and that einsum is just loops and sum, it kind of screams you are doing it the wrong way when using numerical libraries ranging from BLAS to cuda. Loops can also be automatically differentiated. Loops can also be reordered. Loops can be tiled. There are some techniques (polyhedral loop optimization) to do it automatically. Loops can be parallelized and vectorized. Most of the optimization are always the same and due to the memory hierarchy. The problem space is well defined. Candidates solutions can be benchmarked automatically. Some architectures like to have their memory access in a Structures of Arrays, other in Array of Structures, or even Arrays of structures of Arrays. Benchmark depend of sizes in various dimensions. Are you really always going to write the optimizations manually, and do the complex code transforms by hand. This feels like one area ripped for disruption. The other day there was this thread about not writing your own computer language, but here one domain specific language can have a huge impact. Take the naive form and automatically generate the forward, backward, (nth order autodiff) optimized version given an instruction set whether for cpu, gpu or tpu, and for various memory usage profile. Optimizing code, even though it feels quite fun and rewarding to do manually is a job for compilers. Write naive algorithms but smart compilers. Extra kudos if you handle the mathematical function : For example if you want an exponential : write it as an infinite taylor sum. If you want matrix exponential write it as a taylor sum also ( and if your compiler automatically derive the Pade approximation for improved numerical stability you know you are on the right track).

MaxBarraclough 5y ago
> Optimizing code, even though it feels quite fun and rewarding to do manually is a job for compilers. > Write naive algorithms but smart compilers. I'm not sure where you're coming from here. The article's charts show that in this case, naive code gets you less than 2% the performance of an optimised implementation.

d_tr 5y ago

I wrote a matrix multiplication routine on a Haswell a few years ago, for no particular reason. I chose an 8x8 minimum block size since it seemed like the most convenient one, and I had one or two more block sizes for cache utilization IIRC. It was a few tens of lines of code and the performance in GFLOPS was 80% of peak which was actually a pleasant surprise. I think it was because of fmaddps latency. My code would like four but on Haswell it is five cycles. So the other types of instructions (prefetch, broadcast, loops) basically did not have any impact on speed.

Traster 5y ago

I guess that's an thorough way of saying "The problem with building a library is supporting all the use cases of a library, not solving the simple underlying problem". Sure, you've got BLAS level performance on this one CPU, but you've literally hard coded stuff according to the exact cache size of each of your caches and the exact low level ISA extensions your CPU has. I think it would have been far more interesting to take this and then put it on a Threadripper to see "Given I've optimized the hell out of this for one CPU, how portable is it".

d_tr 5y ago
True, and this operation is close to as simple as it gets.
paulsutter 5y ago
The solution is code generation based on CPU parameters, not general portable code
jcranberry 5y ago
Cache sizes can be gotten programatically. I don't see that as an impediment. Would need to be benchmarked of course to see how that affects performance.
- henrydark 5y ago
  Eigen does this
- mkbosmans 5y ago
  Sure, but that doesn't change the fact that you cannot simply plug that into a simple formula to get the e.g. the optimal block size for cache-blocking. To be more precise: you can try to figure out how cache size would affect optimal block size in theory and in practice this gets you pretty decent performance. But to get the last bit of performance out of the hardware you need to benchmark several configurations on the hardware and select the best one based on run time parameters like u-arch and cache size.

throwamon 5y ago

How is Julia doing in that regard?

adgjlsfhk1 5y ago
Really well with LoopVectorization The following code is 4x faster than openblas for small matrices, and 10% faster at 1000x1000. It still loses to MKL at larger sizes, but if you look at Octavian.jl, it's a pure Julia library that beats or is equal to MKL at pretty much all sizes. using LoopVectorization function mul!(C,A,B) @tturbo for i in 1:size(C,1) for j in 1:size(C,2) C[i,j] = zero(eltype(C)) for k in 1:size(A,1) C[i,j] += A[i,k]*B[k,j] end end end return C end

kazinator 5y ago

If you're "using Intel intrinsics", it's not C. It's effectively a form of inline assembly. BLAS-level performance in C means you write ISO C, and the compiler spits out the required instructions.

bigdict 5y ago
Any pointers to docs showing how to get GCC/Clang to emit vectored code based on scalar C source?
- deleted 5y ago
  [deleted]
- gnufx 5y ago
  -ftree-vectorize, but you'll typically want something like -Ofast to allow lax numerics; and -mcpu=native, of course. You also want to ensure you use C "restrict" appropriately, though GCC will multi-version for possible aliasing. That's not the real problem with level 3 BLAS, though, it's taking account of the micro-architecture and matrix dimensions with appropriate blocking and packing. You can get >80% of hand-optimized performance with plain C [1]. Leaving pure C, you should be able to improve that with prefetching, and I should try to figure why that didn't work in my tests. 1. https://fx.srht.site/posts/2022-01-03-pure-c-gemm.html https://fx.srht.site/posts/2022-01-03-pure-c-gemm.html Edit: You probably want -funroll-loops too.
  - moonchild 5y ago
    > -Ofast to allow lax numerics Please don't do this if you want your results to be reproducible.
    - stephencanon 5y ago
      -Ofast doesn't license any transforms that can break GEMM. Everything that it licenses is already used universally by GEMM kernels, because those transforms have no effect on the standard backwards stability bounds for matrix multiplication.

twoodfin 5y ago

This reminds me of a well-read paper some years ago that was titled something like “The Missing 100X”, on the disparity between raw CPU computational throughput gains and measured software performance + an evaluation of techniques to close it. My Google-fu has failed me, however.

deleted 5y ago
[deleted]

nspattak 5y ago

I would be very happy if we lived in a world where a code using SIMD instructions and hand crafted cache-blocking and threaded code was a "naive implementation" but we are not. a "naive implementation" is a simple nested loop. The rest are complex and they become even more complex when one tries to implement them in a way which compile efficiently on multiple CPUs (from different caches to different architectures). Interesting article but the authors' conclusion is completely misguided.

mainland 5y ago

An interesting, older paper that explores how to bridge the "ninja gap" between naive code and hand-optimized code is Can Traditional Programming Bridge the Ninja Performance Gap for Parallel Computing Applications? <https://web.eecs.umich.edu/~msmelyan/papers/isca-2012-paper.pdf https://web.eecs.umich.edu/~msmelyan/papers/isca-2012-paper....>. They reduce the ninja gap for a range of benchmarks from 24x to 1.3x by relying on "smart" compilers and using only basic program transformations---no intrinsics!

gnufx 5y ago

I'm unconvinced by the need for intrinsics in many cases (other than possibly for prefetch). I took out those out of the final version of [1] and got code that actually ran faster on a contemporary processor, targeting the same old SIMD. 1. https://github.com/flame/how-to-optimize-gemm/tree/master/src https://github.com/flame/how-to-optimize-gemm/tree/master/sr...

eternityforest 5y ago

Pls don't actually do this in production though. It won't support other CPUs, nor newer or older opcodes on the same CPU, unless you add that too, but that's just adding fragmentation.