Symbolic Factorization for Sparse Gaussian Elimination with Partial Pivoting

Abstract

Let $Ax = b$ be a large sparse nonsingular system of linear equations to be solved using Gaussian elimination with partial pivoting. The factorization obtained can be expressed in the form $A = P_1 M_1 P_2 M_2 \cdots P_{n - 1} M_{n - 1} U$, where $P_k $ is an elementary permutation matrix reflecting the row interchange that occurs at step k during the factorization, $M_k $ is a unit lower triangular matrix whose kth column contains the multipliers, and U is an upper triangular matrix. Consider the kth step of the elimination. Suppose we replace the structure of row k of the partially reduced matrix by the union of the structures of those rows which are candidates for the pivot row and then perform symbolically Gaussian elimination without partial pivoting. Assume that this is done at each step k, and let $\bar L$ and $\bar U$ denote the resulting lower and upper triangular matrices, respectively. Then the structures of $\bar L$ and $\bar U$, respectively, contain the structures of $\Sigma _{k = 1}^{n - 1} M_k $ and U. This paper describes an algorithm which determines the structures of $\bar L$ and $\bar U$, and sets up an efficient data structure for them. Since the algorithm depends only on the structure of A, the data structure can be created in advance of the actual numerical computation, which can then be performed very efficiently using the fixed storage scheme. Although the data structure is more generous than it needs to be for any specific sequence $P_1 ,P_2 , \cdots ,P_{n - 1} $, experiments indicate that the approach is competitive with conventional methods. Another important point is that the storage scheme is large enough to accommodate the $QR$ factorization of A, so it is also useful in the context of computing a sparse orthogonal decomposition of A. The algorithm is shown to execute in time bounded by $|\bar L| + |\bar U|$, where $|M|$ denotes the number of nonzeros in the matrix M.