# 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
• Not too much I can add to what's already been said but here are a couple of small points. Almost all modern CPUs have multiple cores and a state of the art BLAS will have the option of using all of them particularity for a BLAS3 operation like dgemm. Second, as has already been mentioned, the vector operations like SSE2, SSE3, AVX can substantially improve performance. The more advanced vector instruction sets like AVX are available only on more modern, sophisticated CPUs. – Bill Greene Oct 24 '16 at 0:53

## 3 Answers

Consolidating the comments:

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

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