BLAS-level CPU Performance in 100 Lines of C | Hacker News Reader

BLAS-level CPU Performance in 100 Lines of C

cs.stanford.edu

216 pointshenning62 comments

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".

True, and this operation is close to as simple as it gets.

The solution is code generation based on CPU parameters, not general portable code

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.

Eigen does this

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.

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.)

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.

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.

The core of their multiplication is in

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/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

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.

Any pointers to docs showing how to get GCC/Clang to emit vectored code based on scalar C source?

-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

Edit: You probably want -funroll-loops too.

> -Ofast to allow lax numerics

Please don't do this if you want your results to be reproducible.

-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.

Indeed, the BLIS test suite with generic C kernels passes, not just level 3, though only the -funsafe-math-optimizations of -ffast-math is used by default. (I don't know the specifics for OpenBLAS, but it must be similar.)

At least some of BLAS was defined before there even was an IEEE floating point standard.

Your life will be a lot easier if you decide not to chase bit-level reproducibility when doing floating point calculations. You can try setting CFLAGS to get the same result on all hardware, but that is a lot of effort, especially if you don't want to leave any performance on the table.

And for what gain, just to make it simpler for your tests to compare against a fixed reference? In practice you don't gain anything else by having reproducible calculations because your inputs will have a lot more uncertainty in them than just beyond the 16th digit.

Better to embrace the difference in results between hardware architecture, parallel execution, compiler version, optimization level, etc. This means spending time making your test suite more robust to small variations, but it saves you time investigating failed test runs that are insignificant. And as a bonus, it forces you to understand the numerical properties of your code better.

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 ...

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.

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.

You probably want to look at the kernels themselves: https://github.com/HazyResearch/blocking-tutorial/blob/maste...

I was referring to 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...

That’s the timing harness for the more advanced implementations, I think.

Yes, GCC fully vectorizes (and unrolls) a micro-kernel loop nest.

Each architecture has a "library" of these bizarrely named functions that are mapped by the compiler to CPU-specific instructions. I'd suggest 90%+ of C programmers couldn't name a single one of them, and that's fine, because they have no place in C outside of tight platform-specific optimization work, which very few people have a need for or mandate to actually work on

Intrinsics-heavy code does looks scary because you see tons of weirdly-named functions with no clue as to what they do.

However, you can learn the basics in an hour, and read/write small amounts of code (like in the article) very quickly. The trick that you don't learn all those functions, instead you constantly use the amazing interactive guide:

https://www.intel.com/content/www/us/en/docs/intrinsics-guid...

I guess real pros do memorize the API of course, but I have written many (<1000) lines of intrinsics code like this.

PSA the database of instructions and intrinsics is in this xml file

https://www.intel.com/content/dam/develop/public/us/en/inclu...

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

Due to optimisations similar to those in this article. The asymptotically faster algorithms are not useful in practice.

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).

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.

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.

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.

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.

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).

> 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.

> The article's charts show that in this case, naive code gets you less than 2% the performance of an optimised implementation.

That's not what the chart shows at all. The charts shows what non optimized code compiles to.

Just run the code through any optimizing compiler released in the past decade and you'll find that do the vectorization for you. No need to use non-portable intrinsics.

Even the block optimisation can be done in plain C.

While optimised algorithms are better than naive algorithms, optimising machine code is indeed a compiler task. Some languages tend to do a better job of describing the algorithm than others, but that's a different can worms entirely.

I think that by non-optimized the chart refers to not hand-optimized. The code is likely still compiled via -O3.

In theory the GCC is capable to do all these optimizations, even loop tiling, vectorization and autopar as it indeed has a polyhedral representation of loops. In practice it might be another story of course.

> The code is likely still compiled via -O3.

That doesn't matter. You have to specify the target platform as well.

Just take the naive code and specify -O3 in godbolt.org for example.

Assembly exerpt of the inner loop:

  L5:
        movss   xmm0, DWORD PTR [rax]
        mulss   xmm0, DWORD PTR [rdx]
        add     rax, 4
        add     rdx, rdi
        addss   xmm1, xmm0
        movss   DWORD PTR [rcx], xmm1
        cmp     rsi, rax
        jne     .L5

No AVX in sight.

Now specify -mavx2 -mfma and compare the inner loop:

   .L5:
        vmovss  xmm2, DWORD PTR [rax]
        vfmadd231ss     xmm0, xmm2, DWORD PTR [rdx]
        add     rax, 4
        add     rdx, rdi
        vmovss  DWORD PTR [rcx], xmm0
        cmp     rsi, rax
        jne     .L5

Case in point: without providing the compiler flags used for compilation, charts are misleading at best and outright wrong at worst.

I believe the flags are in the compile.sh script in the repo:

   OMP_NUM_THREADS=8 g++ -Wall -std=c++11 -fopenmp -D_USE_OPENBLAS -O3 -L/usr/lib  -lopenblas   -I /usr/include/openblas -march=native -mavx2 -mfma -funroll-loops compare_blas.cpp  -lopenblas -o matmul

Of course it is possible that different flags were used to benchmark the naive version.

edit: and of course those f*ss suffixes show that the compiler is not actually vectorizing the loops, just using the avx encoding of scalar instructions (plus fma).

edit2: likely -ffast-math is also required for vectorization as it will change the sequence of operations, but after experimenting with it a bit it seems that in this case is still not sufficient.

Literally no-one can write compilers smart enough.

Julia+LoopVectorization is already smart enough for this (it can't yet deal with more complicated transforms like lu or ffts, but it already gives you near optimal matmul code from a triple loop.

Does HN use a GPT-3 detector? This smells like the slightly imbricated (sic) prose it generates. It seems to make sense unless you actually know something this topic.

Since you supposedly know something about this topic, please contribute it.

Okay. Learn how to install and link a BLAS version appropriate for your machine architecture. Do not assume the particular compiler you are currently using has a magical optimizer.

This will future-proof your code and make it more portable, not to mention instill confidence that your code is correct. Of course if you are just playing around with a toy project, this is probably not worth the effort.

Sorry, I should have proofread it. A french-anglicism for nested-loops ("boucles imbriquées") slipped through.

> and if your compiler automatically derive the Pade approximation for improved numerical stability you know you are on the right track

I don't know what you're smoking, but I want some of it.

This is, I think, the gist of array langs like APL. Vectors + tight inner loops are very close to be great.

> Forget about performance for a moment, and just write those imbricated loops.

Right, if we forget about performance a three-nested-loop GEMM is certainly better than any other implementation. But that's completely meaningless and you could find such examples for almost any task. This is almost as saying writing hash functions is easy as long as you forget about collisions and avalanche effects.

> 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.

That is just not true: It is correct that deep learning tasks mainly reduce to convolutions and GEMM. Convolutions can be mapped to contraction kernels which again can be mapped to GEMM. However, Numpy's einsum is a great example how you shouldn't do it as it's notoriously slow. Also, the main task in general tensor contraction functions is to find an optimal contraction ordering rather than executing GEMM (you can have exponential savings if you choose a good contraction sequence). You can read more about this in e.g. https://arxiv.org/abs/2002.01935.

> 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.

This is true in theory and while modern compilers get better and better there is still no compiler framework that could, given an arbitrary linear algebra task outperform a manually optimized BLAS or LAPACK library.

Also, you are mentioning just a very few low-level optimization techniques that can be applied while modern compiler research currently focuses on different aspects as implemented in XLA for example.

> Optimizing code, even though it feels quite fun and rewarding to do manually is a job for compilers. > Write naive algorithms but smart compilers.

Still, don't write bad code and use compilers as an excuse.

The response to this comment is a bit sad.

Yes, "we just needs a sufficiently advanced compiler" is a meme, but 1) not everyone knows that and 2) that's not what this comment is asking for:

> but here one domain specific language can have a huge impact

There's a big difference between a sufficiently advanced compiler for general C code, and a compiler for a domain-specific array processing language, for example.

Have a look at halide: https://halide-lang.org , where you specify the computation to perform and a schedule (loop order, blocking, vectorisation) separately. There's even an automatic scheduler that can find a good schedule given some example array sizes and parameters.

The fftw paper is worth a read, too -- this shows that for a restricted domain it's possible to automatically find code that works as well as or better than hand-tuned algorithms. Generalising this is of course very hard.

For deep learning I think a lot of this is already done, see XLA for example.

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....>.

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!

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/sr...

"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

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.

On the other hand: If you have those 100 lines, it's one library less you have to think about and carry around.

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.

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.

Well, O(10000) = O(1), so it's not really a problem.

:-) yeah I was kind of expecting that comment

How is Julia doing in that regard?

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

cool approach to implement matrix multiplication

Reply on news.ycombinator.com