FLAME-derived LTLT factorization algorithms for skew-symmetric matrices, with fused BLAS-like operations, greatly outperform PFAPACK and Pfaffine

The work systematically derives a family of algorithms for the skew-symmetric LTLT factorization using the FLAME methodology, including unpivoted and pivoted right-looking (Parlett-Reid variant), left-looking (Aasen variant), two-step right-looking (Wimmer variant), three blocked right-looking algorithms, and one blocked left-looking algorithm. Prior algorithms for skew-symmetric LTLT factorization were scattered and mostly adapted from the symmetric case; this work is the first to uniformly derive them from a specification via the FLAME formal workflow, providing correctness proofs and an algorithm family. The paper presents specifications, partitioned matrix expressions, loop invariants, and algorithm pseudocode, with full derivations in a companion technical report; the paper itself states it repeats only brief highlights while remaining self-contained.

The paper identifies new level-2 and level-3 BLAS-like operations required by this factorization, including skew rank2, gen rank2, skew tridiag gemv, skew tridiag rankk, skew tridiag gemm, and skew rank2k. The traditional BLAS interface contains no skew-symmetric operations; this work explicitly lists these operations and where they appear in each algorithm, providing a basis for future BLAS extensions. The paper tabulates each operation and which algorithms use it, and notes that some skew-symmetric operations have recently been added to the reference implementation.

By fusing the skew-symmetric tridiagonal rank-k update (skew tridiag rankk) with BLIS packing and parallelizing level-2 operations and block pivot application, performance improves by more than 5x (unpivoted) and 3x (pivoted) over the unfused version. Prior implementations typically decomposed operations into traditional BLAS calls, incurring extra data movement and workspace; this work leverages BLIS 2.0's custom packing mechanism to avoid these overheads. The paper reports step-by-step optimization experiments on a 2x AMD EPYC 7763 system using 64 cores, showing the performance impact of each optimization step.

The best implementations greatly exceed the performance of the only known prior skew-symmetric LTLT factorization software, PFAPACK and Pfaffine, while matching or exceeding related symmetric factorization software. PFAPACK does not exceed 5.8 GFLOPs compared to a peak of 135 GFLOPs for Pfaffian; the presented implementations force the full factorization mode and thus perform twice as many FLOPs, yet the performance gap is much larger than a factor of two. The paper compares on 64 cores (one socket) with an optimal block size determined from experiments, and notes that PFAPACK and Pfaffine were compiled from source with the same compiler.

This work targets dense skew-symmetric LTLT factorization on shared-memory multi-core CPUs, applicable to scenarios requiring the full L and T factors (such as fast-updating Pfaffians or solving Xv = w); the pivoted versions improve numerical stability. The algorithm family and BLAS-like operation definitions can be modified back for symmetric tridiagonalization, and the fusion and BLAS extension ideas can transfer to other linear algebra and tensor operations.

This is an incomplete reading; specific performance data and curve details in the figures could not be fully captured, so the precise speedup of each optimization step and performance differences across block sizes cannot be verified. Additionally, the parallel efficiency of pivoted algorithms is notably lower than unpivoted ones (the best pivoted algorithm achieves only 24% parallel efficiency at 64 cores), with pivot application wasting cache lines and exhibiting poor TLB reuse due to column-major storage and row swaps; further optimization of this bottleneck remains to be explored. The paper notes that whether a left-looking partial factorization algorithm can perform even fewer operations is an open question, and the formal inclusion of pivoting in FLAME derivation is a future direction.