Achieving numerical accuracy and high performance using recursive tile LU factorization with partial pivoting

Abstract

The LU factorization is an important numerical algorithm for solving systems of linear equations in science and engineering and is a characteristic of many dense linear algebra computations. For example, it has become the de facto numerical algorithm implemented within the LINPACK benchmark to rank the most powerful supercomputers in the world, collected by the TOP500 website. Multicore processors continue to present challenges to the development of fast and robust numerical software due to the increasing levels of hardware parallelism and widening gap between core and memory speeds. In this context, the difficulty in developing new algorithms for the scientific community resides in the combination of two goals: achieving high performance while maintaining the accuracy of the numerical algorithm. This paper proposes a new approach for computing the LU factorization in parallel on multicore architectures, which not only improves the overall performance but also sustains the numerical quality of the standard LU factorization algorithm with partial pivoting. While the update of the trailing submatrix is computationally intensive and highly parallel, the inherently problematic portion of the LU factorization is the panel factorization due to its memory-bound characteristic as well as the atomicity of selecting the appropriate pivots. Our approach uses a parallel fine-grained recursive formulation of the panel factorization step and implements the update of the trailing submatrix with the tile algorithm. Based on conflict-free partitioning of the data and lockless synchronization mechanisms, our implementation lets the overall computation flow naturally without contention. The dynamic runtime system called QUARK is then able to schedule tasks with heterogeneous granularities and to transparently introduce algorithmic lookahead. The performance results of our implementation are competitive compared to the currently available software packages and libraries. For example, it is up to 40% faster when compared to the equivalent Intel MKL routine and up to threefold faster than LAPACK with multithreaded Intel MKL BLAS. HIGH-PERFORMANCE RECURSIVE TILE LU FACTORIZATION 1409Successful high-performance results have already been reported for one-sided factorizations (e.g., QR/LQ, LU, and Cholesky factorizations) and, more recently, for the tridiagonal reduction needed to solve the symmetric eigenvalue problems [3] and the bidiagonal reduction required for the singular value decomposition [4][5][6]. These implementations are based on tile algorithms, which operate on the original matrix using small square regions called tiles; they alleviated the bottlenecks of and block algorithms on column-major storage by bringing the parallelism to the fore, minimizing the synchronization overhead, and relying on dynamic scheduling of fine-grained tasks. The data within the tiles can be contiguously stored in memory (i.e., tile data layout) or left as is, following the standard column-major format, which makes the data layout completely indep...

Recursion based parallelization of exact dense linear algebra routines for Gaussian elimination

“…This optimization used in the computation of the slab factorization improved the computation speed by a factor of 2. This approach shares similarities with the recursive computation of the panel described in [9].…”

Energy Footprint of Advanced Dense Numerical Linear Algebra Using Tile Algorithms on Multicore Architectures

“…We refrain from publishing a comprehensive set of performance charts for the LU factorizations as we have done elsewhere [26]. Figure 1 shows the achieved performance of mixed precision solvers with iterative refinement on both systems for non-symmetric matrices (LU).…”

Porting the PLASMA Numerical Library to the OpenMP Standard

“…Because of the row pivoting, the panel factorization (factorization of a block of columns) cannot be tiled, and the panel has to be dealt with as a whole, which handicaps cache efficiency. In recent years, this problem has been addressed by algorithms that try to achieve good level of cache residency for this operation [12,16]. Many alternative approaches have also been developed [14], but this article focuses on the classic formulations.…”

“…This is in contrast with LAPACK, where one tall panel (block of columns) is eliminated at a time, making it difficult to achieve cache efficiency and apply multithreading. In the course of the PLASMA project, tile algorithms have been developed for a wide range of algorithms, including: Cholesky, LU and QR factorizations [11,14,16], as well as reductions to band forms for solving the singular value problem or the eigenvalue problem [23,31].…”