On the parallel I/O optimality of linear algebra kernels

Abstract

Matrix factorizations are among the most important building blocks of scientific computing. However, state-of-the-art libraries are not communication-optimal, underutilizing current parallel architectures. We present novel algorithms for Cholesky and LU factorizations that utilize an asymptotically communication-optimal 2.5D decomposition. We first establish a theoretical framework for deriving parallel I/O lower bounds for linear algebra kernels, and then utilize its insights to derive Cholesky and LU schedules, both communicating N 3 /(P √ M) elements per processor, where M is the local memory size. The empirical results match our theoretical analysis: our implementations communicate significantly less than Intel MKL, SLATE, and the asymptotically communication-optimal CANDMC and CAPITAL libraries. Our code outperforms these state-of-the-art libraries in almost all tested scenarios, with matrix sizes ranging from 2,048 to 524,288 on up to 512 CPU nodes of the Piz Daint supercomputer, decreasing the time-to-solution by up to three times. Our code is ScaLAPACK-compatible and available as an open-source library.

Pebbles, Graphs, and a Pinch of Combinatorics: Towards Tight I/O Lower Bounds for Statically Analyzable Programs

“…𝑋 −𝑆 (Lemma 2 in Kwasniewski et al [27]). The expression 𝜌 = 𝜒 (𝑋 ) 𝑋 −𝑆 is called the computational intensity.…”

“…To bound the sizes of rectangular subcomputations, we use two lemmas given by Kwasniewski et al [27]: [27]) For statement 𝑆𝑡, given 𝑫, the size of subcomputation H (number of vertices of 𝑆 computed during H ) is bounded by the sizes of the iteration variables' sets 𝐷 𝑡 , 𝑡 = 1, . .…”

I/O-Optimal Algorithms for Symmetric Linear Algebra Kernels

“…Independently, using explicit enumeration of data reuse, Kwasniewski et al [8] obtain a corresponding lower bound for LU factorization: their proof is in a parallel context, but their arguments show that the minimum number of data transfers is lower bounded by 2 √ lower bound, which however makes the implicit assumption that there is no data reuse related to the symmetry of the matrix as discussed in Section 1.…”

“…Regarding Cholesky, Ballard et al [2] reviewed existing parallel distributed algorithms and, based on their lower bound on communication, proved that LAPACK and other block recursive implementations are asymptotically optimal for a carefully selected block size. The work on lower bounds by Kwasniewskiet al [8] leads to the design of parallel distributed 2.5 LU (COnfLUX) and Cholesky (COnfCHOX) algorithms. These algorithms perform a volume of communication per node of 3 √ + O ( 2 ).…”

“…Indeed it is akin to computing the lower half of C = A • B, where A is a × matrix and B is a × matrix. Similarly, a lower bound for Cholesky is derived in [8] under the constraint that it is forbidden to use , when , is available. The lower bound obtained with such a constraint is potentially too large: there may exist algorithms which perform fewer data transfers by making a better use of symmetry.…”