Figure 9.12: Process-Grid Data Distribution of Ax=b. Representation of a concurrent matrix, and distributed-replicated concurrent vectors on a logical process grid. The solution of Ax=b first appears in x, a column-distributed vector, and then is normally ``transposed'' via a global combine to the row-distributed vector y.
We solve the problem Ax=b where A is large, and includes many zero entries. We assume that A is unsymmetric both in sparsity pattern and in numerical values. In general, the matrix A will be computed in a distributed fashion, so we will inherit a distribution of the coefficients of A (cf., Figures 9.12 and 9.13). Following the style of Harwell's MA28 code for unsymmetric sparse matrices, we use a two-phase approach to this solution. There is a first LU factorization called A-mode or ``analyze,'' which builds data structures dynamically, and employs a user-defined pivoting function. The repeated B-mode factorization uses the existing data structures statically to factor a new, similarly structured matrix, with the previous pivoting pattern. B-mode monitors stability with a simple growth factor estimate. In practice, A-mode is repeated whenever instability is detected. The two key contributions of this sparse concurrent solver are reduced communication pivoting, and new data distributions for better overall performance.
Figure 9.13: Example of Process-Grid Data Distribution. An array with block-linear rows (B=2) and scattered columns on a logical process grid. Local arrays are denoted at left by where is the grid position of the process on . Subscripts (i.e., ) are the global (I,J) indices.
Following Van de Velde [Velde:90a], we consider the LU factorization of a real matrix A, . It is well known (e.g., [Golub:89a], pp. 117-118), that for any such matrix A, an LU factorization of the form
exists, where are square, (orthogonal) permutation matrices, and are the unit lower- and upper-triangular factors, respectively. Whereas the pivot sequence is stored (two N-length integer vectors), the permutation matrices are not stored or computed with explicitly. Rearranging, based on the orthogonality of the permutation matrices, We factor A with implicit pivoting (no rows or columns are exchanged explicitly as a result of pivoting). Therefore, we do not store directly, but instead: . Consequently, , , and . The ``unravelling'' of the permutation matrices is accomplished readily (without implication of additional interprocess communication) during the triangular solves.
For the sparse case, performance is more difficult to quantify than for the dense case, but, for example, banded matrices with bandwidth can be factored with work; we expect subcubic complexity in N for reasonably sparse matrices, and strive for subquadratic complexity, for very sparse matrices. The triangular solves can be accomplished in work proportional to the number of entries in the respective triangular matrix L or U. The pivoting strategy is treated as a parameter of the algorithm and is not predetermined. We can consequently treat the pivoting function as an application-dependent function, and sometimes tailor it to special problem structures (cf., Section 7 of [Velde:88a]) for higher performance. As for all sparse solvers, we also seek subquadratic memory requirements in N, attained by storing matrix entries in linked-list fashion, as illustrated in Figure 9.14.
Figure 9.14: Linked-list Entry Structure of Sparse Matrix. A single entry consists of a double-precision value (8 bytes), the local row (i) and column (j) index (2 bytes each), a ``Next Column Pointer'' indicating the next current column entry (fixed j), and a ``Next Row Pointer'' indicating the next current row entry (fixed i), at 4 bytes each. Total: 24 bytes per entry.
For further discussion of LU factorizations and sparse matrices, see [Golub:89a], [Duff:86a].