Skew-Symmetric Matrix Decompositions on Shared-Memory Architectures

Fig. 2 depicts the “basic” right- and left-looking unblocked algorithms, which are equivalent to the algorithms of Parlett and Reid, and of Aasen, respectively. The algorithms as written contain some extra provision for the case in which the first column of $L$ (below the diagonal) is not assumed to be zero as it typically is when performing a full factorization. Actually, as Rozložník pointed out [13], the first column of $L$ may be initialized arbitrarily, leading to distinct but functionally equivalent factorizations. When used to solve the panel sub-problem in a blocked algorithm, this first column will typically contain the last (and usually non-zero) elimination vector from the previous panel. However, whether or not this column should be considered in the sub-problem depends on the structure of the blocked algorithm (see Fig. 3).

The left-looking algorithm performs a matrix-vector multiply each iteration, which is the leading-order cost in terms of FLOPs. While the intermediate product (a portion of $TL^{T}$) can be computed separately in $\mathcal{O}(n)$ time, we denote the entire $LTL^{T}$ matrix-tridiagonal-vector multiplication as a single fused operation gemv-sktri. At iteration $i$, this then incurs $2(n-i-1)(i+1)$ FLOPs and hence a total leading-order cost of $n^{3}/3$ FLOPs. Instead, the right-looking algorithm performs a skew-symmetric rank-2 update (skr2) every iteration $i$ at a cost of $2(n-i-2)^{2}$ FLOPs for a total leading-order cost of $2n^{3}/3$ FLOPs. Because the unblocked algorithm deals only with a “tall-and-skinny” trapezoidal region of $X$ for a blocked sub-problem, the trailing region of the skr2 operation becomes a general rank-2 update (denoted here as ger2).

Wimmer [22] halved the cost of the right-looking algorithm by skipping every other column of $L$, leading to a partial tridiagonal structure of $T$. This structure is sufficient for computing the Pfaffian but not for inversion or solving systems of equations. We derive a roughly similar algorithm as the “two-step” right-looking algorithm in Fig. 4. This algorithm performs a single rank-2 update for every iteration (i.e. every two rows/columns), thus yielding the same $n^{3}/3$ leading-order cost as Wimmer’s partial factorization. However, the two-step algorithm does in fact compute every column of $L$. This example illustrates the potential for reducing computation through operation fusion (in this algorithm, the updates from the two computed columns of $L$ are naturally combined), and importantly, for enabling and identifying such opportunities through an alternative arrangement of the problem definition (the PME, invariant, and iteration conditions) and subsequent formal analysis, rather than through manual post-optimization.

A blocked right-looking algorithm was first proposed by Rozložník et al. [13] and was implemented in LAPACK [26]. This algorithm has significant similarity to our algorithm, depicted in Fig. 3, with the most significant difference being the need to maintain a partial block of $H=LT$ in the former. Instead, our algorithm directly exposes the skew-symmetry of the block update step as an operation with the form of a symmetric rank-$k$ update, but with a skew-symmetric tridiagonal matrix inserted in the middle of the product. We denote this operation as gemmt-sktri due to the fact that pre-combining $H=LT$ or $\tilde{H}=TL^{T}$ leads to a gemmt (general matrix-matrix multiplication updating only the upper or lower half), but with implicit generation of the Hessenberg factor from a tridiagonal-times-matrix multiplication. As $H$ is also not required in the unblocked sub-problem, it can be omitted entirely, reducing workspace and FLOPs. The leading-order cost for both our algorithm and Rozložník’s is only $n^{3}/3$; this is in contrast to the blocked Parlett-Reid algorithm introduced by Wimmer which uses a skew-symmetric rank-$2k$ update (skr2k) and hence requires twice as many FLOPs. All blocked algorithms derived to date utilize the left-looking unblocked (Aasen) algorithm to factorize each panel. When pivoting, this is necessary as right-looking updates are only applied within the panel. Columns pivoted in from outside the panel will then not satisfy the necessary pre-condition. For completeness, we also implement blocked unpivoted algorithms which use a right-looking unblocked algorithm for the sub-problem. However, we opt not to use our new two-step algorithm for blocked sub-problems.

A blocked left-looking algorithm was also derived in [17] and is depicted in Fig. 3. While Yamazaki et al. hint at the possibility of a left-looking blocked Aasen algorithm [26], they do not provide any algorithm nor, to our knowledge, has any implementation been made. Like the right-looking algorithm, it requires $n^{3}/3$ FLOPs, but features a general matrix-matrix product (gemm) to update the next panel of $X$. While considering an intermediate Hessenberg product does perturb the skew-symmetry of the matrix or the efficiency of the gemm, it does require additional storage and memory accesses to perform the tridiagonal-matrix product. We denote the combined general matrix-tridiagonal-matrix product as gemm-sktri. The left-looking algorithm is conventionally considered as lower-performing compared to the right-looking algorithm due to the reduced level of parallelism available in the matrix-matrix product compared to the rank-$k$ update of the right-looking algorithm. Additionally, the left-looking blocked algorithm cannot incorporate partial pivoting as the panel update cannot be applied until pivots are chosen and vice versa.

Finally, we implement the novel blocked right-looking algorithm which starts with a more fine-grained partitioning of the PME and invariant [17]. This algorithm “shifts” the computation by one row/column, resulting in the factorization of column $x_{43}$ (via $L_{44}e_{f}$) during the panel update, and hence allows the trailing update to be computed entirely in the form $X\mathrel{-}=LTL^{T}$ (a fused gemmt-sktri operation). As with the two-step algorithm, an alternative setup of the problem followed by formal analysis exposes additional opportunities for fusion and optimization.

Similarly to other DLA routines, such as those implemented in LAPACK, it is often desirable to reuse the storage of the input matrix for the decomposed form. For an $X=LTL^{T}$ decomposition, the strictly lower or upper $n(n-1)/2$ elements are available, which can exactly hold the unique non-unit/zero elements of $T$ and $L$. For example, the lower part of $L$ (except the leading and presumably zero column) may be stored, “shifted” left by one column, below the sub-diagonal of $X$, with elements of $T$ stored on the sub-diagonal. Because the algorithms here do not require any additional arrays or matrices, they could then be implemented without any external workspace (apart from the pivot vector $p$). However, we note that the updates often explicitly reference the diagonal of $L$, which is implicitly one. For example, in the blocked right-looking algorithm, the update appears,

When stored as described above, the “1” which is in fact $\lambda_{33}$ would refer to the same storage as $x_{32}^{T}e_{l}$, which in turn physically stores $\tau_{32}$. This value cannot simply be temporarily overwritten with an explicit 1, as it is also referenced in the central tridiagonal part of the update. One option is to split the update into separate steps, for example,

Instead, we opt to store the sub-diagonal elements of $T$ in an external vector, and to store explicit 1’s on the sub-diagonal of $X$ as columns of $L$ are computed. This allows the entire update to be computed using a single, fused gemmt-sktri operation. Interestingly, the choice of the storage location of the $L$ and $T$ matrices is expressible as a part of the PME and invariants by imposing additional conditions. Dependency analysis of the resulting logical expressions at the “before” and “after” points can then illuminate when such fusion optimizations can and cannot be performed. Storing $T$ as an external array also allows elements to be accessed contiguously.

While the leading-order cubic cost of the blocked algorithms arises mostly from level-3 BLAS or BLAS-like operations, the unblocked algorithm which operates within blocks works entirely within level-1 and level-2 and, as can be seen in our results, accounts for a considerable portion of the running time. We apply two optimizations to these operations: first, we include OpenMP parallelization of all level-2 operations (including the application of blocks of pivots, see e.g. Fig. 3). As most BLAS libraries do not parallelize level-2 operations, including BLIS, we implement in-house versions, although we use BLIS kernels where applicable. Second, we have implemented “fused” versions of those level-2 operations not appearing in standard BLAS (we count skr2 as a standard operation for the present purposes since it is a minor modification of syr2). These new operations include gemv-sktri and ger2. In the latter case, a speedup approaching $2\times$ may be expected due to accessing the output matrix only once, while in the former case we can only expect a very slight speedup over a naive implementation which first computes $z=Tx$ (with $\mathcal{O}(n)$ memory accesses) and then $y=Az$ ($\mathcal{O}(n^{2})$ accesses) rather than $y=ATx$ (also $\mathcal{O}(n^{2})$ accesses) all at once. Where possible, we apply loop unrolling and jamming (e.g. computing multiple axpys with a single loop) in order to maximize register and prefetch pipeline usage.

Level-3 BLAS operations are typically compute-bound, performing $\mathcal{O}(n^{3})$ memory accesses for only $\mathcal{O}(n^{2})$ memory accesses.¹¹1Actually, any level-3 implementation must perform $\mathcal{O}(n^{3})$ accesses on one matrix (usually the output matrix), although in practice these accesses are carefully prefetched and amortized such that the actual overhead is dominated by accessing the other two matrices in $\mathcal{O}(n^{2})$ time Thus, adding $\mathcal{O}(n^{2})$ overhead, for example to compute $X\mathrel{-}=LTL^{T}$ as $H=TL^{T};X\mathrel{-}=LH$, does not seem like a limitation. However, the skew-symmetric tridiagonal-times-matrix product $H=TL^{T}$ is, as with level-2 operations, entirely memory-bound. In order to fuse these operations together, we take advantage of the flexibility of the BLIS framework, where user-defined “packing” operations can be specified, to create a custom BLAS-like operation. Packing is a typical part of gemm (and other level-3) implementations, and serves to load blocks of the input matrices into a special storage format which provides linear, vectorized access in the inner kernel as well as better TLB reuse. Because packing already necessitates accessing memory in the input matrices, we compute the $TL^{T}$ product on-the-fly and store it into the packed buffer. This custom packing operation has negligible increase in cost over standard packing because, although each element of $TL^{T}$ requires reading two elements of $L$ (due to the skew-symmetric tridiagonal structure), the accesses can take advantage of high temporal locality. The rest of the operation is handled without modification by BLIS, yielding the same overall performance as standard level-3 operations.

A sample implementation of the column-pivoted blocked right-looking algorithm is depicted in Fig. 5. We make use of MArray,[11] which is a header-only DLA library developed by us, similar in spirit to Eigen [8], Armadillo [14], Blitz++ [20, 21], and others. We have extended this library with a range-based partitioning API which expresses the partition–repartition pattern of FLAME. Ranges are strongly typed, with singleton ranges (corresponding to partitions containing only one column or row) typed as scalar integers, and others typed as a C++ class representing the half-open range $[from,to)$. Indexing matrix or vector types with such ranges produces a result which depends on the range types, yielding a sub-matrix, sub-vector, or single scalar entry. Finally, BLAS and BLAS-like operations (besides simple scalar-vector operations which use the expression template functionality of YYY) use custom wrapper functions such that matrix lengths, strides, and other properties can be provided by the relevant objects. Even for unblocked algorithms, which are sensitive to overhead in terms of branches, function calls, and additional memory accesses (register spills, updating structs, etc.), typical compilers (in our case, we use gcc 11.4.0) are capable of optimizing the combination of while loop and range partition–repartition pattern into the same assembly as a plain for loop, with the exception of a small number of spurious conditional moves. As these instructions are spurious (referring to conditions already known), they can be perfectly predicted. This observation matches our profiling data which shows negligible overhead in the ltlt_unb_var function bodies.

In order to investigate the impact of the several optimizations described in Section IV, we measured the performance of a sequence of algorithms in which the optimizations are applied in “steps”:

1. Step 0:

   Only existing BLAS functions (plus skr2) are used from BLIS. Algorithms as in Figs. 2 and 3 are used. Both $L$ and $T$ are stored in $X$ as they are computed and updates involving the diagonal of $L$ are unfused, see Section IV-A.

2. Step 1:

   Fused level-2 operations gemv-sktri and ger2 are introduced, without parallelization.

3. Step 2:

   All level-2 operations as well as block pivoting are parallelized.

4. Step 3:

   $T$ is stored in an external vector and fused updates are used which reference (unit) diagonal elements of $L$.

5. Step 4:

   Fused level-3 operations gemm-sktri and gemmt-sktri are introduced.

6. Step 5:

   The fused blocked right-looking algorithm of Fig. 4 is used.

Optimizations at each step are applied cumulatively. The performance of the blocked right-looking/unblocked left-looking algorithm (with and without pivoting) at each step is presented in Fig. 6, using one socket (64 threads).

Three different algorithmic block sizes are tested: 128 (top row), 256 (middle row), and 512 (bottom row). Successively large block sizes imply more time spent in the unblocked panel factorization, but also more aggregation into level-3 trailing updates. Step 1 has a very small positive impact, as might be expected from the fusion of an essentially level-1 step into the gemv operation. Step 2 has a much larger effect, particularly at large core counts, as approximately 10 cores are required to fully saturate the memory bandwidth. This effect grows with algorithmic block size due to the larger amount of work in level-2 operations. The impact of Step 3 is also very small. At this step, essentially a gemv-sktri is fused with an “adjacent” gemmt-sktri (at the blocked level, fusion at the unblocked level is less important). As the input sub-matrices for the gemv-sktri come only from the current block they typically reside in L2 or L3 cache. The vector to update is typically in main memory, but is only $\mathcal{O}(n)$ in size. This is in contrast to Steps 4 and 5, where now the operations which are fused operate on $\mathcal{O}(n^{2})$ data. The effect of these optimizations is much larger. The effect of fusing the trailing skr2 update (Step 5) is particularly pronounced for small block sizes. Cumulatively, these optimizations increase performance by more than $5\times$ in the unpivoted case, and $3\times$ in the pivoted case.

We measure the parallel scalability of our algorithms by performing a strong scaling experiment from one to all 128 cores of the system for a single $20,\!000\times 20,\!000$ skew-symmetric $X=LTL^{T}$ factorization. The results are presented in Fig. 7, with unpivoted and pivoted blocked algorithms in the top and bottom panels, respectively. All algorithms involve non-trivial memory bound work, primarily in the unblocked sub-problem and in pivoting. Additionally, the unfused blocked right-looking algorithm performs an level-2 skr2 operation which adds considerable overhead (see Section V-A). The impact of these operations can be clearly seen in the parallel scaling, as memory bandwidth only increases with core count up to $\sim 10$ cores, and latency is significantly impacted for inter-CCD and especially inter-socket memory traffic. The unblocked algorithms experience a drop of $\sim 25\%$ in parallel efficiency (which is proportional to per-core performance) from 1 to 16 or 32 threads. Some of this is attributable to a drop in processor frequency due to thermal management. The single core results certainly show an effect of frequency boosting as the base clock single core peak is only 39.2 GFLOPs, while peak at max boost frequency is 56 GFLOPs. Assuming max boost, the fused right-looking algorithm (with or without pivoting) achieves $\sim 85\%$ of peak.

In order to measure this effect, we consider the MKL implementation of gemm as the “speed of light” and present relative performance of our algorithms in Fig. 7c. Indeed, the relative performance for 16+ threads is higher than indicated by panel a), although only slightly. The difference grows with core count, showing a progressive drop in processor frequency, although significant boosts are evident even for 32 or more cores. From panel c), the biggest parallel efficiency drops are at 16 cores, likely due to exceeding a single 8-core cluster (CCD), and then at 64 threads and beyond. Based on profiling data, the main effect is due to increased time in level-2 operations, indicating main memory bandwidth as the limiting factor (inter-cache bandwidth would already be affected at core counts greater than 8 due to inter-CCD communication). In the two socket case, the impact of NUMA remote accesses is likely an additional factor. The scalability of the pivoted algorithms is even more severely affected. In particular, we measure the application of block pivots (see Fig. 3 for example) as the largest overhead. This is primarily because the columns of $L$ to pivot are stored in column-major layout, while rows are swapped. Thus, accesses “waste” most of each cache line loaded and TLB reuse is poor, especially for very large matrices. We attempt to optimize stores by blocking across the row dimension to increase temporal locality. Overall, the best unpivoted algorithm achieves 41% parallel efficiency (53% of gemm) at 64 cores, while the best pivoted algorithm achieves only 24% (29% of gemm). We observe similar poor scalability of existing pivoted factorizations, e.g. pivoted Cholesky (pstrf) and Bunch-Kaufman symmetric tridiagonalization (sytrf) although data is not shown here.

In Fig. 8 we compare the performance of our implementation to related algorithms as well as existing implementations of skew-symmetric tridiagonal factorization. In Fig. 8ab, we compare to common unpivoted and pivoted operations which share computational and algorithmic similarities to $X=LTL^{T}$, but which are expected to be highly optimized given their ubiquity in scientific computing. Panel a) shows a comparison to Cholesky decomposition (potrf). Cholesky decomposition incurs the same number of FLOPs as our algorithms, but is expected to be more efficient due to several factors which we do not take advantage of. The first is recursive decomposition (more than the two levels employed by us), which is necessary to achieve the communication/IO lower bound [5]. The second is the use of level-3 BLAS operations (trsm in the case of Cholesky) to update the trailing portion of blocks of the triangular factor. While the MKL implementation does exceed our performance comfortably, we find that the fused right-looking algorithm with left-looking unblocked updates is competitive with the OpenBLAS Cholesky implementation and much more efficient than that in AOCL. The unfused right-looking algorithm, as well as algorithms which use unblocked right-looking updates experience a significant drop in performance. This is almost exclusively due to the poor parallel performance of skr2. We use a non-uniform thread partitioning scheme, but load imbalance and cache conflicts still significantly impact performance. The left-looking blocked algorithm, even though it does not perform an additional level-2 operation at the blocked level, falls behind the right-looking algorithm. As mentioned above, this is at least in part due to reduced parallelism in block updates, depending on the stage of the algorithm. In panel b), we compare to pivoted Cholesky (pstrf) and column-pivoted $LU$ decomposition (getrf). Both algorithms exhibit similar pivoting requirements and other algorithmic properties, although $LU$ decomposition performs twice as many FLOPs. We see that both of our pivoted algorithms achieve significantly better performance compared to pstrf, but are in turn outperformed by getrf. As getrf is one of the most-used DLA operations it is unsurprisingly highly optimized in all major libraries. Additionally, getrf can make use of recursive blocking and trsm for updates to $U$ (as with Cholesky), and can aggregate most of row swaps due to pivoting while also not requiring complementary column swaps.

In comparison to symmetric tridiagonal factorizations (sytrf and sytrf_aa, Fig. 8c), our algorithms achieve uniformly and significantly higher performance than all implementations, with the exception of MKL’s sytrf for matrix sizes less than $6000\times 6000$. In large part, this difference is not so much attributable to algorithmic deficiencies in these algorithms (for example, the Rozložník algorithm implemented in sytrf_aa should have similar intrinsic performance to our fused blocked right-looking algorithm), but rather due to lack of optimization in the implementation. The difference between the MKL and other sytrf implementation highlights this effect. While we applied several optimizations in our work which might not be leveraged in the lower-performing implementations, the use of a flexible, expressive, and performant C++ interface along with the ability to enable and expose optimization opportunities through formal derivation methods is immediately applicable to many (if not almost all) other DLA routines in LAPACK and beyond. Thus, the use of our techniques should significantly lower the barrier to producing highly-efficient DLA implementations. Finally, in Fig. 8d we compare our algorithms to implementations of Wimmer’s blocked right-looking algorithm from PFAPACK and Pfaffine. These implementations are very similar, with the latter switching from Fortran to C++ and replacing the reference, unoptimized skr2k of PFAPACK with calls to an optimized version, presumably implemented in terms of gemmt. This difference is highly evident from the results where PFAPACK does not exceed 5.8 GFLOPs, compared to a peak of 135 GFLOPs for Pfaffian. Both implementations also perform twice as many FLOPs as our algorithms, thus a reduction in FLOPs plus other optimizations applied here (parallelization of level-2 operations and fusion) accounts for the observed performance difference.

Overall, these results validate the optimizations and implementation strategy that we have applied, as the observed performance is roughly comparable to existing optimized implementations of related operations, and much higher than existing implementations of the specific $X=LTL^{T}$ factorization (symmetric or skew-symmetric).