There are tons of solvers, there're a lot of libraries and SciPy, Fenics, most likely even MATLAB, Octave, FreeFEM++ all of them use LAPACK and BLAS (at least in their most performant flavours), all calling eventually the FORTRAN subroutines written decades ago. And yes, the tricky parts of this is also that how Fortran and C store 2D arrays, https://en.wikipedia.org/wiki/Row-_and_column-major_order.

So no one re-writes their own stuff, unless for academic purposes. I was on the verge of trying to extend the SciPy library to support block matrices 4 years ago for my thesis, but was able to work around it because duck-typing had generalised quite enough to encapsulate even that.

Sparse and Full matrices have a totally different data structure. hence the difference in the solvers. Engineers and Mathematicians that want more performant solving, rely more on the modelling of the problem (physically or mathematically) than improving the solvers.The other option is hardware-assisted optimization. For full solvers, I suppose the AI/ML folk have come up with TPUs, also because they have different tolerances for floating point calculations in classifying an image compared with launching a satellite on Mars.