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.
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.
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
EDIT: the https://ginkgo-project.github.io/ does something very similar for C++, minus the automated dispatch and Apple support. Definitely worth looking into how they do things.
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 )!
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.
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.
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.
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.
In that sense, definitely planning on being readable. That’s the beauty of Zig, you can write highly performant code while still remaining elegant and idiomatic.
Do bear in mind that the entire point of the project is to offload as much work as possible to libraries that have spent far more resources on optimization, so figuring out how the actual computations are getting carried out might still be hard (OpenBLAS, OpenRNG) or impossible (cu*, roc*).
This is what Numpy/SciPy do. It is also the business model for matlab (yuck!).
I think Zig as a toolchain has an edge here. So many times, I’ve seen code rot to death because poor performance, where a few dgemms would have done a big difference, but alas! the dread of getting some cmake versions to give me blas without braking boost … narration stops here to prevent tears. Point is, just being able to zig fetch --save blas would be great.
Numpy/SciPy and also matlab are very successful because they are very useful.
I came across this the other day, big shout out to @adworacz.
I’m well aware - needing to learn MATLAB for uni was part of the inspiration for this project, despite the fact that it can easily be beaten by well-written FORTAN calling a good LAPACK implementation.
Toolchain integration is going to be interesting, given the dependence on third-party libraries. Shipping everything statically as part of the binary is going to mean massive bloat, but we probably want some sort of high-performance fallback in case a user e.g. misconfigures their system.
I’ve read enough of the specs for blas/fortran and cuSPARSEMp to come to the conclusion that, even with dma-buf on linux, operating on user-allocated memory is going to involve prohibitive copy overhead on non-CPU devices. Memory management in general is going to be interesting, but we are not going to be able to simply export the library’s zig functions to get a standard C/Fortran interface without compromising the design of the zig interface. Quite probably, C/Fortran compat. is going to end up as a separate project entirely, if it ever materializes.
(I’ve started working on a POC, but it’s not ready for sharing yet).
Done more research, found out that most of the nvidia high-performance api’s (e.g. cuBLASmp, cuFFTmp) are only available on Linux. This means that the high-performance nvidia backend can use useful features like HMM that are only available on linux without compatibility problems (because nvidia themselves are the root of compatibility issues). On the other hand, this basically means that windows users of nvidia hardware are in a pickle. Amusing as that is, we should either provide a cross-platform vulkan implementation as a fallback, or implement a lower-performance nvidia backend. Probably, the best course of action is to make vulkan first, to support a wider variety of hardware, then make the nvidia one when someone needs it.
That said, I still haven’t finished the CPU backend yet… back to work!
Edit: thought about it more. Maybe an initial implementation with the non-mp variants of the libraries would be a better, more portable first pass, and the multi-gpu implementation can come later.
)!