4.1.1 Partial pivoting

The LAPACK block LU factorization is the main point of reference here, and the LAPACK naming convention is followed. The initial ‘D’ in names denotes double precision routines. The LU factorization of a matrix $A$ has the form

$$\begin{eqnarray}PA=LU,\end{eqnarray}$$

where $L$ is a unit lower triangular matrix, $U$ is an upper triangular matrix, and $P$ is a permutation matrix. The LAPACK algorithm proceeds in the following steps. Initially, a set of $n_{b}$ columns (the panel) is factored and a pivoting pattern is produced (implemented by the DGETF2 routine). Then the elementary transformations, resulting from the panel factorization, are applied in a block fashion to the remaining part of the matrix (the trailing submatrix). This involves swapping up to $n_{b}$ rows of the trailing submatrix (DLASWP), according to the pivoting pattern, application of a triangular solve with multiple right-hand sides to the top $n_{b}$ rows of the trailing submatrix (DTRSM), and finally application of matrix multiplication of the form $A_{ij}\leftarrow A_{ij}-A_{ik}\times A_{kj}$ (DGEMM), where $A_{ik}$ is the panel without the top $n_{b}$ rows, $A_{kj}$ is the top $n_{b}$ rows of the trailing submatrix, and $A_{ij}$ is the trailing submatrix without the top $n_{b}$ rows. Then the procedure is applied repeatedly, descending down the diagonal of the matrix: see Figure 4.1. The block algorithm is described in detail in section 2.6.3 of the book by Demmel (Reference Demmel1997).

Figure 4.1.

The block LU factorization (Level 3 BLAS algorithm of LAPACK).

Instead of the Level 2 BLAS panel factorization (DGETF2) of LAPACK, we use a recursive partial pivoting panel factorization, which uses Level 3 BLAS operations. For a panel of $k$ columns, the algorithm recursively factors the left $k/2$ columns, updates the right $k/2$ columns with DLASWP, DTRSM, and DGEMM, then recursively factors the right $k/2$ columns, and finally applies swaps to the left $k/2$ columns with DLASWP. The recursion terminates at a single column, where it simply searches for the maximum pivot, swaps elements, and scales the column by the pivot. Parallelism is introduced by splitting the panel into $p$ block rows, which are assigned to different processors.

4.1.2 Incremental pivoting

The most performance-limiting aspect of Gaussian elimination with partial pivoting is the panel factorization operation. First, it is an inefficient operation, usually based on a sequence of calls to Level 2 BLAS. Second, it introduces synchronization, by locking an entire panel of the matrix at a time. Therefore, it is desirable to split the panel factorization into a number of smaller, finer-granularity operations, which is the basic premise of the incremental pivoting implementation, also known in the literature as the tile LU factorization.

In this algorithm, instead of factoring the panel one column at a time, the panel is factored one tile at a time. The operation proceeds as follows. First the diagonal tile is factored using the standard LU factorization procedure. Then the factored tile is combined with the tile directly below it and factored. Then the refactored diagonal tile is combined with the next tile and factored again. The algorithm descends down the panel until the bottom of the matrix is reached. At each step, the standard partial pivoting procedure is applied to the tiles being factored. Also, at each step, all the tiles to the right of the panel are updated with the elementary transformations resulting from the panel operations: see Figure 4.2. This way of pivoting is basically the idea of pairwise pivoting applied at the level of tiles, rather than individual elements. The main benefit comes from the fact that updates of the trailing submatrix can proceed alongside panel factorizations, leading to a very efficient parallel execution, where multiple steps of the algorithm are smoothly pipelined.

Figure 4.2.

Incremental LU factorization.

A straightforward implementation of incremental pivoting incurs more floating-point operations than partial pivoting. To reduce this additional cost, a second level of blocking, called inner blocking, is used within each tile (Joffrain, Quintana-Orti and van de Geijn Reference Joffrain, Quintana-Orti, van de Geijn, Dongarra, Madsen and Waśniewski2006, Buttari, Langou, Kurzak and Dongarra Reference Buttari, Langou, Kurzak and Dongarra2009). The total cost to factor an $n\times n$ matrix is

$$\begin{eqnarray}{\textstyle \frac{2}{3}}n^{3}\biggl(1+\frac{i_{b}}{2n_{b}}\biggr)+O(n^{2}),\end{eqnarray}$$

where $n_{b}$ is the tile size and $i_{b}$ is the inner blocking size. A small value for $i_{b}$ decreases the flops but may hurt the efficiency of Level 3 BLAS operations, so $n_{b}$ and $i_{b}$ should be chosen to offer the best compromise between extra operations and BLAS efficiency.

With $n_{b}=1$, incremental pivoting reduces to pairwise pivoting, hence the stability of incremental pivoting is assumed to be the same as pairwise pivoting (Joffrain et al. Reference Joffrain, Quintana-Orti, van de Geijn, Dongarra, Madsen and Waśniewski2006). For pairwise pivoting, the growth factor bound is $4^{n-1}$. Trefethen and Schreiber (Reference Trefethen and Schreiber1990) observed an average growth factor of $n$ for $n\leq 1024$, but more recently, Grigori, Demmel and Xiang (Reference Grigori, Demmel and Xiang2011) observed an average growth factor ${>}~n$ for $n\geq 4096$, leaving the question of stability open for large $n$.

4.1.3 Tournament pivoting

Classic approaches to the panel factorization with partial pivoting communicate asymptotically more than the established lower bounds (Demmel, Grigori, Hoemmen and Langou Reference Demmel, Grigori, Hoemmen and Langou2008a). The basic idea of communication-avoiding LU is to minimize communication by replacing the pivot search performed at each column with a block reduction of all the pivots together. This is done thanks to a new pivoting strategy referred to as tournament pivoting, which performs extra computations and is shown to be stable in practice. Tournament pivoting factors the panel in two steps. First, using a tournament selection, it identifies rows that can be used as good pivots for the factorization of the whole panel. Second, it swaps the selected pivot rows to the top of the panel, and then factors the entire panel without pivoting. With this strategy, the panel is efficiently parallelized, and the communication is provably minimized.

Figure 4.3 presents the first step of tournament pivoting for a panel $W$ using a binary tree for the reduction operation. First, the panel is partitioned into $p_{r}$ blocks, that is, $W=[W_{00},W_{10},\ldots ,W_{p_{r}-1,0}]$, where $W_{ij}$ represents the block owned by thread $i$ at step $j$ of the reduction operation. Figure 4.3 shows an example with $p_{r}=4$.

Figure 4.3.

Example of tournament pivoting with $p_{r}=4$ processors, using a binary tree for the reduction operation.

At the first step of the reduction operation, each thread, $i$, applies Gaussian elimination with partial pivoting to its block, $W_{i0}$; then the resulting permutation matrix $P_{i0}$ is applied to the original unfactored block, $W_{i0}$, and the first $n_{b}$ rows of the permuted block $P_{i0}W_{i0}$ are selected as pivot candidates. These first pivot candidates represent the leaves of the reduction tree. At each node of the tree, the pivot candidates of two child nodes are merged on top of each other, and Gaussian elimination with partial pivoting is applied on the merged block. The resulting permutation matrix is then applied on the original unfactored merged block and the first $n_{b}$ rows are selected as new pivot candidates. By using a binary tree, this step is repeated $\log _{2}p_{r}$ times. The pivots obtained at the root of the tree are then considered to be good pivots for the whole panel. Once these pivots are permuted to the top of the panel, each thread applies Gaussian elimination without partial pivoting to its block, $W_{i0}$.

The example presented in Figure 4.3 uses a binary tree with two tiles reduced together at each level, but any reduction tree can be used, depending on the underlying architecture. The tournament pivoting implementation in PLASMA, used for experiments in this section, reduces four tiles at each level. This number of tiles is chosen because it gives a good ratio of kernel efficiency over one single core, relative to the time spent to perform the factorization of the subset.

For tournament pivoting, the best growth factor bound proved so far is $2^{nH}$, where $H$ is the height of the reduction tree. With a binary tree, $H=\log _{2}p_{r}$, so the bound is $p_{r}^{n}$. However, Grigori et al. (Reference Grigori, Demmel and Xiang2011) have never observed a growth factor as large as $2^{n-1}$, even for pathological cases such as the Wilkinson matrix (Wilkinson Reference Wilkinson1988), and they conjecture that the upper bound is the same as partial pivoting, $2^{n-1}$. They observed an average growth factor of $1.5n^{2/3}$, independent of $p$. Our experiments in Section 4.2.3 confirm that tournament pivoting is as stable as partial pivoting in practice.

4.1.4 Random butterfly transform

As an alternative to pivoting, the partial random butterfly transformation (PRBT) preconditions the matrix as $A_{r}=W^{T}\!AV$, such that, with probability close to 1, pivoting is unnecessary. This technique was proposed by Parker (Reference Parker1995b) and later adapted by Baboulin, Dongarra, Herrmann and Tomov (Reference Baboulin, Dongarra, Herrmann and Tomov2013) to reduce the computational cost of the transformation. An $n\times n$ butterfly matrix is defined as

$$\begin{eqnarray}B^{(n)}=\frac{1}{\sqrt{2}}\left[\begin{array}{@{}cc@{}}R & \phantom{-}S\\ R & -S\end{array}\right]\!,\end{eqnarray}$$

where $R$ and $S$ are random diagonal, non-singular matrices. $W$ and $V$ are recursive butterfly matrices of depth $d$, defined by

$$\begin{eqnarray}W^{(n,d)}=\left[\begin{array}{@{}cccc@{}}B_{1}^{(n/2^{d-1})} & 0 & \ldots \, & 0\\ 0 & \ddots & & \vdots \\ \vdots & & \ddots & 0\\ 0 & \ldots \, & 0 & B_{2^{d-1}}^{(n/2^{d-1})}\\ \end{array}\right]\times \cdots \times B^{(n)}.\end{eqnarray}$$

We use a depth $d=2$, previously found to be sufficient in most cases (Baboulin et al. Reference Baboulin, Dongarra, Herrmann and Tomov2013). Since each $R$ and $S$ is diagonal, $W$ and $V$ can be stored as $n\times d$ arrays. Due to the regular sparsity pattern, multiplying an $m\times n$ matrix by $W$ and $V$ is an efficient, $O(dmn)$ operation.

After applying the PRBT, Gaussian elimination without pivoting is used to obtain the solution, as indicated in Algorithm 1.

While the randomization reduces the need for pivoting, the lack of pivoting can still be unstable, so we use iterative refinement to reduce the potential for instability. The cost of pivoting is thus avoided, at the expense of applying the PRBT and iterative refinement. As with no-pivoting LU, the growth factor remains unbounded. This can easily be seen by letting $A=W^{-T}BV^{-1}$, where the first pivot, $b_{00}$, is zero. However, probabilistically no-pivoting LU (and hence PRBT) has been shown to have an average growth factor of $n^{3/2}$ (Yeung and Chan Reference Yeung and Chan1995).

4.1.5 No pivoting

This implementation of Gaussian elimination completely abandons pivoting. This can be done very rarely in practice without risking serious numerical consequences, or even a complete breakdown of the algorithm if a zero is encountered on the diagonal. In the following experiments, the implementation serves only as a performance baseline. Dropping pivoting increases performance for two reasons. First, the overhead of swapping matrix rows disappears. Second, the level of parallelism dramatically increases, since the panel operations now become parallel, and can also be pipelined with the updates to the trailing submatrix.

4.2.1 Hardware and software

The experiments were run on an Intel system with 16 cores and an AMD system with 48 cores. The Intel system has two sockets with 8-core Intel Sandy Bridge CPUs clocked at 2.6 GHz, with a theoretical peak of $16~\text{cores}\times 2.6~\text{GHz}\times 8~\text{ops per cycle}\simeq 333~\text{Gflops}$ in double precision arithmetic. The AMD system has eight sockets with 6-core AMD Istanbul CPUs clocked at 2.8 GHz, with a theoretical peak of $48~\text{cores}\times 2.8~\text{GHz}\times 4~\text{ops per cycle}\simeq 538~\text{Gflops}$ in double precision arithmetic.

All presented LU codes were built using the PLASMA framework, relying on the CCRB tile layout and the QUARK dynamic scheduler. The GCC 4.1.2 compiler was used for compiling the software stack and Intel MKL (Composer XE 2013) was used to provide an optimized implementation of serial BLAS. On both systems, to avoid important variation in NUMA effects, the PLASMA framework relies on the hwloc library (Broquedis et al. Reference Broquedis, Clet-Ortega, Moreaud, Furmento, Goglin, Mercier, Thibault and Namyst2010) to efficiently bind the threads to cores located in multiple CPU sockets. And all experiments were run with the numactl command to distribute the data in a round-robin fashion over the NUMA locales of the machine. When the number of threads is insufficient to keep all sockets occupied then the memory pages are only bound to the NUMA locales with active PLASMA threads rather than the entire machine.

4.2.2 Performance

We now study the performance of our implementations on square random matrices in double real precision, and compare their performance with that of the LU factorization in MKL (DGETRF). While the absolute performance is different for single, single-complex, and double-complex precisions, we have observed trends similar to double precision for all precisions. In all cases, the matrix is initially in column-major layout, and we include the conversion to and from tiled layout in the time. Since each of our implementations uses a different pivoting strategy, we ran iterative refinement with all the algorithms to achieve the same level of accuracy, for a fair performance comparison. Namely, for MKL, we used the MKL iterative refinement routine DGERFS, while for the tile algorithms, we implemented a tiled iterative refinement with the same stopping criterion as that in DGERFS, that is, the iteration terminates when one of the following three criteria is satisfied.

1. 1 The componentwise backward error, $\max _{i}|r_{i}|/(|A||\widehat{x}|+|b|)_{i}$, is less than or equal to $((n+1)\ast \mathbf{sfmin})/\mathbf{eps}$, where $r=A\widehat{x}-b$ is the residual vector and $\widehat{x}$ is the computed solution, $\mathbf{eps}$ is the relative machine precision, and $\mathbf{sfmin}$ is the smallest value such that $1/\mathbf{sfmin}$ does not overflow.

2. 2 The componentwise backward error is not reduced by half.

3. 3 The number of iterations is equal to ten.

For all of our experiments, iterative refinement converged in less than ten iterations. We observed that even with partial pivoting, it requires a couple of iterations to satisfy this stopping criterion (see Section 4.2.3). Furthermore, in many cases, the MKL version of DGERFS did not scale as well as our implementation. We suspect this is due to the need to compute $|A||\widehat{x}|+|b|$ at each iteration.

Figure 4.4.

Asymptotic performance comparison of LU factorization algorithms on AMD Opteron.

Figure 4.5.

Asymptotic performance comparison of LU factorization algorithms on Intel Sandy Bridge.

The performance of our implementations is sensitive to the tile size $n_{b}$. Hence, for each matrix dimension $n$ on a different number of cores, we studied the performance of each algorithm with the tile sizes of $n_{b}=80$, $160$, $240$, $320$, and $400$. We observed that on $48$ cores of our AMD machine, the performance is especially sensitive to the tile size, and we tried the additional tile sizes of $n_{b}=340$, $360$, and $380$. In addition, the performance of our incremental pivoting is sensitive to the inner blocking size, and we tried using the block sizes of $i_{b}=10$, $20$, and $40$ for both $n_{b}=80$ and $160$; $i_{b}=10$, $20$, $30$, $40$, $60$, and $80$ for $n_{b}=240$; and $i_{b}=10$, $20$, $40$, and $80$ for both $n_{b}=320$ and $400$. Figures 4.4 and 4.5 show the performance obtained using the tile and block sizes that obtained the highest performance. For all algorithms, we compute the effective Gflops using the flops for partial pivoting, given in Table 4.1. For the solve with iterative refinement, we also include a pair of forward and backward substitutions, $2n^{2}$. For the tile algorithms, we included the data layout conversion time as a part of the solution time. We summarize our findings below.

Table 4.1.

Summary of floating-point operations (add, multiply, divide), floating-point comparisons, best known growth factor bound, and average-case growth factor. Differences in flops compared to partial pivoting are shown in bold. References: [1] Trefethen and Schreiber (Reference Trefethen and Schreiber1990); [2] Grigori et al. (Reference Grigori, Demmel and Xiang2011); [3] Buttari et al. (Reference Buttari, Langou, Kurzak and Dongarra2009); [4] Baboulin et al. (Reference Baboulin, Dongarra, Herrmann and Tomov2013); [5] Yeung and Chan (Reference Yeung and Chan1995).

• ∙ PRBT with the default transformation depth of two added only a small overhead over no-pivoting in all the test cases.

• ∙ In comparison to other pivoting strategies, incremental pivoting could exploit a large number of cores more effectively. As a result, when the performance is not dominated by the trailing submatrix updates (e.g. for a small matrix dimension on 48 cores), it obtained performance that is close to that of no-pivoting. However, for a large matrix dimension, due to the extra computation and special kernels required by incremental pivoting to update the trailing submatrix, its performance was lower than that of the partial or tournament pivoting LU that uses Level 3 BLAS DGEMM for their trailing submatrix updates.

• ∙ In comparison to MKL, our partial pivoting and tournament pivoting LU could effectively utilize a larger number of cores, as seen for medium-sized matrices on multiple sockets. This is probably a result of the extra parallelism attained by using a superscalar scheduler, as well as parallel panel factorizations.

• ∙ In this shared memory environment, partial pivoting outperformed tournament pivoting in all cases. However, this observation should not be extrapolated to distributed memory machines, where the communication latency becomes a significant part of the performance. There, the reduced communication in tournament pivoting may be more favourable.

• ∙ MKL performed well for large matrices, where the trailing submatrix update dominates the performance. Note that the combined cost for all the panels is $O(n_{b}n^{2})$, compared to $O(n^{3})$ for the trailing submatrix update. Moreover, on a single socket, MKL outperformed no-pivoting for a sufficiently large matrix dimension. This could be because a tiled implementation loses efficiency due to the smaller BLAS kernels used during its trailing submatrix updates.

• ∙ For small matrices, iterative refinement incurred a significant overhead, which diminishes for large matrices as the $O(n^{3})$ factorization cost dominates the $O(n^{2})$ solve cost. The refinement overhead was particularly large for MKL, which we suspect is due to $|A||\widehat{x}|+|b|$ not being effectively parallelized in MKL.

4.2.3 Accuracy

To study the numerical behaviour of our implementations, we used the synthetic matrices from two recent papers (Grigori et al. Reference Grigori, Demmel and Xiang2011, Baboulin et al. Reference Baboulin, Dongarra, Herrmann and Tomov2013). We tested all these matrices, but here we present only a representative subset that demonstrates the numerical performance of the pivoting strategies. We also conducted the numerical experiments using all the matrix dimensions used in Section 4.2.2, but here we show only the results of $n=30\,000$, which represent the numerical performance trends of the pivoting strategies for all other matrix dimensions. Table 4.2 shows some properties of these test matrices and the stability results of using partial pivoting. In the table, the second to sixth test matrices are from the paper by Grigori et al. (Reference Grigori, Demmel and Xiang2011), where the first two have relatively small condition numbers, while the rest are more ill-conditioned. The last three matrices are from the paper by Baboulin et al. (Reference Baboulin, Dongarra, Herrmann and Tomov2013), where the last test matrix gfpp is one of the pathological matrices that exhibit an exponential growth factor using partial pivoting. Since the condition number of the gfpp matrix increases rapidly with the matrix dimension, we used the matrix dimension of $m=1000$ for our study. Finally, incremental and tournament pivoting exhibit different numerical behaviour using different tile sizes. For the numerical results presented here, we used the tile and block sizes that obtain the best performance on the 16 core Intel Sandy Bridge. All the results are in double real precision.

Table 4.2.

Properties of test matrices and stability results of using partial pivoting $(n=m=30\,000)$.

Figure 4.6.

Numerical accuracy of LU factorization algorithms, showing error before refinement (top point of each line) and after each refinement iteration (subsequent points). Number of iterations is given at the bottom of (a); ‘F’ indicates failure (e.g. overflow).

Figure 4.6(a) shows the componentwise backward errors,

$$\begin{eqnarray}\max _{i}|r_{i}|/(|A||\widehat{x}|+|b|)_{i},\end{eqnarray}$$

at each step of iterative refinement. For these experiments, the right-hand side $b$ is chosen such that the entries of the exact solution $x$ are uniformly distributed random numbers in the range of $[-0.5,0.5]$. We summarize our findings below.

• ∙ For all the test matrices, tournament pivoting obtained initial backward errors comparable to those of partial pivoting.

• ∙ No-pivoting was unstable for five of the test matrices (ris, fiedler, orthog, {-1,1}, and gfpp). For the rest of the test matrices, the initial backward errors of no-pivoting were significantly greater than those of partial pivoting, but were improved after a few refinement iterations.

• ∙ Incremental pivoting failed for the fiedler, orthog, and gfpp matrices. For other test matrices, its backward errors were greater than those of partial pivoting, but were improved to be of the same order as those of partial pivoting after a few refinement iterations. The only exception was with the ris matrix, where the refinement stagnated before reaching an accuracy similar to that of partial pivoting. In each column of the ris matrix, entries with smaller magnitudes are closer to the diagonal (i.e. $a_{ij}=0.5/(n-i-j+1.5)$). As a result, the permutation matrix $P$ of partial pivoting has ones on the anti-diagonal.

• ∙ When no-pivoting was successful, its backward errors were similar to those of PRBT. On the other hand, PRBT was more stable than no-pivoting, being able to obtain small backward errors for the fiedler, {-1,1}, and gfpp matrices.

• ∙ Partial pivoting was not stable for the pathological matrix gfpp. On the other hand, PRBT randomizes the original structure of the matrix and was able to compute the solution to reasonable accuracy. It is also possible to construct pathological test matrices where partial pivoting is unstable while tournament pivoting is stable, and vice versa (Grigori et al. Reference Grigori, Demmel and Xiang2011).

Figure 4.6(b) shows the relative forward error norms of our implementations, which were computed as $\Vert x-\widehat{x}\Vert _{\infty }/\Vert x\Vert _{\infty }$. In the convergence of the forward error norms we observed similar trends to the backward errors. One difference was with the orthog test matrix, where iterative refinement could not adequately improve the forward errors of incremental pivoting and PRBT. Also, even though the backward errors of the ris test matrix were of the order of machine epsilon with partial and tournament pivoting, their relative forward errors were $O(1)$ due to the large condition number.

7.1.1 LAPACK: Bunch–Kaufman algorithm

For solving a symmetric indefinite linear system, LAPACK implements a partitioned LDL$^{\text{T}}$factorization with the Bunch–Kaufman algorithm (Bunch and Kaufman Reference Bunch and Kaufman1977, Anderson and Dongarra Reference Anderson and Dongarra1989) that computes

(7.1)$$\begin{eqnarray}PAP^{T}=LDL^{T},\end{eqnarray}$$

where $D$ is a block-diagonal matrix with either $1\times 1$ or $2\times 2$ diagonal blocks. This algorithm is backward stable, subject to growth factors (Ashcraft, Grimes and Lewis Reference Ashcraft, Grimes and Lewis1998), and performs about $\frac{1}{3}n^{3}+O(n^{2})$ flops; see Algorithm 5.

To select the pivot at each step of the panel factorization, the algorithm scans two columns of the trailing submatrix, and depending on the numerical values of the scanned matrix entries, it uses either a $1\times 1$ or $2\times 2$ pivot. This results in synchronizations and irregular data access (due to the symmetric matrix storage), which have become expensive on modern computers. As a result, in many cases, the implementation of the algorithm cannot obtain high performance or exploit parallelism.

7.1.2 PLASMA: random butterfly transformation

PLASMA provides a set of dense linear algebra routines based on tiled algorithms that break a given algorithm into fine-grained computational tasks that operate on small square submatrices called tiles. Since each tile is stored contiguously in memory and fits in a local cache memory, this algorithm can take advantage of the hierarchical memory architecture used in modern computing hardware. Furthermore, by dynamically scheduling the tasks as soon as all of their dependencies are resolved, PLASMA can exploit a fine-grained parallelism and utilize a large number of cores.

Randomized algorithms are gaining popularity in linear algebra since they can often exploit more parallelism than corresponding deterministic algorithms (Parker Reference Parker1995a). To solve a symmetric indefinite linear system, PLASMAextends a randomization technique developed for the LU factorization (Baboulin et al. Reference Baboulin, Dongarra, Herrmann and Tomov2013) to symmetric indefinite systems (Becker, Baboulin and Dongarra Reference Becker, Baboulin, Dongarra, Wyrzykowski, Dongarra, Karczewski and Wasniewski2012, Baboulin, Becker and Dongarra Reference Baboulin, Becker and Dongarra2012). That is, it uses a multiplicative preconditioner by means of random matrices called recursive butterfly matrices $U_{k}$:

$$\begin{eqnarray}(U_{d}^{T}U_{d-1}^{T}\ldots U_{1}^{T})A(U_{1}U_{2}\ldots U_{d}),\end{eqnarray}$$

where

$$\begin{eqnarray}U_{k}=\left(\begin{array}{@{}ccc@{}}B_{1} & & \\ & \ddots & \\ & & B_{2^{k-1}}\end{array}\right)\!,\quad B_{k}=\frac{1}{\sqrt{2}}\left(\begin{array}{@{}cc@{}}R_{k} & S_{k}\\ R_{k} & -S_{k}\end{array}\right)\!,\end{eqnarray}$$

and $R_{k}$ and $S_{k}$ are random diagonal matrices. As a result, the original matrix $A$ is transformed into a matrix that is sufficiently random so that, with probability close to one, pivoting is not needed.

This random butterfly transformation (RBT) requires only $2dn^{2}+O(n)$ flops compared to the $\frac{1}{3}n^{3}+O(n^{2})$ flops required for the factorization. In practice, $d=2$ achieves satisfying accuracy, and this is the default setup used in PLASMA. Since $A$ is factorized without pivoting after RBT, it allows a scalable factorization of a dense symmetric indefinite matrix.

The main drawback of this method is reliability. There is no theory demonstrating its deterministic stability. Though it may fail in certain cases, numerical tests showed that with iterative refinement, it is reliable for many test cases, including pathological ones (Becker et al. Reference Becker, Baboulin, Dongarra, Wyrzykowski, Dongarra, Karczewski and Wasniewski2012). Furthermore, since iterative refinement is required, the failure of the method can be signalled without extra computation.

7.1.3 PLASMA: Aasen’s algorithm

To solve a symmetric indefinite linear system, Aasen’s algorithm (Aasen Reference Aasen1971) factorizes $A$ into an LTL$^{\text{T}}$decomposition of the form

(7.2)$$\begin{eqnarray}PAP^{T}=LTL^{T},\end{eqnarray}$$

where $P$ is a permutation matrix, $L$ is a unit lower triangular matrix, and $T$ is a symmetric tridiagonal matrix; see Algorithm 6. Aasen’s algorithm takes advantage of the symmetry in $A$ and performs $\frac{1}{3}n^{3}+O(n^{2})$ flops, which is a half of the flops needed for an LU factorization of $A$, and the same as the number of the flops required by the Bunch–Kaufman algorithm. Furthermore, it is backward stable and subject to a growth factor. Once the factorization is computed, the solution $x$ is computed by successively solving the linear systems with the matrices $L$, $T$, and $L^{T}$.

To exploit the memory hierarchy on a modern computer, a partitioned version of Aasen’s algorithm was proposed (Rozloz̆ník, Shklarski and Toledo Reference Rozloz̆ník, Shklarski and Toledo2011). This algorithm first factorizes a panel in a left-looking fashion, and then uses Level 3 BLAS operations to update the trailing submatrix in a right-looking way. Compared to a standard column-wise algorithm, this partitioned algorithm slightly increases the operation count, performing

$$\begin{eqnarray}\frac{1}{3}\biggl(1+\frac{1}{n_{b}}\biggr)n^{3}+O(n^{2}n_{b})\end{eqnarray}$$

flops with a block size of $n_{b}$. However, Level 3 BLAS can be used to perform most of these flops; i.e.,

$$\begin{eqnarray}\frac{1}{3}\biggl(1+\frac{1}{n_{b}}\biggr)n^{3}\end{eqnarray}$$

flops. Since the Level 3 BLAS operations have higher ratios of flop counts over communication volumes than the Level 2 BLAS or Level 1 BLAS operations do, this partitioned algorithm is shown to significantly shorten the factorization time on modern computers, where data transfer is much more expensive than floating-point operations (Rozloz̆ník et al. Reference Rozloz̆ník, Shklarski and Toledo2011). A serial implementation of the partitioned Aasen’s algorithm is shown to be as efficient as the Bunch–Kaufman algorithm of LAPACK on a single core (Rozloz̆ník et al. Reference Rozloz̆ník, Shklarski and Toledo2011). However, the panel factorization is still based on Level 1 BLAS and Level 2 BLAS operations. As a result, this panel factorization often obtains only a small fraction of the peak performance on modern computers, and could become the bottleneck, especially in a parallel implementation.

The blocked version of Aasen’s algorithm was proposed to avoid this bottleneck at the panel factorization (Ballard et al. Reference Ballard, Becker, Demmel, Dongarra, Druinsky, Peled, Schwartz, Toledo and Yamazaki2014). It computes an LTL$^{\text{T}}$factorization of $A$, where $T$ is a banded matrix (instead of tridiagonal) with its half-bandwidth being equal to the block size $n_{b}$, and then uses a banded matrix solver to compute the solution. In this blocked algorithm, each panel can be factorized using an existing LU factorization algorithm, such as recursive LU (Castaldo and Whaley Reference Castaldo and Whaley2010, Gustavson Reference Gustavson1997, Dongarra, Faverge, Ltaief and Luszczek Reference Dongarra, Faverge, Ltaief and Luszczek2011, Toledo Reference Toledo1997) and communication-avoiding LU (CALU) (Grigori et al. Reference Grigori, Demmel and Xiang2011, Demmel, Grigori, Hoemmen and Langou Reference Demmel, Grigori, Hoemmen and Langou2012). In comparison with the panel factorization algorithm used in the partitioned Aasen’s algorithm, these LU factorization algorithms reduce communication, and hence are expected to speed up the whole factorization process.

Figure 7.1.

Residual norms of blocked Aasen’s algorithm on a random matrix.

Figure 7.2.

Parallel performance on random matrix with $n_{b}=250$ on up to eight $6$-core $2.8$ MHz AMD Opteron processors.

We implemented this blocked Aasen’s algorithm on multicore architectures, and analysed its parallel performance (Ballard et al. Reference Ballard, Becker, Demmel, Dongarra, Druinsky, Peled, Schwartz, Toledo and Yamazaki2013). Our implementation follows the framework of PLASMA and uses a dynamic scheduler called QUARK. To efficiently utilize a large number of cores in parallel, all the existing factorization routines in PLASMA update the trailing submatrix in a right-looking fashion. Hence, our implementation was not only the first parallel implementation of the blocked Aasen’s algorithm, but it was also the first implementation of a left-looking algorithm in PLASMA. In order to fully exploit the limited parallelism in this left-looking algorithm, we studied several performance enhancing techniques (e.g. parallel reduction to update the panel, recursive LU and CALU for panel factorization, and parallel symmetric pivoting). Our performance results on up to 48 AMD Opteron processors demonstrate that a left-looking algorithm can be implemented efficiently on a multicore architecture. In addition, our numerical results show that, in comparison with the widely used stable algorithm (the Bunch–Kaufman algorithm of LAPACK), our implementation loses only one or two digits in the computed residual norms when a recursive LU is used on a panel: see Figure 7.1. Figure 7.2 compares the performance of different algorithms.

7.2.1 Bunch–Kaufman algorithm

Our first implementation of the Bunch–Kaufman algorithm is based on a hybrid CPU/GPU programming paradigm where the block column (commonly referred to as a panel) is factorized on the CPUs (e.g. using the threaded MKL), while the trailing submatrix is updated on the GPU. This is often an effective programming paradigm for many of the LAPACK subroutines because the panel factorization is based on Level 1 BLAS or Level 2 BLAS, which can be efficiently implemented on the CPUs, while Level 3 BLAS is used for the submatrix updates, which exhibit high data parallelism and can be efficiently implemented on the GPU. Unfortunately, at each step of the panel factorization, the Bunch–Kaufman algorithm may select the pivot from the trailing submatrix. Hence, although copying the panel from the GPU to the CPU can be overlapped with the update of the rest of the trailing submatrix on the GPU, the look-ahead – a standard optimization technique to overlap the panel factorization on the CPUs with the trailing submatrix update on the GPU – is prohibited. In addition, when the pivot column is on the GPU, this leads to the expensive data transfer between the GPU and the CPU at each step of the factorization. To avoid this expensive data transfer, our second implementation performs the entire factorization on the GPU. Though the CPUs may be more efficient performing the Level 1 BLAS and Level 2 BLAS based panel factorization, this implementation often obtains higher performance by avoiding the expensive data transfer.

When the entire factorization is implemented on the GPU, to select a pivot at each step of the Bunch–Kaufman algorithm, up to two columns of the trailing submatrix must be scanned – the current column and the one with an index corresponding to the row index of the element with the maximum modulus in the first column. This not only leads to the expensive global reduce on the GPU, but also to irregular data accesses since only the lower or upper triangular part of the submatrix is stored. This makes it difficult to obtain high performance on the GPU. In the next two sections, we describe two other algorithms (i.e. communication-avoiding and randomization algorithms) that aim to reduce this bottleneck.

7.2.2 Random butterfly transformation

The application of a depth 1 butterfly matrix is performed using a CUDA kernel where the computed part of the matrix $A$ is split into blocks. For each of these blocks, the corresponding part of the matrix $U$ is stored in the shared memory to improve the memory access performance. Matrix $U$ is small enough to fit into the shared memory due to its packed storage.

To compute the LDL$^{\text{T}}$ factorization of the randomized matrix without pivoting, we implemented a block variant of no-pivoting LDL factorization on multicore CPUs with a GPU. In our implementation, the matrix is first copied to the GPU, then the CPU is used to compute the LDL$^{\text{T}}$ factorization of the diagonal block. Once the resulting LDL$^{\text{T}}$ factors of the diagonal block are copied back to the GPU, the corresponding off-diagonal blocks of the $L$-factor are computed by the triangular solve on the GPU. Finally, we update each block column of the trailing submatrix, calling a matrix–matrix multiply on the GPU.

Figure 7.3.

Performance of symmetric factorization on two 8-core Intel Xeon X5680 CPUs with an NVIDIA Kepler K20 GPU (double precision).

7.2.3 Aasen’s algorithm

Though the GPU has a greater memory bandwidth than the CPU, the memory accesses are still expensive compared to the arithmetic operations. Hence, our implementation is based on the CA Aasen’s algorithm. Though this algorithm performs most of the flops using Level 3 BLAS, most of the operations are on the submatrices of the block size $n_{b}$. In order to exploit parallelism between the small BLAS calls, we use GPU streams extensively. In addition, for an efficient application of the symmetric pivots after the panel factorization, we apply the pivots in two steps. The first step copies all the columns of the trailing submatrix, which need to be swapped, into an $n\times 2n_{b}$ workspace. Here, because of the symmetry, the $k$th block column consists of the blocks in the $k$th block row and those in the $k$th block column. Then, in the second step, we copy the columns of the workspace back to a block column of the submatrix after the column pivoting is applied. The same pivoting strategy is used to exploit the parallelism on multicore CPUs (Ballard et al. Reference Ballard, Becker, Demmel, Dongarra, Druinsky, Peled, Schwartz, Toledo and Yamazaki2013). We tested using the LU factorization with partial pivoting as the panel factorization, either using the threaded MKL on the CPUs or using its native GPU implementation of MAGMA on GPU. Though the Level 1 BLAS and Level 2 BLAS based panel factorization may be more efficient on the CPUs, the second approach avoids the expensive data transfer to copy the panel from the GPU to the CPU. Figure 7.3 compares the performance of different factorization algorithms.

8.2.1 The classical reduction to bidiagonal condensed form

Due to its high computational complexity of $O({\textstyle \frac{8}{3}}n^{3})$ (for square matrices) and interdependent data access patterns, the bidiagonal reduction phase is the most challenging stage to develop and optimize – both algorithmically and from the implementation standpoint. The classical approach from LAPACK (Anderson et al. Reference Anderson, Bai, Bischof, Blackford, Demmel, Dongarra, Ducroz, Greenbaum, Hammarling, McKenney and Sorensen1999) is denoted by ‘one-stage algorithm’, whereby the Householder transformations are grouped and applied in a blocked fashion to directly reduce the dense matrix to bidiagonal form. The one-stage reduction to bidiagonal form suffers from lack of efficiency. To understand it better, one has to focus on the two computational procedures (panel and update) that are repeated until the matrix is reduced. The reduction proceeds by steps of size $n_{b}$. Below we give the detailed cost of step $i$, as a cost for the panel and a cost for the update.

• ∙ The panel factorization computes the similarity transformations (Householder reflectors) to introduce zeros to the entries below the subdiagonal within a single block of $n_{b}$ columns. The factorization of every column is dominated primarily by two matrix–vector products with the trailing matrix. Thus, this computation is critical as the entire trailing submatrix needs to be loaded into memory, and very few floating-point operations are executed on all of the transferred data. As the memory bandwidth becomes more limited and does not scale with the number of cores, the panel factorization step is not expected to scale for large matrices that do not fit in cache. The cost of a panel is $4n_{b}l^{2}+{\rm\Theta}(n)$, where $l$ is the size of the trailing matrix. For simplicity, we omit ${\rm\Theta}(n)$ and round up the cost of the panel by the cost of the matrix–vector product.

• ∙ The trailing submatrix update consists of applying the Householder reflectors, generated during the panel factorization, to the trailing matrix from both the left and right side using compute-bound operations that utilize Level 3 BLAS, which have a high data re-use ratio allowing for highly tuned implementations. In addition, Level 3 BLAS lend themselves to parallelization due to their inherent data locality properties, which can be exploited to achieve high performance rates for large amounts of operations, with most of it completely independent and thus perfectly suited to multicore processors. Unfortunately, each panel factorization must be synchronized with the corresponding update of the trailing submatrix, which prevents asynchronous execution and overlap of memory-bound and compute-bound steps that could potentially alleviate the effects of the former. The computational cost of this procedure is

  $$\begin{eqnarray}A_{i+n_{b}:n,i+n_{b}:n}\leftarrow A_{i+n_{b}:n,i+n_{b}:n}-V\times Y^{T}-X\times U^{T},\end{eqnarray}$$
  where $V$ and $Y$ are the Householder reflectors computed during the panel phase, $X$ and $Y$ are two rectangular matrices needed for the update, and also computed during the panel phase. This update phase can be performed by two matrix–matrix products using the GEMM routine, and its cost is $2\times 2n_{b}k^{2}$, where $k$ is the size of the trailing matrix at step $i$.

Figure 8.1.

The percentage of the time spent in each phase of the SVD solver using the classical one-stage approach to compute the bidiagonal form.

For all steps ($n/n_{b}$), the trailing matrix size varies from $n$ to $n_{b}$ by steps of size $n_{b}$, where $l$ varies from $n$ to $n_{b}$ and $k$ varies from $(n-n_{b})$ to $2n_{b}$. Thus the total cost for the $n/n_{b}$ steps is

(8.1)$$\begin{eqnarray}{\approx}4n_{b}\mathop{\sum }_{n_{b}}^{n/n_{b}}l^{2}+4n_{b}\mathop{\sum }_{2n_{b}}^{\frac{n-n_{b}}{n_{b}}}k^{2}\approx \frac{4}{3}n_{\mathtt{GEMV}}^{3}+\frac{4}{3}n_{\text{Level 3}}^{3}\approx \frac{8}{3}n^{3}.\end{eqnarray}$$

Following the equations above, we derive the maximum performance $P_{\max }$ that can be reached by any of these reduction algorithms. In particular, for large matrix sizes $n$,

$$\begin{eqnarray}P_{\max }=\frac{\text{number of operations}}{\text{minimum time }t_{\min }},\end{eqnarray}$$

and thus $P_{\max }$ is expressed as

(8.2)$$\begin{eqnarray}\frac{\frac{8}{3}n^{3}}{\frac{4}{3}n^{3}\ast \frac{1}{P_{\mathtt{GEMV}}}+\frac{4}{3}n^{3}\ast \frac{1}{P_{\text{Level 3}}}}=\frac{2\ast P_{\text{Level 3}}\ast P_{\mathtt{GEMV}}}{P_{\text{Level 3}}+P_{\mathtt{GEMV}}}\approx 2P_{\mathtt{GEMV}}\end{eqnarray}$$

when $P_{\text{Level 3}}\gg P_{\mathtt{GEMV}}$. Thus, equation (8.2) proves the well-known low-performance behaviour of the classical bidiagonal reduction algorithm, and for that it is considered the most time-consuming phase. Figure 8.1 shows the percentage of the total time for each of the four components of the SVD solver using the standard one-stage reduction approach when all the singular vectors are computed. The figure makes it clear that the reduction to the bidiagonal form requires more than 70% of the total time for the case when all the singular vectors are computed. In addition, when only singular values are needed, the reduction requires about 90% of the total computing time. Because of the expense of the reduction step, renewed research has focused on improving this step, resulting in a novel technique based on a ‘two-stage’ algorithm (Lang Reference Lang1999, Bientinesi et al. Reference Bientinesi, Igual, Kressner and Quintana-Ortí2010, Haidar et al. Reference Haidar, Ltaief and Dongarra2011, Haidar et al. Reference Haidar, Ltaief, Luszczek and Dongarra2012, Ltaief et al. Reference Ltaief, Luszczek, Haidar, Dongarra, Wyrzykowski, Dongarra, Karczewski and Wasniewski2012, Haidar, Kurzak and Luszczek Reference Haidar, Kurzak and Luszczek2013a), where a first stage uses Level 3 BLAS operations to reduce $A_{s}$ to band form, and a second stage further reduces the matrix to the proper bidiagonal form. We developed our singular value solver based on the ‘two-stage’ techniques, which allow us to exploit more parallelism and use recent hardware more efficiently.

Figure 8.2.

Two-stage technique for the reduction phase.

8.2.2 Two-stage reduction

The two-stage reduction is designed to overcome the limitations of the one-stage approach, which relies heavily on memory-bound operations. It also increases the use of compute-intensive operations that benefit from the increase in CPU core count. Many algorithms have been studied extensively in the context of the symmetric eigenvalue problems (Bischof, Lang and Sun Reference Bischof, Lang and Sun2000, Bischof and Van Loan Reference Bischof and Van Loan1987), and, more recently, tile algorithms have achieved good performance (Luszczek, Ltaief and Dongarra Reference Luszczek, Ltaief and Dongarra2011, Haidar et al. Reference Haidar, Ltaief and Dongarra2011). The idea behind them is to split the original one-stage approach into a compute-intensive phase (first stage) and a memory-bound phase (second or bulge chasing stage) as represented in Figure 8.2. The first stage reduces the original general dense matrix to a band form (either upper or lower), and the second stage reduces the band form to the canonical bidiagonal form (again, either upper or lower). The two-stage approach is used in our implementation, and exhibits some similarities to what has been developed for the symmetric eigenvalue problem (Haidar et al. Reference Haidar, Ltaief and Dongarra2011). To put our work into a proper perspective, we start by briefly describing the first stage (reduction from full dense to band), and then we explain in more detail the reduction from band to bidiagonal form. We will also talk about the scheduling techniques on which our implementation relies.

First stage: compute-intensive and efficient kernels.

The first stage applies a sequence of blocked Householder transformations to reduce the general dense matrix to either an upper or a lower band matrix. This stage uses compute-intensive matrix-multiply kernels that eliminate the memory-bound matrix–vector product from the one-stage panel factorization. It also illustrates a beneficial data access pattern through the use of compute-intensive operations based on Level 3 BLAS for large portions of the code (Bischof et al. Reference Bischof, Lang and Sun2000, Dongarra, Sorensen and Hammarling Reference Dongarra, Sorensen and Hammarling1989, Gansterer, Kvasnicka and Ueberhuber Reference Gansterer, Kvasnicka and Ueberhuber1999, Haidar et al. Reference Haidar, Tomov, Dongarra, Solcà and T.2014). Moreover, this stage is made highly parallel because of the tile algorithm formulation (Agullo, Hadri, Ltaief and Dongarra Reference Agullo, Hadri, Ltaief and Dongarra2009b), which brings the parallelism to the fore. Conceptually and physically, the matrix is split into $nt\times nt$ tiles ($nt=n/nb$, $nb$ is the size of the tile and the resulting matrix bandwidth). The data within a tile is stored contiguously in memory. The algorithm then proceeds as a collection of interdependent tasks that operate on the tile data layout. Figure 8.3 highlights the execution breakdown during the second step of the first stage of the reduction. A QR factorization is computed for the tile $A_{2,2}$ (the solid grey tile). After this QR is finished, a set of independent tasks is released and they can all be executed in parallel. All tiles $A_{2,\bullet }$ (the black tiles of Figure 8.3(a)) can be updated by applying the Householder transformations that are generated by the QR factorization of $A_{2,2}$. Also, all the tiles $A_{\bullet ,2}$ (the grey diamond pattern tiles of Figure 8.3(a)) can also be independently annihilated one after another with the R factor of $A_{2,2}$. After each tile $A_{i,2}$ is annihilated (e.g. the black checkerboard pattern tile of Figure 8.3(a)), a set of parallel tasks may be launched to update all the tiles of the block row $i$ (the stripe pattern tiles of Figure 8.3(a)). Moreover, when $A_{2,3}$ is updated, then an LQ factorization is performed for this tile (the solid grey tile of Figure 8.3(b)). Similarly to the QR process, after LQ, all the tiles in the third column of tiles, $A_{3:nt,3}$ (the grey diamond pattern tiles of Figure 8.3(b)), can now be independently updated by the Householder vectors from the LQ factorization, provided that they have been updated with the transformation from the QR factorization. Similarly, all the tiles $A_{2,4:nt}$ (the black checkerboard pattern tiles of Figure 8.3(b)) can also be annihilated. Likewise, each annihilation of $A_{1,i}$ (the black checkerboard pattern tile of Figure 8.3(b)) enables a set of tasks to update the block column $i$ (the stripe pattern tiles of Figure 8.3(b)). We note that this requires implementations of new computational kernels to be able to operate on the new data structures. The details of the implementation of this stage are provided elsewhere (Haidar et al. Reference Haidar, Ltaief and Dongarra2011, Luszczek et al. Reference Luszczek, Ltaief and Dongarra2011). This process of interleaving the QR and the LQ factorizations at each step repeats until the end, and, as a result, we obtain a band matrix with a bandwidth of size $nb$. As we mentioned above, the tile formulation of the algorithm resulted in the creation of a large number of parallel tasks. These tasks are organized into a directed acyclic graph (DAG) (Buttari et al. Reference Buttari, Dongarra, Kurzak, Langou, Luszczek, Tomov, Kågström, Elmroth, Dongarra and Wasniewski2006, Chan, Quintana-Orti, Quintana-Orti and van de Geijn Reference Chan, Quintana-Orti, Quintana-Orti and van de Geijn2007), with the nodes representing the computational tasks and the edges – the data dependencies between them. Thus, restructuring the linear algebra algorithms as a sequence of tasks that operate on tiles of data removes the fork–join bottleneck that befits the LAPACK-style implementations. It also avoids idle time for individual cores, while, at the same time, increasing the data locality for each core. Moreover, in order to increase the efficiency of our algorithm, we improved our dynamic scheduler by developing two main features that played an important role in targeting high-performance execution rates. Two auxiliary options were added to the scheduler that allow us to label some of the tasks PRIORITY and LOCALITY. For example, when the QR factorization of $A_{1,1}$ and a sequence of QR updates ensues, it is better to increase the priority of all the tasks that modify the tile $A_{1,2}$, in such a way that the LQ process of $A_{1,2}$ starts as soon as possible. As a result, this technique will increase the number of parallel tasks and will aid the interleaving of the tasks from both the QR and the LQ update sequences. As described below, this has a big impact on data locality. The second task label that we developed is the LOCALITY flag. It is used for the update of any tile of $A_{2:n,2:n}$. It makes it possible for $A_{2,2}$ to be updated by the QR’s Householder process and by the LQ Householder process. Hence, our scheduler has an opportunity to let the same core update $A_{2,2}$ by the two transformations, one after the other. This will result in the tile data being loaded from the main memory only once.

Figure 8.3.

Kernel execution of the BRD algorithm during the first stage.

Second stage: cache-friendly computational kernels.

The band form is further reduced to the final condensed form using the bulge chasing technique. This procedure annihilates the extra off-diagonal elements by chasing the created fill-in elements down to the bottom right side of the matrix using successive orthogonal transformations at each sweep. This stage involves memory-bound operations and requires the band matrix to be accessed from multiple disjoint locations. In other words, there is an accumulation of substantial latency overhead each time different portions of the matrix are loaded into cache memory, which is not compensated for by the low execution rate of the actual computations (the so-called surface-to-volume effect). To overcome these critical limitations, we developed a bulge chasing algorithm, very similar to the novel bulge chasing techniques for symmetric eigenvalue problems (reduction from band to bidiagonal) developed in Haidar et al. (Reference Haidar, Ltaief and Dongarra2011), but we differ from it in using a column-wise elimination instead of an element-wise elimination. In addition, we also differ by developing our kernels to deal with general matrices instead of symmetric matrices. When the singular vectors need to be computed, the most problematic aspect of the standard procedure is the element-wise elimination (Haidar et al. Reference Haidar, Ltaief and Dongarra2011). Such an implementation is very suitable when only singular value is required, but is limited when singular vectors are required. In particular, it generates element-wise Householder reflectors, and thus the update of the singular vectors by these reflectors becomes the bottleneck as it is based on Level 1 BLAS operations. Our modification adds a small amount of extra work, but it allows the use of the Level 3 BLAS kernels to compute the transformations or to apply them in the form of the orthogonal matrix $U_{2}$ and $V_{2}$ – the result of computation in this phase. Moreover, we designed our algorithm to extensively use cache-friendly kernels combined with fine-grained, memory-aware tasks in an out-of-order scheduling technique which considerably enhances data locality.

Figure 8.4.

Kernel execution of the BRD algorithm during the second stage.

The bulge chasing algorithm consists of a succession of three new kernels designed to increase cache re-use. The idea is to load a block of data in the actual cache memory and to apply all the possible computation to it before unloading it. The first kernel, called xGBCW1 and illustrated in Figure 8.4(a), manipulates the black checkerboard pattern block of data. It triggers the beginning of each sweep by annihilating the extra non-zero entries within a single row, then applies the computed elementary Householder reflector from the right within this block. Hence, it subsequently generates triangular bulges as shown in Figure 8.4(a) (the black block). Note that this triangular bulge must be annihilated eventually in order to avoid the excessive growth of the fill-in structure. A classical implementation will eliminate the whole triangular bulge. However, for an appropriate study of the bulge chasing procedure, let us remark that the elimination of the row $i+1$ (the sweep $i+1$), at the next step, creates a triangular bulge which will overlap this one-by-one column shift to the right and one row to the bottom, as shown in Figure 8.4(d), where the reader can see that the lower triangular portion of the grey checkerboard pattern block (the bulge created in sweep $i+1$) overlaps with the lower triangular portion of the black checkerboard pattern block (corresponding to the bulges created by the previous sweep $i$). As a result, we can reduce the computational cost, and instead of eliminating the whole triangular bulge created for sweep $i$, we only eliminate the non-overlapped region of it – its first column. The remaining columns can be delayed to the upcoming annihilation sweeps. In this way, we can avoid the growth of the bulges and reduce the extra cost accrued when the whole bulge is eliminated. Moreover, we designed a cache-friendly kernel that takes advantage of the fact that the created bulge (the black checkerboard pattern block) remains in the cache and therefore it directly eliminates its first column and applies the corresponding left update to the remaining column of the black checkerboard pattern block. The second kernel, xGBCW2, loads the next block and it applies the necessary left updates derived from the previous kernel. This will eventually generate triangular bulges as shown in Figure 8.4(b). Hence, this kernel will also annihilate the first row of the created bulge and update its black stripe pattern block from the right. Finally, the third kernel, xGBCW3, loads the next block (the third black checkerboard pattern block of Figure 8.4(c)) and it continues applying, from the right, the transformations of the previous kernel 2. Like kernel 1, it also creates a bulge which is removed, and the black checkerboard pattern block is updated correspondingly from the left, similar to the process undertaken by kernel 1. Accordingly, the annihilation of each sweep can be described as a single call to kernel 1 followed by repetitive calls to a cycle of kernel 2 and kernel 3.

The implementation of this stage is done by using either a dynamic or a static runtime environment, whose detailed implementation we developed and discussed in Haidar et al. (Reference Haidar, Kurzak and Luszczek2013a, Reference Haidar, Tomov, Dongarra, Solcà and T.2014). This stage is, in our opinion, one of the main challenges for algorithms as it is difficult to track the data dependencies. The annihilation of the subsequent sweeps will generate computational tasks, which will partially overlap within the data used by the tasks of the previous sweeps (see Figure 8.4(d)) – the main challenge of dependency tracking. We have used our data translation layer (DTL) and functional dependencies (Luszczek et al. Reference Luszczek, Ltaief and Dongarra2011, Haidar et al. Reference Haidar, Ltaief and Dongarra2011) to handle the dependencies and to provide crucial information to the runtime to achieve the correct scheduling. As mentioned above, the amount of data involved with each task is of size $nb\times nb$, but is handled efficiently because our kernels increase cache re-use. We also developed our scheduling technique to minimize the memory traffic and increase the efficiency of our algorithm in such a way that the subsequent tasks that involve the same region of data will be executed by the same thread. As an example, the first task $T_{1}^{(2)}$ of the annihilation of sweep 2, which operates on the grey stripe pattern block of Figure 8.4(d), overlaps the region of the data that has been used by the first task $T_{1}^{(1)}$ of the sweep 1 (black checkerboard pattern block). However, it is obvious to execute task $T_{1}^{(2)}$ by the same thread that executed $T_{1}^{(1)}$. To ensure maximum re-use, we force the scheduler to distribute the tasks according to their data location. The main goal of our data-based scheduling is to define fine-grained local tasks that increase cache re-use and minimize communication. For small matrix sizes, we prefer the use of a subsequent set of threads (the number of threads that fit the matrix into its fast memory) while leaving the other threads working on different portions of the code rather than using all the available resources. The implementation of this phase has been well optimized. It has been observed that it takes between 5% and 10% of the global time of the reduction from dense to bidiagonal.

8.2.3 The bidiagonal singular solver

A bidiagonal singular solver computes the spectral decomposition of a bidiagonal matrix $B$ such that

(8.3)$$\begin{eqnarray}B=\widetilde{U}{\rm\Sigma}\widetilde{V}^{H}\quad \text{with }\widetilde{U}\widetilde{U}^{H}=I\text{ and }\widetilde{V}\widetilde{V}^{H}=I,\end{eqnarray}$$

where $\widetilde{U}$ and $\widetilde{V}^{H}$ are the singular vectors and ${\rm\Sigma}$ the singular values of $B$. There are two main implementations of the singular solver, one based on the QR algorithm (Golub and Kahan Reference Golub and Kahan1965), and one based on the divide and conquer algorithm (Jessup and Sorensen Reference Jessup and Sorensen1989, Gu and Eisenstat Reference Gu and Eisenstat1995).

8.2.4 The singular vector computation

The classical one-stage approach reduces the dense matrix $A$ to bidiagonal $B$ and computes its singular values and vectors. The singular values of $A$ are the same as the ${\rm\Sigma}$ computed for the bidiagonal matrix $B$. The singular vectors of $A$ ($U$ and $V^{H}$) are computed from the left and right singular vectors of $B$ by applying the same orthogonal matrices $Q$ and $P^{H}$ that were used in the reduction to condensed form, respectively, that is,

(8.4)$$\begin{eqnarray}\displaystyle \begin{array}{@{}rclcl@{}}U\ & =\ & Q\widetilde{U}\ & =\ & (I-V_{1}T_{1}V_{1}^{H})\widetilde{U},\\ \text{and}\quad V^{H}\ & =\ & \widetilde{V}^{H}P_{1}^{H}\ & =\ & \widetilde{V}^{H}(I-W_{1}Tr_{1}^{H}W_{1}^{H}),\end{array} & & \displaystyle\end{eqnarray}$$

where ($V_{1},T_{1}$ and $W_{1},Tr_{1}$) represent the left and right Householder reflectors generated during the reduction to the bidiagonal form, respectively. From this representation, either $Q$ and $P^{H}$ can be formed explicitly using DORGBR; or we can multiply by $Q$$P^{H}$ in an implicit fashion using DORMBR, which is less expensive. In either case, applying $Q$ and $P^{H}$ becomes a series of DGEMM operations. Algorithm 7 describes the back transformation of the one-stage reduction to the left singular vectors. The same will be applied to the right vectors. The cost of this phase is equal to $2n^{3}$ for each singular vector. Since all of its computation is based on Level 3 BLAS, the performance upper bound of this phase is considered to be on a par with the performance of the Level 3 BLAS, which is usually a large fraction of the peak of the machine.

In the case of the two-stage approach, the first stage reduces the original general dense matrix $A$ to a general band matrix by applying a two-sided transformation to $A$ such that $Q_{1}^{H}AP_{1}=A_{\text{band}}$. Similarly, the second stage – bulge chasing – reduces the band matrix $A_{\text{band}}$ to the bidiagonal form by applying the transformation from both the left and right side to $A_{\text{band}}$ such that $Q_{2}^{H}A_{\text{band}}P_{2}=B$. Let us now denote the SVD of $A$ by $A=U{\rm\Sigma}V^{H}$, where $U$ and $V^{H}$ are, respectively, the left and right singular vectors of $A$, and the diagonal entries of ${\rm\Sigma}$ hold the singular values. These singular vectors must be multiplied by both of $Q_{\ast }$ and $P_{\ast }$, according to

(8.5)$$\begin{eqnarray}\begin{array}{@{}rclcl@{}}U\ & =\ & Q_{1}Q_{2}\widetilde{U}\ & =\ & (I-V_{1}T_{1}V_{1}^{H})(I-V_{2}T_{2}V_{2}^{H})\widetilde{U},\\ \text{and}\quad V^{H}\ & =\ & \widetilde{V}^{H}P_{2}^{H}P_{1}^{H}\ & =\ & \widetilde{V}^{T}(I-W_{2}Tr_{2}^{H}W_{2}^{H})(I-W_{1}Tr_{1}^{H}W_{1}^{H}),\end{array}\end{eqnarray}$$

where ($V_{1},T_{1}$ and $W_{1},Tr_{1}$) and ($V_{2},T_{2}$ and $W_{2},Tr_{2}$) represent the left and right Householder reflectors generated during the first and second stages of the reduction to the bidiagonal form. It is clear that the two-stage approach introduces a non-trivial amount of extra computation – the application of $Q_{2}$ and $P_{2}^{H}$ – for the case when the singular vectors are needed. Thus, one of the crucial procedures of the two-stage algorithm, when optimizing for performance, is the update of the singular vectors by the Householder transformations that were generated during the two stages of the reduction to the bidiagonal form. Obviously, the implementation of this scheme is not as straightforward as simply parallelizing a loop. In particular, because of complications of the bulge chasing mechanism, the order of generated task dependencies is quite intricate. Due to the geometric shape of non-zero entries during the second stage of the reduction from $A_{\text{band}}$ to bidiagonal, called bulge chasing, the application of $Q_{2}$ and $P_{2}$ has to follow the order mentioned in Figure 8.5, where each application consists of the diamond shape of the $V$ applied to $U$. Figure 8.5 shows the Householder reflectors for both stages. More details of the optimized technique to apply $Q_{2}$ and $P_{2}$ can be found in Haidar, Luszczek and Dongarra (Reference Haidar, Luszczek and Dongarra2014b). The algorithm for the update of $U$ using $Q_{2}$ is described in Algorithm 7, and the cost of this step can be summarized by a

$$\begin{eqnarray}2\biggl(1+\frac{i_{b}}{n_{b}}\biggr)n^{3}\end{eqnarray}$$

operation using the Level 3 BLAS. Similarly, the cost of applying $P_{2}^{H}$ can be computed in the same manner.

Figure 8.5.

(a) Tiling of $V_{1}$. (b) Blocking used to apply $V_{2}$.

The application of $Q_{1}$ can be done easily and efficiently using our tile algorithm (Haidar et al. Reference Haidar, Luszczek and Dongarra2014b) and is described in Algorithm 9. The parallelism in this step comes from two sources: first, the singular vector matrix is viewed as a set of independent tiles to be updated; second, the parallelism can also be extracted by applying $Q_{1}$ as a set of independent Householder reflectors, row-wise, meaning that the only constraint to follow is their column-wise order (e.g., once the Householder reflector of block$(i,j)$ has been applied, that of block$(i,j-1)$ can be applied). As a result, the design of the tile algorithm generates a large number of independent tasks that can be applied in an asynchronous manner using either a static or a dynamic scheduler. The $V_{1}$ are stored in a tile fashion, as shown in Figure 8.5, to increase data locality. The total cost of this step is approximately $2n^{3}$. Consequently, the total cost of updating the singular vectors when using the two-stage techniques is

$$\begin{eqnarray}2\biggl(1+\frac{i_{b}}{n_{b}}\biggr)n^{3}+2n^{3}\end{eqnarray}$$

for the left singular vectors and the same for the right singular vectors. Compared to the one-stage approach, we can observe that the two-stage technique has an extra cost of the order of

$$\begin{eqnarray}2\biggl(1+\frac{i_{b}}{n_{b}}\biggr)n^{3}\end{eqnarray}$$

operations for the left singular vector and the same for the right singular vector.

Hence, when only the singular values are computed, our new two-stage implementation is expected to be largely advantageous when compared to the classic algorithm (e.g., it can reach about $7\times$ speedup). However, when computing the left and right singular vectors, our two-stage approach has an extra cost of

$$\begin{eqnarray}2\times 2\biggl(1+\frac{i_{b}}{n_{b}}\biggr)n^{3}\end{eqnarray}$$

operations from the back-transformation of the second stage. This will clearly affect the roughly $8\times$ observed speedup for the reduction to the bidiagonal form. However, since the extra cost is computed by the Level 3 BLAS operations, one can expect that we can still get speedup over the classical approach – particularly for accelerated architectures. Our experiments described in the next section show that we can reach still reach up to $2\times$ speedup on a multicore architecture when both the left and right singular vectors are computed. Also, one can expect that if a percentage of the singular vectors is needed, the extra cost is dramatically decreased, and a speedup of about $5\times$ was observed.

8.2.5 Experimental results

This section presents the performance comparisons of our tile algorithm for two-stage SVD against state-of-the-art numerical linear algebra libraries.

Performance results.

The following experiments illustrate the superior efficiency and the scalability of our proposed SVD solver with respect to the state-of-the-art optimized vendor numerical linear algebra libraries. Each graph below presents comparison curves that we will describe. We performed an extensive study with a large number of experiments on two different machines (one with a large number of cores, system A, and another with small number of cores, system B) in order to give the reader as much information as possible. We computed the SVD decomposition, where we either compute only the singular values, or both singular values and vectors, with two different algorithms (DGESVD, DGESDD), varying the size of the matrices from $2000$ to $26\,000$ using our two systems described above. Since many applications only need a portion of the singular vectors, we also present results using DGESDD (the divide and conquer SVD solver) to compute 20% of the singular vectors. In addition, in order to make our experiments self-contained, we report the result of the improvements that our two-stage implementation brought to the reduction to bidiagonal form, compared against the one-stage approach from the state-of-the-art numerical linear algebra libraries.

Figure 8.6.

The speedup obtained by our implementation of DGESDD versus its counterpart from the Intel MKL library on 48 AMD cores of system A.

Figure 8.7.

The speedup obtained by our implementation of DGESDD versus its counterpart from the Intel MKL library on 16 Sandy Bridge cores of system B.

In particular, Figures 8.6 and 8.7 show speedups and efficiencies for computing the SVD using the divide and conquer technique. For each speedup curve (representing a routine), we depict the ratio of the runtime between the routine from Intel’s MKL library and its counterpart from our implementation within the same computing environment. These results show four types of behaviour.

Let us first comment on the reduction to bidiagonal form (DGEBRD: triangles) and on the SVD decomposition when only the singular values are computed (DGESDD NO Vectors: squares). The speedups shown are remarkable, our implementation asymptotically achieves more than $8\times$ speedup on the 48 cores of system A and more than $4\times$ speedup on the 16 cores of system B. This was not unexpected, and the results obtained here confirm the importance of the considerations discussed in our study when designing high-performance libraries. The performance gain is due to the tile algorithm that we developed and to its efficient implementation of the first stage (reduction to band), which is the compute-intensive stage, and from the design of both the second stage (bulge chasing) and the orthogonal transformation update, which maps both the algorithm and the data to the hardware using cache-friendly kernels and scheduling based on increasing data locality. Our new algorithm scales as the matrix sizes increase and asymptotically achieves a perfect speedup, despite the side effects of running on a NUMA system.

In contrast to the standard approach where the reduction to bidiagonal dominates the overall time (as described above in Figure 8.1), the two-stage reduction to bidiagonal consists of less than 20% of the overall time. Thus, when only singular values or a small portion of the vectors are needed, our approach is particularly favourable. One of the most fruitful advantages of our two-stage SVD algorithm is the attractive speedup shown when a portion of the singular vectors is computed. The performance obtained by DGESDD, when 20% of the singular vectors are required, shows more than a $4\times$ increase when using the 48 cores of system A (circles in Figure 8.6), and can get more than $2\times$ speedup when using the 16 cores of system B (circles in Figure 8.7). The underlying results are straightforward. The computing portion of the singular vectors significantly minimizes the overhead of the extra cost (applying $Q_{2}$ and $P_{2}$) introduced by the two-stage approach, and therefore the speedup observed here is relatively high.

Moreover, we also evaluated the performance and speedup of our algorithm when all the singular vectors are computed (diamonds in Figures 8.6 and 8.7). We can observe that our algorithm performs consistently better than the state-of-the-art optimized MKL routine DGESDD, with the exception of small matrices. Our implementation is around $3\times$ faster when using a large number of cores (48 cores of system A) and is around $1.5\times$ faster on a smaller number of cores (16 cores of system B).

Our fourth observation is related to the scalability of our implementation. To explain it briefly, we performed experiments with only 12 cores of system A, and compared it with the results from the full 48 cores. Our solver accomplished very good scalability: the execution time on 48 cores is about $3.5\times$ faster than the one obtained on 12 cores. It also appears from the speedup illustrated in our figures that the scalability increases with respect to the number of cores, and for that reason we believe that this implementation is suitable for large distributed and shared memory machines.

Figure 8.8.

The speedup obtained by our implementation of DGESVD versus its counterpart from the Intel MKL library on 48 AMD cores of system A.

Figure 8.9.

The speedup obtained by our implementation of DGESVD versus its counterpart from the Intel MKL library on 16 Sandy Bridge cores of system B.

Since we prefer to cover all the available SVD solvers and not be limited by the divide and conquer technique, we also performed the same set of experiments using the implicit zero-shift QR algorithm (DGESVD) as another solver variant. The speedup results obtained on both of our machines, system A and system B, are presented in Figures 8.8 and 8.9, respectively. Here we note that we could not compute a set of the singular vectors because the routine which computes the SVD decomposition of the bidiagonal matrix (DBDSQR) does not provide this option, and requires the transformation matrices $Q$ and $P^{H}$ to be provided explicitly as input. However, experiments with either all the singular vectors or none of the singular vectors were reported. The speedup of the DGEBRD routine has been reported in Figures 8.8 and 8.9, for the sake of completeness.

The performance shown in Figures 8.8 and 8.9 is very close to that of the divide and conquer solver. The performance drops slightly compared to Figures 8.6 and 8.7, because the DBDSQR solver used here requires the matrix $Q$ and $P$ to be computed explicitly as input, and so it updates the singular vectors $L_{s}$ and $R_{s}$ internally. Thus it does not benefit from all the optimizations we implemented. However, the explicit generation of the matrices $Q$ and $P$ take advantage of all the techniques described in this section. As a consequence, the speedup curve trend in Figures 8.8 and 8.9 slightly decreases compared to that of Figures 8.6 and 8.7. Therefore, reaching a twofold speedup (Figure 8.8) is worth the effort.

Finally, we demonstrate that our algorithm is very efficient and can achieve more than a twofold speedup over the well-known state-of-the-art optimized libraries. It is especially suitable when only the singular values – or when a portion of the singular vectors – are needed. Results show a $4\times$ to $10\times$ speedup. We believe that this achievement makes our algorithm a very good candidate for current and next generation computing hardware.

8.3.1 The classical reduction to tridiagonal condensed form

The classical approach (LAPACK algorithms) to reducing a matrix to tridiagonal form is to use a one-stage algorithm (Golub and Van Loan Reference Golub and Van Loan1996). Similar to the one-sided factorizations (Cholesky, LU, QR), the two-sided factorizations are split into a panel factorization and a trailing matrix update. Unlike the one-sided factorizations, however, the panel factorization requires computing the Level 2 BLAS symmetric matrix–vector product with the entire trailing matrix. This requires loading the entire trailing matrix into memory, which incurs a significant amount of memory-bound operations. It creates data dependencies and produces artificial synchronization points between the panel factorization and the trailing submatrix update steps. This prevents the use of standard techniques to increase the computational intensity (e.g. look-ahead), which are used extensively in the one-sided LU, QR, and Cholesky factorizations. Let us compute the cost of the algorithm. The reduction proceeds by steps of size $n_{b}$, where each step consists of the cost of the panel and the cost of the update.

• ∙ The panel is of size $n_{b}$ columns. The factorization of every column is dominated primarily by one symmetric matrix–vector product with the trailing matrix. Thus the cost of a panel is

  $$\begin{eqnarray}2n_{b}l^{2}+{\rm\Theta}(n),\end{eqnarray}$$
  where $l$ is the size of the trailing matrix. For simplicity, we omit ${\rm\Theta}(n)$ and round up the cost of the panel by the cost of the matrix–vector product.

• ∙ The update of the trailing matrix consists of applying the Householder reflectors generated during the panel factorization to the trailing matrix from both the left and right side according to

  $$\begin{eqnarray}A_{i+n_{b}:n,i+n_{b}:n}\leftarrow A_{i+n_{b}:n,i+n_{b}:n}-V\times W^{T}-W\times V^{T},\end{eqnarray}$$
  where $V$ and $W$ were computed during the panel phase. This Level 3 BLAS operation is computed by the SYR2K routine and its cost is $2n_{b}k^{2}$, where $k=n-in_{b}$ is the size of the trailing matrix at step $i$.

For all steps ($n/n_{b}$), the trailing matrix size varies from $n$ to $n_{b}$, by steps of size $n_{b}$, where $l$ varies from $n$ to $n_{b}$ and $k$ varies from $(n-n_{b})$ to $2n_{b}$. Thus, the total cost for the $n/n_{b}$ steps is

(8.6)$$\begin{eqnarray}\text{flops}\approx 2n_{b}\mathop{\sum }_{n_{b}}^{n/n_{b}}l^{2}+2n_{b}\mathop{\sum }_{2n_{b}}^{\frac{n-n_{b}}{n_{b}}}k^{2}\approx \frac{2}{3}n_{\mathtt{SYMV}}^{3}+\frac{2}{3}n_{\mathtt{SYR2K}}^{3}\approx \frac{4}{3}n^{3}.\end{eqnarray}$$

According to the equations above, we derive the maximum performance $P_{\max }$ that can be reached by such an algorithm. In particular, for large matrix sizes $n$, $P_{\max }$ is expressed as

(8.7)$$\begin{eqnarray}\displaystyle P_{\max } & = & \displaystyle \frac{\text{number of operations}}{\text{minimum time }t_{\min }}\nonumber\\ \displaystyle & = & \displaystyle \frac{\frac{4}{3}n^{3}}{t_{\min }({\textstyle \frac{2}{3}}n^{3}\text{ flops in }\mathtt{SYMV})+t_{\min }(\frac{2}{3}n^{3}\text{ flops in }\mathtt{SYR2K})}\nonumber\\ \displaystyle & = & \displaystyle \frac{\frac{4}{3}n^{3}}{\frac{2}{3}n^{3}\ast \frac{1}{P_{\mathtt{SYMV}}}+\frac{2}{3}n^{3}\ast \frac{1}{P_{\text{Level 3}}}}=\frac{2\ast P_{\text{Level 3}}\ast P_{\mathtt{SYMV}}}{P_{\text{Level 3}}+P_{\mathtt{SYMV}}}\nonumber\\ \displaystyle & \leq & \displaystyle 2P_{\mathtt{SYMV}}\quad \text{when }P_{\text{Level 3}}\gg P_{\mathtt{SYMV}}.\end{eqnarray}$$

Thus, from equation (8.7) we can understand the well-known low-performance behaviour of the classical tridiagonal reduction algorithm. In practice, it is the most time-consuming phase (about 70% of the time when all eigenvectors are requested, and about 90% of the time when only eigenvalues are needed). Figure 8.10 depicts the percentage of the time spent in each phase of the eigensolver using the classical one-stage approach to compute the tridiagonal form and eigenvectors.

Figure 8.10.

The percentage of the time spent in each phase of the eigensolver using the classical one-stage approach for the tridiagonal reduction.

8.3.2 Two-stage reduction

Similarly to the singular value decomposition, we can take advantage of the two-stage approach for the reduction to tridiagonal form. The two-stage reduction is designed to increase the utilization of compute-intensive operations. Many algorithms have been investigated using this two-stage approach. The idea is to split the original one-stage approach into a compute-intensive phase (first stage) and a memory-bound phase (second or ‘bulge chasing’ stage). In this section we will cover the description for the symmetric case. The first stage reduces the original symmetric dense matrix to a symmetric band form, while the second stage reduces from band to tridiagonal form, as depicted in Figure 8.11. The same algorithm used for reducing the general matrix $A$ to band form in the singular value decomposition case (Section 8.2.2) can be adapted to the symmetric case, where the symmetric matrix $A$ is reduced to tridiagonal form.

Figure 8.11.

Two-stage technique for the reduction phase.

Second stage: cache-friendly computational kernels.

The band matrix $A_{\text{band}}$ is further reduced to the tridiagonal form $T$ using the bulge chasing technique, similar to the one described in Section 8.2.2, but using kernels designed for a symmetric case. This procedure annihilates the extra off-diagonal elements by chasing the created fill-in elements down to the bottom right side of the matrix using successive orthogonal similarity transformations. Each annihilation of the $n_{b}$ non-zero element below the off-diagonal of the band matrix is called a sweep, as described in Figure 8.12. This stage involves memory-bound operations and requires the band matrix to be accessed from multiple disjoint locations. In other words, there is an accumulation of substantial latency overhead each time different portions of the matrix are loaded into cache memory, which is not compensated for by the low execution rate of the actual computations (the so-called surface-to-volume effect). To overcome these critical limitations, we used a bulge chasing algorithm to use cache-friendly kernels, combined with fine-grained memory-aware tasks, in an out-of-order scheduling technique, which considerably enhances data locality. We refer the reader to Haidar et al. (Reference Haidar, Ltaief and Dongarra2011, Reference Haidar, Tomov, Dongarra, Solcà and T.2014) for a detailed description of the technique.

Figure 8.12.

Two-stage technique for the reduction phase.

8.3.3 The tridiagonal eigensolver

The tridiagonal eigensolver is used to compute eigenpairs of the tridiagonal matrix:

(8.8)$$\begin{eqnarray}T=Z{\rm\Lambda}Z^{T},\end{eqnarray}$$

where $Z$ is the matrix of orthogonal eigenvectors of $T$, and ${\rm\Lambda}$ is the diagonal matrix encoding the eigenvalues. Four algorithms are available: QR iterations, bisection and inverse iteration (BI), divide and conquer (D&C), and multiple relatively robust representations (MRRR). The first two, QR and BI, are presented in Demmel (Reference Demmel1997), but a performance comparison made by Demmel, Marques, Parlett and Vomel (Reference Demmel, Marques, Parlett and Vomel2008b) compared LAPACK algorithms and concluded that D&C and MRRR are the fastest available solvers. However, while D&C requires a larger extra workspace, MRRR is less accurate. Accuracy is a fundamental parameter, because the tridiagonal eigensolver is known to be the part of the overall symmetric eigensolver where accuracy can be lost. D&C is more robust than MRRR, which can fail to provide an accurate solution in some cases. In theory, MRRR is a ${\rm\Theta}(n^{2})$ algorithm, whereas D&C is between ${\rm\Theta}(n^{2})$ and ${\rm\Theta}(n^{3})$, depending on the matrix properties. In many real-life applications, D&C is often less than cubic, while MRRR seems to be slower than expected due to the number of floating divisions and the cost of the iterative process. The main asset of MRRR is that a subset computation is possible, reducing the complexity to ${\rm\Theta}(nk)$ for computing $k$ eigenpairs. Such an option was not included within the classical D&C implementations, or it only reduced the cost of the updating phase of the last step of the algorithm.

8.3.4 The eigenvector computation

After the reduction to condensed form, the eigensolver finds the eigenvalues ${\rm\Lambda}$ and eigenvectors $Z$ of $T$. The eigenvalues are the same as for the original matrix $A$. To find the eigenvectors of the original matrix $A$, the eigenvectors $Z$ of $T$ need to be back-transformed by applying the same orthogonal matrices (e.g. $Q$ for the one-stage approach or $Q_{1}$ and $Q_{2}$ for the two-stage approach) that were used in the reduction to condensed form. For this purpose, the block Householder transformations $Q_{i}=I-V_{i}T_{i}V_{i}^{T}$ are used. For the one-stage approach, there is only one $Q$ that has been used in the reduction phase, and its multiplication by $Z$ can be formed explicitly using DORGTR; or we can multiply by $Q$ in an implicit fashion using DLARFB following the procedure described in Algorithm 7. In either case, applying $Q$ becomes a series of DGEMM operations. Since all of its computation is based on Level 3 BLAS, the performance upper bound of this phase is considered to be on a par with the performance of the Level 3 BLAS, which is usually a large fraction of the machine’s peak, and – for this phase – is considered to be a small percentage of the global time of the eigensolver. Figure 8.10 illustrates the percentage of time occupied by the one-stage algorithm’s three phases, and we can see that this phase consumes less than 10% of the total time.

Figure 8.13.

Speedup of our new two-stage implementation over the classic one-stage approach when both eigenvalues and eigenvectors are computed on 28 Intel E5-2697 cores.

Below, we discuss the application of the Householder reflectors generated from the two stages of the reduction to tridiagonal form. The first stage reduces the original symmetric and/or Hermitian matrix $A$ to a band matrix by applying a two-sided transformation to $A$ such that $A=Q_{1}A_{b}Q_{1}^{H}$. Similarly, the second stage, called bulge chasing, reduces the band matrix $B$ to the tridiagonal form by applying the transformation from both the left and right side of $A_{b}$ such that $A_{b}=Q_{2}TQ_{2}^{H}$. Thus, when the eigenvector matrix $X$ of $A$ is requested, the eigenvector matrix $Z$ that is produced by the eigensolver needs to be updated from the left by the two Householder reflectors generated during the reduction phase, according to the formula

(8.9)$$\begin{eqnarray}X=Q_{1}Q_{2}Z=(I-V_{1}T_{1}V_{1}^{H})(I-V_{2}T_{2}V_{2}^{H})Z,\end{eqnarray}$$

where ($V_{1},T_{1}$) and ($V_{2},T_{2}$) represent the Householder reflectors generated during the first and second reduction stages, respectively. The application of the $V_{2}$ reflectors is not straightforward, while the application of the $V_{1}$ reflectors has to follow the tile fashion style. The same technique used for the singular value decomposition is used here as well. Thus, Algorithms 8 and 9 describe the procedure of applying $Q_{2}$ and $Q_{1}$, respectively. More details of the optimized technique can be found elsewhere (Haidar et al. Reference Haidar, Luszczek and Dongarra2014b). Both multiplications are based on Level 3 BLAS and are thus considered to perform very well on recent hardware, which alleviates the extra cost of the application of $Q_{2}$. The overall speedup and efficiency of this technique is demonstrated in the experimental results section below.

Figure 8.14.

Speedup of our new two-stage implementation over the classic one-stage approach when only eigenvalues are computed using 28 Intel Xeon E5-2697 cores.

8.3.5 Experimental results

This section presents the performance comparisons of our tile algorithm for a two-stage symmetric eigenvalue solver against state-of-the-art numerical linear algebra libraries. To verify our findings in practice, we have performed the experiments on our two-socket multicore machine featuring a 14-core Intel Xeon E5-2697 processor in each socket running at 2.6 GHz, with 128 GB of main memory and 35 MB Level 2 cache. Figure 8.13 illustrates the speedup obtained by our two-stage approach over the classical one-stage implementation. The speedup curve shows that our two-stage algorithm outperforms the classical one-stage implementation by more than $2\times$ when all eigenvalues and eigenvectors are needed.

When only the eigenvalues are needed, there is no extra cost related to updating the vectors, and thus we can expect a large speedup similar to that of the singular value decomposition. Figure 8.14 shows the obtained speedup from our experiments for eigenvalues only. The speedup obtained here is between $5.5\times$ and $6\times$ faster when using the two-stage approach. In other words, the speedup varies depending on the portion of eigenvectors requested. For example, if 20% of the eigenvectors are needed, the last merge of our implementation of the divide and conquer algorithm can be modified to apply the update on only 20% of the vectors, thereby reducing the cost of the eigensolver by about 50%. The back transformation cost will be reduced by 80% for both the one-stage and two-stage approaches. The observed speedup in our experiment is about $4\times$ to $5\times$.

8.4.1 Hessenberg reduction

Initially, the matrix $A$ is reduced to an upper Hessenberg matrix $H$ using an orthogonal similarity transformation, $H=Q_{1}^{T}AQ_{1}$. The similarity transform preserves the eigenvalues, so $A$ and $H$ have the same eigenvalues, while using orthogonal matrices ensures good stability. Using elementary transformations based on an LU decomposition is also possible, which incurs half the floating-point operations as orthogonal matrices, but at reduced stability (Kabir, Haidar, Tomov and Dongarra Reference Kabir, Haidar, Tomov and Dongarra2015). The purpose of reducing the matrix to Hessenberg form is to reduce the operations in the second phase, QR iteration. From the Hessenberg form, only one sub-diagonal needs to be eliminated via QR iteration to achieve triangular Schur form.

Algorithm 11 presents a high-level view of the Hessenberg reduction. The reduction proceeds by introducing zeros below the sub-diagonal in each column $a_{j}$ using a Householder transformation,

$$\begin{eqnarray}H_{(j)}=I-{\it\tau}_{j}v_{j}v_{j}^{T}.\end{eqnarray}$$

As in the QR factorization in Section 5, instead of applying the Householder updates one by one using inefficient Level 2 BLAS operations, a blocked Householder transformation is used, where at step $j$,

$$\begin{eqnarray}Q_{(j)}=H_{(1)}H_{(2)}\cdots H_{(j)}=I-V_{(j)}T_{(j)}V_{(j)}^{T}\end{eqnarray}$$

with

$$\begin{eqnarray}\displaystyle V_{(j)} & = & \displaystyle \left[\begin{array}{@{}c@{}}v_{1},\ldots ,v_{j}\end{array}\right],\nonumber\\ \displaystyle T_{(j)} & = & \displaystyle \left[\begin{array}{@{}cc@{}}T_{(j-1)} & -{\it\tau}_{j}T_{(j-1)}V_{(j-1)}^{T}v_{j}\\ 0 & {\it\tau}_{j}\end{array}\right].\nonumber\end{eqnarray}$$

The update to the trailing matrix of $A$ is then

(8.10)$$\begin{eqnarray}\hat{A}=Q_{(j)}^{H}AQ_{(j)}=(I-V_{(j)}T_{(j)}^{H}V_{(j)}^{H})(A-Y_{(j)}T_{(j)}V_{(j)}^{H}),\end{eqnarray}$$

where

$$\begin{eqnarray}Y_{(j)}=AV_{(j)}.\end{eqnarray}$$

While we delay updating the trailing matrix until accumulating a block of $n_{b}$ Householder transformations, we must immediately update the next column of the panel with all previous transformations before we can compute its Householder transformation. For column $a_{j+1}$, this requires computing column $j+1$ of $Q_{(j)}^{H}AQ_{(j)}$. We can do this without computing the entire product, but it does require computing $y_{j}=Av_{j}$.

This matrix–vector product $Av_{j}$ is a memory-bound, Level 2 BLAS operation that becomes a major bottleneck in the algorithm, limiting the maximum performance to $P_{\max }\leq 5P_{\text{GEMV}}$. To achieve higher performance, we utilize GPUs to accelerate both GEMV and GEMM operations. Because the GEMV is memory-bound, it benefits from the increased memory bandwidth available on the GPU, while the GEMM benefits from the high computational rate of the GPU.

The GPU implementation is shown in Algorithm 12; differences from the CPU-only algorithm occur in lines 3, 6, 7, 8, 16 and 17. The entire matrix is stored on the GPU. The current panel is copied to the CPU to be factored there. During the panel factorization, the vector $v_{j}$ is sent to the GPU, where $y_{j}=Av_{j}$ is computed, and the resulting $y_{j}$ sent back to the CPU. After each panel factorization, the trailing matrix is updated on the GPU.

A multi-GPU implementation is also available. In this case, the matrix is distributed in a one-dimensional block-cyclic fashion across the GPUs. The GEMV is done independently on each GPU, with partial results sent to the CPU where a final summation is done. For instance, to compute $y=Av$ using three GPUs,

$$\begin{eqnarray}y=Av=\left[\begin{array}{@{}ccc@{}}A_{1} & A_{2} & A_{3}\end{array}\right]\left[\begin{array}{@{}c@{}}v_{1}\\ v_{2}\\ v_{3}\end{array}\right]=A_{1}v_{1}+A_{2}v_{2}+A_{3}v_{3}.\end{eqnarray}$$

Each GPU computes its part, $A_{i}v_{i}$, which are then sent to the CPU to be summed for the final result, $y$. After the panel factorization, the $V$, $Y$, and $T$ matrices are replicated on all the GPUs, then each GPU updates its portion of the trailing matrix. Figure 8.16 shows the performance of the multi-GPU version compared to Intel MKL running on two 8-core CPUs. A single GPU achieves a $3.3\times$ speedup over the CPU-only version, while three GPUs achieve up to a $9.2\times$ speedup. This shows good parallel speedup, with three GPUs being $2.8\times$ faster than a single GPU. The $P_{\max }<5P_{\text{GEMV}}$ bound is also shown in Figure 8.16 for both CPU and GPU. The GPU GEMV achieves 45.6 Gflops, while the CPU GEMV achieves up to 14.7 Gflops. On the GPU, the Hessenberg achieves 77% of the bound, while on the CPU, it achieves 72% of the bound.

Figure 8.16.

Execution time of Hessenberg reduction (DGEHRD) using two 8-core CPUs and up to three NVIDIA Kepler K40 GPUs. Dashed lines indicate performance bound, $P_{\max }=5P_{\text{GEMV}}$, for CPU and single GPU implementations.

8.4.2 Eigenvector computation

When eigenvectors are desired, the third phase computes eigenvectors of the triangular Schur form $S$, then back-transforms them to eigenvectors of the original matrix $A$. In LAPACK, this phase is implemented in the TREVC (triangular eigenvector computation) routine. We will assume only right eigenvectors are desired, but the computation of left eigenvectors is similar and amenable to the same techniques as described here. After the Hessenberg and QR iteration phases, the diagonal entries of $S$ are the eigenvalues ${\it\lambda}_{k}$ of $A$. To determine the corresponding eigenvectors, we solve $Sz_{k}={\it\lambda}_{k}z_{k}$ by considering the decomposition (Golub and Van Loan Reference Golub and Van Loan1996)

(8.11)$$\begin{eqnarray}\left[\begin{array}{@{}ccc@{}}S_{11} & u & S_{13}\\ 0 & {\it\lambda}_{k} & v^{T}\\ 0 & 0 & S_{33}\\ \end{array}\right]\left[\begin{array}{@{}c@{}}\hat{z}\\ 1\\ 0\\ \end{array}\right]={\it\lambda}_{k}\left[\begin{array}{@{}c@{}}\hat{z}\\ 1\\ 0\\ \end{array}\right],\end{eqnarray}$$

which yields $(S_{11}-{\it\lambda}_{k}I)\hat{z}=-u$. Thus computing each eigenvector $z_{k}$ of $S$ involves a $k-1\times k-1$ triangular solve, for $k=2,\ldots ,n$. Each solve has a slightly different $S$ matrix, with the diagonal modified by subtracting ${\it\lambda}_{k}$. The resulting eigenvector $z_{k}$ of $S$ must then be back-transformed by multiplying with the $Q$ formed in the Hessenberg and QR iteration phases to get the eigenvector $x_{k}=Qz_{k}$ of the original matrix $A$.

Note that if two eigenvalues, ${\it\lambda}_{k}$ and ${\it\lambda}_{j}$ ($k>j$), are identical, then $S_{11}-{\it\lambda}_{k}I$ is singular. More generally, $S_{11}-{\it\lambda}_{k}I$ can be badly conditioned. Therefore, instead of using the standard BLAS triangular solver (TRSV), a specialized triangular solver (LATRS) is used, which scales columns to protect against overflow, and can generate a consistent solution for a singular matrix.

This method works in complex arithmetic, but the case in real arithmetic is more complicated. For a real matrix $A$, the eigenvalues can still be complex, coming in conjugate pairs, ${\it\lambda}_{k}$ and $\bar{{\it\lambda}}_{k}$. The eigenvectors are likewise conjugate pairs, $z_{k}$ and $\bar{z}_{k}$. In real arithmetic, the closest that QR iteration can come to triangular Schur form is quasi-triangular real Schur form, which has a $2\times 2$ diagonal block for each conjugate pair of eigenvalues. A specialized quasi-triangular solver is required, which factors each $2\times 2$ diagonal block, and protects against overflow and dealing with singular matrices. In LAPACK this solver is implemented inline as part of the DTREVC routine.

8.4.3 Blocking back-transformation

Our first step to improve the eigenvector computation is to block the $n$GEMV operations for the back-transformation into $n/n_{b}$GEMM operations, where $n_{b}$ is the block size. This requires two $n\times n_{b}$ workspaces: one for the vectors $z_{k}$, the second for the back-transformed vectors $x_{k}$ before copying to the output $V$.

Pseudocode for the blocked back-transformation is shown in Algorithm 13 along with the parallel solver described in Section 8.4.4. For each block, we loop over $n_{b}$ columns, performing a triangular solve for each column and storing the resulting eigenvectors $z_{k}$ in workspace $Z$. After filling up $n_{b}$ columns of $Z$, a single GEMM back-transforms all $n_{b}$ vectors, storing the result in workspace $\tilde{Z}$. The vectors are then normalized and copied to $V$. On input, the matrix $V=Q$. Recall from equation (8.11) that the bottom $n-k$ rows of eigenvector $z_{k}$ are 0, so the last $n-k$ columns of $Q$ are not needed for the GEMM. Therefore, we start from $k=n$ and work down to $k=1$, writing each block of eigenvectors to $V$ over columns of $Q$ after they are no longer needed.

The real case is similar, but has the minor complication that complex conjugate pairs of eigenvalues will generate conjugate pairs of eigenvectors, $z_{k}=a+bi$ and $\bar{z}_{k}=a-bi$, which are stored as two columns, $a$ and $b$, in $Z$. When the first eigenvalue of each pair is encountered, both columns are computed; then the next eigenvalue (its conjugate) is skipped. Once $n_{b}-1$ columns are processed, if the next eigenvector is complex it must be delayed until the next block.

8.4.4 Multithreading triangular solver

After blocking the back-transform, the triangular solver remains a major bottleneck because it is not parallelized. Recall that the triangular matrix being solved is different for each eigenvector – the diagonal is modified by subtracting ${\it\lambda}_{k}$. This prevents blocking and solving multiple eigenvectors together using a Level 3 BLAS DTRSM operation.

In the complex case, LAPACK’s ZTREVC uses a safe triangular solver, ZLATRS. Unlike the standard ZTRSV BLAS routine, ZLATRS uses column scaling to avoid numerical instability, and handles singular triangular matrices. Therefore, to avoid jeopardizing the accuracy or stability of the eigensolver, we continue to rely on ZLATRS instead of the optimized, multithreaded ZTRSV. However, processing $S$ column by column prevents parallelizing individual calls to ZLATRS, as is typically done for BLAS functions. Instead, we observe that multiple triangular solves could occur in parallel. One obstacle is that a different ${\it\lambda}_{k}$ is subtracted from the diagonal in each case, modifying $S$ in memory. Our solution is to write a modified routine, ZLATRSD (triangular solve with modified diagonal), which takes both the original unmodified $S_{11}$ and the ${\it\lambda}_{k}$ to subtract from the diagonal. The subtraction is done as the diagonal elements are used, without modifying $S$ in memory. This allows us to pass the same $S$ to each ZLATRSD call and hence solve multiple eigenvectors in parallel, one in each thread.

As previously mentioned, the real case requires a special quasi-triangular solver to solve each $2\times 2$ diagonal block. In the original LAPACK code, this quasi-triangular solver is embedded in the DTREVC routine. To support multithreading, we refactor it into a new routine, DLAQTRSD, a quasi-triangular solver with modified diagonal. Unlike the complex case, instead of passing ${\it\lambda}_{k}$ separately, DLAQTRSD computes it directly from the diagonal block of $S$. If ${\it\lambda}_{k}$ is real, DLAQTRSD computes a single real eigenvector. If ${\it\lambda}_{k}$ is one of a complex conjugate pair, DLAQTRSD computes a complex eigenvector as two real vectors.

To deal with multithreading, we use a thread pool design pattern. As shown in Algorithm 13, the main thread inserts LATRSD tasks into a queue. Worker threads pull tasks out of the queue and execute them. For this application, there are no dependencies to be tracked between the triangular solves. After a block of $n_{b}$ vectors has been computed, we back-transform them with a GEMM. We could call a multithreaded GEMM, as available in MKL, but to simplify thread management, we found it more suitable to use the same thread pool for the GEMM as for LATRSD. For $p$ threads, the GEMM is split into $p$ tasks, each task multiplying a single block row of $Q$ with $Z$. After the GEMM, the next block of $n_{b}$ vectors is computed. Within each thread, the BLAS calls are single-threaded.

8.4.5 Using accelerators

To further accelerate the eigenvector computation, we observe that the triangular solves and the back-transformation GEMM can be done in parallel. In particular, the GEMM can be done on a GPU while the CPU performs the triangular solves. Data transfers can also be done asynchronously and overlapped with the triangular solves. To facilitate this asynchronous computation, we double-buffer the CPU workspace $Z$, using it to store results from LATRSD, then swap with $\tilde{Z}$, which is used to send data to the GPU while the next block of LATRSD solves are performed with $Z$. The difference from Algorithm 13 is shown in Algorithm 14.

9.1.1 Direct methods

A common approach to the solution of linear systems, either dense or sparse, is to perform the LU factorization of the coefficient matrix using Gaussian elimination. First, the coefficient matrix $A$ is factored into the product of a lower triangular matrix $L$ and an upper triangular matrix $U$. Partial row pivoting is, in general, used to improve numerical stability resulting in a factorization $PA=LU$, where $P$ is a permutation matrix. The solution for the system is achieved by first solving $Ly=Pb$ (forward substitution) and then solving $Ux=y$ (backward substitution). Due to round-off errors, the computed solution $x$ carries a numerical error magnified by the condition number of the coefficient matrix $A$.

In order to improve the computed solution, we can apply an iterative process which produces a correction to the computed solution at each iteration, which then yields the method that is commonly known as the iterative refinement algorithm. As Demmel (Reference Demmel1997) points out, the non-linearity of the round-off errors makes iterative refinement equivalent to Newton’s method applied to the function $f(x)=b-Ax$. Provided that the system is not too ill-conditioned, the algorithm produces a solution correct for the working precision. Iterative refinement in double–double precision is a fairly well-understood concept and was analysed by Wilkinson (Reference Wilkinson1963), Moler (Reference Moler1967) and Stewart (Reference Stewart1973).

The algorithm can be modified to use a mixed precision approach. The factorization $PA=LU$, the solution of the triangular systems $Ly=Pb$, and $Ux=y$ are computed using single precision arithmetic. The residual calculation and the update of the solution are computed using double precision arithmetic and the original double precision coefficients: see Algorithm 15. The most computationally expensive operation, the factorization of the coefficient matrix $A$, is performed using single precision arithmetic and takes advantage of its higher speed. The only operations that must be executed in double precision are the residual calculation and the update of the solution (denoted by ${\it\varepsilon}_{d}$ in Algorithm 15). We observe that the only operation with computational complexity of $O(n^{3})$ is handled in single precision, while all operations performed in double precision are of at most $O(n^{2})$ complexity. The coefficient matrix $A$ is converted to single precision for the LU factorization and the resulting factors are stored in single precision while the initial coefficient matrix $A$ needs to be kept in memory. Therefore, one drawback of the following approach is that it uses 50% more memory than the standard double precision algorithm.

The method in Algorithm 15 can offer significant improvements for the solution of a sparse linear system in many cases.

1. 1 Single precision computation is significantly faster than double precision computation.

2. 2 The iterative refinement procedure converges in a small number of steps.

3. 3 The cost of each iteration is small compared to the cost of the system factorization. If the cost of each iteration is too high, then a low number of iterations will result in a performance loss with respect to the full double precision solver. In the sparse case, for a fixed matrix size, both the cost of the system factorization and the cost of the iterative refinement step may vary substantially depending on the number of non-zeros and the matrix sparsity structure. In the dense case, results are more predictable.

Note that the choice of the stopping criterion in the iterative refinement process is critical. Formulas for the error computed at each step of Algorithm 15 can be obtained, for instance, in Oettli and Prager (Reference Oettli and Prager1964) and Demmel et al. (Reference Demmel, Hida, Kahan, Li, Mukherjee and Riedy2006).

9.1.2 Iterative methods

Direct methods are usually a very robust tool for the solution of sparse linear systems. However, they suffer from fill-in, which results in high memory requirements, long execution time, and non-optimal scalability in parallel environments. To overcome these limitations, various pivot reordering techniques are commonly applied to minimize the amount of generated fill-in and to enable better exploitation of parallelism. Still, there are cases where direct methods pose too high a memory requirement or deliver poor performance. Iterative methods are a valid alternative even though they are less robust and have less predictable behaviour. Iterative methods do not require more memory than is needed for the original coefficient matrix. Furthermore, time to solution can be better than that of direct methods if convergence is achieved in relatively few iterations (Barrett et al. Reference Barrett, Berry, Chan, Demmel, Donato, Dongarra, Eijkhout, Pozo, Romine and van der Vorst1994a, Saad Reference Saad2003).

In the context of iterative methods, the refinement method outlined in Algorithm 15 can be represented as

(9.2)$$\begin{eqnarray}x_{i+1}=x_{i}+M(b-Ax_{i}),\end{eqnarray}$$

where $M$ is $(LU)^{-1}P$. Note that in this case $M$ approximates $A^{-1}$, which makes it known, in the area of linear algebra and numerical analysis, as a preconditioner of $A$. The preconditioner should approximate $A^{-1}$, and the quality of the approximation determines the convergence properties of the iterative solver. In general, a preconditioner is intended to improve the robustness and the efficiency of iterative methods. Note that (9.2) can also be interpreted as a Richardson method’s iteration in solving $MAx=Mb$, which is called left preconditioning. An alternative is to use right preconditioning, whereby the original problem $Ax=b$ is transformed into a problem of solving

$$\begin{eqnarray}AMu=b,\quad x=Mu\end{eqnarray}$$

iteratively. Later on, we will use right preconditioning for mixed precision iterative methods.

$M$ needs to be easy to compute, apply, and store to guarantee the overall efficiency. Note that these requirements were addressed in the mixed precision direct methods above by replacing $M$ (coming from LU factorization of $A$ followed by matrix inversion) with its single precision representation so that arithmetic operations can be performed more efficiently on it. Here, however, we go two steps further. We not only replace $M$ with an inner loop, which is an incomplete iterative solver working in single precision arithmetic (Turner and Walker Reference Turner and Walker1992), but we also replace the outer loop with a more sophisticated iterative method (e.g. based on Krylov subspace).

Note that replacing $M$ with an iterative method leads to nesting of two iterative methods. Variations of this type of nesting, also known in the literature as an inner–outer iteration, have been studied, both theoretically and computationally (Golub and Ye Reference Golub and Ye2000, Saad Reference Saad1993, Simoncini and Szyld Reference Simoncini and Szyld2003, Axelsson and Vassilevski Reference Axelsson and Vassilevski1991, Notay Reference Notay2000, Vuik Reference Vuik1995, van den Eshof, Sleijpen and van Gijzen Reference van den Eshof, Sleijpen and van Gijzen2005). The general appeal of these methods is that the computational speedup hinges on the inner solver’s ability to use an approximation of the original matrix $A$ that is fast to apply. In our case, we use single precision matrix–vector products as a fast approximation of the double precision operator in the inner iterative solver. Moreover, even if no faster matrix–vector product is available, speedup can often be observed due to improved convergence (see e.g. Simoncini and Szyld Reference Simoncini and Szyld2003, who explain the possible benefits of FGMRES–GMRES over restarted GMRES).

To illustrate the above concepts, we demonstrate an inner–outer non-symmetric iterative solver in mixed precision. The solver is based on restarted generalized minimal residual (GMRES). In particular, consider Algorithm 16, where the outer loop uses flexible GMRES (FGMRES: Saad Reference Saad1993, Reference Saad2003) and the inner loop uses GMRES in single precision arithmetic (GMRES$_{\text{SP}}$). FGMRES, being a minor modification of standard GMRES, is meant to accommodate non-constant preconditioners. Note that in our case this non-constant preconditioner is GMRES$_{\text{SP}}$. The resulting method is denoted by FGMRES($m_{\text{out}}$)–GMRES$_{\text{SP}}$($m_{\text{in}}$), where $m_{\text{in}}$ is the restart for the inner loop and $m_{\text{out}}$ for the outer FGMRES. Algorithm 16 checks for convergence at every $m_{\text{out}}$ outer iterations, but an actual implementation can check for convergence at every inner iteration with simple tricks, for almost no computational cost.

The potential benefits of FGMRES compared to GMRES are becoming better understood (Simoncini and Szyld Reference Simoncini and Szyld2003). Numerical experiments confirm improvements in speed, robustness, and sometimes memory requirements for these methods (Buttari et al. Reference Buttari, Dongarra, Kurzak, Luszczek and Tomov2008a, Baboulin et al. Reference Baboulin, Buttari, Dongarra, Kurzak, Langou, Langou, Luszczek and Tomov2009). The memory requirements for the method are the matrix $A$ in CRS format, the non-zero matrix coefficients in single precision, $2\,m_{\text{out}}$ number of vectors in double precision, and $m_{\text{in}}$ number of vectors in single precision.

The generalized conjugate residuals (GCR) method (Vuik Reference Vuik1995, van der Vorst and Vuik Reference van der Vorst and Vuik1994) is a possible replacement for FGMRES as the outer iterative solver. It is not yet well understood whether one should choose GCR or FGMRES.

As in the dense case, the choice of the stopping criterion in the iterative refinement process is critical. In the sparse case, formulas for the errors can be computed following the work of Arioli, Demmel and Duff (Reference Arioli, Demmel and Duff1989).

9.2.1 Direct methods

The mixed precision iterative refinement solvers for dense matrices are available in LAPACK (Anderson et al. Reference Anderson, Bai, Bischof, Blackford, Demmel, Dongarra, Ducroz, Greenbaum, Hammarling, McKenney and Sorensen1999). For the non-symmetric case, step 1 in Algorithm 15 is implemented by means of the SGETRF subroutine, steps 2, 3 and 5, 6 with the SGETRS subroutine, step 4 with the DGEMM subroutine, and step 7 with the DAXPY subroutine. For the symmetric case the SGETRF, SGETRS, and DGEMM subroutines were replaced by the SPOTRF, SPOTRS, and DSYMM subroutines, respectively. Further details of these implementations can be found in Langou et al. (Reference Langou, Langou, Luszczek, Kurzak, Buttari and Dongarra2006) and Buttari et al. (Reference Buttari, Dongarra, Langou, Langou, Luszczek and Kurzak2007).

As already mentioned, iterative refinement solvers require $1.5\times$ as much memory as a regular double precision solver, because the mixed precision iterative refinement solvers need to store both the single precision and the double precision versions of the coefficient matrix at the same time. This applies to both dense and sparse matrices.

Numerical results show that on older single-core architectures – such as the AMD Opteron, Intel Woodcrest, or IBM PowerPC – the mixed precision iterative solver could provide a speedup of up to $1.8\times$ for the non-symmetric solver and $1.5\times$ for the symmetric solver for sufficiently large problem sizes (Baboulin et al. Reference Baboulin, Buttari, Dongarra, Kurzak, Langou, Langou, Luszczek and Tomov2009). For small problem sizes, the cost of even a few iterative refinement iterations is high compared to the cost of the factorization, and thus the mixed precision iterative solver is less efficient than the double precision solver (Baboulin et al. Reference Baboulin, Buttari, Dongarra, Kurzak, Langou, Langou, Luszczek and Tomov2009).

The mixed precision iterative refinement solvers were ported and optimized for multicore architectures in the PLASMA library. On the STI Cell BE architecture, for example, where the peak for single precision operations was $14\times$ higher than the peak for double precision operations, the mixed precision solver performed up to $7\times$ faster than the double precision peak in the non-symmetric case and $11\times$ faster for the symmetric positive definite case. Implementation details for this case can be found in Kurzak and Dongarra (Reference Kurzak and Dongarra2007) and Kurzak, Buttari and Dongarra (Reference Kurzak, Buttari and Dongarra2008).

Figure 9.1.

Mixed precision, iterative refinement method from MAGMA 1.6.2 for the solution of dense linear systems on the NVIDIA’s TITAN X GPU.

Parallel implementations of Algorithm 15 have been available for GPUs since MAGMA’s 0.2 release (Tomov, Nath, Du and Dongarra Reference Tomov, Nath, Du and Dongarra2009). Figure 9.1 illustrates the performance of MAGMA 1.6.2 on a Titan X GPU. Due to the $32\times$ difference between the speed of FP32 and FP64 on this GPU, the mixed precision solver performs up to $26\times$ faster than the double precision peak. Implementation details can be obtained from the MAGMA sources.Footnote ⁶ Results on older GPUs, including the use of multiple GPUs, can be found in Tomov, Nath, Ltaief and Dongarra (Reference Tomov, Nath, Ltaief and Dongarra2010b).

For sparse matrices, the methodology can be applied using sparse direct methods: with either multifrontal (Duff and Reid Reference Duff and Reid1983) or supernodal (Ashcraft et al. Reference Ashcraft, Grimes, Lewis, Peyton and Simon1987) factorization approaches. For results on supernodal solvers see Buttari et al. (Reference Buttari, Dongarra, Kurzak, Luszczek and Tomov2008a). For multifrontal solvers, there are a number of freely available packages. The use of MUMPS (Amestoy, Duff and L’Excellent Reference Amestoy, Duff and L’Excellent2000, Amestoy, Duff, L’Excellent and Koster Reference Amestoy, Duff, L’Excellent and Koster2001, Amestoy, Guermouche, L’Excellent and Pralet Reference Amestoy, Guermouche, L’Excellent and Pralet2006), for example, comprises three separate steps.

1. 1 System analysis. In this phase the system sparsity structure is analysed in order to estimate the element fill-in, which provides an estimate of the memory that will be allocated in the following steps. Also, pivoting is performed based on the structure of $A+A^{T}$, ignoring numerical values. Only integer operations are performed at this step.

2. 2 Matrix factorization. In this phase the $PA=LU$ factorization is performed. Computationally, this is the most expensive step of the system solution.

3. 3 System solution. The system is solved in two steps, $Ly=Pb$ and $Ux=y$.

The first and second phases correspond to step 1 in Algorithm 15, while the third phase corresponds to steps 2, 3 and 5, 6.

Table 9.1.

Performance improvements for direct sparse methods when going from a full double precision solve (reference time) to a mixed precision solve.

Table 9.1 shows the speedup of the mixed precision iterative refinement approach over the double precision approach for sparse direct methods, using the matrices described in Table 9.2.

Table 9.2.

Test matrices for sparse mixed precision, iterative refinement solution methods.

9.2.2 Iterative methods

Similar to the case of sparse direct solvers, we demonstrate the numerical performance of Algorithm 16 on the matrices from Table 9.2.

Figure 9.2.

Mixed precision iterative refinement with FGMRES–GMRES$_{\text{SP}}$ from Algorithm 16 versus FGMRES–GMRES$_{\text{DP}}$.

Figure 9.2 shows the performance ratio of the mixed precision inner–outer FGMRES–GMRES$_{\text{SP}}$ versus the full double precision inner–outer FGMRES–GMRES$_{\text{DP}}$. In other words, we compare two inner–outer algorithms that are virtually the same. The only difference is that their inner loop’s incomplete solvers are performed in single and double precision arithmetic.

For further details see Baboulin et al. (Reference Baboulin, Buttari, Dongarra, Kurzak, Langou, Langou, Luszczek and Tomov2009), which includes an experiment showing that inner–outer-type iterative methods may be very competitive compared to their original counterparts.

10.4.1 Methodology based on batched BLAS

In a batched problem solution methodology that is based on batched BLAS, there are many small dense matrices that must be factorized simultaneously (as illustrated in Figure 10.1). This means that all the matrices will be processed simultaneously by the same kernel. However, each matrix problem is still solved independently, identified by a unique batch ID. We follow this model in our batched implementations and developed the following set of new batched CUDA kernels.

• ∙ Cholesky panel. This provides the batched equivalent of LAPACK’s DPOTF2 routine. At step $j$ of a panel of size $(m,nb)$, the column vector $A(j:m,j)$ must be computed. This requires a dot-product using row $A(j,1:j)$ to update element $A(j,j)$, followed by a DGEMV$A(j+1,1)A(j,1:j)=A(j+1:m,j)$, and finally a DSCAL on column $A(j+1:m,j)$. This routine involves two Level 1 BLAS calls (dot and scal), as well as a Level 2 BLAS DGEMV. Since there are $nb$ steps, these routines are called $nb$ times, and thus one can expect the performance to depend on the performances of Level 2 and Level 1 BLAS operations. Hence, it is a slow, memory-bound algorithm. We used shared memory to load both row $A(j,1:j)$ and column $A(j+1:m,j)$ to re-use them, and wrote a customized batched DGEMV kernel to read and write these vectors from/into the shared memory.

• ∙ LU panel. This provides the batched equivalent of LAPACK’s DGETF2 routine to factorize panels of size $m\times nb$ at each step of the batched LU factorizations. It consists of three Level 1 BLAS calls (idamax, DSWAP and DSCAL) and one Level 2 BLAS call (DGER). The DGETF2 procedure proceeds as follows. Find the maximum element of the $i$th column, then swap the $i$th row with the row owning the maximum, and scale the $i$th column. To achieve higher performance and minimize the effect on the Level 1 BLAS operation, we implemented a tree reduction to find the maximum pivot where all the threads contribute to find the pivot. Since it is the same column that is used to find the max, then scaled, we load it into the shared memory. These are the only data that we can re-use within one step.

• ∙ QR panel. This provides the batched equivalent of LAPACK’s DGEQR2 routine to perform the Householder panel factorizations. It consists of $nb$ steps where each step calls a sequence of the DLARFG and the DLARF routines. At every step (to compute one column), the DLARFG involves a norm computation followed by a DSCAL that uses the results of the norm computation in addition to some underflow/overflow checking. The norm computation is a sum reduce and thus a synchronization step. To accelerate it, we implemented a two-layer tree reduction where, for sizes larger than 32, all 32 threads of a warp progress required to do a tree reduction similar to the MPI_REDUCE operation – and the latest 32 element – are reduced by only one thread. Another optimization is to allow more than one thread block to execute the DLARFG kernel, which means that the kernel needs to be split in two – one for norm and one for scaling – in order to guarantee the synchronization. Customized batched implementations of both DLARFG and the DLARF have been developed.

• ∙ Trailing matrix updates. Trailing matrix updates are mainly Level 3 BLAS operations. However, for small matrices, it might be difficult to extract performance from very small Level 3 BLAS kernels. The DGEMM is the best Level 3 BLAS kernel: it is GPU-friendly, highly optimized, and achieves the highest performance among BLAS. Thus, high performance can be achieved if we redesign our update kernels to be represented by DGEMMs. For Cholesky, the update consists of the DSYRK routine. It performs a rank-$nb$ update on either the lower or upper portion of $A_{22}$. Since cuBLAS does not provide a batched implementation of this routine, we implemented our own. It is based on a sequence of customized DGEMMs in order to extract the best possible performance. The trailing matrix update for Gaussian elimination (LU) is composed of three routines: the DLASWP, which swaps the rows on the left and right of the panel in consideration, followed by the DTRSM to update $A_{12}\leftarrow L_{11}^{-1}A_{12}$, and finally a DGEMM for the update $A_{22}\leftarrow A_{22}-A_{21}L_{11}^{-1}A_{12}$. The swap (or pivoting) is required to improve the numerical stability of Gaussian elimination. However, pivoting can be a performance killer for matrices stored in column-major format because rows in that case are not stored continuously in memory, and thus cannot be read en masse. Indeed, a factorization stored in column-major format can be $2\times$ slower (depending on hardware and problem sizes) than implementations that transpose the matrix in order to use a row-major storage format internally (Volkov and Demmel Reference Volkov and Demmel2008). Nevertheless, experiments showed that this conversion is too expensive for batched problems. Moreover, the swapping operations are serial, that is row by row, which limits the parallelism. To minimize this penalty, we propose a new implementation that emphasizes a parallel swap and allows coalescent read/write. We also developed a batched DTRSM. It loads the small $nb\times nb$$L_{11}$ block into shared memory, inverts it with the DTRTRI routine, and then the $A_{12}$ update is accomplished by a DGEMM. Generally, computing the inverse of a matrix may compromise the numerical stability of the method, but since $A_{11}$ results from numerically stable LU with partial pivoting, and its size is just $nb\times nb$, or in our case $32\times 32$, we do not have this problem (DuCroz, Dongarra and Higham Reference Ducroz, Dongarra and Higham1992). For the Householder QR decomposition the update operation is used by the DLARFB routine. We implemented a batched DLARFB that is composed of three calls to the batched DGEMM:

  $$\begin{eqnarray}A_{22}\leftarrow (I-VT^{H}V^{H})A_{22}\equiv (I-A_{21}T^{H}A_{21}^{H})A_{22}.\end{eqnarray}$$

10.5.1 Parallel swapping

Profiling the batched LU reveals that more than 60% of the time is spent in the swapping routine. Figure 10.2 shows the execution trace of the batched LU for $2000$ matrices of size $512$. We can observe on the top trace that the classical DLASWP kernel is the most time-consuming part of the algorithm. The swapping consists of $nb$ successive interchanges of two rows of the matrices. The main reason that this kernel is the most time-consuming is because the $nb$ row interchanges are performed in a sequential order, and the data of a row are not coalescent, thus the thread warps do not read/write it in parallel. It is clear that the main bottleneck here is the memory access. Indeed, slow memory access compared to high computing capability has been a persistent problem for both CPUs and GPUs. CPUs alleviate the effect of long latency operations and bandwidth limitations by using hierarchical caches.

Figure 10.2.

Execution trace of the batched LU factorization using either classical swap (a) or our new parallel swap (b).

Accelerators, on the other hand, in addition to hierarchical memories, use thread level parallelism (TLP), where threads are grouped into warps and multiple warps are assigned for execution on the same SMX unit. The idea is that when a warp issues an access to the device memory, it stalls until the memory returns a value, while the accelerator’s scheduler switches to another warp. In this way, even if some warps stall, others can execute, keeping functional units busy while resolving data dependencies, branch penalties, and long latency memory requests. In order to overcome the bottleneck of swapping, we propose to modify the kernel to apply all $nb$ row swaps in parallel. This modification will also allow the coalescent write back of the top $nb$ rows of the matrix. Note that the first $nb$ rows are those used by the DTRSM kernel that is applied right after the DLASWP, so one optimization is to use shared memory to load a chunk of the $nb$ rows, and apply the DLASWP followed by the DTRSM. We changed the algorithm to generate two pivot vectors. The first vector stores the indices of the final destination (e.g. row indices) of the top $nb$ rows of the panel. The second stores the row indices of the $nb$ rows to swap out and bring into the top $nb$ rows of the panel. Figure 10.2 depicts the execution trace (bottom) when using our parallel DLASWP kernel. The experiment shows that this reduces the time spent in the kernel from 60% to around 10% of the total elapsed time. Note that the colours between the top and the bottom traces do not match each other; this is because the NVIDIA profiler automatically chooses the tracing colour. As a result, the performance gain obtained is about $1.8\times$ as shown in Figure 10.3. We report each of the proposed optimizations for the LU factorization in Figure 10.3, but we would like to mention that the percentage of improvement obtained for the Cholesky and QR factorization is similar, and to simplify we reported the LU factorization only. Note that starting from this version we were able to be faster than the cuBLAS implementation of the batched LU factorization.

Figure 10.3.

Performance in Gflops of the different versions of our batched LU factorizations compared to the cuBLAS implementation for different matrix sizes, where $m=n$.

10.5.2 Recursive nested blocking

The panel factorizations described in Section 10.4.1 factorize the $nb$ columns one after another, similarly to the LAPACK algorithm. At each of the $nb$ steps, either a rank-one update is required to update the vectors to the right of the factorized column $i$ (this operation is done by the DGER kernel for LU and the DLARF kernel for QR), or a left-looking update of column $i$ by the columns on its left is required, before factorizing it (this operation is done by DGEMV for the Cholesky factorization). Since we cannot load the entire panel into the shared memory of the GPU, the columns to the right (in the cases of LU and QR) or to the left (in the case of Cholesky) are loaded back and forth from the main memory at every step. Thus, one can expect that this is the most time-consuming part of the panel factorization. A detailed analysis using the profiler reveals that the DGER kernel requires more than $80$% of the panel time and around $40$% of the total LU factorization time. Similarly for the QR decomposition, the DLARF kernel used inside the panel computation needs $65$% of the panel time and $33$% of the total QR factorization time. Likewise, the DGEMV kernel used within the Cholesky panel computation needs around $91$% of the panel time and $30$% of the total Cholesky factorization time. These routines’ inefficient behaviour is also due to the memory access. To overcome that bottleneck, we propose to improve the efficiency of the panel and to reduce the memory access by using a recursive level of blocking technique: see Figures 10.4 and 10.5. In principle, the panel can be blocked recursively until it is a single element. However, in practice, $2$–$3$ blocked levels are sufficient to achieve high performance. The above routines must be optimized for each blocked level, which complicates the implementation. More than a $30$% boost in performance is obtained with this optimization, as demonstrated in Figure 10.3 for the LU factorization. The same trend has been observed for both the Cholesky and the QR factorization.

Figure 10.4.

Recursive nested blocking.

Figure 10.5.

Execution trace of the batched LU factorization using either classical GETF2 (a) or our recursive GETF2 (b).

10.5.3 Trading extra computation for higher performance

The challenge discussed here is as follows. For batched problems there is a need to minimize the use of low-performance kernels on the GPU even if they are Level 3 BLAS. For the Cholesky factorization, this concerns the DSYRK routine that is used to update the trailing matrix. The performance of DSYRK is important to the overall performance, since it takes a big part of the runtime. We implemented the batched DSYRK routine as a sequence of DGEMM routines, each of size $M=m,N=K=nb$. In order to exclusively utilize the DGEMM kernel, our implementation writes both the lower and the upper portion of the $nb\times nb$ diagonal blocks of the trailing matrix. This results in $nb^{3}$ extra operations for the diagonal block. However, since $nb$ is small (e.g. $nb=32$), these extra operations can be considered free. In practice the extra operation allows us to use DGEMM and thus achieve higher performance than the one that touches the lower/upper portion of the $nb\times nb$ diagonal blocks. Tests show that our implementation of DSYRK is $2\times$ faster than the DGEMM kernel for the same matrix size. This shows that our DSYRK is very well optimized for it to reach the performance of DGEMM (which is $2\times$ slower since it computes twice as many flops).

We applied the same technique in the DLARFB routine used by the QR decomposition. The QR trailing matrix update uses the DLARFB routine to perform

$$\begin{eqnarray}A_{22}=(I-VT^{H}V^{H})A_{22}=(I-A_{21}T^{H}A_{21}^{H})A_{22}.\end{eqnarray}$$

The upper triangle of $V$ is zero with ones on the diagonal. $A_{21}$, which stores $V$ in its lower triangular part and $R$ (part of the upper $A$) in its upper triangular part, is available in the classical DLARFB. Therefore, the above is computed using DTRMM for the upper part of $A_{21}$ and DGEMM for the lower part. Also, the $T$ matrix is an upper triangular and therefore the classical DLARFB implementation uses DTRMM to perform the multiplication with $T$. Thus, if one can guarantee that the lower portion of $T$ is filled with zeros and the upper portion of $V$ is filled zeros and ones on the diagonal, the DTRMM can be replaced by DGEMM. With that in mind, we implemented a batched DLARFB that uses three DGEMM kernels by initializing the lower portion of $T$ with zeros, and filling up the upper portion of $V$ with zeros and ones on the diagonal. Note that this brings $3nb^{3}$ extra operations, but again, the overall time spent in the new DLARFB update using the extra computation is around 10% less than DTRMM.

Similarly to DSYRK and DLARFB, we implemented the batched DTRSM (that solves $AX=B$) by inverting the small $nb\times nb$ block $A$ and using DGEMM to get the final result $X=A^{-1}B$.

10.5.4 Block recursive DLARFT algorithm

The DLARFT is used to compute the upper triangular matrix $T$ that is needed by the QR factorization in order to update either the trailing matrix or the right-hand side of the recursive portion of the QR panel. The classical LAPACK approach computes $T$ column by column in a loop over the $nb$ columns as described in Algorithm 19. Such an implementation takes up to $50$% of the total QR factorization time. This is due to the fact that the kernels needed – DGEMV and DTRMV – require implementations where threads go through the matrix in different directions (horizontal versus vertical, respectively). An analysis of the mathematical formula of computing $T$ allowed us to redesign the algorithm to use Level 3 BLAS and to increase the data re-use by putting the column of $T$ in shared memory. One can observe that the loop can be split into two loops – one for DGEMV and one for DTRMV. The DGEMV loop that computes each column of $\widehat{T}$ can be replaced with one DGEMM to compute all the columns of $\widehat{T}$ if the triangular upper portion of $A$ is zero and the diagonal is made of ones. This is already needed for our implementation’s trailing matrix update in order to use DGEMM in the DLARFB, and thus can be exploited here as well. For the DTRMV phase, we load the $T$ matrix into shared memory as this allows all threads to read/write from/into shared memory during the $nb$ steps of the loop. The redesign of this routine is depicted in Algorithm 20. Since we developed a recursive blocking algorithm, we have to compute the $T$ matrix for every level of the recursion. Nevertheless, the analysis of Algorithm 20 lets us conclude that the portion of the $T$s computed in the lower recursion level is the same as the diagonal blocks of the $T$ of the upper level (white triangular blocks in Figure 10.6), and thus we can avoid their (re-)computation. For that, we modified Algorithm 20 in order to compute either the whole $T$ or the upper rectangular portion that is missed (dark/pale grey portions in Figure 10.6).