Classical and detailed treatments of triangular solves that also cover sparse issues are given in the papers Brayton et al. (1970), Gilbert & Peierls (1988), and Gilbert (1994). For reducibility theory that is closely connected to the general theory of matrices, see Brualdi & Ryser (1991), which includes, for example, a proof of Theorem 3.4.

Algorithm 3.6 for computing strongly connected components of a digraph is introduced in Tarjan (1972); see also Sharir (1981) and Duff & Reid (1978) for an early implementation.

For identifying supervariables, Algorithm 3.7 follows Reid & Scott (1999), but see also Ashcraft (1995) and Hogg & Scott (2013a) (the latter presents an efficient variant that employs a stack). The approximate block partitioning of Section 3.5.2 is from the paper by Saad (2003a), which also describes some modifications of the basic approach; more sophisticated schemes with overlapping blocks are given in Fritzsche et al. (2013).

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Chapter 4 Sparse Cholesky Solver: The Symbolic Phase**

*The modern view of numerical linear algebra as being to a large extent the study and systematic use of matrix decompositions has certainly been influenced by Cholesky's posthumously published work – Benzi (2017).*

This chapter focuses on the symbolic phase of a sparse Cholesky solver. The sparsity pattern S{*A*} of the symmetric positive definite (SPD) matrix *A* is used to determine the nonzero structure of the Cholesky factor *L* without computing the numerical values of the nonzeros. The subsequent numerical factorization is discussed in the next chapter. Because the symbolic phase works only with S{*A*} (the values of the entries of *A* are not considered), it is also used for symmetric indefinite matrices and sometimes within LU factorizations of symmetrically structured nonsymmetric problems. It is implicitly assumed that all the diagonal entries of *A* are included in S{*A*} (even if they are zero). During the factorization phase, it may be necessary to amend the data structures to allow for indefiniteness. This makes the factorization of indefinite matrices potentially more expensive and more complex; this is considered further in Chapter 7.

A fundamental difference between dense and sparse Cholesky factorizations is that, in the latter, each column of *L* depends on only a subset of the previous columns. The elimination tree is a data structure that describes the dependencies among the columns of *A* during its factorization. A key result that assists in the understanding of sparse Cholesky factorizations is that the sparsity pattern of column *j* of *L* is the union of the pattern of column *j* of the lower triangular part of *A* and the patterns of the children of *j* in the elimination tree; this is shown in Section 4.3. Furthermore, the fact that disjoint parts of the elimination tree can be factored independently offers the potential for high-level tree-based parallelism that does not exist for dense matrices.

#### **4.1 Column Replication Principle**

We begin by looking at how the sparsity pattern of a computed column of *L* influences the patterns of subsequent Schur complements. From (3.2), the Schur complement *S(k)* can be written as ⎛⎜⎝⎞⎟⎠

This on suusequent \*\*scınır\*\* compnenens.\*\* From (5.2), une \*\*scınır\*\* in be written as\*\*

$$S^{(k)} = A\_{k:n,k:n} - \sum\_{j=1}^{k-1} \begin{pmatrix} l\_{kj} \\ \vdots \\ l\_{nj} \end{pmatrix} (l\_{kj} \dots l\_{nj}) \,. \tag{4.1}$$

Consider column *j* of *L* (1 ≤ *j* ≤ *k* − 1), and let *lij* = 0 for some *i>j* . The involvement of *lij* in the outer product in (4.1) allows the following observation.

**Observation 4.1** *For any i>j* ≥ 1 *such that lij* = 0

$$\mathcal{S}\{L\_{i:n,j}\} \subseteq \mathcal{S}\{L\_{i:n,i}\}.\tag{4.2}$$

*This is called the* **column replication principle** *because the pattern of column j of L (rows i to n) is replicated in the pattern of column i of L.*

Denote the row index of the first subdiagonal nonzero entry in column *j* of *L* by *parent (j )*, that is,

$$parent(j) = \min\{i \mid i > j \text{ and } l\_{lj} \neq 0\}.\tag{4.3}$$

If there is no such entry, set *parent (j )* = 0. The row index *parent (parent (j ))* is denoted by *parent*2*(j )*, and so on. Applying column replication recursively implies the sparsity pattern of column *j* of *L* is replicated in that of column *parent (j )*, which in turn is replicated in the pattern of column *parent*2*(j )*, and so on. This is illustrated in Figure 4.1. Here *j* = 1, and because the first subdiagonal entry in column 1 is in row 3, *parent (*1*)* = 3. Likewise, *parent (*3*)* = *parent*2*(*1*)* = 5.

**Figure 4.1** An illustration of column replication. On the left are the entries in *L* before step 1 of a Cholesky factorization (that is, the entries in the lower triangular part of *A*); in the centre, we show the replication of the nonzeros from column 1 in the pattern of column *parent (*1*)* = 3 (red entries *f* ); on the right, we show the subsequent replication in column *parent*2*(*1*)* = 5.

The following result shows that, provided *A* is irreducible, the mapping *parent (j )* has nonzero values given by (4.3) for all *j<n*.

**Theorem 4.1 (Liu 1986)** *If A is SPD and irreducible, then in each column j (*1 ≤ *j<n) of its Cholesky factor L there exists an entry lij* = 0 *with i>j .*

*Proof* From Parter's rule, each step of the Cholesky factorization corresponds to adding new edges into the graph of the corresponding Schur complement. If *A* is irreducible, then the graphs corresponding to the Schur complements are connected. Consequently, for any vertex *j* (1 ≤ *j<n*) in any of these graphs, there is at least one vertex *i* with *i>j* to which *j* is connected. This corresponds to the nonzero entry in column *j* of *L*.

With the convention *parent*1*(j )* = *parent (j )*, the next theorem shows that if entry *lij* of *L* is nonzero, then *parentt (j )* = *i* for some *t* ≥ 1, and there is an entry in row *i* of *L* in each of the columns in the replication sequence *j, parent*1*(j ), parent*2*(j ), . . . , parentt (j ).*

**Theorem 4.2 (Liu 1990; George 1998)** *Let A be SPD, and let L be its Cholesky factor. If lij* = 0 *for some j<i* ≤ *n, then there exists t* ≥ 1 *such that parentt (j )* = *i and lik* = 0 *for k* = *j , parent*1*(j ), parent*2*(j ), . . . , parentt (j ).*

*Proof* If *i* = *parent*1*(j )*, the result is immediate. Otherwise, there exists an index *k*, *j<k<i* of a subdiagonal entry in column *j* of *L* such that *k* = *parent*1*(j )*. Column replication implies *lik* = 0. Applying an inductive argument to *lik*, the result follows after a finite number of steps.

If there is a sequence of nonzeros in a row of *L*, it is natural to ask where the sequence begins. It is straightforward to see if there is no *k* ≥ 1 such that *aik* = 0, no replication of nonzeros can start in row *i*. The main result on the replication of nonzeros of *A* is summarized as Theorem 4.3.

**Theorem 4.3 (Liu 1986)** *Let A be SPD, and let L be its Cholesky factor. If aij* = 0 *for some* 1 ≤ *j<i* ≤ *n, then there is a filled entry lij* = 0 *if and only if there exist k<j and t* ≥ 1 *such that aik* = 0 *and parentt (k)* = *j .*

#### **4.2 Elimination Trees**

The discussion of column replication is significantly simplified using elimination trees. The **elimination tree** (or **etree**) T *(A)* (or simply T ) of an SPD matrix has vertices 1*,* 2*,...,n* and an edge between each pair *(j, parent (j ))*, where *parent (j )* is given by (4.3); *j* is a root vertex of the tree if *parent (j )* = 0. The edges of T are considered to be directed from a child to its parent, that is,

$$\mathcal{E}(\mathcal{T}) = \{ (j \longrightarrow i) \mid i = parent(j) \}.$$

**Figure 4.2** An illustration of a sparse matrix *A* with a symmetric sparsity pattern and its elimination tree T *(A)*. The root vertex is 8. The filled entries in S{*L* + *LT* } are denoted by *f* .

If T has a single component, then the root vertex is *n*. Despite the terminology, the elimination tree need not be connected and in general is a **forest**. For simplicity, in our discussions, we assume T has a single component, and we say that T is described by the vector *parent*.

An example of a matrix and its elimination tree is given in Figure 4.2. Here and elsewhere, following conventional notation, directional arrows are omitted from the tree plot.

Concepts such as child, leaf, ancestor, and descendant vertices introduced in Section 2.3 for directed rooted trees can be applied to T . Additionally, *anc*T {*j* } and *desc*T {*j* } are defined to be the sets of ancestors and descendants of vertex *j* in T . We denote by T *(j )* the **subtree** of T induced by *j* and *desc*T {*j )*; *j* is the root vertex of T *(j )*. The **size** |T *(j )*| is the number of vertices in the subtree. A **pruned subtree** of T *(j )* is the connected subgraph induced by *j* and a subset of *desc*T {*j )*. That is, for any vertex *i* in a pruned subtree of T *(j )*, all the ancestors of *i* belong to the pruned subtree. A pruned subtree of T shares the mapping *parent* with T .

The following observation is straightforward.

**Observation 4.2** *If i* ∈ *anc*T {*j* } *for some j* = *i, then i>j .*

The connection between the mapping *parent* and the sets of ancestors and descendants is emphasized by the next observation.

**Observation 4.3** *If i and j are vertices of the elimination tree* T *with j<i* ≤ *n, then*

*i*∈ *anc*T {*j* } *if and only if j* ∈ *desc*T {*i*} *if and only if parentt (j )*= *i for some t* ≥1*.*

The results in Section 4.1 can be expressed using rooted trees. Consider, for example, Theorem 4.2. Instead of stating that there exists *t* ≥ 1 such that *parentt (j )* = *i*, we can write that *i* ∈ *anc*T {*j* }. Rewriting Theorem 4.3 as the following corollary provides a clear characterization of the sparsity patterns of the rows of *L*.

**Figure 4.3** The row subtree T*r(*5*)* of the elimination tree T from Figure 4.2 (left). Vertex 3 has been pruned because *a*35 = 0. The row subtree T*r(*8*)* (right) differs from T = T *(A)* because vertex 1 has been pruned (*a*18 = 0).

**Corollary 4.4 (Liu 1986)** *Consider the elimination tree* T *and the Cholesky factor L of A. If i and j are vertices of* T *with j<i* ≤ *n and aij* = 0*, then lij* = 0 *if and only if there exists k<j such that j* ∈ *anc*T *(k) and aik* = 0*.*

The subtree of T with vertices that correspond to the nonzeros of row *i* of *L* is called the *i*-th **row subtree** and is denoted by T*r(i)*. Formally, it is a pruned subtree of T induced by the union of the vertex set

$$\{i\} \cup \{k \mid a\_{ik} \neq 0 \text{ and } k < i\}$$

with all vertices on the directed paths in T from *k* to *i*, that is, with all their ancestors from T*r(i)*. The root vertex is *i*, and the leaf vertices are a subset of the column indices in the *i*-th row of the lower triangular part of *A*. Figure 4.3 illustrates row subtrees for the matrix and elimination tree from Figure 4.2. Note that row subtrees are connected subgraphs of T , even if T is not connected. If T can be found without determining the pattern of *L*, then T*r(i)* can be used to derive the sparsity pattern of row *i* of *L*, without having to store each entry explicitly.

Theorem 4.5 characterizes the ancestors of a given vertex *j* using paths in G*(A)*. The proof helps clarify the relationship between T and paths in G*(A)*.

**Theorem 4.5 (Schreiber 1982; Liu 1986)** *If i and j are vertices in the elimination tree* T *with j<i* ≤ *n, then i* ∈ *anc*T {*j* } *if and only if there exists a path*

$$j \xleftarrow[{\mathcal{G}(A)}\_{\{1,\ldots,i\}}] \ i. \tag{4.4}$$

*Proof* Assume *i* ∈ *anc*T {*j* }. Then there is a path *j* T ⇒ *i* of length *l* ≥ 1. Each edge of this path belongs to G*(L)* and corresponds either to an edge in G*(A)* or to a fill-path in G*(A)*. Connecting these paths together gives (4.4).

Conversely, if the path (4.4) exists, induction on its length can be used to prove the result. If the path is of length 1, then the result holds because *i* and *j* are connected in G*(A)* by an edge. Consequently, from Theorem 4.2, *i* is an ancestor of *j* . Now assume that the result is true for all paths of length less than *l* (*l >* 1), and consider a path of length *l*. Let *m* be the largest vertex on this path. If *m<j* , then (4.4) is a fill-path connecting *i* and *j* and, therefore, *i* ∈ *anc*T {*j* }. Otherwise, for *m* ≥ *j* , the assumption implies *i* ∈ *anc*T {*m*}∪{*m*} and *m* ∈ *anc*T {*j* }∪{*j* }, that is, *i* ∈ *anc*T {*j* }.

Given a vertex *j* in T , the following corollary indicates how to find *parent (j )* (if it exists). If the set of ancestors of *j* is non-empty, then the lowest numbered one is its parent.

**Corollary 4.6 (Liu 1986, 1990)** *Vertex i is the parent of vertex j in* T *if and only if i is the lowest numbered vertex satisfying j<i* ≤ *n for which there is a path (4.4).*

The existence of (4.4) is equivalent to requiring *i* and *j* belong to the same component of the graph G*(A*1:*i,*1:*i)* corresponding to the *i* × *i* principal leading submatrix *A*1:*i,*1:*i* of *A*. Figure 4.4 depicts G*(A)* for the matrix *A* given in Figure 4.2. Consider vertex 4. Its set of ancestors for which paths from Theorem 4.5 exist comprises vertices 5, 6, and 8. Vertex 7 is not an ancestor of 4 because there is no path from 7 to 4 in the graph G*(A*1:7*,*1:7*)*. Among the ancestors of 4, vertex 5 fulfils the condition from Corollary 4.6 and is thus the parent of 4.

T = T *(A)* can be constructed by stepwise extensions of the elimination trees of the principal leading submatrices of *A*. Assume we have T *(A*1:*i*−1*,*1:*i*−1*)* and we want to construct T *(A*1:*i,*1:*i)*. Initialize T *(A*1:*i,*1:*i)* = T *(A*1:*i*−1*,*1:*i*−1*)*. If there are no entries in row *i* of *A* to the left of the diagonal, then there is nothing to do, and only an isolated vertex *i* is added. Otherwise, *i* is the root of the row subtree T*r(i)* and an ancestor of some vertex *j* in T . The ancestors *k* of *j* with *k<i* are in T *(A*1:*i*−1*,*1:*i*−1*)*. Because row subtrees are connected subgraphs of T , a directed path in T *(A*1:*i,*1:*i)* with *parentt (j )* = *i* exists for some *t* ≥ 1. The search for this path starts from *j root* = *j* and continues, while *parent (j root)* = 0 and *parent (j root)* = *i*, using a sequence of assignments *j root* = *parent (j root)*. It terminates once *parent (j root)* = *i* or *i* is found to have already been added when

**Figure 4.4** The graph G*(A)* of the matrix from Figure 4.2 illustrating Theorem 4.5 and Corollary 4.6.

tracing the path from another entry *j* in row *i*. The construction of T is summarized in Algorithm 4.1.

The most expensive part of Algorithm 4.1 is the **while** loop that searches for subtree roots. Because the directed path from *j* to its root *parentt (j )* is unique, shortcuts can be incorporated; this is called **path compression**. Having found a directed path from *j* to *k*, subsequent searches can be made more efficient by introducing a vector *ancestor* and setting *ancestor(j )* = *k*. The modified algorithm is outlined in Algorithm 4.2. It maintains two structures using the current values of *parent* and *ancestor*. The tree described by *ancestor* is termed the **virtual tree**.

Figure 4.5 shows a matrix for which path compression makes constructing T significantly more efficient. For this example, T is determined by the mapping *parent (*6*)* = 0; *parent (i)* = *i* + 1 for *i* = 1*,...,* 5. The complexity of Algorithm 4.1 is *O(n*2*)*, but for this example the complexity of Algorithm 4.2 is *O(n)*. Formally, the complexity of Algorithm 4.2 is *O(nz(A)*log2*(n))*, where *nz(A)* is the number of nonzeros of *A*, but the logarithmic factor is rarely reached. Additional modifications can reduce the theoretical complexity to *O(nz(A) g(nz(A), n))*, where *g(nz(A), n)* is a very slowly increasing function called the functional inverse of Ackermann's function. This means that, in practice, the complexity of constructing T , and hence of obtaining an implicit representation of S{*L*}, is close to linear in *nz(A)* (which in general is much smaller than *nz(L)*).

**ALGORITHM 4.2 Construction of an elimination tree using path compression Input:** *A* with a symmetric sparsity pattern and its undirected graph G. **Output:** Elimination tree T described by the vector *parent*.

1: **for** *i* = 1 : *n* **do** Loop over the rows of *A* 2: *parent (i)* = 0, *ancestor(i)* = 0 Initialisation 3: **for** *j* ∈ *adj*G{*i*} and *j<i* **do** Loop over the below diagonal entries in row *i* 4: *j root* = *j* 5: **while** *ancestor(j root)* = 0 and *ancestor(j root)* = *i* **do** 6: *l* = *ancestor(j root)* 7: *ancestor(j root)* = *i* Path compression to accelerate future searches 8: *j root* = *l* 9: **end while** 10: **if** *ancestor(j root)* = 0 **then** 11: *ancestor(j root)* = *i* and *parent (j root)* = *i* 12: **end if** 13: **end for** 14: **end for**

```
⎛
⎜⎜⎜⎝
 ∗∗∗∗∗∗
 ∗ ∗
 ∗ ∗
 ∗ ∗
 ∗ ∗
 ∗ ∗
         ⎞
         ⎟⎟⎟⎠
```
**Figure 4.5** A sparse matrix for which computing the elimination tree using Algorithm 4.2 is much more efficient than using Algorithm 4.1.

The following simple theorem states that there is no edge in G*(L* + *LT )* between vertices belonging to subtrees of T with different vertex sets. If there was such an edge *(s, t)*, then from Theorem 4.2, one of the vertices *s* and *t* must be an ancestor of the other, which is a contradiction. The importance of this result is that it implies that for any such pairs of vertices the corresponding column sparsity patterns in *L* can be computed in parallel.

**Theorem 4.7 (Liu 1990)** *Consider the elimination tree* T *and the Cholesky factor L of A. Let* T *(i) and* T *(j ) be two vertex-disjoint subtrees of* T *. Then for all s* ∈ T *(i) and t* ∈ T *(j ), the entry lst of L is zero.*

#### **4.3 Sparsity Pattern of** *L*

The explicit structure of *L* is not always required; sometimes only the numbers of nonzeros in each row and column of *L* are needed. For example, when comparing the amount of fill-in in the factors for different initial orderings of *A*, allocating factor storage, finding relaxed supernodes (see Section 4.6), and determining load balance and synchronization events in parallel factorizations.

Let *rowL*{*i*} denote the sparsity pattern of the off-diagonal part of row *i* of *L*, that is,

$$row\_L\{i\} = \mathcal{S}\{L\_{i,1:i-1}\} = \{j \mid j < i, \ l\_{lj} \neq 0\}, \quad 1 \le i \le n.$$

The number of entries in *L* is

$$\begin{aligned} \left| \begin{array}{l} \left( \left( \left( \right)^{n} \right)^{-1} \right)^{-1} \right| & \left( \left( \left( \right)^{n} \right)^{-1} \right) \\\\ \text{is} \\\\ \left| nz(L) = \sum\_{i=1}^{n} \left| row\_{L}\{i\} \right| + n \dots \end{aligned} $$

Corollary 4.4 implies *rowL*{*i*} is given by the vertices of the row subtree T*r(i)*. This suggests Algorithm 4.3. Here the vector *mark* is used to flag vertices so as to avoid including them more than once within a row subtree. The complexity of the algorithm is *O(nz(L))*.

#### **ALGORITHM 4.3 Computation of the row sparsity patterns of the Cholesky factor** *L*

**Input:** *A* with a symmetric sparsity pattern, its undirected graph G and elimination tree T described by the vector *parent*.

**Output:** Row sparsity patterns *rowL*{*i*} of the Cholesky factor *L* of *A* (1 ≤ *i* ≤ *n*).

1: **for** *i* = 1 : *n* **do** Loop over the rows of *A* 2: *rowL*{*i*}=∅ Initialisation 3: *mark(i)* = *i* 4: **for** *k* ∈ *adj*G{*i*} and *k<i* **do** Loop over the below diagonal entries in row *i* 5: *j* = *k* 6: **while** *mark(j )* = *i* **do** Column *j* not yet encountered in row *i* 7: *mark(j )* = *i* Flag *j* as encountered in row *i* 8: *rowL*{*i*} = *rowL*{*i*}∪{*j* } Add *j* to the sparsity pattern of row *i* 9: *j* = *parent (j )* Move up the elimination tree 10: **end while** 11: **end for** 12: **end for**

**Figure 4.6** An illustration of the sparsity pattern of *A* and its graph G*(A)* (left) and the sparsity pattern of the corresponding skeleton matrix *A*− and graph G*(A*−*)* (right). The entries in *A* and edges of G*(A)* that do not belong to the skeleton matrix and graph are depicted in red.

Efficiency can be improved by employing the **skeleton graph** G*(A*−*)* that is obtained from G*(A)* by removing every edge *(i, j )* for which *j<i* and *j* is not a leaf vertex of T*r(i)*. G*(A*−*)* is the smallest subgraph of G*(A)* with the same filled graph as G*(A)*. The corresponding matrix is the **skeleton matrix**. An example is given in Figure 4.6. The complexity of constructing the elimination tree using the skeleton matrix and its graph G*(A*−*)* is *O(nz(A*−*) g(nz(A*−*), n))*, where *nz(A*−*)* is the number of entries in the skeleton matrix. Because *nz(A*−*)* is often significantly smaller than *nz(A)*, an implementation that processes G*(A*−*)* rather than G*(A)* can be substantially faster.

Analogously to the row sparsity patterns, let *colL*{*j* } denote the sparsity pattern of the off-diagonal part of column *j* of *L*, that is,

$$col\_L\{j\} = \mathcal{S}(L\_{j+1:n,j}) = \{i \mid i > j, \ l\_{lj} \neq 0\}, \quad 1 \le j \le n.$$

The column replication principle can be written as

$$col\_L\{j\} \subseteq col\_L\{parent(j)\} \cup parent(j).$$

Theorem 4.8 describes *colL*{*j* } using the vertices of the subtree T *(j )*.

**Theorem 4.8 (George & Liu 1980c, 1981)** *The column sparsity pattern colL*{*j* } *of the Cholesky factor L of the matrix A is equal to the adjacency set of vertices of the subtree* T *(j ) in* G*(A), that is,*

**Figure 4.7** Two topological orderings of an elimination tree.

$$col\_L\{j\} = adj\_{\mathcal{G}(A)}\{\mathcal{T}(j)\}.\tag{4.5}$$

*Proof* If *i* ∈ *colL*{*j* }, then *j* ∈ *rowL*{*i*}, and Theorem 4.3 implies *j* ∈ *anc*T {*k*} for some *k* such that *aik* = 0. That is, *i* ∈ *adj*G{T *(j )*}. Conversely, *i* ∈ *adj*G{T *(j )*} implies that in row *i* the entry in column *j* of *L* is nonzero. Thus, *j* ∈ *rowL*{*i*}, and hence, *i* ∈ *colL*{*j* }.

Algorithm 4.3 can be used to compute the column counts and the column sparsity patterns because when *j* is added to *rowL*{*i*} at line 8, *i* can be added to *colL*{*j* }. This does not generally obtain the column sparsity patterns sequentially. To derive an approach that does compute them sequentially, rewrite (4.5) as follows: ⎛⎝⎞⎠

$$\begin{aligned} & \text{is not generally obtain the column sparsity patterns sequentially.}\\ & \text{each that does compute them sequentially, rewrite } (4.5) \text{ as follows:}\\ & \begin{pmatrix} \\ \end{pmatrix} = \left( \underset{\{k \mid k \in \overline{\mathcal{T}}(j) \mid \{j\}\}}{\text{adj}} \underset{\{k \mid \}}{\text{col}} \right) \end{aligned}$$

Using the column replication, this can be significantly simplified ⎝⎠

$$\text{g the column relation, this can be significantly simplified}$$

$$col\_L\{j\} = \left(adj\_{\mathcal{G}(A)}\{j\} \bigcup\_{\{k \mid j = parent(k)\}} col\_L\{k\}\right) \tag{4.6}$$

This is used to obtain Algorithm 4.4, which constructs the sparsity pattern of each column *j* of *L* as the union of the sparsity pattern of column *j* of *A* (*adj*G*(A)*{*j* }) and the patterns of the children of *j* in T *(A)*. Here *child*{*j* } denotes the set of children of *j* . Because any child *k* of *j* satisfies *k<j* , the *j* -th outer step has the information needed to compute the sparsity pattern described by (4.6). Observe that T *(A)* does not need to be input.

### **ALGORITHM 4.4 Determining the sparsity patterns of each column of** *L* **Input:** *A* with symmetric sparsity pattern and its undirected graph G. **Output:** Column sparsity patterns *colL*{*j* } of the Cholesky factor *L* of *A* (1 ≤ *j* ≤ *n*).

```
1: for j = 1 : n do  Loop over the columns of L
2: child{j }=∅  Initialisation
3: colL{j } = adjG{j }\{1,...,j − 1}
4: for k ∈ child{j } do  Unifying child structures in (4.6)
5: colL{j } = colL{j } ∪ colL{k}\{j }
6: end for
7: if colL{j }
= ∅ then
8: l = min{i | i ∈ colL{j }}
9: child{l} = child{l}∪{j }  Parent of j detected using Corollary 4.6
10: end if
11: end for
```
### **4.4 Topological Orderings**

The outer loop in Algorithm 4.4 does not have to be performed in the strict order *j* = 1*,...,n*. What is necessary is that for each step *j* , the column sparsity pattern for each child of *j* has already been computed. An ordering of the vertices in a tree (and, more generally, in a DAG) is a topological ordering if, for all *i* and *j* , *j* ∈ *desc*T {*i*} implies *j<i* (Section 2.2). Observation 4.2 confirms that the ordering of vertices in the elimination tree T is a topological ordering. A new topological ordering of T defines a relabelling of its vertices corresponding to a symmetric permutation of *A*. This is illustrated in Figure 4.7. The sparsity patterns of the Cholesky factors of *A* and *P APT* can be different, but the following result shows that the amount of fill-in is the same.

**Theorem 4.9 (Liu 1990)** *Let* S{*A*} *be symmetric. If P is the permutation matrix corresponding to a topological ordering of the elimination tree* T *of A, then the filled graphs of A and P APT are isomorphic.*

There are many topological orderings of T . One class is obtained using the depthfirst search given by Algorithm 2.1. This searches all the components of T starting at their root vertices. In this case, once vertex *i* has been visited, all the vertices of the subtree T *(i)* are visited immediately after *i* and *i* is labelled as the last vertex of T *(i)*. A topological ordering of T is a **postordering** if the vertex set of any subtree T *(i)* (*i* = 1*,...,n*) is a contiguous sublist of 1*,...,n*. Unless additional rules on how vertices are selected are imposed, a postordering is generally not unique, as demonstrated in Figure 4.8. One possible postordering is defined in Algorithm 2.1. In this case, there is some freedom in the depth-first search to choose from the vertices that have not been visited, resulting in different postorderings.

**Figure 4.8** An example to illustrate the non-uniqueness of postorderings of an elimination tree.

#### **4.5 Leaf Vertices of Row Subtrees**

Leaf vertices of row subtrees play a key role in graph algorithms related to sparse Cholesky factorizations. They can be used to find the skeleton matrix described in Section 4.3, and they are important in parallel processing based on fundamental supernodes (see Section 4.6.1). Theorem 4.10 describes the relation between standard subtrees of T and row subtrees obtained by pruning (Section 4.2). This pruning is determined by the leaf vertices of row subtrees.

**Theorem 4.10 (Liu 1986)** *Let the elimination tree* T *of A be postordered. Let the column indices of the nonzeros in the strictly lower triangular part of row i of A be c*1*,...,cs with s* ≥ 1 *and* 0 *< c*1 *<...< cs < i. Then ct is a leaf vertex of the row subtree* T*r(i) if and only if*

$$t = 1 \quad or \ \ (1 < t \le s \quad \text{and} \ c\_{t-1} \notin \mathcal{T}(c\_t)).$$

*Proof c*1 is always a leaf vertex of T*r(i)*. If this is not the case, then there exists a directed path from some vertex *k, k* = *c*1 to *i* via *c*1 such that *k* ∈ T*r(i)* and *aik* = 0. Postordering of T implies *k<c*1. This is a contradiction because *c*1 is the index of the first nonzero in row *i*.

Consider now *t >* 1. Assume that *ct*−1 ∈ T *(ct)* and that *ct* is a leaf vertex of T*r(i)*. Row replication (Theorem 4.2) implies any *k* ∈ *anc*T {*ct*−1}∪{*ct*−1} such that *ct*−1 ≤ *k<i* satisfies *lik* = 0. Because T is postordered, *ct*−1 ≤ *k* ≤ *ct* , and there is at least one *k<ct* satisfying this inequality. It follows that *k* = *ct*−1. Because *k* belongs to T*r(i)*, *ct* cannot be a leaf vertex of T*r(i)*, which is a contradiction.

Conversely, assume for *t >* 1 that *ct*−1 ∈ T *(ct)* and *ct* is not a leaf vertex of T*r(i)*. From the second part of the assumption and the fact that *ct* ∈ T*r(i)*, it follows that there is at least one leaf vertex *k<i* of T*r(i)* from which there is a directed path to *i* via *ct* . Thus *k<ct* . From the definition of the postordering of T , all vertices *l* with *k<l* ≤ *ct* are vertices of T *(ct)*. Vertex *ct*−1 must be among them and *ct*−1 ∈ T *(ct)*. This contradiction completes the proof. **ALGORITHM 4.5 Find the sizes of subtrees** T *(i)* **of** T **Input:** Elimination tree T described by the vector *parent*. **Output:** Subtree sizes |T *(i)*| (1 ≤ *i* ≤ *n*).

1: |T *(*1 : *n)*| = 1 2: **for** *i* = 1 : *n* − 1 **do** 3: *k* = *parent (i)* 4: |T *(k)*|=|T *(k)*|+|T *(i)*| 5: **end for**

**Corollary 4.11 (Liu 1986)** *Under the assumptions of Theorem 4.10, ct is a leaf vertex of* T*r(i) if and only if*

$$t = 1 \text{ or } \begin{array}{c} \text{or} \ \vert \; 1 < t \le s \end{array} \text{ and } \begin{array}{c} c\_{l-1} < c\_l - \vert \mathcal{T}(c\_l) \vert + 1 \text{)}.$$

Subtree sizes can be computed using Algorithm 4.5. Correctness of Algorithm 4.5 is guaranteed because *parent* defines a topological ordering of T .

Theorem 4.12 relaxes the condition that the entries in the rows of *A* are sorted by increasing column indices. This allows the leaf vertices of the row subtrees to be determined by columns.

**Theorem 4.12 (Liu et al. 1993)** *Consider the elimination tree* T *of A. Vertex j is a leaf vertex of some row subtree of* T *if and only if there exists i* ∈ *adj*G*(A)*{*j* }*, j<i* ≤ *n, such that i* ∈ *adj*G*(A)*{*k*} *for all k* ∈ T *(j )* \ {*j* }*.*

*Proof* Assume that for some *i* ∈ *anc*T {*j* } vertex *j* is a leaf vertex of T*r(i)*. That is, *i* ∈ *adj*G*(A)*{*j* }, *i>j* . Suppose there exists *k* in T *(j )*\{*j* } such that *i* ∈ *adj*G*(A)*{*k*}. Then all the ancestors of *k, k* ≤ *i*, in particular *j* , belong to T*r(i)* and *j* cannot be a leaf vertex of T*r(i)*. This is a contradiction.

Conversely, assume that *j* is not a leaf vertex of any row subtree of T and that there exists *i* ∈ *adj*G*(A)*{*j* }, *j<i* ≤ *n*, such that *i* ∈ *adj*G*(A)*{*k*} for all *k* ∈ T *(j )* \ {*j* }. Because *j* is not a leaf vertex of any such T*r(i)*, Theorem 4.3 implies that there exists *k* ∈ T *(j )* \ {*j* } such that *aik* = 0, which gives a contradiction and completes the proof.

To find leaf vertices of row subtrees of T , Algorithm 4.6 uses a marking scheme based on Theorem 4.12 and exploits the postordering of T . The auxiliary vector *prev*\_*nonz* stores the column indices of the most recently encountered entries in the rows of the strictly lower triangular part of *A*.

#### **4.6 Supernodes and the Assembly Tree**

Because of column replication, the columns of *L* generally become denser as the Cholesky factorization proceeds. Exploiting this density can significantly enhance

## **ALGORITHM 4.6 Find leaf vertices of row subtrees of** T

**Input:** *A* with a symmetric sparsity pattern and a corresponding postordered elimination tree T .

**Output:** Logical vector *isleaf* with entries set to true for leaf vertices of row subtrees.

1: *isleaf (*1 : *n)* = *f alse*, *prev*\_*nonz(*1 : *n)* = 0 2: Compute |T *(*1 : *n)*| Use Algorithm 4.5 3: **for** *j* = 1 : *n* **do** Loop over the columns of *A* 4: **for** *i* such that *i>j* and *aij* = 0 **do** Row index in strictly lower triangular part of *A* 5: *k* = *prev*\_*nonz(i)* Column index of most recently seen entry in row *i* 6: **if** *k<j* − |T *(j )*| + 1 **then** 7: *isleaf (j )* = *true j* is a leaf vertex by Corollary 4.11 8: **end if** 9: *prev*\_*nonz(i)* = *j* Flag *j* as the most recently seen entry in row *i* 10: **end for** 11: **end for**

the performance of the numerical factorization in terms of both computation time and memory requirements. For this, we require the concept of supernodes. The idea is to group together columns with the same sparsity structure, so that they can be treated as a dense matrix for storage and computation. Let 1 ≤ *s,t* ≤ *n* with *s* + *t* − 1 ≤ *n*. A set of contiguously numbered columns of *L* with indices *S* = {*s, s* + 1*,...,s* + *t* − 1} is a **supernode** of *L* if

$$col\_L\{\mathbf{s}\} \cup \{\mathbf{s}\} = col\_L\{\mathbf{s} + t - 1\} \cup \{\mathbf{s}, \dots, \mathbf{s} + t - 1\},\tag{4.7}$$

and *S* cannot be extended for *s >* 1 by adding *s* − 1 or for *s* + *t* − 1 *< n* by adding *s* + *t*. Because *S* cannot be extended, it is a **maximal** subset of column indices. In graph terminology, a supernode is a **maximal clique** of contiguous vertices of G*(L* + *LT )*. A supernode may contain a single vertex. Figure 4.9 illustrates the supernodes in a Cholesky factor of order 8.

The **supernodal elimination** or **assembly tree** is defined to be the reduction of the elimination tree that contains only supernodes. Each vertex of the elimination tree is associated with one elimination, and a single integer (the index of its parent) is needed. Associated with each vertex of the assembly tree is an index list of the row indices of the nonzeros in the columns of the supernode. These implicitly define the sparsity pattern of *L*. An example that demonstrates the difference between the elimination and assembly trees is given in Figure 4.10. Here the elimination tree is postordered, and there are 5 supernodes: {1*,* 2}, 3, 4, 5, {6*,* 7*,* 8*,* 9}. For supernode 1 that comprises columns 1 and 2, the row index list is {1*,* 2*,* 8*,* 9}.

$$L = \begin{array}{ccccc} & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \\ & & 2 & \begin{pmatrix} \* & & & & & & \\ \* & & \* & & & & & \\ & & \* & & & & & \\ & & & \* & & & & \\ & & \* & \* & \* & & & \\ \* & \* & \* & \* & \* & \* & \\ & & & \* & \* & \* & \* & \\ & & & \* & \* & \* & \* & \* \\ & & & \* & \* & \* & \* & \* & \* \\ \end{array} \end{array}$$

**Figure 4.9** An example to illustrate supernodes in *L*. The first supernode comprises columns 1 and 2, the second columns 3 and 4, and the third columns 5–8.

**Figure 4.10** A sparse matrix and its postordered elimination tree (left) and postordered assembly tree (right). Filled entries in S{*L* + *LT* } are denoted by *f* . For the assembly tree, the vertices are in red and the index lists associated with each vertex are given.

Supernodes can be characterized by the following result on the column counts of *L*, from which we see that supernodes can be found using column counts rather than the column sparsity patterns that appear in (4.7).

**Theorem 4.13 (Liu et al. 1993)** *The set of columns of L with indices S* = {*s, s* + 1*,...,s* + *t* − 1} *is a supernode of L if and only if it is a maximal set of contiguous columns such that s* + *i* − 1 *is a child of s* + *i for i* = 1*,...,t* − 1 *and*

$$\left|col\_L\{\mathbf{s}\}\right| = \left|col\_L\{\mathbf{s} + t - 1\}\right| + t - 1. \tag{4.8}$$

*Proof* Let *S* be a supernode. For *i, j* ∈ *S* with *i>j* , we have *i* ∈ *colL*{*j* }. This implies that in the postordered elimination tree the vertex *i* = *j* + 1 is the parent of *j* for *j* = *s,... , s* + *t* − 2. Moreover, from Observation 4.2, for any *i, j* ∈ *S* with *i>j* , *i* ∈ *colL*{*j* } implies *colL*{*j* }\{1*,...,i*} ⊆ *colL*{*i*}*.* Therefore,

$$|\operatorname{col}\_{L}\{\mathbf{s} + i\}| \ge |\operatorname{col}\_{L}\{\mathbf{s} + i - 1\}| - 1, \quad i = 1, \ldots, t - 1,\tag{4.9}$$

with equality if and only if

$$col\_L\{\mathbf{s} + i\} = col\_L\{\mathbf{s} + i - 1\} \mid \{\mathbf{s} + i\},$$

that is, if *S* is a supernode.

Conversely, assume *S* is a maximal set of contiguous columns such that, for *i* = 1*,...,t* − 1, *s* + *i* − 1 is a child of *s* + *i* and *S* satisfies (4.8). Because of column replication, such a sequence of parent and child vertices must satisfy (4.9) with equality if and only if (4.7) is satisfied. It follows that *S* is a supernode.

Supernodes enhance the efficiency of sparse factorizations and sparse triangular solves because they enable floating-point operations to be performed on dense submatrices rather than on individual nonzeros, thus improving memory hierarchy utilization and allowing the use of highly efficient dense linear algebra kernels (such as Level 3 BLAS kernels). Because the rows and columns of a supernode have a common sparsity structure, this only needs to be stored once, reducing indirect addressing. Supernodes help to increase the granularity of tasks, which is useful for improving the computation to overhead ratio in a parallel implementation. Fill-in results in supernodes near the root of the assembly tree often being much larger than those close to the leaf vertices.

Observe that the columns within a supernode are numbered consecutively, but they can be numbered within the supernode in any order without changing the number of nonzeros in *L* (assuming the corresponding rows are permuted symmetrically). On some architectures, particularly those using GPUs, this freedom can be exploited to improve the factorization efficiency. Essentially, it is desirable to order the columns within a supernode such that the entries of *L* form fewer but less fragmented dense blocks.

Some applications, such as power grid analysis, in which the basis of the linear system is not a finite element or finite difference discretization of a physical domain, can lead to sparse matrices that incur very little fill-in during factorization. The supernodes can then be very small, and the costs associated with identifying them may not be offset by the increase in performance resulting from the potential for block operations. However, as supernodes can offer such significant performance gains, it can be advantageous to merge (small) supernodes that have similar (but not exactly the same) nonzero patterns, despite this increasing the overall fill-in and operation count. This process is termed **supernode amalgamation**, and the resultant nodes are often referred to as **relaxed supernode**.

#### *4.6.1 Fundamental Supernodes*

In practice, fundamental supernodes are easier to work with in the numerical factorization. Let 1 ≤ *s,t* ≤ *n* with *s* + *t* − 1 ≤ *n*. A maximal set of contiguously numbered columns of *L* with indices *S* = {*s, s* +1*,...,s* +*t* −1} is a **fundamental supernode** if for any successive pair *i* −1 and *i* in the list, *i* −1 is the only child of *i* in T and *colL*{*i*}∪{*i*} = *colL*{*i* − 1}. *s* is termed the starting vertex. An example is given in Figure 4.11. The difference between the sets of supernodes and fundamental supernodes is normally not large, with the latter having (slightly) more blocks in the resulting partitioning of *L*. Note that fundamental supernodes are independent of the choice of the postordering of T . Theorem 4.14 describes the relationship between fundamental supernodes and the leaf vertices of row subtrees of T . In particular, it characterizes starting vertices of the fundamental supernodes. The leaf vertices of T are trivially starting vertices of fundamental supernodes. But, possibly surprisingly, so too are the leaf vertices of row subtrees.

**Theorem 4.14 (Liu et al. 1993)** *Assume* T *is postordered. Vertex s is the starting vertex of a fundamental supernode if and only if it has at least two child vertices in* T *or it is a leaf vertex of a row subtree of* T *.*

*Proof* If *s* has at least two child vertices then, from the definition of a fundamental supernode, it must be the starting vertex of a fundamental supernode. Assume that, for some *i>s*, *s* is a leaf vertex of T*r(i)*. If *s* is also a leaf vertex of T , then *s* is a starting vertex of a supernode. The remaining case is *s* having only one child. Because T is postordered, this child must be *s* − 1. Theorem 4.3 then implies *ais* = 0 and *ai,s*−1 = 0*,* that is, *i* ∈ *colL*{*s*} and *i /*∈ *colL*{*s* − 1}. It follows that

$$
\mathcal{S}\{L\_{s-1:n,s-1}\} \subsetneq \mathcal{S}\{L\_{s:n,s}\} \cup \{s-1\},
$$

and vertices *s* and *s* − 1 cannot belong to the same supernode. Hence, *s* is the starting vertex of a new fundamental supernode.

**Figure 4.11** A matrix *A* and its postordered elimination tree T for which the set of supernodes {1*,* 2} and {3*,* 4*,* 5*,* 6} and the set of fundamental supernodes {1*,* 2}*,*{3*,* 4} and {5*,* 6} are different. The filled entries in S{*L* + *LT* } are denoted by *f* .

Conversely, assume that *s* is the starting vertex of a fundamental supernode *S*. If *s* has no child vertices or at least two child vertices, the result follows. If *s* has exactly one child vertex, postordering implies this child is *s* − 1. Because *S* is maximal, there exists *i* such that *i* ∈ *colL*{*s* − 1} and *i* ∈ *colL*{*s*} (otherwise *S* could be extended by adding *s* − 1). Hence, *s* is a leaf vertex of T*r(i)*.

Because fundamental supernodes are characterized by their starting vertices, they can be found by modifying Algorithm 4.6 to incorporate marking leaf vertices of the row subtrees and vertices with at least two child vertices. Once the elimination tree has been computed, the complexity is *O(n*+*nz(A))*. The computation can be made even more efficient by using the skeleton graph G*(A*−*)*.

#### **4.7 Notes and References**

The excellent monographs by Tewarson (1973), George & Liu (1981), and Davis (2006) represent milestones in the development of contemporary symbolic factorization algorithms and their implementation. A complementary way to follow many of the developments is by looking at the early software (and accompanying user documentation), such as YSMP (Eisenstat et al., 1982) and SPARSPAK (George & Ng, 1984). In addition, there are several influential survey articles focusing on sparse Cholesky algorithms and emphasizing the crucial role of the elimination tree, for example, Liu (1990), George (1998); see also Bollhöfer & Schenk (2006), Hogg & Scott (2013a) and the more recent comprehensive survey of Davis et al. (2016). The latter provides a general overview of much of the research related to sparse direct methods and includes pointers to many specialized references.

There are a large number of journal articles that provide a fuller understanding of the theory and algorithms employed in symbolic factorizations. Schreiber (1982) defines the elimination tree of a sparse symmetric matrix. The seminal paper of Liu (1986) describes elimination tree construction, while for an extensive overview of the roles of elimination trees and topological orderings as well as the determination of the column sparsity patterns of the factor *L*, we refer to Liu (1990). If only row and column counts of *L* are needed, the fastest known algorithms are described in Gilbert et al. (1994). This paper also refers to another admirable paper of Liu et al. (1993) that describes the efficient computation of fundamental supernodes based on the leaf vertices of row subtrees of the elimination tree.

A key driver behind research into efficient (in terms of time and memory) sparse Cholesky algorithms has always been the development of computational codes. Many currently available packages implement not only sparse Cholesky factorizations but also more general LDLT factorizations of sparse symmetric indefinite matrices. The software is necessarily highly sophisticated and is therefore generally accompanied by technical reports and/or journal publications that explain the data structures and choices that were made in the algorithm and software design as well as providing details of the different options that are offered (examples include Duff (2004), Reid & Scott (2009), Hogg et al. (2010)).

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Chapter 5 Sparse Cholesky Solver: The Factorization Phase**

*The adoption of Cholesky's method owes not a little to the publicity given to it shortly after the end of World War II by British mathematicians and computer pioneers, including Alan Turing, Leslie Fox, Jim Wilkinson, and especially John Todd – Benzi (2017).*

*Achieving high performance for sparse direct solvers in general, and sparse Cholesky factorization, in particular, is a very well researched topic – Rennich et al. (2016)*

Having considered the symbolic phase of a sparse Cholesky solver in the previous chapter, the focus of this chapter is the subsequent numerical factorization phase. If *A* is a symmetric positive definite (SPD) matrix, then it is factorizable (strongly regular) and (in exact arithmetic) its Cholesky factorization *A* = *LLT* exists. LDLT factorizations of general symmetric indefinite matrices are considered in Chapter 7.

### **5.1 Dense Cholesky Factorizations**

Because efficient implementations of sparse Cholesky factorizations rely heavily on exploiting dense blocks, we first consider algorithms for the Cholesky factorization of dense matrices that can be applied to such blocks. Algorithm 5.1 is a basic leftlooking algorithm. It is an in-place algorithm because *L* can overwrite the lower triangular part of *A* (thus reducing memory requirements if *A* is no longer required).

Writing *A* in the block form (1.2), the computation can be reorganized to give Algorithm 5.2. This allows the exploitation of Level 3 BLAS for the computationally intensive components (dense matrix-matrix multiplies and dense triangular solves). Here *A* has *nb* block columns, which are referred to as **panels**. Step 6 can be performed using Algorithm 5.1.

Algorithms 5.1 and 5.2 are left-looking. This means that the updates are not applied immediately. Instead, all updates from previous (block) columns are applied together to the current (block) column before it is factorized. In a right-looking

#### **ALGORITHM 5.1 In-place dense left-looking Cholesky factorization Input:** Dense SPD matrix *A*.

**Output:** Factor *L* such that *A* = *LLT* .

1: **for** *j* = 1 : *n* **do** 2: *Lj* :*n,j* = *Aj* :*n,j* Only the lower triangular part of *A* is required 3: **for** *k* = 1 : *j* − 1 **do** 4: *Lj* :*n,j* = *Lj* :*n,j* − *Lj* :*n,k ljk* Update column *j* using previous columns 5: **end for** 6: *ljj* = *(ljj )*1*/*2 Overwrite the diagonal entry with its square root 7: *Lj*+1:*n,j* = *Lj*+1:*n,j / ljj* Scale off-diagonal entries in column *j* 8: **end for**

#### **ALGORITHM 5.2 In-place dense left-looking panel Cholesky factorization**

**Input:** Dense SPD matrix *A* in the form (1.2) with *nb* panels. **Output:** Factor *L* such that *A* = *LLT* .

approach (Algorithm 5.3), outer product updates are applied to the part of the matrix that has not yet been factored as they are generated.

The large panel updates can be split into operations involving only blocks. This is shown in Algorithm 5.4 for the right-looking approach.

The panel and block descriptions of the factorization enable efficient parallelization. The three main block operations, which are called tasks, are **factorize***(j b)*, **solve***(ib, j b)*, and **update***(ib, j b, kb)*. There are the following dependencies between the tasks.

**factorize***(j b)* depends on **update***(j b, kb, j b)* for all *kb*=1*,...,jb* − 1. **solve***(ib, j b)* depends on **update***(ib, kb, j b)* for all *kb*=1*,...,jb* − 1, and **factorize***(j b)*.

**update***(ib, j b, kb)* depends on **solve***(ib, kb)*, **solve***(j b, kb)*.

**ALGORITHM 5.3 In-place dense right-looking panel Cholesky factorization Input:** Dense SPD matrix *A* in the form (1.2) with *nb* panels. **Output:** Factor *L* such that *A* = *LLT* . 1: **for** *j b* = 1 : *nb* **do**

2: *Lj b*:*nb,j b* = *Aj b*:*nb,j b* 3: **end for** 4: **for** *j b* = 1 : *nb* **do** 5: Compute in-place factorization of *Lj b,j b* Overwrite *Lj b,j b* with its Cholesky factor 6: *Lj b*+1:*nb,j b* = *Lj b*+1:*nb,j b L*−*T j b,j b* Dense triangular solve 7: **for** *kb* = *j b* + 1 : *nb* **do** 8: *Lkb*:*nb,kb* = *Lkb*:*nb,kb* − *Lkb*:*nb,j b LT kb,j b* 9: **end for** 10: **end for**

### **ALGORITHM 5.4 In-place dense right-looking block Cholesky factorization Input:** Dense SPD matrix *A* in the form (1.2) with *nb* × *nb* blocks. **Output:** Factor *L* such that *A* = *LLT* .

1: **for** *j b* = 1 : *nb* **do** 2: *Lj b*:*nb,j b* = *Aj b*:*nb,j b* 3: **end for** 4: **for** *j b* = 1 : *nb* **do** 5: Compute in-place factorization of *Lj b,j b* Task **factorize***(j b)* 6: **for** *ib* = *j b* + 1 : *nb* **do** 7: *Lib,j b* = *Lib,j b L*−*T j b,j b* Task **solve***(ib, j b)* 8: **for** *kb* = *j b* + 1 : *ib* **do** 9: *Lib,kb* = *Lib,kb* − *Lib,j b LT kb,j b* Task **update***(ib, j b, kb)* 10: **end for** 11: **end for** 12: **end for**

A dependency graph can be used to schedule the tasks. Its vertices correspond to tasks and dependencies between tasks are represented as directed edges. The result is a directed acyclic graph (DAG). A task is ready for execution if and only if all tasks with incoming edges to it have completed. DAG-driven linear algebra uses either a static or dynamic schedule based on these graphs to implement the tasks in a parallel environment. In practice, it is not necessary to explicitly compute the

⎝

⎠

task DAG: it can be constructed on-the-fly taking into account the dependencies. The task DAG allows a lot of flexibility in the order in which tasks are carried out: the left- and right-looking approaches correspond to particular restricted orderings of the tasks.

#### **5.2 Introduction to Sparse Cholesky Factorizations**

There are several classes of algorithms that implement sparse Cholesky factorizations. Their major differences relate to how they schedule the computations. This affects the use of dense kernels, the amount of memory required during the factorization as well as the potential for parallel implementations. As in the dense case, the factorization is split into tasks that involve computations on and between dense submatrices and the precedence relations among them can be captured by a task graph.

We start by extending the dense Cholesky factorizations to the sparse case in a straightforward way. In practice, it is essential for efficiency to exploit the supervariables of *A* and the supernodes of *L*. Thus, while for simplicity of the descriptions and notation, we refer to rows and columns of *A* and *L*, these typically represent block rows and block columns and, as in the above discussion of the dense block factorization algorithm, the entries of *A* and *L* are then submatrices. ⎛⎞⎛⎞

The entries of *L* satisfy the relationship ⎝⎠

$$\begin{aligned} \text{block factorization algorithm, the entries of } A \text{ and } L \text{ are then submanifolds.}\\ \text{The entries of } L \text{ satisfy the relationship} \\ L\_{j+1:n,j} &= \left( A\_{j+1:n,j} - \sum\_{k=1}^{j-1} L\_{j+1:n,k} l\_{jk} \right) / l\_{jj} \quad \text{with} \quad l\_{jj} = \left( a\_{jj} - \sum\_{k=1}^{j-1} l\_{jk}^2 \right)^{1/2}, \end{aligned}$$

and from this we deduce the following result.

**Theorem 5.1 (Liu 1990)** *The numerical values of the entries in column j>k of L depend on the numerical values in column k of L if and only if ljk* = 0*.*

The theoretical background of the previous chapter based on the elimination tree T enables the dependencies in Theorem 5.1 to be searched for efficiently. In particular, T allows the row (or column) counts of *L* to be computed and they can be used to allocate storage for *L*. It can also be used to find supernodes and the resulting (block) elimination tree can then be employed to determine the (block) column structure of *L*. In practice, it can be beneficial to split large supernodes into smaller panels to better conform to computer caches.

Algorithms 5.5 and 5.6 are simplified sparse left- and right-looking Cholesky factorization algorithms that are straightforward sparse variants of Algorithms 5.1 and 5.4, respectively (the latter with *nb* = *n*, that is, without considering blocks). Here, we assume that the sparsity pattern of *L* has already been determined in the symbolic phase and static storage formats based, for example, on compressed columns and/or rows are used.

## **ALGORITHM 5.5 Simplified sparse left-looking Cholesky factorization Input:** SPD matrix *A* and sparsity pattern S{*L*}. **Output:** Factor *L* such that *A* = *LLT* . 1: *lij* = *aij* for all *(i, j )* ∈ S{*L*} Filled entries in *L* are initialised to zero 2: **for** *j* = 1 : *n* **do** 3: **for** *k* ∈ {*k<j* | *ljk* = 0} **do** 4: **for** *i* ∈ {*i* ≥ *j* | *lik* = 0} **do** 5: *lij* = *lij* − *likljk* 6: **end for** 7: **end for** 8: *ljj* = *(ljj )*1*/*2

```
9: for i ∈ {i>j | lij
= 0} do
10: lij = lij / ljj
11: end for
```
12: **end for**

#### **ALGORITHM 5.6 Simplified sparse right-looking Cholesky factorization**

**Input:** SPD matrix *A* and sparsity pattern S{*L*}. **Output:** Factor *L* such that *A* = *LLT* .

1: *lij* = *aij* for all *(i, j )* ∈ S{*L*} Filled entries in *L* are initialised to zero 2: **for** *j* = 1 : *n* **do** 3: *ljj* = *(ljj )*1*/*2 4: **for** *i* ∈ {*i>j* | *lij* = 0} **do** 5: *lij* = *lij / ljj* 6: **end for** 7: **for** *k* ∈ {*k>j* | *lkj* = 0} **do** 8: **for** *i* ∈ {*i* ≥ *k* | *lij* = 0} **do** 9: *lik* = *lik* − *lij lkj* 10: **end for** 11: **end for** 12: **end for**

An alternative for sparse matrices held in row-wise format is to compute *L* one row at a time. This is sometimes called an **up-looking** factorization because rows 1 to *i* − 1 are employed to compute row *i* (*i >* 1). The approach is asymptotically optimal in the work performed and for highly sparse matrices it is potentially extremely efficient because the entries of *A* are used in the natural order in which they are stored. However, it is difficult to incorporate high level BLAS.

The following relation holds for the *i*-th row of *L*

$$L\_{l,1:l-1}^T = L\_{1:l-1,1:l-1}^{-1} A\_{1:l-1,l} \quad \text{with} \quad l\_{ll}^2 = a\_{ll} - L\_{l,1:l-1} L\_{l,1:l-1}^T.$$

The application of *L*−1 1:*i*−1*,*1:*i*−1 can be implemented by solving the triangular system

$$L\_{1:l-1,1:l-1} \mathbf{y} = A\_{1:l-1,l} \mathbf{y}$$

and setting *LT i,*1:*i*−1 = *y*. The following result can be used to determine the sparsity pattern of *y*.

**Theorem 5.2 (Gilbert 1994)** *Consider a sparse lower triangular matrix L and the DAG* G*(LT ) with vertex set* {1*,* 2*,...,n*} *and edge set* {*(j* −→ *i)* | *lij* = 0}*. The sparsity pattern* S{*y*} *of the solution y of the system Ly* = *b is the set of all vertices reachable in* G*(LT ) from* S{*b*}*.*

*Proof* From Algorithm 3.4 and assuming the non-cancellation assumption, we see that (a) if *bi* = 0, then *yi* = 0 and (b) if for some *j<i*, *yj* = 0 and *lij* = 0, then *yi* = 0. These two conditions can be expressed as a graph transversal problem in G*(LT )*. (a) adds all vertices in S{*b*} to the set of visited vertices and (b) states that if vertex *j* has been visited, then all its neighbours in G*(LT )* are added to the set of visited vertices. Thus S{*y*} = R*each(*S{*b*}*)* ∪ S{*b*}.

Figure 5.1 illustrates the sparsity patterns of a lower triangular matrix *L* and vector *b* together with G*(LT )*. The vertices that are reachable from S{*b*}={2*,* 4} are 5 and 6 and thus S{*y*}={2*,* 4*,* 5*,* 6}.

Algorithm 5.7 outlines a sparse row Cholesky factorization that is based on the repeated solution of triangular linear systems. Theorem 5.2 can be used to determine the sparsity pattern of row *i* at Step 3, that is, by finding all the vertices that are reachable in G*(LT* 1:*j*−1*,*1:*j*−1*)* from the set {*i* | *aij* = 0*,i<j* }. A depth-first search

**Figure 5.1** An example to illustrate *L*, *b* and G*(LT )*.

of G*(LT* 1:*j*−1*,*1:*j*−1*)* determines the vertices in the row sparsity patterns in topological order, and performing the numerical solves in that order correctly preserves the numerical dependencies. Alternatively, because nonzeros of *Li,*1:*i*−1 correspond to the vertices in the *i*-th row subtree T*r(i)* that are not equal to *i*, another option is to find the row subtrees using T *(A)*.

#### **5.3 Supernodal Sparse Cholesky Factorizations**

The simplified schemes form the basis of sophisticated supernodal algorithms that are designed to be efficient in parallel computational environments. Consider the right-looking variant and recall that a supernode consists of one or more consecutive columns of *L* with the same sparsity pattern. These nonzeros are stored as a dense trapezoidal matrix (only the lower triangular part of the block on the diagonal needs to be stored and the rows of zeros in the columns of the supernode are not held). This is termed a **nodal matrix** (see Figure 5.2).

Once a supernode is ready to be factorized, a dense Cholesky factorization of the block on the diagonal of the nodal matrix is performed (one of the approaches of Section 5.1 can be used). Then a triangular solve is performed with the computed factor and the rectangular part of the nodal matrix. The next step is to iterate over ancestors of the supernode in the assembly tree. For each parent, the rows of the current supernode corresponding to the parent's columns are identified, and then the outer product of those rows and the part of the supernode below those columns formed (update operations). The resulting matrix can be held in a temporary buffer. The rows and columns of this buffer are matched against indices of the ancestors and are added to them in a sparse scatter operation. For efficiency, the updates may use panels so that the temporary buffer remains in cache.

**Figure 5.2** An illustration of a supernode (left), the corresponding nodal matrix (centre), and the nodal matrix with two panels (right). The shaded lower triangular part of the block on the diagonal and the shaded block rows are treated as dense.

#### *5.3.1 DAG-Based Approach*

The DAG-based approach can also be extended to the sparse case. Each nodal matrix is subdivided into blocks. The factorization is split into tasks in which a single block is revised. The key difference compared to the dense case is that it is necessary to distinguish between two types of update operations: **update\_internal** performs the update between blocks in the same nodal matrix and **update\_between** performs the update when the blocks belong to different nodal matrices. Thus the sparse Cholesky factorization is split into the following tasks; the first two are illustrated in Figure 5.3. In this example, the nodal matrix has two block columns that do not contain the same number of columns.

Again, the tasks are partially ordered and a task DAG is used to capture the dependencies. For example, the updating of a block of a nodal matrix from a block

**Figure 5.3** An illustration of a blocked nodal matrix with two block columns. The first block on the diagonal is triangular and the second one is trapezoidal. The task **factorize\_block** is illustrated on the left and in the centre; the task **solve\_block** is illustrated on the right.

column of *L* that is associated with a descendant of the supernode has to wait until all the relevant rows of the block column are available. At each stage of the factorization, tasks will be executing (in parallel) while others are held (in a stack or pool of tasks) ready for execution.

#### **5.4 Multifrontal Method**

The **multifrontal** method is an alternative way to compute a sparse Cholesky factorization. To discuss this popular approach, we use the following result that determines which rows and columns influence particular Schur complements using the terminology of the elimination tree.

**Theorem 5.3 (Liu 1990)** *Let A be SPD and let* T *be its elimination tree. The numerical values of entries in column k of the Cholesky factor L of A only affect the numerical values of entries in column i of L for i* ∈ *anc*T {*k*} *(*1 ≤ *k<i* ≤ *n* − 1*).* ⎛⎞

*Proof* From (4.1), setting *S(*1*)* = *A*, for *k* ≥ 2 the *(n* − *k* + 1*)* × *(n* − *k* + 1*)* Schur complement *S(k)* can be expressed as ⎜⎝⎟⎠

$$\begin{aligned} \text{for } & \text{from } (\star, 1), \text{ secure } \mathcal{S}^{(k)} = \mathcal{A}, \text{ for } \kappa \ge 2 \text{ and } (\mu - \kappa + 1) \times (\mu - \kappa + 1) \text{ such} \\\\ \text{emplement } \mathcal{S}^{(k)} \text{ can be expressed as} \\\\ \mathcal{S}^{(k)} = S\_{k:n, k:n}^{(k-1)} - \begin{pmatrix} l\_{k,k-1} \\ \vdots \\ l\_{n,k-1} \end{pmatrix} \begin{pmatrix} l\_{k,k-1} \dots \ l\_{n,k-1} \end{pmatrix} = S\_{k:n, k:n}^{(k-1)} - L\_{k:n, k-1} L\_{k:n, k-1}^T. \end{aligned} \tag{5.1}$$

Theorem 4.2 implies that all nonzero off-diagonal entries *lik* in column *k* of *L* explicitly used in the update (5.1) are such that *i* ∈ *anc*T {*k*}. Considering the Cholesky factorization as a sequence of Schur complement updates, only columns *i* with *i* ∈ *anc*T {*k*} can be influenced numerically by the Schur complement update in the *k*-th step of the factorization, and the result follows.

The computation of subsequent Schur complements by adding individual updates as in (5.1) is straightforward; the multifrontal method employs further modifications and enhancements of this basic concept. First, because the vertices of T are topologically ordered, the order in which the updates are applied progresses up the tree from the leaf vertices to the root vertex. This allows the computation of *S(k)* to be rewritten as *S(k)* = *Ak*:*n,k*:*n* −

$$S^{(k)} = A\_{k:n,k:n} - \sum\_{j \in \mathcal{T}(k) \backslash \{k\}} L\_{k:n,j} L\_{k:n,j}^T,$$

emphasizing the role of T . In place of Schur complements, the multifrontal method uses frontal matrices connected to subtrees of T . Assume *k, k*1*,...,kr* are the row indices of the nonzeros in column *k* of *L*. The **frontal matrix** *Fk* of the *k*-th subtree T *(k)* of T is the dense *(r* + 1*)* × *(r* + 1*)* matrix defined by ⎛⎜⎜⎜⎞⎟⎟⎟⎛⎜⎜⎜⎞⎟⎟⎟

$$\begin{aligned} \text{(i) of } \mathcal{T} \text{ is the dense } (r+1) \times (r+1) \text{ matrix defined by} \\ F\_k = \begin{pmatrix} a\_{kk} & a\_{kk\_1} \dots a\_{kk\_r} \\ a\_{k\_1k} & 0 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ a\_{k\_rk} & 0 & \dots & 0 \end{pmatrix} - \sum\_{j \in \mathcal{T}(k) \backslash \{k\}} \begin{pmatrix} l\_{kj} \\ l\_{k\_1j} \\ \vdots \\ l\_{k\_rf} \end{pmatrix} \begin{pmatrix} l\_{kj} \ l\_{k\_1j} \dots \ l\_{k\_rf} \end{pmatrix} .\end{aligned} \tag{5.2}$$

One step of the Cholesky factorization of *Fk* can be written as ⎜⎜⎜⎟⎟⎟⎝⎠⎜⎜⎜

$$F\_k = \begin{pmatrix} l\_{kk} & 0 & \dots & 0 \\ l\_{k1k} & & & \\ \vdots & & I & \\ l\_{k,k} & & & \end{pmatrix} \begin{pmatrix} 1 & & \\ & & \\ & & V\_k \end{pmatrix} \begin{pmatrix} l\_{kk} \ l\_{k,k} \ \dots \ l\_{k,k} \\ 0 \\ \vdots & & I \\ 0 \end{pmatrix} \tag{5.3}$$

⎟

⎟

⎟

$$\mathbf{u} = \begin{pmatrix} l\_{kk} \\ l\_{k\_1k} \\ \vdots \\ l\_{k\_rk} \end{pmatrix} \begin{pmatrix} l\_{kk} \ l\_{k\_1k} \ \dots \ l\_{k\_rk} \end{pmatrix} + \begin{pmatrix} 0 \\ & \\ & V\_k \end{pmatrix},\tag{5.4}$$

where *Vk* is termed a **generated element** (it is also sometimes called an **update matrix** or a **contribution block**). The name "generated element" is because the multifrontal method has its origins in the simpler **frontal method**, which uses a single frontal matrix. The frontal method was originally proposed for problems arising in finite element problems to avoid the need to explicitly construct the system matrix *A*; it was later generalized to non-element problems. It works with a single frontal matrix and has less scope for parallelisation compared to the multifrontal method; it is no longer widely used.

Equating the last *r* rows and columns in (5.2) and (5.4) yields ⎜⎝⎟⎠

⎝

⎛

⎞

$$\text{rows and columns in } (5.2) \text{ and } (5.4) \text{ yields}$$

$$V\_k = -\sum\_{j \in \mathcal{T}(k)} \begin{pmatrix} l\_{k\_1 j} \\ \vdots \\ l\_{k\_r j} \end{pmatrix} \begin{pmatrix} l\_{k\_1 j} & \dots & l\_{k\_r j} \end{pmatrix} . \tag{5.5}$$

Assume that *cj* (*j* = 1*,...,s*) are the children of *k* in T . The set T *(k)* \ {*k*} is the union of disjoint sets of vertices in the subtrees T *(cj )*. Each of these subtrees is represented in the overall update by the generated element (5.5). Thus, *Fk* can be written in an recursive form using the generated elements of the children of *k* as follows ⎛⎜⎜⎜⎞⎟⎟⎟

⎠

$$F\_k = \begin{pmatrix} a\_{kk} \ a\_{kk\_1} \dots a\_{kk\_r} \\ a\_{k1k} \ 0 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ a\_{k,k} \ 0 & \dots & 0 \end{pmatrix} \oplus V\_{c\_1} \oplus \dots \oplus V\_{c\_r}.\tag{5.6}$$

Here, the operation ←→ denotes the addition of matrices that have row and column indices belonging to subsets of the same set of indices (in this case, *k, k*1*,...,kr*); entries that have the same row and column indices are summed. This is referred to as the **extend-add operator**.

Adding a row and column of *A* and the generated elements into a frontal matrix is called the **assembly**. A variable is **fully summed** if it is not involved in any rows and columns of *A* that have still to be assembled or in a generated element. Once a variable is fully summed, it can be eliminated. A key feature of the multifrontal method is that the frontal matrices and the generated elements are compressed and stored without zero rows and columns as small dense matrices. Integer arrays are used to maintain a mapping of the local contiguous indices of the frontal matrices to the global indices of *A* and its factors. Symmetry allows only the lower triangular part of these matrices to be held. Algorithm 5.8 outlines the basic multifrontal method.

### **ALGORITHM 5.8 Basic multifrontal Cholesky factorization Input:** SPD matrix *A* and its elimination tree. **Output:** Factor *L* such that *A* = *LLT* .

**ALGORITHM 5.9 Multifrontal Cholesky factorization using the assembly tree Input:** SPD matrix *A* and its assembly tree. **Output:** Factor *L* such that *A* = *LLT* .

We have the following observation.

**Observation 5.1** *Each generated element Vk is used only once to contribute to a frontal matrix Fparent (k). Furthermore, the index list for the frontal matrix Fk is the set of row indices of the nonzeros in column k of the Cholesky factor L.*

In practical implementations, efficiency is improved by using the assembly tree (Section 4.6) because it allows more than one elimination to be performed at once. This is outlined in Algorithm 5.9. Here *kb* is used to index the frontal matrix on the *kb*-th step (1 ≤ *kb* ≤ *nsup*).

As an example, consider the matrix and its assembly tree given in Figure 4.10. The *nsup* = 5 supernodes are {1*,* 2}*,* 3*,* 4*,* 5*,*{6*,* 7*,* 8*,* 9} and so variables 1 and 2 can be eliminated together on the first step. Assembling rows/columns 1 and 2 of the original matrix, the frontal matrix *F*1 and generated element *V*1 have the structure ⎛⎞

⎟

$$\begin{array}{ccccccccc} & & & & & \\ & & & & & \\ & & & 1 & 2 & 8 & 9 \\ & & & & & 1 & 2 & 8 \\ F\_1 = & 2 & \begin{pmatrix} \* & \* & \* & \* & \* \\ \* & \* & \* & \* & \* \\ \* & \* & \* & & \\ \* & \* & \* & & \end{pmatrix}, & V\_1 = & \begin{matrix} 8 & 9 \\ 5 & -1 \end{matrix} \\ \end{array}$$

where *f* denotes fill-in entries (only the lower triangular entries are stored in practice). Similarly, ⎛⎞⎠

$$\begin{array}{ccccc} 3 & 4 & 8\\ 3 & \left( \* & \* & \* \\ F\_{2} = 4 & \begin{pmatrix} \* & \* & \* \\ \* & \* & \* \\ \* & \* & \* \end{pmatrix} \right) , \quad V\_{2} = \frac{4}{8} \begin{pmatrix} \* & \* \\ \* & \* \end{pmatrix} . \end{array}$$

The frontal matrix *F*3 and generated element *V*3 are given by

⎝

⎜

#### 5.5 Parallelism Within Sparse Cholesky Factorizations 85 ⎛⎞

⎝

⎛

⎝

⎜

*F*3 = 478 4 ∗∗∗ 7 ∗ ∗ 8 ∗ ∗ ←→ *V*2*, V*3 = 7 8 7 ∗ *f* 8 *f* ∗ *.*

⎞

⎠

⎠

Then

$$F\_4 = \begin{array}{cccc} & \mathbf{s} & \mathbf{7} & \mathbf{8} \\ & \mathbf{3} & \begin{pmatrix} \* & \* & \* \\ \* & \* & \* \\ \* & \* & \* \end{pmatrix}, & V\_4 = \begin{array}{ccc} & \mathbf{7} & \mathbf{s} \\ & \mathbf{8} & \begin{pmatrix} \* & f \\ f & \* \end{pmatrix} . \end{array}$$

⎞

⎟

and, finally, with *kb* = 5 we have ⎛

$$F\_5 = \begin{pmatrix} 6 & 7 & 8 & 9 \\ \* & \* & \* & \* \\ 8 & \begin{pmatrix} \* & \* & \* \\ \* & & \* \\ \* & \* & \* \end{pmatrix} & \bigoplus {\*} V\_4 \xleftrightarrow{\bigoplus} V\_3 \xleftrightarrow{\bigoplus} V\_1... $$

An important implementation detail is how and where to store the generated elements. The partial factorization of *Fkb* at supernode *kb* can be performed once the partial factorizations at all the vertices belonging to the subtree of the assembly tree with root vertex *kb* are complete. If the vertices of the assembly tree are ordered using a depth-first search, the generated elements required at each stage are the most recently computed ones amongst those that have not yet been assembled. This makes it convenient to use a stack. This affects the order in which the variables are eliminated but in exact arithmetic, the results are identical.

Nevertheless, the memory demands of the multifrontal method can be very large. Not only is it dependent on the initial ordering of *A* but the ordering of the children of a vertex in the assembly tree can significantly affect the required stack size. Some implementations target limiting stack storage requirements. An attractive feature of the multifrontal method is that the generated elements can be held using auxiliary storage (in files on disk) to restrict the in-core memory requirements, allowing larger problems to be solved than would otherwise be possible.

#### **5.5 Parallelism Within Sparse Cholesky Factorizations**

Sparse Cholesky factorizations use supernodes and task graphs (the assembly tree for the multifrontal method) to control the computation. The number of rows and columns in a supernode typically increases away from the leaf vertices and towards the root of the task graph because a supernode accumulates fill-in from its ancestors in the task graph. As a result, tasks that are relatively close to the root tend to have more work associated with them. On the other hand, the width of the task graph shrinks close to the root. In other words, a typical task graph for sparse matrix factorization tends to have a large number of small independent tasks close to the leaf vertices, but a small number of large tasks close to the root. An ideal parallelization strategy that would match the characteristics of the problem is as follows. Initially, assign the relatively plentiful independent tasks at or near the leaf vertices to parallel threads or processes. This is called **task** or **tree level** parallelism; it is influenced by the ordering of *A*. As tasks complete, other tasks become available and are scheduled similarly. This continues while there are enough independent tasks to keep all the threads or processes busy. When the number of available parallel tasks becomes too small, the only way to keep the latter busy is to assign more than one to a task. This is termed **node level** parallelism. The number of threads or processes working on individual tasks should increase as the number of parallel tasks decreases. Eventually, all threads or processes are available to work on the root task. The computation corresponding to the root task is equivalent to factoring a dense matrix of the size of the root supernode.

The multifrontal method is often the formulation of choice for highly parallel implementations of sparse matrix factorizations. This is because of its natural data locality (most of the work of the factorization is performed in the dense frontal matrices) and the ease of synchronization that it permits. In general, each supernode is updated by multiple other supernodes and it can potentially update many other supernodes during the course of the factorization. If implemented naively, all these updates may require excessive locking and synchronization in a shared-memory environment or generate excessive message-traffic in a distributed environment. In the multifrontal method, the updates are accumulated and channelled along the paths from the leaf vertices of the assembly tree to its root vertex. This gives a manageable structure to the potentially haphazard interaction among the tasks.

In Section 1.2.4, bit compatibility was discussed. While different orderings of the children of a vertex in the assembly tree do not affect the total number of floatingpoint operations that are performed in the multifrontal method, in finite-precision arithmetic changing the order of the assemblies into the frontal matrices can lead to slightly different results. Given that the number of children is typically small and that large matrices can be partitioned such that summations can be safely performed in parallel, the overhead in the multifrontal method of enforcing a defined order of the summation is relatively small. By contrast, in the supernodal approach, for each data block a number of matrices equal to the block dependencies are summed. Given the relatively large numbers (several thousand) for many nodes, an enforced order may be detrimental to efficiency.

#### **5.6 Notes and References**

Exploiting panels and blocks in both left- and right-looking Cholesky factorization algorithms is extremely important. The development of sparse supernodal factorizations for uniprocessors and multiprocessors in the 1990s is discussed by Ng & Peyton (1993a,b); Rothberg & Gupta (1993) presents an early comparison of various types of block Cholesky factorizations. PaStiX of Hénon et al. (2002) is a parallel left-looking supernodal solver that is primarily designed for positive definite systems. Rotkin & Toledo (2004) introduce a hybrid left-looking/right-looking algorithm and Rozin & Toledo (2005) show that no sparse numerical factorization is uniformly better than the others. An up-looking approach, which is fast in practice for very sparse matrices, is employed in the widely used CHOLMOD solver of Chen et al. (2008). The package HSL\_MA87 implements a sparse DAG-based Cholesky factorization for shared-memory architectures; further details of the approach can be found in Hogg et al. (2010).

The multifrontal algorithm has its origins in the simpler frontal method of Irons (1970), which was developed by the civil engineering community from the 1960s onwards to solve the linear systems that arise within finite element methods. At a time when the main memory of even the most powerful computers was extremely limited, the frontal method was heavily influenced by the need to minimize the memory requirements of the linear solver. It was initially designed for SPD banded linear systems and was subsequently extended to nonsymmetric problems by Hood (1976) and to the symmetric indefinite case by Reid (1981); Duff (1984) generalizes the approach to non-element problems. The frontal method proceeds by alternating the assembly of the finite elements into a single dense frontal matrix with the elimination and update of variables. Once variables have been eliminated they are no longer needed during the factorization and so they are removed from the frontal matrix and stored elsewhere (for example, not in main memory but on an external disk) until needed during the solve phase. This frees up space to accommodate the next element to be assembled. Because the frontal method does not use the assembly tree, the frontal matrix can be much larger than those in the multifrontal method, leading to higher operation counts but also allowing the use of BLAS with larger block sizes. Efficient implementations were developed up until the late 1990s. For example, by Duff & Scott (1996, 1999), who provide a package MA62 for SPD problems in element form that employs a single array of length *n*, exploits Level 3 BLAS, and holds the computed factors on disk; a coarse-grained parallel version is also available, see Duff & Scott (1994) and Scott (2001).

The frontal method and the work of Speelpenning (1978) on the so-called generalized element method led to the development by Duff & Reid (1983) of the multifrontal method for solving general symmetric systems (including systems in element form). A detailed matrix-based explanation is given in Liu (1992). The method is implemented in some of the most important sparse direct solvers. The MUMPS (2022) package, which has been actively developed over many years, provides a state-of-the-art distributed memory general-purpose multifrontal solver that uses shared-memory parallelism within each MPI process. Other important parallel multifrontal solvers are HSL\_MA97 (Hogg & Scott, 2013b) and WSMP (2020), while the serial package MA57 of Duff (2004) (which superseded the original and perhaps most well-known multifrontal solver MA27 of Duff & Reid, (1983)) remains very popular. An attractive feature of HSL\_MA97 is that it computes bit-compatible solutions. HSL\_MA77 of Reid & Scott (2009) is designed to minimize memory requirements by allowing the factors and the multifrontal stack to be efficiently held outside of main memory (an option that is also offered by MUMPS). In common with earlier frontal solvers, HSL\_MA77 allows the user to input the system matrix in element form (that is, *A* is not explicitly assembled for problems coming from finite element applications but is input one element at a time).

The use of GPUs is well-suited to a multifrontal or supernodal factorization because these approaches rely on regular block computations within dense submatrices. Implementing the multifrontal method (including for symmetric indefinite matrices) on GPU architectures is discussed in Hogg et al. (2016), while Lacoste et al. (2012) and Rennich et al. (2016) present GPU-accelerated supernodal factorizations. Discussion of the use of GPUs within direct solvers is included in the comprehensive survey of Davis et al. (2016).

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Chapter 6 Sparse LU Factorizations**

*The closer one looks, the more subtle and remarkable Gaussian elimination appears – Trefethen (1985)*

*Gaussian elimination is living mathematics. It has mutated successfully for the last two hundred years to meet changing social needs – Grcar (2011)*

This chapter considers the LU factorization of a general nonsymmetric nonsingular sparse matrix *A*. In practice, numerical pivoting for stability and/or ordering of *A* to limit fill-in in the factors is often needed and the computed factorization is then of a permuted matrix *P AQ*. Pivoting is discussed in Chapter 7 and ordering algorithms in Chapter 8.

## **6.1 Sparse LU Factorizations and Their Graph Models**

In Chapter 4, graphs were used to describe structural changes during a sparse Cholesky factorization. In particular, the elimination tree was shown to play a key role and, in the previous chapter, the use of DAGs was discussed. For general matrices, there are a number of ways that graphs can be employed.

## *6.1.1 Use of Elimination DAGs*

The first graph model uses the elimination DAGs associated with *L* and *U* that were defined in (2.1)–(2.2). The following observation, which is illustrated in Figure 6.1, generalizes Observation 4.1 to nonsymmetric matrices.

**Observation 6.1** *If i>j and uj i* = 0*, then the* **column replication principle** *states*

**Figure 6.1** An illustration of the column and row replication principles of sparse LU factorizations. The matrix *A* is on the left. In the centre, we show in red the filled entries in *L* resulting from the replication of the first column in the second column because *u*12 = 0. On the right, we show in blue the filled entries in *U* resulting from the replication of the second row in the third row because *l*32 = 0. Other filled entries resulting from subsequent steps of the factorization are denoted in black.

$$
\mathcal{S}\{L\_{i:n,j}\} \subseteq \mathcal{S}\{L\_{i:n,i}\},
$$

*that is, the pattern of column j of L (rows i to n) is replicated in the pattern of column i of L. Analogously, if i>j and lij* = 0*, then the* **row replication principle** *states*

$$
\mathcal{S}\{U\_{j,i:n}\} \subseteq \mathcal{S}\{U\_{i,i:n}\},
$$

*that is, the pattern of row j of U (columns i to n) is replicated in the pattern of row i of U.*

Algorithm 6.1 outlines a basic sparse LU factorization. Here it is assumed that *A* is factorizable so that pivoting is not needed. The remainder of this chapter looks at techniques that can be used to develop the approach into an efficient one.

The following theorem formulates the recursive column replication and the replication of nonzeros along rows of *L* using directed paths in G*(U )*; an analogous result holds for the rows of *U* and directed paths in G*(LT )*.

**Theorem 6.1 (Gilbert & Liu 1993)** *Assume that for some k<j there is a directed path k* G*(U )* ⇒ *j . Then*

$$\mathcal{S}\{L\_{j:n,k}\} \subseteq \mathcal{S}\{L\_{j:n,j}\}.\tag{6.1}$$

*Moreover, if lik* = 0 *for some i>j , then lis* = 0 *for all vertices s on this path.*

The next two theorems generalize Theorem 4.3 to *A* being a general nonsymmetric matrix.

**Theorem 6.2 (Gilbert & Liu 1993)** *If aij* = 0 *and i>j , then there is a filled entry lij* = 0 *if and only if there exists k<j such that aik* = 0 *and there is a directed path k* G*(U )* ⇒ *j .*

#### **ALGORITHM 6.1 Basic sparse LU factorization**

**Input:** Nonsymmetric and factorizable matrix *A* = *LA* + *DA* + *UA*. **Output:** LU factorization *A* = *LU*.

1: *L* = *I* + *LA* Identity plus strictly lower triangular part of *A* 2: *U* = *DA* + *UA* Diagonal plus strictly upper triangular part of *A* 3: **for** *k* = 1 : *n* − 1 **do** 4: **for** *i* ∈ {*i>k* | *lik* = 0} **do** 5: *lik* = *lik/ukk* 6: *Ui,i*:*n* = *Ui,i*:*n* − *Uk,i*:*nlik* Update row *i* of *U* 7: **end for** 8: **for** *j* ∈ {*j>k* | *ukj* = 0} **do** 9: *Lj*+1:*n,j* = *Lj*+1:*n,j* − *Lj*+1:*n,kukj* Update column *j* of *L* 10: **end for** 11: **end for**

**Theorem 6.3 (Gilbert & Liu 1993)** *If aij* = 0 *and i<j , then there is a filled entry uij* = 0 *if and only if there exists k<i such that akj* = 0 *and there is a directed path k* G*(LT )* ⇒ *i.*

Theorems 6.2 and 6.3 are demonstrated in Figure 6.2. Consider the directed path 1 → 3 → 5 → 6 in G*(U )*. Existence of this path implies the fill-in in *L*, first in

**Figure 6.2** The sparsity patterns of *A* (left) and *L*+*U* (right) together with the graphs G*(A)* (left), G*(LT )* (centre) and G*(U )* (right). The filled entries are denoted by *f* and the corresponding edges are the red dashed lines.

**Figure 6.3** Example to show the transitive reduction of a DAG. G is on the left, its transitive reduction G0 is in the centre, and one possible G that is equireachable with G is on the right.

column 3, then in columns 5 and 6. Similarly, the directed path 2 → 4 → 5 → 6 in G*(LT )* implies fill-in at positions *(*4*,* 7*)*, *(*5*,* 7*)* and *(*6*,* 7*)* in *U*.

### *6.1.2 Transitive Reduction and Equireachability*

To employ G*(LT )* and G*(U )* in efficient algorithms, they need to be simplified. One possibility is to use transitive reductions that are sparser and preserve reachability within the graphs. A subgraph G0 = *(*V*,* E0*)* is a **transitive reduction** of G = *(*V*,* E*)* if the following conditions hold:

A transitive reduction is unique for a DAG, as shown in the following theorem and illustrated in Figure 6.3.

**Theorem 6.4 (Aho et al. 1972)** *Let* G *be a DAG. The transitive reduction* G0 *of* G *is unique and is the subgraph that has an edge for every path in* G *and has no proper subgraph with this property.*

If S{*A*} is symmetric, then, as illustrated in Figure 6.4, the role of the transitive reduction is played by the elimination tree.

**Theorem 6.5 (Liu 1990; Eisenstat & Liu 2005a)** *If A is symmetrically structured, then the transitive reduction of the DAG* G*(LT ) (=* G*(U )) is the elimination tree* T *(A).*

Obtaining the exact transitive reduction of a DAG can be expensive. Instead, approximate reductions that drop the minimality condition may be computed. A directed graph G with the same vertex set as G that satisfies condition *(T* 1*)* is said

**Figure 6.4** The sparsity patterns of *L*+*U* of a symmetrically structured *A* together with the DAG G*(LT )* (left) and the elimination tree T *(A)* (right). The filled entries are denoted by *f* and the corresponding edges are the red dashed lines. It is straightforward to see that T *(A)* is obtained as the transitive reduction of G*(LT )*.

to be **equireachable** with G. The next result is a simplification of Theorem 6.1; an analogous result holds for the sparsity patterns of the rows of *U*.

**Theorem 6.6 (Gilbert & Liu 1993)** *Assume* G *is equireachable with* G*(U ) and for some k<j there is a directed path k* G ⇒ *j . Then (6.1) holds. Moreover, if lik* = 0 *for some i>j , then lis* = 0 *for all vertices s on the directed path.*

Equireachability enables sparse triangular linear systems to be solved more efficiently. In Chapter 5, Theorem 5.2 describes how to obtain the sparsity pattern J of the solution of a lower triangular system using paths in G*(LT )*. This graph can be replaced by any graph that is equireachable with G*(LT )*. Equireachability also allows Theorems 6.2 and 6.3 to be rewritten using paths in a graph G that is equireachable with G.

**Theorem 6.7 (Gilbert & Liu 1993)** *If aij* = 0 *and i>j , then there is a filled entry lij* = 0 *if and only if there exists k<j such that aik* = 0 *and a directed path k* G *(U )* ⇒ *j , where* G *(U ) is equireachable with* G*(U ).*

**Theorem 6.8 (Gilbert & Liu 1993)** *If aij* = 0 *and i<j , then there is a filled entry uij* = 0 *if and only if there exists k<i such that akj* = 0 *and a directed path k* G *(LT )* ⇒ *i, where* G *(LT ) is equireachable with* G*(LT ).*

Figure 6.5 depicts G*(U )* and G *(U )* for the matrix in Figure 6.2.

A description of the sparsity patterns of the columns of *L* can be obtained from the Schur complement (3.2) as follows:

**Figure 6.5** The DAG G*(U )* for the matrix from Figure 6.2 (left) and G *(U )* which is equireachable with G*(U )* (right). S{*Lj* :*n,j* } = S{*Aj* :*n,j* }

$$\mathcal{S}\{L\_{j;n,j}\} = \mathcal{S}\{A\_{j;n,j}\} \bigcup\_{k$$

Theorem 6.7 implies that not all the terms in this union are needed to obtain S{*Lj* :*n,j* }. This result is given in Theorem 6.9, which shows how S{*L*} can be computed by columns if G *(U )* that is equireachable with G*(U )* is known.

**Theorem 6.9 (Gilbert & Liu 1993)** *If* G *(U ) is equireachable with* G*(U ), then*

S{*Lj* :*n,j* } = S{*Aj* :*n,j* } *(k*→*j )*∈E*(*G *(U ))* S{*Lj* :*n,k*}*,* 1 ≤ *j* ≤ *n.* (6.2)

*Proof* Consider an edge *(k* → *j )* in G*(U )* but not in G *(U )*. Repeatedly applying (6.1) along the directed path *k* G *(U )* ⇒ *j* , we see that *Lj* :*n,k* is contained in the right-hand side of (6.2) and therefore S{*Lj* :*n,j* } is contained in the right-hand side of (6.2). Because the right-hand side of (6.2) is trivially contained in the left-hand side, the result follows.

An analogous result holds for the rows of *U*.

**Theorem 6.10 (Gilbert & Liu 1993)** *If* G *(L) is equireachable with* G*(L), then*

$$\begin{aligned} \text{plus result holds for the rows of } U. \\ \text{6.10 (GIibert \& \text{Lui } 1993) \quad If \mathcal{G}'(L) \text{ is equireacable with } \mathcal{G}(L) \\ \mathcal{S}\{U\_{l,l:n}\} = \mathcal{S}\{A\_{l,l:n}\} \underbrace{\bigcup\_{(k \to i) \in \mathcal{S}(\mathcal{G}'(L^T))} \mathcal{S}\{U\_{k,l:n}\}}\_{(k \to i) \in \mathcal{S}(\mathcal{G}'(L^T))}, \quad 1 \le i \le n. \end{aligned}$$

As an example of Theorem 6.9, consider the matrix in Figure 6.2. Because *(*3 → 5*)* is the only edge of G *(U )* in the union on the right-hand side of (6.2), S{*L*5:7*,*5} is given by

$$\mathcal{S}\{L\_{\ $?7,\$ }\} = \mathcal{S}\{A\_{\ $?7,\$ }\} \cup \mathcal{S}\{L\_{\ $?7,\$ }\}.$$

We can see this from the graph G *(U )* in Figure 6.5 (top right).

#### *6.1.3 Symbolic LU Factorizations Using DAGs*

Factorization by bordering can be used to obtain S{*L*} by rows and S{*U*} by columns. Assume the sparsity patterns of the first *k* − 1 rows of *L* and the first *k* − 1 columns of *U* (1 *< k* ≤ *n*) have been computed. At step *k*, the factors satisfy *Ak,*1:*k*−1 *akk* 0 *ukk*

$$A\_{1:k,1:k} = \begin{pmatrix} A\_{1:k-1,1:k-1} \ A\_{1:k-1,k} \\ A\_{k,1:k-1} \end{pmatrix} = \begin{pmatrix} L\_{1:k-1,1:k-1} & 0 \\ L\_{k,1:k-1} & 1 \end{pmatrix} \begin{pmatrix} U\_{1:k-1,1:k-1} \ U\_{1:k-1,k} \\ 0 \end{pmatrix} . \tag{6.3}$$

Equating terms for the *(*2*,* 1*)* block, row *k* of *L* satisfies

$$L\_{k,1:k-1}U\_{1:k-1,1:k-1} = A\_{k,1:k-1},$$

or, equivalently, if *y* denotes the off-diagonal part of the column *k* of *LT* , then it is the solution of the lower triangular system

$$U\_{1:k-1,1:k-1}^T \mathcal{Y} = A\_{k,1:k-1}^T$$

From Theorem 5.2, the sparsity pattern of *y* is the set of all vertices reachable in the DAG G*(U*1:*k*−1*,*1:*k*−1*)* (or in a graph that is equireachable with it) from the nonzeros in *Ak,*1:*k*−1. Similarly, equating terms in (6.3) for the *(*1*,* 2*)* block, column *k* of *U* satisfies

$$L\_{1:k-1,1:k-1}U\_{1:k-1,k} = A\_{1:k-1,k}\cdot\bar{z}$$

Again, its sparsity pattern can be determined using Theorem 5.2 and the DAG G*(LT* 1:*k*−1*,*1:*k*−1*)*. The diagonal entry *ukk* is then computed as *akk* −*Lk,*1:*k*−1*U*1:*k*−1*,k*. This shows that determining the sparsity patterns of *L* and *U* and computing their numerical values is coupled: computation of the factors needs be mutually interleaved because computing part of one requires information from a part of the other.

#### *6.1.4 Graph Pruning*

Consider the matrices in Figure 6.6. The one in the centre is the same as the one on the left except that the entries in positions *(*4*,* 6*)* and *(*6*,* 4*)* have been removed (that is, pruned). Both matrices have the same sets of reachable vertices in G*(LT )* and G*(U )*. This suggests how to find G *(LT )* and G *(U )* that are equireachable with G*(LT )* and G*(U )*, respectively.

**Theorem 6.11 (Eisenstat & Liu 1992)** *If for some j<s both lsj* = 0 *and uj s* = 0*, then there are no edges (j* → *k) with k>s in the transitive reductions of* G*(U ) and* G*(LT ).*

**Figure 6.6** An example of symmetric pruning. On the left is S{*L*+*U*}. In the centre is the reduced sparsity pattern obtained by symmetric pruning. On the right is the reduced sparsity pattern that results from symmetric path pruning.

*Proof* Let *(j* → *k)* be an edge of G*(U )*, that is, *ujk* = 0. Because *lsj* = 0 and *ujk* = 0 implies that *usk* = 0, there is a path *j* → *s* → *k* in G*(U )* and the edge *(j* → *k)* does not belong to the transitive reduction of G*(U )*. The result for G*(LT )* can be seen analogously.

This theorem implies that if for some *s >* 1 there are edges

$$j \xrightarrow{\mathcal{G}(L^{\overline{r}})} s \quad \text{and} \quad j \xrightarrow{\mathcal{G}(U)} s,$$

then all edges *(j* → *k)* in G*(U )* and G*(LT )* with *k>s* can be pruned. The resulting DAGs G *(U )* and G *(LT )* have fewer edges and are equireachable with G*(U )* and G*(LT )*, respectively. The removal of redundant edges based on Theorem 6.11 is called **symmetric pruning**.

There are other ways to perform pruning. For example, if for some *s >* 1 there are paths

$$j \xrightarrow{\mathcal{G}(L^{T})} s \quad \text{and} \quad j \xrightarrow{\mathcal{G}(U)} s,$$

then for all *k>s* **symmetric path pruning** removes the edges *(j* → *k)* from G*(U )* and G*(LT )*. Consider again Figure 6.6. In the centre is the sparsity pattern after symmetric pruning and on the right is the reduced sparsity pattern that results from symmetric path pruning. The edge *(*1 → 6*)* is not required in G *(LT )* or G *(U )* because there are paths

$$1 \xrightarrow{\mathcal{G}(L^{\top})} 2 \xrightarrow{\mathcal{G}(L^{\top})} 4 \xrightarrow{\mathcal{G}(L^{\top})} \mathfrak{S} \xrightarrow{\mathcal{G}(L^{\top})} 6 \quad \text{and} \quad 1 \xrightarrow{\mathcal{G}(U)} \mathfrak{S} \xrightarrow{\mathcal{G}(U)} 6.$$

**Figure 6.7** An example of the sparsity pattern of a nonsymmetric matrix *A* (left), S{*L* + *U*} with filled entries denoted by *f* (right) and its elimination tree.

#### *6.1.5 Elimination Trees for Nonsymmetric Matrices*

The elimination DAGs G*(L)* and G*(U )* can be combined into a single structure called the **nonsymmetric elimination tree** in which edges are replaced by paths. This can be advantageous because it is more compact. From (4.3), if S{*A*} is symmetric, then its elimination tree is defined in terms of the mapping

$$parent(j) = \min\{i \mid i > j \text{ and } l\_{ij} \neq 0\}.$$

The condition *lij* = 0 is equivalent to *i* G*(L)* −−−→ *j* G*(LT )* −−−→ *i*. In the nonsymmetric case, the definition can be generalized using directed paths

$$parent(j) = \min\{i \mid i > j \text{ and } i \stackrel{\mathcal{G}(L)}{\right\}} \stackrel{\mathcal{G}(U)}{\longrightarrow} j \stackrel{\mathcal{G}(U)}{\longrightarrow} i\text{)}.\tag{6.4}$$

This is illustrated in Figure 6.7. Vertices 6, 8, and 10 are the only ones with cycles of the form

$$i \xrightarrow{\mathcal{G}(L)} 2 \xrightarrow{\mathcal{G}(U)} i,$$

namely,

## **ALGORITHM 6.2 Basic computation of the elimination tree for nonsymmetric** *A* **Input:** Digraph G*(A)*.

**Output:** The elimination tree given by the mapping *parent*.

1: *parent (*1 : *n)* = 0 2: **for** *i* = 1 : *n* **do** 3: Find the vertex set V*C* of the strong component of G*(A*1:*i,*1:*i)* that contains *i* 4: **for** *j* ∈ V*C* \ {*i*} **do** 5: **if** *parent (j )* = 0 **then** 6: *parent (j )* = *i* 7: **end if** 8: **end for** 9: *parent (i)* = 0

10: **end for**

## 6 G*(L)* −−−→ 2 G*(U )* −−−→ 5 G*(U )* −−−→ 6*,* 8 G*(L)* −−−→ 2 G*(U )* −−−→ 8 and 10 G*(L)* −−−→ 6 G*(L)* −−−→2 G*(U )* −−−→10*.*

In this example, *parent (*2*)* = 6.

Theorem 6.12, which can be regarded as a generalization of Corollary 4.6, shows how the elimination tree for nonsymmetric *A* can be constructed.

**Theorem 6.12 (Eisenstat & Liu 2005a)** *Let A be a nonsymmetric matrix. i* = *parent (j ) if and only if i>j and i is the smallest vertex that belongs to the same strong component of* G*(A*1:*i,*1:*i) as vertex j .*

This result is employed in Algorithm 6.2. The complexity of finding the strong components of a digraph with *m* edges and *n* vertices is *O(n* + *m)*. Hence, the complexity of Algorithm 6.2 is *O(nz(A) n)*. More sophisticated approaches with complexity *O(nz(A)*log *n)* exist.

To illustrate Algorithm 6.2, consider the matrix and its elimination tree depicted in Figure 6.7. The main loop sets the first nonzero value in the array *parent* when *i* = 3 because this is the first *i* for which the set V*C* \ {*i*} is non empty; it is equal to {1} and thus *parent (*1*)* = *i* = 3. For *i* = 4, the vertex set {1*,* 3*,* 4} forms a strong component of G*(A*1:4*,*1:4*)* and so *parent (*3*)* = 4. For *i* = 5, the single vertex {5} is a strong component of G*(A*1:5*,*1:5*)* and, therefore, 5 is not a parent of any other vertex (it is a leaf vertex). G*(A*1:6*,*1:6*)* has two strong components with vertex sets {1*,* 3*,* 4} and {2*,* 5*,* 6}. *i* = 6 belongs to the second of these and thus the algorithm sets *parent (j )* = *i* = 6 for *j* = 2 and 5.

An attractive idea for constructing S{*L* + *U*} and subsequently computing the LU factorization is based on using the **column elimination tree** T *(AT A)*.

**Theorem 6.13 (George & Ng 1985; Grigori et al. 2009)** *Assume all the diagonal entries of A are nonzero and let L L T be the Cholesky factorization of AT A. Then for any row permutation matrix P such that P A* = *LU the following holds:*

**Figure 6.8** The sparsity patterns of *A* and *L* + *U* (top) and of *AT A* and *L* + *L T* , where *AT A* = *L L T* (bottom). Filled entries are denoted by *f* . The corresponding elimination trees are also given.

$$
\mathcal{S}\{L+U\} \subseteq \mathcal{S}\{\widehat{L} + \widehat{L}^T\}.
$$

An important feature of Theorem 6.13 is that it holds for *any* row permutation matrix *P* applied to *A*. This allows partial pivoting (Section 3.1.2) to be used. The following result states that T *(AT A)* represents the potential dependencies among the columns in an LU factorization and that for strong Hall matrices no tighter prediction is possible from the sparsity structure of *A*.

**Theorem 6.14 (Gilbert & Ng 1993)** *If P A* = *LU is any factorization of A with partial pivoting, then the following hold.*

Figure 6.8 illustrates the differences in the sparsity patterns of *A* and *AT A* and of their factors; the corresponding elimination trees are also given. This reveals a potential problem with the column elimination tree: S{*AT A*} can have significantly more entries than S{*L* + *U*}. An extreme example is when *A* has one or more dense rows because *AT A* is then fully dense.

### *6.1.6 Supernodes in LU Factorizations*

Supernodes group together columns of the factors with the same nonzero structure, allowing them to be treated as a dense submatrix for storage and computation. When solving SPD systems, supernodes can be determined during the symbolic phase. For nonsymmetric matrices, supernodes are harder to characterize. The need to incorporate pivoting means it may not be possible to predict the sparsity structures of the factors before the numerical factorization and they must be identified on-thefly. While there are several possible ways to define supernodes, the simplest (which is widely used in practice) follows the symmetric case and defines a supernode to be a set of contiguously numbered columns of *L* with the triangular diagonal block treated as dense and the columns as having the same structure below the diagonal block.

In a Cholesky solver, fundamental supernodes (Section 4.6.1) are made contiguous by symmetrically permuting the matrix according to a postordering of its elimination tree; this does not change the sparsity of the Cholesky factor. For nonsymmetric *A*, before the numerical factorization, T *(AT A)* can be constructed and the columns of *A* then permuted according to its postordering to bring together supernodes. The following result extends Theorem 4.9.

**Theorem 6.15 (Li 1996)** *Let A have column elimination tree* T *(AT A). Let p be a permutation vector such that if pi is an ancestor of pj in* T *(AT A), then i>j . Let P be the permutation matrix corresponding to p and let A* = *P APT . Then* T *(A T A) is isomorphic to* T *(AT A); in particular, relabelling each vertex i of* T *(A T A) as pi yields* T *(AT A). If, in addition, A* = *L U is an LU factorization without pivoting then PT LP and PT UP are lower triangular and upper triangular matrices, respectively, so that A* = *(PT LP )(P T UP) is also an LU factorization.*

In practice, for many matrices the average size of a supernode is only 2 or 3 columns and many comprise a single column. Larger artificial supernodes may be created by merging vertex *j* with its parent vertex *i* in T *(AT A)* if the subtree rooted at *i* has fewer than some chosen number of vertices.

#### **6.2 LU Multifrontal Method**

The multifrontal method (Section 5.4) can be generalized to nonsymmetric *A* by modifying the definitions of the frontal matrices and generated elements to conform to an LU factorization. But natural generalizations to rectangular frontal and generated element matrices do not simultaneously satisfy the statements from Observation 5.1. These statements can be rewritten for the LU factorization as follows.

(a) Each generated element *Vj* is used only once to contribute to a frontal matrix.

(b) The row and column index lists for the rectangular frontal matrix *Fj* correspond to the nonzeros in column *Lj* :*n,j* and nonzeros in row *Uj,j* :*n*, respectively.

These conditions cannot both hold. An approach that satisfies (a) can be based on the sparsity pattern of S{*A* + *AT* } and storing some explicit zeros if S{*A*} is not symmetric. It is then analogous to the symmetric multifrontal method. In this case, although the frontal and generated elements may be numerically nonsymmetric, they are square and structurally symmetric. This approach performs well if S{*A*} is close to symmetric, that is, the symmetry index of *A* is close to unity.

An approach that satisfies (b) and not necessarily (a) splits the generated elements into smaller ones that are embedded into further rectangular frontal matrices. We illustrate this using the example from Figure 6.7, that is, ⎛⎞

⎟

*,*

⎜

 ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 2 3 4 5 6 7 8 9 10 1 ∗ ∗ 2 ∗ ∗ *f* ∗ ∗∗ 3 ∗ ∗∗ 4 ∗∗ ∗ ∗ 5 ∗∗∗ ∗ *f f* 6 ∗ *fff* ∗ *ff f* 7 ∗ ∗ *f f* 8 ∗ ∗ *ffff* ∗ *f f* 9 ∗ ∗ 10 ∗ *f* ∗ *f f* ∗ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

where ∗ are entries in *A* and filled entries in *L* + *U* are denoted by *f* . Taking the entries in the first row and column, the sparsity patterns of the first frontal matrix and the corresponding generated element are ⎛⎜⎞⎟⎛⎞

$$F\_1 = \begin{matrix} 1 & 3 \\ 2 & \* \\ 3 & \* \\ 8 & \* \end{matrix}, \quad V\_1 = \begin{matrix} 3 \\ 2 \\ \* \\ 8 \end{matrix}$$

To construct *F*2 that satisfies (b) we can only use part of *V*1. From the row and column replication principles, because *a*13 = 0, the sparsity pattern of column 1 is replicated in that of column 3 of the factors. While the entry in position *(*2*,* 3*)* belongs to *F*2, because of the row replication of the sparsity pattern of the first row in that of the second row, the remaining entries contribute to *F*3 and so we split *V*1 into two as follows

\*\*Lemma:\*\*,\*\* une normam \*\*cours comuncave to 1 y una \*\*cours\*\*
$$\begin{array}{ll} \mathfrak{3} & \mathfrak{3} \\ \mathrm{V}\_{1}^{2} = 2 & \left( \begin{array}{c} f \\ \end{array} \right), & \mathrm{V}\_{1}^{3} = \frac{3}{8} \begin{pmatrix} \* \\ f \end{pmatrix}, & \mathrm{V}\_{1} = \mathrm{V}\_{1}^{2} \oplus \mathrm{V}\_{1}^{3}, \end{array}$$

where ←→ is the extend-add operator and *V* 2 1 and *V* 3 1 contribute to *F*2 and *F*3, respectively. Then *F*2 and the corresponding generated element *V*2 are ⎛⎞⎛⎞

⎠

⎛

⎞

⎞

⎟

⎝

$$\begin{array}{ccccccccc} & 2 & 5 & 8 & 10 & & & 2 & 3 & 5 & 8 & 10 \\ & 2 & \left(\* & \* & \* & \* & \\ & \* & & & & & \oplus & V\_1^2 = 6 \\ & & \* & & & & 8 & \left(\* & & & \end{\*}\right) & \bigoplus \, V\_1^2 = \begin{array}{ccccc} & & 3 & 5 & 8 & 10 \\ \* & & & & & & \\ \* & & & & & & \end{array} \end{array}$$

Consider the following splitting of *V*2

⎠

$$\mathbb{V}\_2 = \overset{\mathfrak{G}}{\underset{\mathfrak{S}}{6}} \begin{pmatrix} f \\ f \end{pmatrix} \oplus \overset{\mathfrak{G}}{\underset{\mathfrak{S}}{6}} \begin{pmatrix} f \\ f \end{pmatrix} \oplus \overset{\mathfrak{G}}{\underset{\mathfrak{S}}{6}} \begin{pmatrix} f & f \\ \* & f \end{pmatrix} \equiv V\_2^3 \rightsquigarrow V\_2^5 \oplus V\_2^6.$$

⎛

⎜

The next frontal matrix is

$$\begin{array}{ccccc} & & 3 & 4 & & & 4\\ & 3 & 4 & & & & 4\\ F\_3 = \frac{3}{4} \left( \begin{matrix} \* & \*\\ \* & \* \end{matrix} \right) \oplus \begin{matrix} \* & \* \end{matrix} \oplus \begin{array}{c} & & & 4\\ \* & \* & \* \end{matrix} & & & & & 4\\ \oplus \begin{pmatrix} f & \*\\ f & \* \end{pmatrix}, & V\_3 = 6 \begin{pmatrix} \*\\ f\\ f \end{pmatrix} \cdot \\\ \end{array}$$

The subsequent steps can be described in a similar way.

Theorem 6.16 expresses the nested relationship between the nonsymmetric multifrontal method and the nonsymmetric elimination tree.

**Theorem 6.16 (Eisenstat & Liu 2005b)** *Assume A is a general nonsymmetric matrix and t* = *parent (k) in* T *(A). Then*

$$
\mathcal{S}\{L\_{t\wedge n,k}\} \subseteq \mathcal{S}\{L\_{t\wedge n,l}\} \quad \text{and} \quad \mathcal{S}\{U\_{k,t\wedge n}\} \subseteq \mathcal{S}\{U\_{l,t\wedge n}\}.
$$

*Proof* Because *t* is the parent of *k*, by definition *t* G*(L)* ⇒ *k* G*(U )* ⇒ *t*. If *uij* = 0, then a multiple of column *i* is added to column *j* during the LU factorization. Thus, by a simple induction argument, for each *j* on the path *k* G*(U )* ⇒ *t*, we must have S{*Lj* :*n,k*} ⊆ S{*Lj* :*n,j* }*.* In particular, this holds for column *t*. The second part follows by a similar argument using the path *t* G*(L)* ⇒ *k*.

This result shows that the parent relationship in the nonsymmetric elimination tree guarantees that both row and column replications can be applied at the same time. Hence all entries of the submatrices of the generated element *Vk* with indices greater than or equal to *parent (k)* can be added to *Vparent (k)* using the operation ←→ . To illustrate this, consider again the 10 × 10 example above for which *parent (*1*)* = 3. Theorem 6.16 guarantees that *V*1 can be embedded into *F*3 because S{*L*3:*n,*1} ⊆ S{*L*3:*n,*3} and S{*U*1*,*3:*n*} ⊆ S{*U*3*,*3:*n*}.

⎝

#### **6.3 Preprocessing Sparse Matrices**

We now turn our attention to preprocessing techniques that can help in computing an LU factorization. In particular, we consider when *A* does not have a full transversal (that is, it has one or more zeros on the diagonal). For numerical stability and to reduce the number of permutations required during the factorization, it can be useful to permute *A* before the factorization begins to put large nonzero entries on the diagonal. As an example, consider the matrix *A* in Figure 6.9. It has *a*22 = 0 and we want to know whether it can be permuted so that all the diagonal entries are nonzero. This question and its answer can be formulated in terms of matchings and bipartite graphs.

#### *6.3.1 Bipartite Graphs and Matchings*

Given a graph G = *(*V*,* E*)*, an edge subset M ⊆ E is called a **matching** (or assignment) if no two edges in M are incident to the same vertex. In matrix terms, a matching corresponds to a set of nonzero entries with no two belonging to the same row or column. A vertex is matched if there is an edge in the matching incident on the vertex, and is unmatched (or free) otherwise. The **cardinality** of a matching is the number of edges in it. A **maximum cardinality matching** (or **maximum matching**) is a matching of maximum cardinality. A matching is **perfect** if all the vertices are matched.

A **bipartite graph** is an undirected graph whose vertices can be partitioned into two disjoint sets such that no two vertices within the same set are adjacent, that is, each set is an **independent set**. Let the *n*×*n* matrix *A* have entries {*aij* }. Associated with *A* is a bipartite graph defined as a triple G*b* = *(*V*row,* V*col,* E*)*, where the row vertex set V*row* = {*i* |*aij* = 0} and the column vertex set V*col* = {*j* |*aij* = 0} correspond to the rows and columns of *A* and there is an (undirected) edge *(i, j )* ∈ E if and only if *aij* = 0. This is illustrated in Figure 6.9. We use prime to distinguish between the independent set of row vertices and the independent set of column vertices, that is, *i* denotes a row vertex and *i* denotes a column vertex.

If *A* is structurally nonsingular, a matching M in G*b* is perfect if it has cardinality *n*. A perfect matching defines an *n* × *n* permutation matrix *Q* with entries *qij* given by

$$q\_{ij} = \begin{cases} 1, & \text{if } (j, i') \in \mathcal{M}, \\ 0, & \text{otherwise}. \end{cases}$$

Both *QA* and *AQ* have the matching entries on the (zero-free) diagonal. *Q* and the column permuted matrix *AQ* for the example in Figure 6.9 are given in Figure 6.10.

**Figure 6.9** A sparse matrix and its bipartite graph G*b* (left). The matched matrix entries are in blue and edges that belong to a perfect matching in G*b* are given by the blue dashed lines (right). Note that the perfect matching is not unique (an alternative is in Figure 6.11).

Q = ⎛ ⎜⎜⎜⎜⎜⎝ 123456 1 1 2 1 3 1 4 1 5 1 6 1 ⎞ ⎟⎟⎟⎟⎟⎠ AQ = ⎛ ⎜⎜⎜⎜⎜⎝ 3- 1- 4- 2- 5- 6- 1 ∗ ∗ 2 ∗ ∗ ∗ 3 ∗ ∗ 4 ∗ ∗ 5 ∗ ∗ 6 ∗ ∗ ∗ ⎞ ⎟⎟⎟⎟⎟⎠

**Figure 6.10** The permutation matrix *Q* and the column permuted matrix *AQ* corresponding to the matrix in Figure 6.9. The matched entries are on the diagonal of *AQ*.

#### *6.3.2 Augmenting Paths*

If a perfect matching exists, it can be found using augmenting paths. A path P in a graph is an ordered set of edges in which successive edges are incident to the same vertex. P is called an M**-alternating path** if the edges of P are alternately in M and not in M. An M-alternating path is an M**-augmenting path** in G*b* if it connects an unmatched column vertex with an unmatched row vertex. Note that the length of an M-augmenting path is an odd integer.

## **ALGORITHM 6.3 Maximum matching algorithm Input:** An undirected graph. **Output:** Output maximum matching. 1: Find an initial matching M For example, M = ∅ 2: **while** there exists a M-augmenting path P **do**

**Figure 6.11** An illustration of the search for a perfect matching using augmenting paths. On the left, the dashed lines represent a matching with cardinality 5. In the centre, the blue line is an augmenting path with end vertices 2 and 2 . On the right is the perfect matching with cardinality 6 that is obtained using the augmenting path.

Let M and P be subsets of E and define the symmetric difference

$$
\mathcal{M} \oplus \mathcal{P} := (\mathcal{M} \backslash \mathcal{P}) \cup (\mathcal{P} \backslash \mathcal{M}),
$$

that is, the set of edges that belongs to either M or P but not to both. If M is a matching and P is an M-augmenting path, then M ⊕ P is a matching with cardinality |M|+1. Growing the matching in this way is called augmenting along P. The next result shows that augmenting paths can be used to find a maximum matching (Algorithm 6.3).

**Theorem 6.17 (Berge 1957)** *A matching* M *in an undirected graph is a maximum matching if and only if there is no* M*-augmenting path*

Figure 6.11 demonstrates the procedure. On the left is a bipartite graph with a matching with cardinality 5. In the centre, an augmenting path 2 ⇒ 3 ⇒ 3 ⇒ 4 ⇒ 4 ⇒ 2 is shown. Augmenting the matching along this path, the cardinality of the matching increases to 6 and M ⊕ P is a perfect matching.

### *6.3.3 Weighted Matchings*

While the maximum matching algorithm finds a permutation of *A* such that the permuted matrix has nonzero diagonal entries, there are more sophisticated variations that aim to ensure the absolute values of the diagonal entries of the permuted matrix (or their product) are in some sense large. This can increase the likelihood that the permuted matrix is strongly regular and reduce the need for partial pivoting during the LU factorization. The core problem is as follows: given an *n* × *n* matrix *A*, find a matching of the rows to the columns such that the product of the matched entries is maximized. That is, find a permutation vector *q* that maximizes

$$\prod\_{l=1}^{n} |a\_{lq\_l}|.\tag{6.5}$$

Define a matrix *C* corresponding to *A* with entries *cij* ≥ 0 as follows:

$$c\_{ij'} = \begin{cases} \log(\max\_l |a\_{lj'}|) - \log |a\_{lj'}|, & \text{if } a\_{lj'} \neq 0\\ \infty, & \text{otherwise.} \end{cases} \tag{6.6}$$

It is straightforward to see that finding a *q* that solves (6.5) is equivalent to finding a *q* that minimizes

$$\sum\_{l=1}^{n} |c\_{lq\_l}|,\tag{6.7}$$

which is equivalent to finding a minimum weight perfect matching in an edge weighted bipartite graph. This is a well-studied problem and is known as the bipartite weighted matching or linear sum assignment problem.

If G*b* = *(*V*row,* V*col,* E*)* is the bipartite graph associated with *A* then let G*b(C)* = *(*V*row,* V*col,* E*)* be the corresponding weighted bipartite graph in which each edge *(i, j )* ∈ E has a weight *cij* ≥ 0. The weight (or cost) of a matching M in G*b(C)*, denoted by *csum(*M*)*, is the sum of its edge weights; i.e. *csum(*M*)* =

$$csum(\mathcal{M}) = \sum\_{(i,j') \in \mathcal{M}} c\_{ij'}.$$

A perfect matching M in G*b(C)* is said to be a **minimum weight perfect matching** if it has smallest possible weight, i.e. *csum(*M*)* ≤ *csum(*M*)* for all possible perfect matchings M.

The key concept for finding a minimum weight perfect matching is that of a **shortest augmenting path**. An M-augmenting path P starting at an unmatched column vertex is called **shortest** if

$$
\operatorname{csum}(\mathcal{M} \oplus \mathcal{P}) \le \operatorname{csum}(\mathcal{M} \oplus \widehat{\mathcal{P}}).
$$

for all other possible M-augmenting paths P starting at the same column vertex. A matching M*e* is **extreme** if and only if there exist *ui* and *vj* (which are termed **dual variables**) satisfying

$$\begin{cases} c\_{lj'} = u\_l + v\_{j'}, & \text{if } (i, j') \in \mathcal{M}\_{\mathfrak{e}}, \\ c\_{lj'} \ge u\_l + v\_{j'}, & \text{otherwise.} \end{cases} \tag{6.8}$$

This is employed by the MC64 algorithm. The dual variables will be important when we discuss scaling sparse matrices in Section 7.4.2. The MC64 algorithm is outlined here as Algorithm 6.4. It starts with a feasible solution and corresponding extreme matching and then proceeds to iteratively increase its cardinality by one by constructing a sequence of shortest augmenting paths until a perfect extreme matching is found. The algorithm can be made more efficient if a large initial extreme matching can be found. For example, Step 3 can be replaced by setting *ui* = min{*cij* | *j* ∈ S{*Ai,*1:*n*}} for *i* ∈ V*row* and *vj* = min{*cij* − *ui*| *i* ∈ S{*A*1:*n,j* }} for *j* ∈ V*col*. In Step 4, an initial extreme matching can be determined from the edges for which *cij* − *ui* − *vj* = 0.

There are a number of potential problems with the MC64 algorithm. First, the runtime is hard to predict and depends on the initial ordering of *A*. Second, it is a serial algorithm and as such it can represent a significant fraction of the total factorization time of a direct solver. Because the complexity of Step 6 of Algorithm 6.4 is *O((n* + *nz(A))*log *n)* and the complexity of Step 7 is *O(n)* and of Step 8 is *O(n*+*nz(A)*, MC64 has a worst-case complexity of *O(n(n*+*nz(A))*log *n)*. In practice, this bound is not achieved and the algorithm is widely used.

#### **ALGORITHM 6.4 Outline of the MC64 algorithm**

**Input:** Matrix *A*.

**Output:** A matching M and dual variables *ui*, *vj* .

solution

4: Set M = {*(i, j*

*)*| *ui* + *vj* } Initial extreme matching

#### *6.3.4 Dulmage-Mendelsohn Decompositions*

The importance of preordering *A* to block triangular form was discussed in Section 3.4. The **Dulmage-Mendelsohn decomposition** is based on matchings and is a generalization of the block triangular form. It provides a precise characterization of structurally rank deficient matrices and it can be used to reduce the work required for an LU factorization. It comprises row and column permutations *P* and *Q* such that ⎛⎝⎞⎠

$$
\begin{array}{c}
\mathcal{C}\_1 \quad \mathcal{C}\_2 \quad \mathcal{C}\_3 \\
\mathcal{R}^A \mathcal{Q} = \begin{array}{c}
\mathcal{R}\_1 \\
\mathcal{R}\_2 \\
\mathcal{R}\_3
\end{array} \begin{pmatrix}
A\_1 & A\_4 & A\_6 \\
0 & A\_2 & A\_5 \\
0 & 0 & A\_3
\end{pmatrix} .
\end{array} \tag{6.9}
$$

Here *A*1 is an *m*1 × *n*1 underdetermined matrix (*m*1 *< n*1 or *m*1 = *n*1 = 0), *A*2 is an *m*2 × *m*2 square matrix and *A*3 is an *m*3 × *n*3 overdetermined matrix (*m*3 *> n*3 or *m*3 = *n*3 = 0). It can be shown that *AT* 1 and *A*3 are strong Hall matrices but *A*2 need not be a strong Hall matrix, in which case it can be permuted to block upper triangular form.

If row and column sets R and C form a maximum matching of *A*, then R1 and R2 are subsets of R and |R3 ∩ R| = *n*3, and C2 and C3 are subsets of C and |C1 ∩ C| = *m*1. An example decomposition for a 10 × 10 matrix is given in Figure 6.12. Here R = {1*,* 2*,...,* 9} and C = {2*,* 3*,...,* 10}.

The **coarse Dulmage-Mendelsohn decomposition** orders the unmatched columns as the first columns in *P AQ* and orders the unmatched rows as the last rows in *P AQ.* If *A* is square and has a perfect matching then its coarse decomposition has only the matrix *A*2; otherwise, both *A*1 and *A*3 are present. The coarse decomposition is computed by first finding a maximum matching. Assuming it is not a perfect matching, the rows in *A*1 are determined by performing depth-first searches from the unmatched columns to find all of the row vertices that

.

**Figure 6.12** An example of a coarse Dulmage-Mendelsohn decomposition. The blue entries belong to the maximum matching. *m*1 = 3, *m*2 = 4, *m*3 = 3, *n*1 = 4, *n*2 = 4, *n*3 = 2. Column 1 and row 10 are unmatched.

are reachable from the unmatched columns via alternating augmenting paths. The columns in *A*1 are defined to be the union of the set of unmatched columns and the set of columns matched with the rows in *A*1. Similarly, the columns in *A*3 are determined by performing depth-first searches from the unmatched rows to find all of the column vertices that are reachable from the unmatched rows via alternating augmenting paths. The rows in *A*3 are defined to be the union of the set of rows matched to the columns in *A*3 and the set of unmatched rows.

It may be possible to further permute the matrix to obtain the **fine Dulmage-Mendelsohn decomposition**. The fine Dulmage-Mendelsohn decomposition computes *P* and *Q* such that *A*1 and *A*3 are block diagonal matrices in which each diagonal block is irreducible, and *A*2 is block upper triangular with strongly connected (square) diagonal blocks. Once the coarse decomposition has been computed, *A*1 and *A*3 are searched to find any irreducible blocks and the permutations required to place these on the diagonals of *A*1 and *A*3 are computed. Finally, following Section 3.4, strongly connected components in *A*2 are found and a permutation is formed to reduce *A*2 to block upper triangular form (with the strongly connected components lying on the diagonal). If *A* is reducible and nonsingular, the fine Dulmage-Mendelsohn decomposition can be used to solve the linear system *Ax* = *b* using block back-substitution.

#### **6.4 Notes and References**

Early theoretical results related to sparse LU factorizations can be found in Rose & Tarjan (1978), which extends the systematic understanding of the symbolic elimination introduced in Rose et al. (1976). A key paper that influenced the discussion and development of both the theory and algorithms for predicting sparsity structures in LU factorizations is Gilbert (1994) (first available in 1986 as a Cornell technical report). As the primary and still very useful resource on transitive reduction, we refer to Aho et al. (1972); Gilbert & Liu (1993) extend the concept of an elimination tree to study sparse LU factorizations of nonsymmetric matrices and present theoretical concepts based on DAGs; see also the parallel counterpart in Grigori et al. (2007). Ways to simplify symbolic factorizations and prune DAGs are discussed in Eisenstat & Liu (1992, 1993a). An elegant treatment of both the theoretical and practical aspects of LU factorizations based on DAGs and the nonsymmetric elimination tree (including pruning and pivoting) is given in a series of three papers by Eisenstat & Liu (2005a,b, 2007).

Partial pivoting within the sparse column LU factorization is introduced in Gilbert & Peierls (1988). This paper influenced not only further developments in sparse LU factorizations but also the development of incomplete factorizations. Partial pivoting based on the column elimination tree is first discussed in George & Ng (1985); see also Gilbert & Ng (1993) and Li (1996) for further use of column elimination trees. Further research on exactness of structural predictions is presented by Grigori et al. (2009).