Exact Methods

Most library routines, for example those in the NAG (Numerical Algorithms Group, n.d.) or LaPack (Lapack Numerical Library, n.d.) libraries are based on some variation of Gaussian elimination. The standard procedure is to call a first function which performs an $\bss{L}\bss{U}$ decomposition of the matrix $\bss{A}$,

$$\bss{A} = \bss{L}\bss{U}$$ (3.10)

where $\bss{L}$ and $\bss{U}$ are lower and upper triangular matrices respectively, followed by a function which performs the operations

$$\bss{Y} = \bss{L}^{-1}\bss{B}$$ (3.11)

$$\bss{X} = \bss{U}^{-1}\bss{Y}$$ (3.12)

on each column of $\bss{B}$. The procedure is sometimes described as Gaussian Elimination. A common variation on this procedure is partial pivoting. This is a trick to avoid numerical instability in the Gaussian Elimination (or $\bss{L}\bss{U}$ decomposition) by sorting the rows of $\bss{A}$ to avoid dividing by small numbers. There are several aspects of this procedure which should be noted:

• The $\bss{L}\bss{U}$ decomposition is usually done in place, so that the matrix $\bss{A}$ is overwritten by the matrices $\bss{L}$ and $\bss{U}$.
• The matrix $\bss{B}$ is often overwritten by the solution $\bss{X}$.
• A well written $\bss{L}\bss{U}$ decomposition routine should be able to spot when $\bss{A}$ is singular and return a flag to tell you so.
• Often the $\bss{L}\bss{U}$ decomposition routine will also return the determinant of $\bss{A}$.
• Conventionally the lower triangular matrix $\bss{L}$ is defined such that all its diagonal elements are unity. This is what makes it possible to replace $\bss{A}$ with both $\bss{L}$ and $\bss{U}$.
• Some libraries provide separate routines for the Gaussian Elimination and the Back-Substitution steps. Often the Back-Substitution must be run separately for each column of $\bss{B}$.
• Routines are provided for a wide range of different types of matrices. The symmetry of the matrix is not often used however.

The time taken to solve $N$ equations in $N$ unknowns is generally proportional to $N^3$ for the Gaussian-Elimination ( $\bss{L}\bss{U}$-decomposition) step. The Back-Substitution step goes as $N^2$ for each column of $\bss{B}$. As before this may be modified by parallelism.