Also, the major point is that BLAS has little to no role played here. Algorithms which just hit BLAS are very suboptimal already. There's a tearing step which r

Also, the major point is that BLAS has little to no role played here. Algorithms which just hit BLAS are very suboptimal already. There's a tearing step which reduces the problem to many subproblems which is then more optimally handled by pure Julia numerical linear algebra libraries which greatly outperform OpenBLAS in the regime they are in: https://github.com/YingboMa/RecursiveFactorization.jl#performance https://github.com/YingboMa/RecursiveFactorization.jl#perfor... And there are hooks in the differential equation solvers to not use OpenBLAS in many cases for this reason: https://github.com/SciML/DiffEqBase.jl/blob/master/src/linear_nonlinear.jl#L76-L77 https://github.com/SciML/DiffEqBase.jl/blob/master/src/linea... Instead what this comes out to is more of a deconstructed KLU, except instead of parsing to a single sparse linear solve you can do semi-independent nonlinear solves which are then spawning parallel jobs of small semi-dense linear solves which are handled by these pure Julia linear algebra libraries. And that's only a small fraction of the details. But at the end of the day, if someone is thinking "BLAS", they are already about an order of magnitude behind on speed. The algorithms to do this effectively are much more complex than that.