Algorithm

LU decomposition is performed by Crout's algorithm, described in Sec. 2.3 of [1]; the routines in this module are essentially re-implementations of the Numerical Recipes routines ludcmp() and lubksub(). For large $N$, their operation counts are approximately $N^3/3$ and $N^2$, respectively.

One difference between ludcmp() in [1] and the routines LALSLUDecomp() and LALDLUDecomp() in this module is the way in which singular matrices are handled. In [1], there is a distinction between between a manifestly singular matrix (where an entire row of the matrix is zero) and a numerically singular matrix (if a diagonal element in the decomposed matrix turns out to be zero). In the former case, they raise an error signal; in the latter, they replace the offending element with a ``tiny'' but nonzero number and continue. This treatment does not strike the present author as satisfactory.

Instead, the routines LALSLUDecomp() and LALDLUDecomp() will always return successfully, even with a singular matrix, but will not ``adjust away'' a numerical singularity. Instead, they will signal the presence of the singularity in two ways: First, they will set the permutation sign *sgn to zero; second, they will set all elements of the *indx vector yo zero. This ensures that routines computing the determinant (whose sign depends on *sgn) will correctly give a zero determinant, while the meaningless *indx provides a simple sanity check for routines such as LALSLUBackSub() and LALDLUBackSub() that attempt to invert the linear system. Note that the returned value of *matrix will be meaningless garbage.