# Beating typical BLAS libraries matrix multiplication performance

A dull matrix multiplication algorithm where we use the formula

$$C_{ij}=\sum_{k}A_{ik}B_{kj}$$

By literally following this in 3 loops we'll get a very slow program, because we don't utilize processor's vectorization power. An example of this slow implementation is:

T comp;
for (int i = 0; i < lhs.rows(); i++)
{
for (int j = 0; j < rhs.columns(); j++)
{
comp = 0;
for (int k = 0; k < lhs.columns(); k++)
{
comp += lhs.at(i,k)*rhs.at(k,j);
}
result.at(i,j) = comp;
}
}


My question is: If we consider vectorization, are we going to beat typical libraries like OpenBLAS?.

Why? Because I want to specialize my matrix multiplications so much that BLAS doesn't provide what I need.

If we assume that these matrices are square matrices of side-length $n$, and are stored in column-major order and contiguous in memory (e.g., std::vector, with $n^2$ elements). Consider the following code snippet:

transpose(lhs); //apply in-place transpose
for (int i = 0; i < lhs.columns(); i++)
{
for (int j = 0; j < rhs.columns(); j++)
{
result.at(i,j) =
std::transform(
lhs.begin()+i*lhs.rows(),
lhs.begin()+(i+1)*lhs.rows(),
rhs.begin()+j*rhs.rows(),
std::multiplies<T>);
}
}
transpose(lhs);


Now this doesn't only consider vectorization, but even uses standard algorithms to give the compiler full freedom to do all the optimizations it wants. So my question here is: Will this beat good BLAS libraries? What else should I do to make such code super-efficient?

• Trust the people who've been doing it for more than twenty years? Why do you think a good BLAS library (e.g., OpenBLAS, MKL) doesn't make use of the same ideas -- and then some: for most current architectures, the major kernels are handwritten in architecture-specfic assembly language? – Christian Clason Oct 23 '16 at 20:29
• @ChristianClason Actually if you look into the implementation of dgemm in LAPACK (yes, in netlib's LAPACK package, there's a BLAS implementation), you'll see it's just 3 loops... that's it! That's why I'm surprised. I'm interested in knowing what else they're doing with more information than just "they do it for 30 years, they know better than you". – The Quantum Physicist Oct 23 '16 at 20:46
• @TheQuantumPhysicist: The netlib implementation is just meant as a simple reference implementation -- it's not really at all optimised. In my experience, a highly-optimised, architecture-specific library (i.e. Intel's MKL, gotoBLAS, OpenBLAS, etc) is often many, many times faster. These libraries employ a number of strategies to achieve such efficiency (blocking for cache reuse, architecture-specific vector instructions, etc). You would really have to know what you were doing to even compete with them. – Darren Engwirda Oct 23 '16 at 21:26
• What @Darren said -- nobody uses the LAPACK implementation (which is really just a bare-bones instantiation of the interface, which is the main point). You yourself mentioned OpenBLAS -- have you looked at their implementation(s) of dgemm? – Christian Clason Oct 23 '16 at 22:55
• Also, don't guess -- benchmark! You have your implementation and OpenBLAS; if yours beats the competition for the specific matrices your interested in, good on you! – Christian Clason Oct 23 '16 at 23:07

No, you are very unlikely to beat a typical BLAS library such as Intel's MKL, AMD's Math Core Library, or OpenBLAS.1 These not only use vectorization, but also (at least for the major functions) use kernels that are hand-written in architecture-specific assembly language in order to optimally exploit available vector extensions (SSE, AVX), multiple cores, and cache reuse. For example, OpenBLAS has

dgemm_kernel_16x2_haswell.S
dgemm_kernel_4x4_haswell.S
dgemm_kernel_4x8_haswell.S
dgemm_kernel_4x8_sandy.S
dgemm_kernel_6x4_piledriver.S
dgemm_kernel_8x2_bulldozer.S
dgemm_kernel_8x2_piledriver.S


and that's only for the x86_64 architecture -- not only a different implementation for each instruction set, but also for different register blocks (meaning that different kernels can be used depending on the matrix size). Of course, where possible, they make use of (also optimized) BLAS2 and BLAS1 operations. That's not something a compiler could do automatically.

1. The LAPACK implementation of BLAS available from netlib does not really count -- it's a reference implementation and should rather be considered as a (functional) specification of the interface. (See: http://dx.doi.org/10.1145/355841.355847, http://dx.doi.org/10.1145/42288.42291 and http://dx.doi.org/10.1145/77626.79170)

If we consider vectorization, are we going to beat typical libraries like OpenBLAS?

I do not fully agree with the other answer and would like to say kind of. It is true that libraries like Intel (R) MKL have handwritten kernels which use a lot of optimizations. On the other hand, you cannot compile Intel (R) MKL on your own machine, so it has to detect whether AVX instructions can be used at runtime (if this is possible at all).

That is why I believe that you actually can beat these BLAS implementations using C++ template libraries which make heavy use of SIMD and OpenMP for multithreading.

A good example for that is the Blaze library. The Bitbucket repository also has a benchmark page where they also compare BLAS level 3 routines. For very large matrices Blaze and Intel (R) MKL are almost the same in speed (probably memory limited) but for smaller matrices Blaze beats MKL. It is even more obvious for the BLAS level 2 routines.

• Do you know why the blaze benchmarks you linked to don't include OpenBLAS? – Kirill Jun 27 '17 at 2:48
• @Kirill No idea. Looks like the developers put the main focus on how it compares to other template libraries. – Henri Menke Jun 27 '17 at 2:56
• "On the other hand, you cannot compile Intel (R) MKL on your own machine, so it has to detect whether AVX instructions can be used at runtime (if this is possible at all)." You definitely can detect wich instructions are available on the host machine, see cpuid.See also CPU dispatching topics 1, 2, ... – Catree Mar 4 '18 at 4:57

As mentioned, netlib BLAS is not at all optimized, but it is definetly the "refblas". Using IKML, ACML, OpenBLAS or "your vendor" BLAS, you are (somehow) assured, that the results of the operation of the optimized BLAS is equal to the "refblas" up to a known error. Take into care that: vendors (intel, amd, nvidia, ...) try hard to offer an implmentation that is superior on their own platform. They have professionals who care for accuracy (e.g. avoid the Pentium FDIV bug) and know details on platform strengths and weaknesses and consequently tune for performance. The most specialized "vendor" BLAS for your platform is very likely to be the best you can ever get for that platform. OpenBLAS is optimized for several popular contemporary (not future!) platforms. If your code is even more platform-agnostic, ATLAS might be a good choice. Never try to implement a BLAS function if you don't have to.

This is not an answer, but a reference to explore the topic.

Here there is an article by Higham about the matrix multiplication (~ 1990).

Title: Exploiting Fast Matrix Multiplication Within the Level 3 BLAS