Sparse matrix factorization is a black art and is best done with some existing out-of-the-box solution.<p>The problem is that they&#x27;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.<p>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 &quot;just don&#x27;t.&quot;<p>EDIT: The real reason is that FORTRAN is less apt to copy and more apt to reference than C++.