Is there a version of this that supports sparse LU solves?

I’m trying to make sense of this question.

GEMMs are dense O(N^3) work operations that have roughly the same access pattern and data reuse properties across all matrices. Of course, I’m simplifying things a lot here; tall-skinny and short-fat patterns are much harder to get performance out of but the spirit of the approach is the same as big square matrices.

Sparse LU solves have a different character. There is nowhere near O(N^3) work. You typically expect something closer to O(N^2) but getting performance out of these operations is notoriously difficult because it depends a lot on the sparsity pattern of the linear system. Making matters worse is that you may commonly have a sparse A that factorises to dense L and/or U matrices.