How does BLAS get such extreme performance?

Question

Out of curiosity I decided to benchmark my own matrix multiplication function versus the BLAS implementation... I was to say the least surprised at the result:

Custom Implementation, 10 trials of 1000x1000 matrix multiplication:

Took: 15.76542 seconds.

BLAS Implementation, 10 trials of 1000x1000 matrix multiplication:

Took: 1.32432 seconds.

This is using single precision floating point numbers.

My Implementation:

template<class ValT>
void mmult(const ValT* A, int ADim1, int ADim2, const ValT* B, int BDim1, int BDim2, ValT* C)
{
 if ( ADim2!=BDim1 )
 throw std::runtime_error("Error sizes off");
 memset((void*)C,0,sizeof(ValT)*ADim1*BDim2);
 int cc2,cc1,cr1;
 for ( cc2=0 ; cc2<BDim2 ; ++cc2 )
 for ( cc1=0 ; cc1<ADim2 ; ++cc1 )
 for ( cr1=0 ; cr1<ADim1 ; ++cr1 )
 C[cc2*ADim2+cr1] += A[cc1*ADim1+cr1]*B[cc2*BDim1+cc1];
}

I have two questions:

1. Given that a matrix-matrix multiplication say: nxm * mxn requires n*n*m multiplications, so in the case above 1000^3 or 1e9 operations. How is it possible on my 2.6Ghz processor for BLAS to do 10*1e9 operations in 1.32 seconds? Even if multiplcations were a single operation and there was nothing else being done, it should take ~4 seconds.
2. Why is my implementation so much slower?

• dmckee --- ex-moderator kitten: BLAS has been optimized up one side and down the other by specialist in the field. I assume it is taking advantage of the SIMD floating point unit on your chip and playing lots of tricks to improve the caching behavior as well...
• DeusAduro: Still how do you do 1E10 operations on a 2.63E9 cycles/second processor in 1.3 seconds?
• dmckee --- ex-moderator kitten: Multiple execution units, pipe-lining, and Single Instruction Multiple Data ((SIMD) which means doing the same operation on more than one pair of operands at the same time). Some compilers can target the SIMD units on common chips but you just about always have to explicitly turn in on, and it helps to know how it all works (en.wikipedia.org/wiki/SIMD). Insuring against cache misses is almost certainly the hard part.
• MSalters: Assumption is wrong. There are better algorithms known, see Wikipedia.
• Michael C. Lehn: @DeusAduro: In my answer for How to write a matrix matrix product that can compete with Eigen? I posted a small example on how to implement a cache efficient matrix-matrix product.

Answer 1

A good starting point is the great book The Science of Programming Matrix Computations by Robert A. van de Geijn and Enrique S. Quintana-Ortí. They provide a free download version.

BLAS is divided into three levels:

• Level 1 defines a set of linear algebra functions that operate on vectors only. These functions benefit from vectorization (e.g. from using SIMD such as SSE).

• Level 2 functions are matrix-vector operations, e.g. some matrix-vector product. These functions could be implemented in terms of Level 1 functions. However, you can boost the performance of these functions if you can provide a dedicated implementation that makes use of some multiprocessor architecture with shared memory.

• Level 3 functions are operations like the matrix-matrix product. Again you could implement them in terms of Level 2 functions. But Level 3 functions perform O(N^3) operations on O(N^2) data. So if your platform has a cache hierarchy then you can boost performance if you provide a dedicated implementation that is cache optimized/cache friendly. This is nicely described in the book. The main boost of Level 3 functions comes from cache optimization. This boost significantly exceeds the second boost from parallelism and other hardware optimizations.

By the way, most (or even all) of the high performance BLAS implementations are NOT implemented in Fortran. ATLAS is implemented in C. GotoBLAS/OpenBLAS is implemented in C and its performance-critical parts in Assembler. Only the reference implementation of BLAS is implemented in Fortran. However, all these BLAS implementations provide a Fortran interface such that it can be linked against LAPACK (LAPACK gains all its performance from BLAS).

Optimized compilers play a minor role in this respect (and for GotoBLAS/OpenBLAS the compiler does not matter at all).

IMHO no BLAS implementation uses algorithms like the Coppersmith–Winograd algorithm or the Strassen algorithm. The likely reasons are:

• Maybe it's not possible to provide a cache-optimized implementation of these algorithms (i.e. you would lose more than you would win)
• These algorithms are numerically not stable. As BLAS is the computational kernel of LAPACK this is a no-go.
• Although these algorithms have a nice time complexity on paper, the Big O notation hides a large constant, so it only starts to become viable for extremely large matrices.

Edit/Update:

The new and groundbreaking papers for this topic are the BLIS papers. They are exceptionally well written. For my lecture "Software Basics for High Performance Computing" I implemented the matrix-matrix product following their paper. Actually I implemented several variants of the matrix-matrix product. The simplest variant is entirely written in plain C and has less than 450 lines of code. All the other variants merely optimize the loops

 for (l=0; l<MR*NR; ++l) {
 AB[l] = 0;
 }
 for (l=0; l<kc; ++l) {
 for (j=0; j<NR; ++j) {
 for (i=0; i<MR; ++i) {
 AB[i+j*MR] += A[i]*B[j];
 }
 }
 A += MR;
 B += NR;
 }

The overall performance of the matrix-matrix product only depends on these loops. About 99.9% of the time is spent here. In the other variants I used intrinsics and assembler code to improve the performance. You can see the tutorial going through all the variants here:

ulmBLAS: Tutorial on GEMM (Matrix-Matrix Product)

Together with the BLIS papers it becomes fairly easy to understand how libraries like Intel MKL can gain such performance. And why it does not matter whether you use row or column major storage!

The final benchmarks are here (we called our project ulmBLAS):

Benchmarks for ulmBLAS, BLIS, MKL, openBLAS and Eigen

Another Edit/Update:

I also wrote some tutorials on how BLAS is used for numerical linear algebra problems like solving a system of linear equations:

High Performance LU Factorization

(This LU factorization is for example used by Matlab for solving a system of linear equations.)

(削除) I hope to find time (削除ここまで) to extend the tutorial to describe and demonstrate how to realise a highly scalable parallel implementation of the LU factorization like in PLASMA.

Ok, here you go: Coding a Cache Optimized Parallel LU Factorization

P.S.: I also did make some experiments on improving the performance of uBLAS. It actually is pretty simple to boost (yeah, play on words :) ) the performance of uBLAS:

Experiments on uBLAS.

Here a similar project with BLAZE:

Experiments on BLAZE.

• Ahmed Fasih: New link to "Benchmarks for ulmBLAS, BLIS, MKL, openBLAS and Eigen": apfel.mathematik.uni-ulm.de/~lehn/ulmBLAS/#toc3
• Aurélien Pierre: most of the links are dead
• Alex Shpilkin: A PDF of TSoPMC can be found on the page of the author, at cs.utexas.edu/users/rvdg/tmp/TSoPMC.pdf
• Nihar Karve: Although the Coppersmith-Winograd algorithm has a nice time complexity on paper, the Big O notation hides a very large constant, so it only starts to become viable for ridiculously large matrices.
• Michael C. Lehn: @ÖöTiib Also, it is a common misunderstanding that matrices have to be stored row major in C/C++. Every high performance implementation of BLAS (defined through a Fortan reference implementation) is written in C or C++ or assembly.

Answer 2

So first of all BLAS is just an interface of about 50 functions. There are many competing implementations of the interface.

Firstly I will mention things that are largely unrelated:

• Fortran vs C, makes no difference
• Advanced matrix algorithms such as Strassen, implementations dont use them as they dont help in practice

Most implementations break each operation into small-dimension matrix or vector operations in the more or less obvious way. For example a large 1000x1000 matrix multiplication may broken into a sequence of 50x50 matrix multiplications.

These fixed-size small-dimension operations (called kernels) are hardcoded in CPU-specific assembly code using several CPU features of their target:

• SIMD-style instructions
• Instruction Level Parallelism
• Cache-awareness

Furthermore these kernels can be executed in parallel with respect to each other using multiple threads (CPU cores), in the typical map-reduce design pattern.

Take a look at ATLAS which is the most commonly used open source BLAS implementation. It has many different competing kernels, and during the ATLAS library build process it runs a competition among them (some are even parameterized, so the same kernel can have different settings). It tries different configurations and then selects the best for the particular target system.

(Tip: That is why if you are using ATLAS you are better off building and tuning the library by hand for your particular machine then using a prebuilt one.)

• Robert van de Geijn: ATLAS is no longer the most commonly used open source BLAS implementation. It has been surpassed by OpenBLAS (a fork of the GotoBLAS) and BLIS (a refactoring of the GotoBLAS).
• Andrew Tomazos: @ulaff.net: That maybe. This was written 6 years ago. I think the fastest BLAS implementation currently (on Intel of course) is Intel MKL, but it isn't open source.
• ramslök: I agree with the spirit of your answer. Here's an academic link, but it shows that some have used Strassen-type/Winograd-type algorithms to have real world speedups ics.uci.edu/~paolo/FastMM/FMM-Reference/reference.html

Answer 3

First, there are more efficient algorithms for matrix multiplication than the one you're using.

Second, your CPU can do much more than one instruction at a time.

Your CPU executes 3-4 instructions per cycle, and if the SIMD units are used, each instruction processes 4 floats or 2 doubles. (of course this figure isn't accurate either, as the CPU can typically only process one SIMD instruction per cycle)

Third, your code is far from optimal:

• You're using raw pointers, which means that the compiler has to assume they may alias. There are compiler-specific keywords or flags you can specify to tell the compiler that they don't alias. Alternatively, you should use other types than raw pointers, which take care of the problem.
• You're thrashing the cache by performing a naive traversal of each row/column of the input matrices. You can use blocking to perform as much work as possible on a smaller block of the matrix, which fits in the CPU cache, before moving on to the next block.
• For purely numerical tasks, Fortran is pretty much unbeatable, and C++ takes a lot of coaxing to get up to a similar speed. It can be done, and there are a few libraries demonstrating it (typically using expression templates), but it's not trivial, and it doesn't just happen.

• DeusAduro: Thanks, I've added restrict correct code as per Justicle's suggestion, didn't see much improvement, I like the blockwise idea. Out of curiosity, without knowing the CPU's cache size how would one right optimal code?
• Stack Overflow is garbage: You don't. To get optimal code, you need to know the CPU's cache size. Of course the downside to this is that you're effectively hardcoding your code for best performance on one family of CPU's.
• Peter Cordes: At least the inner loop here avoids strided loads. It looks like this is written for one matrix already being transposed. That's why it's "only" one order of magnitude slower than BLAS! But yeah, it's still thrashing because of the lack of cache-blocking. Are you sure Fortran would help much? I think all you'd gain here is that restrict (no aliasing) is the default, unlike in C / C++. (And unfortunately ISO C++ doesn't have a restrict keyword, so you have to use __restrict__ on compilers that provide it as an extension).

Answer 4

I don't know specfically about BLAS implementation but there are more efficient alogorithms for Matrix Multiplication that has better than O(n3) complexity. A well know one is Strassen Algorithm

• Michael C. Lehn: The Strassen Algorithm is not used in numerics for two reasons: 1) It is not stable. 2) You save some computations but that comes with the price that you can exploit cache hierarchies. In practice you even loose performance.
• Jianyu Huang: For the practical implementation of Strassen Algorithm tightly built upon BLAS library source code, there is a recent publication: "Strassen Algorithm Reloaded" in SC16, which achieves higher performance than BLAS, even for the problem size 1000x1000.

Answer 5

Most arguments to the second question -- assembler, splitting into blocks etc. (but not the less than N^3 algorithms, they are really overdeveloped) -- play a role. But the low velocity of your algorithm is caused essentially by matrix size and the unfortunate arrangement of the three nested loops. Your matrices are so large that they do not fit at once in cache memory. You can rearrange the loops such that as much as possible will be done on a row in cache, this way dramatically reducing cache refreshes (BTW splitting into small blocks has an analogue effect, best if loops over the blocks are arranged similarly). A model implementation for square matrices follows. On my computer its time consumption was about 1:10 compared to the standard implementation (as yours). In other words: never program a matrix multiplication along the "row times column" scheme that we learned in school. After having rearranged the loops, more improvements are obtained by unrolling loops, assembler code etc.

 void vector(int m, double ** a, double ** b, double ** c) {
 int i, j, k;
 for (i=0; i<m; i++) {
 double * ci = c[i];
 for (k=0; k<m; k++) ci[k] = 0.;
 for (j=0; j<m; j++) {
 double aij = a[i][j];
 double * bj = b[j];
 for (k=0; k<m; k++) ci[k] += aij*bj[k];
 }
 }
 }

One more remark: This implementation is even better on my computer than replacing all by the BLAS routine cblas_dgemm (try it on your computer!). But much faster (1:4) is calling dgemm_ of the Fortran library directly. I think this routine is in fact not Fortran but assembler code (I do not know what is in the library, I don't have the sources). Totally unclear to me is why cblas_dgemm is not as fast since to my knowledge it is merely a wrapper for dgemm_.

• jcsahnwaldt Reinstate Monica: Much faster than my initial naive array-based code. Cool! But... still roughly 10 times slower than Julia's built-in BLAS matrix multiplication, even if I turn off Julia's OpenMP multi-threading. Tested with matrix sizes from 1000x1000 to 5000x5000 on my M2 MacBook Air.

Answer 6

With respect to the original code in MM multiply, memory reference for most operation is the main cause of bad performance. Memory is running at 100-1000 times slower than cache.

Most of speed up comes from employing loop optimization techniques for this triple loop function in MM multiply. Two main loop optimization techniques are used; unrolling and blocking. With respect to unrolling, we unroll the outer two most loops and block it for data reuse in cache. Outer loop unrolling helps optimize data-access temporally by reducing the number of memory references to same data at different times during the entire operation. Blocking the loop index at specific number, helps with retaining the data in cache. You can choose to optimize for L2 cache or L3 cache.

https://en.wikipedia.org/wiki/Loop_nest_optimization

Answer 7

This is a realistic speed up. For an example of what can be done with SIMD assembler over C++ code, see some example iPhone matrix functions - these were over 8x faster than the C version, and aren't even "optimized" assembly - there's no pipe-lining yet and there is unnecessary stack operations.

Also your code is not "restrict correct" - how does the compiler know that when it modifies C, it isn't modifying A and B?

• DeusAduro: Sure if you called the function like mmult(A..., A..., A); you would certainly not get the expected result. Again though I wasn't trying to beat/re-implement BLAS, just seeing how fast it really is, so error checking wasn't in mind, just the basic functionality.
• Justicle: Sorry, to be clear, what I'm saying is that if you put "restrict" on your pointers, you'd get much faster code. This is because every time you modifiy C, the compiler doesn't have to reload A and B - dramatically speeding up the inner loop. If you don't believe me, check the disassembly.
• Michael Burr: @DeusAduro: This isn't error checking - it's possible that the compiler is unable to optimize the accesses to the B[] array in the inner loop because it might not be able to figure out that the A and C pointers never alias the B array. If there were aliasing it would be possible for the value in the B array to change while the inner loop is executing. Hoisting the access to the B[] value out of the inner loop and putting it in a local variable might enable the compiler to avoid continual accesses to B[].
• DeusAduro: Hmmm, so I tried first using the '__restrict' keyword in VS 2008, applied to A, B, and C. This showed no change in the result. However moving the access to B, from the innermost loop to the loop outside improved the time by ~10%.
• Justicle: Sorry, I'm not sure about VC, but with GCC you need to enable -fstrict-aliasing. There's also better explanation of "restrict" here: cellperformance.beyond3d.com/articles/2006/05/…

Answer 8

I find the other answers somewhat lacking. It doesn't matter that much what BLAS implementations are doing in detail, but it matters what OP's code is not doing: Taking advantage of SIMD instruction sets like SSE (AVX was not yet available in 2009), taking advantage of cache and taking advantage of compiler optimization. 3 1000x1000 single precision matrices fit only in the very largest caches of 2009, and without caching, all data has to be loaded from memory. According to cpu-monkey.com, a Core 2 Quad Q9650 can achieve about 17 GB/s memory bandwidth and 10 1000x1000 matrix multiplications could require up to 10 * 1e9 memory accesses, or 40 GB for single precision floats. Even without any additional time for computation and assuming perfect utilization of RAM bandwidth, the required time would be about 2.3 seconds.

By working on smaller blocks, like BLAS implementations do, you can fit these blocks into cache and make use of cache's much higher bandwidth (probably at least 32 bytes per clock cycle for level 1 cache even in 2009). If you do it well enough that memory bandwidth doesn't bottleneck you any more, you are bottlenecked by computational performance. In 2009, you would be able to use SSE4.1, which should achieve a throughput of 4 single precision FLOPs per clock cycle using mulps and addps. At OP's 2.6 GHz that makes 10.4 GFLOP/s. To execute 20 * 1e9 (not 10 * 1e9, because you need both additions and multiplications) operations in 1.32 seconds means >15 GFLOP/s, so the BLAS implementation uses some trick to get around that limit and that trick is probably to interleave multiplications and additions so that the CPU can schedule them to different executions units (not sure how many a 2009 CPU had, but Sandy Bridge had 6 per core, so 2+ per core seems reasonable), so that additions and multiplications can be executed simultaneously. That raises the limit to 20.8 GFLOP/s, of which the BLAS implementation achieves 73%.

It is easy to improve OP's code to come closer to those 73%. Simply compiling with -O3 -march=core2 gives me about 5 seconds on an 3.7 GHz i7-8700K, compared to over 40 seconds without optimization. Slightly reordering the computation and unrolling the loops improves that to 3.5 seconds. This slightly improved code achieves about 33% of the performance of the highly optimized OpenBLAS code on my system, which in turn achieves about 84% of my system's theoretical maximum. OP's original code achieves <10% of a highly optimized implementation, so it probably was compiled without the highest compiler optimization and without SSE support. Of course compiler optimization has come a long way since 2009, so I would have to rerun the experiment with a 2009 compiler for a more accurate answer.

Answer 9

For many reasons.

First, Fortran compilers are highly optimized, and the language allows them to be as such. C and C++ are very loose in terms of array handling (e.g. the case of pointers referring to the same memory area). This means that the compiler cannot know in advance what to do, and is forced to create generic code. In Fortran, your cases are more streamlined, and the compiler has better control of what happens, allowing him to optimize more (e.g. using registers).

Another thing is that Fortran store stuff columnwise, while C stores data row-wise. I havent' checked your code, but be careful of how you perform the product. In C you must scan row wise: this way you scan your array along contiguous memory, reducing the cache misses. Cache miss is the first source of inefficiency.

Third, it depends of the blas implementation you are using. Some implementations might be written in assembler, and optimized for the specific processor you are using. The netlib version is written in fortran 77.

Also, you are doing a lot of operations, most of them repeated and redundant. All those multiplications to obtain the index are detrimental for the performance. I don't really know how this is done in BLAS, but there are a lot of tricks to prevent expensive operations.

For example, you could rework your code this way

template<class ValT>
void mmult(const ValT* A, int ADim1, int ADim2, const ValT* B, int BDim1, int BDim2, ValT* C)
{
if ( ADim2!=BDim1 ) throw std::runtime_error("Error sizes off");
memset((void*)C,0,sizeof(ValT)*ADim1*BDim2);
int cc2,cc1,cr1, a1,a2,a3;
for ( cc2=0 ; cc2<BDim2 ; ++cc2 ) {
 a1 = cc2*ADim2;
 a3 = cc2*BDim1
 for ( cc1=0 ; cc1<ADim2 ; ++cc1 ) {
 a2=cc1*ADim1;
 ValT b = B[a3+cc1];
 for ( cr1=0 ; cr1<ADim1 ; ++cr1 ) {
 C[a1+cr1] += A[a2+cr1]*b;
 }
 }
 }
} 

Try it, I am sure you will save something.

On you #1 question, the reason is that matrix multiplication scales as O(n^3) if you use a trivial algorithm. There are algorithms that scale much better.

• Andrew Tomazos: This answer is completely wrong sorry. BLAS implementations are not written in fortran. The performence-critical code is written in assembly, and the most common ones these days are written in C above that. Also BLAS specifies the row/column order as part of the interface, and implementations can handle any combination.
• Michael C. Lehn: Yes, this answer is completely wrong. Unfortunately it is full of common non-sense, e.g. the claim BLAS was faster because of Fortran. Having 20(!) positive ratings is a bad thing. Now this non-sense even spreads further because of the popularity of Stackoverflow!
• Andrew Tomazos: I think you're confusing the unoptimized reference implementation with production implementations. The reference implementation is just for specifying the interface and behaviour of the library, and was written in Fortran for historical reasons. It isn't for production use. In production people use optimized implementations that exhibit the same behaviour as the reference implementation. I've studied the internals of ATLAS (which backs Octave - Linux "MATLAB") which I can confirm first hand is written in C/ASM internally. The commercial implementations are almost certainly as well.
• Andrew Tomazos: @KyleKanos: Yes, here is the source of ATLAS: sourceforge.net/projects/math-atlas/files/Stable/3.10.1 As far as I know it is the most commonly used open source portable BLAS implementation. It is written in C/ASM. High performance CPU manufacturers like Intel, also provide BLAS implementations especially optimized for their chips. I guarantee at the low level parts of Intels library are are written in (duuh) x86 assembly, and I'm pretty sure the mid-level parts would be written in C or C++.
• Andrew Tomazos: @KyleKanos: You're confused. Netlib BLAS is the reference implementation. The reference implementation is much slower than optimized implementations (see performance comparison). When someone says they are using netlib BLAS on a cluster, it doesn't mean they are actually using the netlib reference implementation. That would be just silly. It just means they are using a lib with the same interface as the netlib blas.