Beating typical BLAS libraries matrix multiplication performance

Question

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?

Answer 1

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

Answer 2

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.

Answer 3

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.

Answer 4

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