Almost no part of the algorithm as specified in the blog post is GPU friendly, but it can be improved.
The "symbolic" phase is very hard to port to the GPU in the first place and almost always happens on the CPU side. But even then it's hard (but not impossible) even to parallelize it. This is generally OK because it is usually cheaper than the numeric phase anyways. But it can be trouble on workstations where the CPU exists primarily as a device for shoveling data into many GPUs, you don't want to stall their kernel pipelines with your symbolic analysis.
For the numeric phase you need to introduce the concept of either a "front" or a "supernode." This is a technique where multiple columns get batched together and you get a new elimination tree in terms of those batches of columns rather than individual columns. This turns a lot of irregular memory access (gather/scatter style updates) into densely addressed memory patterns, often just calls into level 3 BLAS or LAPACK.
There are some sparsity patterns for which the above supernode/multifrontal approach does not work very well, but for many practical cases like PDE simulations it does. The sparsity patterns which defeat it usually spread nonzeros out all over the matrix rather than containing them within a specific band. What this causes in practice is when you materialize a supernode's dense matrix, most of its numeric entries will just be 0s, so you create a lot of extra work that serves no productive purpose.
I have no submission for this but I joined the hype in my own way by optimizing the training loop. These tiny models are not really well suited to frameworks like pytorch, and with highly patient AI agents we can now just inline the whole thing into C++ just to see what happens, which I do below:
Basically I just don't want to hear about "the state of SIMD in Rust" unless it is about dramatic improvement in autovectorization in the rust compiler.
80%-90% or so of real life vectorization can be achieved in C or C++ just by writing code in a way that it can be autovectorized. Intrinsics get you the rest of the way on harder code. Autovectorization is essentially a solved problem for the vast majority of floating point code.
Not so with rust, because of a dogmatic approach to floating point arithmetic that assumes bitwise reproducibility is the "right" answer for everyone (actually, it's the right answer to almost nobody) to the point of not even allowing a user to even flag on these optimizations. and once you get to the point of writing intrinsics you have to handwrite code for every new architecture when autovectorizers could have gotten you 80%-90% of the way there with a single source and often this is just enough.
the contention with the above is that if a user needs SIMD they can just use some SIMD API and make their intention more clear. this is essentially an argument that we should handwrite intrinsics. well guess what. I'm a programmer and I use compilers because they _do this for me_ and indeed are able to do so very easily in C or C++ when I instruct it that I'm ok with with reordering operations and other "accuracy impacting" optimizations.
The huge joke on us is that these optimizations generally have the effect of _improving_ accuracy because it will reduce the number of rounding steps either by simply reducing the number of operations or by using fused multiply adds which round only once.
There are two sides to numerical linear algebra. The first is the "linear algebra" part, which is very mathematically sophisticated and the language you choose to represent these concepts is not so important as your understanding. pencil and paper is the ideal place to prove out understanding of this.
The "numerical" part is a minefield because it will take all your math and demolish it. just about every theoretical result you proved out above will hold not true and require extra-special-handholding in code to retain _some_ utility.
As such I think a language which enables you to go as fast as possible from an idea to seeing if it crosses the numerical minefield unscathed is the one to use, and these days that is python. It is just so fast to test a concept out, get immediate feedback in the form of plotting or just plain dumb logging if you like, and you can nearly instantly share this with someone even if you're on ARM +linux & they are Intel+windows
The most problematic issue with python&numpy, as it relates to learning _numerical_ side of linear algebra, is making sure you haven't unintentionally promoted a floating point precision somewhere (for example: if your key claim is that an algorithm works entirely in a working precision of 32 bits but python silently promoted your key accumulator to 64 bits, you might get a misleaing idea of how effective the algorithm was) but these promotions don't happen in a vacuum and if you understand how the language works they won't happen.
edit: & I have worked professionally with fortran for a long time, having known some of the committee members and BLAS working group. so I have no particular bias against the language
The "symbolic" phase is very hard to port to the GPU in the first place and almost always happens on the CPU side. But even then it's hard (but not impossible) even to parallelize it. This is generally OK because it is usually cheaper than the numeric phase anyways. But it can be trouble on workstations where the CPU exists primarily as a device for shoveling data into many GPUs, you don't want to stall their kernel pipelines with your symbolic analysis.
For the numeric phase you need to introduce the concept of either a "front" or a "supernode." This is a technique where multiple columns get batched together and you get a new elimination tree in terms of those batches of columns rather than individual columns. This turns a lot of irregular memory access (gather/scatter style updates) into densely addressed memory patterns, often just calls into level 3 BLAS or LAPACK.
There are some sparsity patterns for which the above supernode/multifrontal approach does not work very well, but for many practical cases like PDE simulations it does. The sparsity patterns which defeat it usually spread nonzeros out all over the matrix rather than containing them within a specific band. What this causes in practice is when you materialize a supernode's dense matrix, most of its numeric entries will just be 0s, so you create a lot of extra work that serves no productive purpose.