## 2 Sparse Direct Algorithms
### An Introduction to the Factorize Phase  
The factorization can be performed in many different ways, depending on the order in which matrix entries are accessed and/or updated. Possible variants include left-looking (fan-in), right-looking (fan-out) and multifrontal algorithms. The (supernodal) right-looking variant computes a (block) row and column of \(L\) at each step and uses them to immediately update all rows and columns in the part of the matrix that has not yet been factored. In the (supernodal) left-looking variant, the updates are not applied immediately; instead, before a (block) column \(k\) is eliminated, all updates from previous columns of \(L\) are applied together to the (block) column \(k\) of \(A\). In the multifrontal method [25; 26], the updates are accumulated; they are propagated from a descendant column \(j\) to an ancestor column \(k\) via all intermediate nodes on the elimination tree path from \(j\) to \(k\). Two of the packages in this study (HSLMA77 and HSLMA97) implement multifrontal methods, while HSLMA86 for indefinite systemsand HSL_MA87 for positive-definite systems use a supernodal left-right looking approach. In Sections 3 and 4, we discuss the two approaches in more detail and, in particular, consider how parallelism can be exploited. For symmetric matrices that are positive definite, the factorization phase can use the chosen pivot order without modification. Moreover, the data structures determined by the analyse phase can be static throughout the factorization phase. For symmetric indefinite problems, using the pivot order from the analyse phase may be unstable or impossible because of (near) zero diagonal entries. Thus, in general, it is necessary to modify the pivot sequence to ensure numerical stability. If symmetry is to be maintained, \(1\times 1\) and \(2\times 2\) pivoting must be performed, and the resulting factorization takes the form  
\[PAP^{T}=LDL^{T}\]  
where \(D\) is a block diagonal matrix with \(1\times 1\) and \(2\times 2\) blocks. The stability of the factorization of symmetric indefinite systems was considered in detail in the paper by [PERSON], [PERSON] and [PERSON] [27]. They showed that bounding the size of the entries of \(L\), together with a backward stable scheme for solving \(2\times 2\) linear systems, suffices to show backward stability for the entire solution process. They found that the widely used strategy of [PERSON] and [PERSON] [28] does not have this property, whereas the threshold pivoting technique first used by [PERSON] and [PERSON] [25] in their original multifrontal solver does. [PERSON] and [PERSON] choose the pivots one-by-one, with the aim of limiting the size of the entries \(l_{i,j}\) in \(L\)  
\[|l_{i,j}|<u^{-1} \tag{1}\]  
where the threshold \(u\) is a user-set value in the range 0 \(<u\leq\) 1.0. Suppose that \(q\) denotes the number of rows and columns of \(D\) found so far. Let \(a_{i,j}\), with \(i>q\) and \(j>q\), denote an entry of the matrix after it has been updated by all the permutations and pivot operations so far. For a \(1\times 1\) pivot in column \(j\) = \(q+1\), the requirement for inequality (1) corresponds to the stability threshold test  
\[|a_{q+1,q+1}|>u\max_{i>q+1}|a_{i,q+1}|. \tag{2}\]  
However, it is not always possible to select a valid \(1\times 1\) pivot. It is sufficient to use \(2\times 2\) pivots in these cases (see, for example, Section 4.4 of [9]). Following [PERSON] _et al_. [29], the corresponding stability test for a \(2\times 2\) pivot is  
\[\left| \begin{pmatrix}a_{q+1,q+1}&a_{q+1,q+2}\\ a_{q+1,q+2}&a_{q+2,q+2}\end{pmatrix} ^{-1}\right| \begin{pmatrix}\max_{i>q+2}|a_ {i,q+1}|\\ \max_{i>q+2}|a_{i,q+2}|\end{pmatrix} < \begin{pmatrix}u^{-1}\\ u^{-1}\end{pmatrix}  \tag{3}\]  
where the absolute value notation for a matrix refers to the matrix of corresponding absolute values. The choice of \(u\) controls the balance between stability and sparsity: in general, the smaller \(u\) is, the sparser \(L\) will be, while the larger \(u\) is, the more stable the factorization will be.