How do I do a plain old "A=LU" decomposition of sparse matrices?

If you just need the LU factorization to solve a system of equations, you can certainly use it with the permutation term. Indeed, you can solve Ax=b with just x = lu(A) \ b and not worry about how it is implemented.

For example, if you look at the algorithm 2.4 in that paper (where you actually use the LU factorizations), it appears that you could implement w = P_F^t U_F^{-1} L_F^{-1} P_F u_F with something like:

w[p_F] = LU_F \ u_F[p_F]

where LU_F is your lu(...) factorization of P_F H_F P_F^t (= L_F U_F in algorithm 2.3) and p_F is the ordering corresponding to their permutation matrix P_F (which you also don’t explicitly construct).

In general, there is a certain amount of translation that has to take place when converting a mathematical algorithm into code. The example I always give my students is that whenever you see x = A^{-1} b, you should read it as “solve Ax=b by the fastest available method”, which typically does not involve computing A^{-1} explicitly.

8 Likes