A Brief Survey of Direct Linear Solvers
  J. Demmel, P. Koev, and X. Li

Many of the most effective methods in this book, in particular those using shift-and-invert, require the solution of linear systems of the form $(A-\sigma I)x=b$ or $(A - \sigma I)^*x=b$, where $\sigma$ is a shift, usually chosen to be an approximate eigenvalue of $A$. Much of the time in these methods can be spent solving these systems, so it is important to use the best methods available.

There are two basic approaches to solving $(A-\sigma I)x=b$: direct and iterative. Direct methods typically use variations on Gaussian elimination. They are effective even when $\sigma$ is close to an eigenvalue and $A-\sigma I$ is nearly singular, which is when iterative methods have most trouble converging. Section 10.4 discusses iterative solvers and when they can be effectively used in some eigenvalue algorithms, such as the Jacobi-Davidson method. This section considers only direct methods.

A direct method for solving $(A-\sigma I)x=b$ will consist of two steps:

1. Compute a factorization of $A-\sigma I$. This is the most expensive part.
2. Use the factorization to solve $(A-\sigma I)x=b$ or $(A - \sigma I)^*x=b$.

Step 1 usually costs much more than step 2 (e.g., $O(n^3)$ vs. $O(n^2)$ in the dense case). The factorization can be reused to solve $(A-\sigma I)x=b$ for many different right-hand sides $b$. The factorization only needs to be recomputed when the shift $\sigma$ is changed. Fortunately, most iterative methods solve many linear systems with different right-hand sides and the same shift, so step 1 is performed many fewer times than step 2.

The choice of algorithm for steps 1 and 2 depends on both the mathematical structure of $A-\sigma I$ and the storage structure of $A$. By mathematical structure we most often mean whether $A-\sigma I$ is Hermitian or non-Hermitian, and definite or indefinite.

Hermitian and definite:

In this case we compute the Cholesky factorization
$A - \sigma I = \pm LL^T$, where $L$ is a lower triangular matrix. The choice of sign depends on whether $A-\sigma I$ is positive definite ($+$) or negative definite ($-$). In the sparse case, we may instead compute $A - \sigma I = \pm PLL^*P^*$, where $P$ is a permutation matrix.

Hermitian and indefinite:

In this case we compute the Hermitian indefinite factorization $A - \sigma I = PLDL^*P^*$, where $L$ is lower triangular, $P$ is a permutation matrix, and $D$ is block diagonal with 1 by 1 or 2 by 2 blocks, or perhaps tridiagonal.

Non-Hermitian:

In this case we compute the triangular factorization $A - \sigma I = P_1 LU P_2$, where $L$ is lower triangular, $U$ is upper triangular, and $P_1$ and $P_2$ are permutations (sometimes one $P_i$ is omitted).

There are many other structures and factorizations as well, such as when $A$ is Toeplitz ($a_{i,j}$ depends only on $i-j$).

By storage structure we mean being dense, banded, sparse (in one of the formats described in §10.1), or structured (such as storing the first row and column of a Toeplitz matrix, which determines the other entries). A sparse non-Hermitian matrix may also be stored in a sparse Hermitian data structure, as in § 10.1.

There are specialized algorithms and software for many combinations of these mathematical and storage structures, and choosing the right algorithm can make orders of magnitude difference in solution time for large matrices. In this section we summarize the best algorithms and software available for these problems. It is often the case that the best algorithm in a research paper is not available as well designed and easily accessible software, and we shall concentrate on available software. We recommend software that is well maintained and in the public domain, or available at low cost, unless none is available.

We consider four cases: methods for dense matrices, methods for band matrices, methods for sparse matrices, and methods for structured matrices, such as Toeplitz matrices.