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 transposefor (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?
Beating typical BLAS libraries matrix multiplication performance
optimization;c++;matrix;performance;blas
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.1These 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.Sand 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) 