Algorithm

A linear system of equations is written in matrix terms as:

$$\mathsf{M}^a{}_b \mathsf{x}^b = \mathsf{v}^a \;,$$ (27.7)

where $\mathsf{M}^a{}_b$ is a known matrix, $\mathsf{v}^a$ a known vector, and $\mathsf{x}^b$ an unknown vector that we wish to solve for. A standard method for solving this is the method of LU decomposition, based on the fact that any non-singular square matrix $\mathsf{M}^a{}_b$ can be decomposed into the product of a lower-triangular matrix $\mathsf{L}^a{}_b$ and an upper-triangular matrix $\mathsf{U}^a{}_b$:

$$\mathsf{M}^a{}_b = \mathsf{L}^a{}_c \mathsf{U}^c{}_b = \lef... ...& \vdots \0 & 0 & \cdots & U^N{}_N \end{array}\right] \;.$$ (27.8)

The linear system can then be broken up as $\mathsf{L}^a{}_b \mathsf{y}^b = \mathsf{v}^a$ and $\mathsf{U}^b{}_c \mathsf{x}^c = \mathsf{y}^b$, where these two equations are trivially invertible:

$$y^i = v^i - \sum_{j=1}^i-1 L^i{}_j y^j \;,\qquad i=1,2,\ldots,N \;,$$ (27.9)

$$x^i = \frac{1}{U^i{}_i}\left( y^i - \sum_{j=i+1}^N U^i{}_j x_j \right) \;,\qquad i=N,N-1,\ldots,1 \;,$$ (27.10)

where the calculations are arranged so that the computation of a given $y^i$ or $x^i$ requires only those values that have been computed already. This process of solving the linear system is called backsubstitution.

The determinant of $\mathsf{M}^a{}_b$ is simply the product of the diagonal elements of the decomposed matrix: $\vert\mathsf{M}^a{}_b\vert=\prod_{i=1}^N U^i{}_i$. The inverse matrix $(\mathsf{M}^{-1}){}^a{}_b$ can be computed by performing a column-by-column backsubstitution of the identity matrix.

The routines in DetInverse.c simply call the routines in DetInverseInternal.c, first using LALSLUDecomp() or LALDLUDecomp() to perform an LU decomposition of the matrix, then either computing the determinant from the diagonal elements, or using LALSLUBackSub() or LALDLUBackSub() to determine the inverse by back-substitution of basis vectors. The routines that compute the determinant will also handle any ``singular matrix'' error code returned by the LU decomposition routines, returning zero as the determinant.

Since the diagonal components $L^i{}_i$ are conventionally assigned to 1, they do not need to be explicitly stored. Therefore we can store both matrices $\mathsf{L}^a{}_b$ and $\mathsf{U}^a{}_b$ ``in-place'', in the same memory block used for the input matrix $\mathsf{M}^a{}_b$. This the procedure taken by LALSLUDecomp() and LALDLUDecomp(); hence on return the routines in this module will leave the input *matrix in this decomposed state. However, these routines also permute the rows of the input matrix, and the information on this permutation is not returned by the routines in DetInverse.c, so the information in *matrix will be irretrievably mangled. If you want to do further work with the LU-decomposed matrix, call the routines in DetInverseInternal.c directly.

Computing the determinant is dominated by the cost of doing the LU decomposition, or of order $N^3/3$ operations. Computing the inverse requires an additional $N$ back-substitutions, but since most of the vector elements are zero, this reduces the average cost of each back-substitution from $N^2$ to $2N^2/3$, for a total operation count of $N^3$.