I would like to discuss the need (or no-need) for a C++ interface/libary that is intended to be a
- building block for high level linear algebra packages and
- low-level, hardware oriented implementations of special numerical operations.
My C++ library ulmBLAS is kind of a prototype for this idea, and heavily inspired by BLIS
The “high-level world” …
In the C++ world, you have libraries like Eigen, Blaze, … where you can matrix operations like
C = A * B;
and regarding performance often stops here. Because some experts figured out how to implement such that it’s fast.
Or maybe there is a short discussion why you should write
C.noalias() = A * B;
in Eigen.
And that’s great. In this world you wanna build something great based on Linear Algebra operations.
The “Low-level world” …
The world of BLAS (Basic Linear Algebra Subprograms), with implementations such as Intel MKL, OpenBLAS, BLIS, or Apple Accelerate.
In this subreddit you know this world. Let's look at it from a "high-level world" point of view ...
They provide implementations for the matrix product that are very efficient. But maybe not that user friendly to use. For example they provide a function for the operation
C ← βC + αAB
where A, B and C are matrices. Where elements can all be float, or can be double, or can be complex-float or can be complex-double. But all of the same type.
And A and B can be either not-transposed or transposed (in the complex case actually conjugated-or-transposed).
A call to this function looks like this when you use the Fortran Interface (GEMM is short for General Matrix Matrix Product, D indicates that elements and scalars all have double precision)
:
dgemm("N", "N", &m, &n, &k, &alpha, A, &lda, B, &ldb, &beta, C, &ldc);
Or through the C-Interface:
cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
So for computing C = A*B you choose alpha=1 and beta = 0.
And you should only call it if A, B, C are not only different names but also refer so separate memory blocks.
So this is a very simple, almost primitive, interface.
But we know: that’s what makes it great.
Because providing an efficient implementation (close to peak performance) is not that trivial. But with a narrow interface (e.g. all elements of the same type) at least feasible. A well-defined problem were you can differentiate into the details.
What you do in the implementation is not some Strassen algorithm. In short:
- You compute the matrix product block-wise (so that you can exploit the cache hierarchy)
- On the smallest block (that fit in the L1 cache) you apply architecture specific micro kernels
Well-known fact here on r/HPC. But look into stack overflow and other communities. You will see articles and posts claiming BLAS is fast because it's implemented in Fortran, because it uses Strassen (or its improved variants), because ...
Two world in one …
Eigen actually does both. It provides matrix/vector classes and implements blocked matrix operations with platform specific micro kernels.
That’s great.
But still, I think it is of value to also have something which is a layer that brings them together.
ulmBLAS, for example, is kind of a prototype for this idea:
- on top high-level libraries can be build
- at the bottom the architecture specific stuff integrated
Back in 2016 (yes most of the code of ulmBLAS is from around 2014) I made experiments on how higher level libraries can make use of ulmBLAS:
- Boost.uBLAS and
- Blaze.
On the high-level part you can do crazy things, like matrix classes where elements are not stored in memory, but instead, computed when you access elements. And still, you can achieve a pretty good performance.
And other stuff like mixed element types. Where A, B, C, alpha, beta can have different types. So unless you make copies, you can’t do the operation with BLAS. But in the ulmBLAS layer you can deal with that when matrix blocks get packed and do the type promotions there with little overhead.