Sparse matrix factorization is a black art and is best done with some existing out-of-the-box solution. The problem is that they're all these FORTRAN packages

Sparse matrix factorization is a black art and is best done with some existing out-of-the-box solution. The problem is that they're all these FORTRAN packages that one person wrote, and then everybody, in the name of human progress, reused and reused ever since -- akin to bootstrapping human intelligence. Nobody ever wants to write the Lanczos bidiagonalization algorithm with partial reorthogonalization from scratch, and when I considered writing my own C++ implementation to bring human progress to the 21st century, my academic advisor told me "just don't." EDIT: The real reason is that FORTRAN is less apt to copy and more apt to reference than C++.