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?
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". $\endgroup$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. $\endgroup$dgemm
? $\endgroup$