Fool's errand - building the fastest numerical computation library

[sfn]

I may or may not end up figuring out how to write my own domain-specific high-performance physics simulators in the near future, but the possibility got me thinking about how to do the underlying calculations in the fastest possible way for any given hardware, from my old ThinkPad to Multi-Processor Multi-GPU racks.

Now, I may be a fool but I’m not an idiot, each hardware vendor writes their own set of libraries for these sorts of workloads and I don’t think to outperform them. The idea would be to write a wrapper library that presents a simple, unified interface to each of these vendor implementations, choosing which one to use when multiple are available (for example, on a box with an intel CPU and Nvidia GPU, choosing what problems to execute where and how).

Here are the libraries I’ve found so far for the sort of computations one might need, not sure if I’ve missed anything important:

Calculation type	x86 CPU	arm CPU	Apple M-series	Intel GPU	NVIDIA GPU	AMD GPU
BLAS/(Sca)LAPACK	oneMKL (ex intel)	APL	A. Accelerate*	oneMKL (SYCL flavour)	cuSOLVERMp	rocSOLVER
FFT	oneMKL	APL	A. Accelerate	oneMKL	cuFFTMp	rocFFT
(P)RNG	oneMKL(?)	openRNG (part of APL)	openRNG(?)	oneMKL	cuRAND	rocRAND
Sparse	oneMKL	APL	A. Accelerate	oneMKL	cuSPARSE/cuSPARSELt	rocSPARSE/rocSPARSELt

An eventual simulator also needs a good non-linear solver and a DAE solver. If these ever come to light, I think they should be separate libraries implemented on top of this computational core.

Acronyms:

• A. = Apple
• APL = arm performance libraries
• MPS = (apple) metal performance shaders

Other notes

• This project intends to adhere to Zig praxis and philosophy: follow the style guide, accept external alloc/io where possible, no AI for writing code, ecc.
• If we can get this library to provide compatible C/Fortran interfaces for each component (e.g. use this library as a drop-in replacement for existing software) that would be great, but performance comes first
• RISC-V cpu’s are excluded from lack of hardware and knowledge, but I see no reason to add them once software/hardware comes out
• I’m aware of the OpenMathLib, they might be a good generic fallback option but they seem to be slightly outperformed by vendor implementations across the board. Also, no device (GPU) capability.
• On Apple M-series chips, matrix multiplication of certain large matrices might be further sped up with shared memory and MPS, see Apple vs. Oranges: Evaluating the Apple Silicon M-Series SoCs for HPC Performance and Efficiency
• Different sparse implementations (for example cuSPARSE vs cuSPARSELt) need to be selected depending on the sparsity ratio of the matrix. Sparse routines might therefore want to take a sparseness value to aid in deciding how to dispatch them, which can either be estimated by the user in advance or calculated.
• FFTW would not, in fact, appear to be the fastest Fourier transform, but I’m not entirely convinced yet: Benchmark source 1 Benchmark source 2
• Need better benchmarks for openRNG vs oneMKL on x86 CPU’s and alternatives to openRNG for Apple hardware
• Given the focus on performance, we need a consistent, scalable benchmark for each component in general
• A suggested approach is to dynamically load the vendor libraries that are installed on the running system, and fall back to open-source generic alternatives (I.e. OpenBLAS, OpenRNG) when we can’t find anything better
• Something similar is kind-of provided by GitHub - uxlfoundation/oneMath: oneAPI Math Library (oneMath) · GitHub, but a) it’s C++ b) no dedicated Apple support and c) requires the user to handle dispatching (choosing which device to use for what).

Feedback appreciated as to the best way to go about doing this (i.e. tell me why this is stupid and I’m wasting my time )!

[smj-edison]

I’ll admit I’m not super familiar with this space, but HPC is definitely on my radar to learn. One promising piece prior work is the Eigen library in C++ if you’re not already aware of it.

[pzittlau]

First of I like the spirit but I don’t quite know if Zig is a real advantage here because most of them are written half in assembly anyway while the other half being macros. The only ones I know of that are trying to be readable are BLIS and ulmBlas which is more trying to be educational.

If you want for it to also abstract over accelerators and maybe even distribution and cluster configuration you could also take a look at sycl.

A nice thing that I would like if it is reasonably readable.

[sfn]

Zig isn’t so much an advantage as an alternative to C/C++ (these libraries provide C interfaces), mostly personal preference (although I suspect the I/O interface over io_uring might come in handy should I ever need to write a task scheduler).

I’ve worked with SYCL in the past (check out the amazing https://adaptivecpp.github.io/ project), but I don’t think the abstraction is work the overhead for this use-case - vendor libraries already handle that for us, and still provide the tools we need for memory management.

What do you mean by readable? Good quality code? An interface that calls functions triangularMatrixMultiply instead of trmm? I must admit, a lot of the older libraries for this sort of workload seem to have been written to a character limit.

[pzittlau]

@sfn

What do you mean by readable? Good quality code? An interface that calls functions triangularMatrixMultiply instead of trmm? I must admit, a lot of the older libraries for this sort of workload seem to have been written to a character limit.

I don’t think that stuff like trmm, gemm, or axpy are really problematic. I would just think that generics likely aren’t done via prefixes. Having a triangularMatrixMultiply that just defers to trmm would likely make it a bit more usable for people that haven’t really worked with blas but isn’t that necessary I think.

For readability I was more talking about the implementation. If I look into something like openBLAS I see a lot of different ifdefs enabling or disabling single lines which is quite common in old C and on it’s own okay, but in addition to all the macros it get’s hard to reason about, at least for me.
Of course, I also wouldn’t really like to have a lot less speed just to make it readable but I think we now know better then to have a soup of conditional compilation.