This dissertation presents an automated system to generate highly efficient, platform-adapted implementations of sparse matrix kernels. We show that conventional implementations of important sparse kernels like sparse matrix-vector multiply (SpMV) have historically run at 10% or less of peak machine speed on cache-based superscalar architectures. Our implementations of SpMV, automatically tuned using a methodology based on empirical-search, can by contrast achieve up to 31% of peak machine speed, and can be up to 4× faster. 
Given a matrix, kernel, and machine; our approach to selecting a fast implementation consists of two steps: (1) we identify and generate a space of reasonable implementations, and then (2) search this space for the fastest one using a combination of heuristic models and actual experiments (i.e., running and timing the code). We build on the SPARSITY system for generating highly-tuned implementations of the SpMV kernel y ← y + Ax, where A is a sparse matrix and x, y are dense vectors. We extend SPARSITY to support tuning for a variety of common non-zero patterns arising in practice, and for additional kernels like sparse triangular solve (SpTS) and computation of ATA·x (or AAT·x) and A ρ·x. 
We develop new models to compute, for particular data structures and kernels, the best absolute performance (e.g., Mflop/s) we might expect on a given matrix and machine. These performance upper bounds account for the cost of memory operations at all levels of the memory hierarchy, but assume ideal instruction scheduling and low-level tuning. We evaluate our performance with respect to such bounds, finding that the generated and tuned implementations of SpMV and SpTS achieve up to 75% of the performance bound. This finding places limits on the effectiveness of additional low-level tuning (e.g., better instruction selection and scheduling). (Abstract shortened by UMI.)

/pdf/automatic-performance-tuning-of-sparse-matrix-kernels-76ql9zn9l5.pdf

Automatic performance tuning of sparse matrix kernels

We consider the problem of adaptive subspace tracking, when the rank of the subspace we seek to estimate is assumed to be known. Starting from the data projection method (DPM), which constitutes a simple and reliable means for adaptively estimating and tracking subspaces, we develop a fast and numerically robust implementation of DPM, which comes at a considerably lower computational cost. Most existing schemes track subspaces corresponding either to the largest or to the smallest singular values, while our DPM version can switch from one subspace type to the other with a simple change of sign of its single parameter. The proposed algorithm provides orthonormal vector estimates of the subspace basis that are numerically stable since they do not accumulate roundoff errors. In fact, our scheme constitutes the first numerically stable, low complexity, algorithm for tracking subspaces corresponding to the smallest singular values. Regarding convergence towards orthonormality our scheme exhibits the fastest speed among all other subspace tracking algorithms of similar complexity.

Fast and Stable Subspace Tracking

Achieving peak performance from the computational kernels that dominate application performance often requires extensive machine-dependent tuning by hand. Automatic tuning systems have emerged in response, and they typically operate by (1) generating a large number of possible, reasonable implementations of a kernel, and (2) selecting the fastest implementation by a combination of heuristic modeling, heuristic pruning, and empirical search (i.e. actually running the code). This paper presents quantitative data that motivate the development of such a search-based system, using dense matrix multiply as a case study. The statistical distributions of performance within spaces of reasonable implementations, when observed on a variety of hardware platforms, lead us to pose and address two general problems which arise during the search process. First, we develop a heuristic for stopping an exhaustive compile-time search early if a near-optimal implementation is found. Secondly, we show how to construct run-time decision rules, based on run-time inputs, for selecting from among a subset of the best implementations when the space of inputs can be described by continuously varying features. We address both problems by using statistical modeling techniques that exploit the large amount of performance data collected during the search. We demonstrate these methods on actual performance data collected by the PHiPAC tuning system for dense matrix multiply. We close with a survey of recent projects that use or otherwise advocate an empirical search-based approach to code generation and algorithm selection, whether at the level of computational kernels, compiler and run-time systems, or problem-solving environments. Collectively, these efforts suggest a number of possible software architectures for constructing platform-adapted libraries and applications.

/pdf/statistical-models-for-empirical-search-based-performance-59v7lo1l9i.pdf

Statistical Models for Empirical Search-Based Performance Tuning

The performance of both serial and parallel implementations of matrix multiplication is highly sensitive to memory system behavior. False sharing and cache conflicts cause traditional column-major or row-major array layouts to incur high variability in memory system performance as matrix size varies. This paper investigates the use of recursive array layouts to improve performance and reduce variability. Previous work on recursive matrix multiplication is extended to examine several recursive array layouts and three recursive algorithms: standard matrix multiplication and the more complex algorithms of Strassen (1969) and Winograd. While recursive layouts significantly outperform traditional layouts (reducing execution times by a factor of 1.2-2.5) for the standard algorithm, they offer little improvement for Strassen's and Winograd's algorithms. For a purely sequential implementation, it is possible to reorder computation to conserve memory space and improve performance between 10 percent and 20 percent. Carrying the recursive layout down to the level of individual matrix elements is shown to be counterproductive; a combination of recursive layouts down to canonically ordered matrix tiles instead yields higher performance. Five recursive layouts with successively increasing complexity of address computation are evaluated and it is shown that addressing overheads can be kept in control even for the most computationally demanding of these layouts.

Recursive array layouts and fast matrix multiplication

We study the high-performance implementation of the inversion of a Symmetric Positive Definite (SPD) matrix on architectures ranging from sequential processors to Symmetric MultiProcessors to distributed memory parallel computers. This inversion is traditionally accomplished in three “sweeps”: a Cholesky factorization of the SPD matrix, the inversion of the resulting triangular matrix, and finally the multiplication of the inverted triangular matrix by its own transpose. We state different algorithms for each of these sweeps as well as algorithms that compute the result in a single sweep. One algorithm outperforms the current ScaLAPACK implementation by 20-30 percent due to improved load-balance on a distributed memory architecture.

/pdf/families-of-algorithms-related-to-the-inversion-of-a-cx7sd8hvav.pdf

Families of algorithms related to the inversion of a Symmetric Positive Definite matrix

Recursion leads to automatic variable blocking for dense linear algebra algorithms. The recursion transforms LAPACK level-2 algorithms into level3 codes. For this and other reasons recursion usually speeds up the algorithms.

Recursion provides a new, easy and very successful way of programming numerical linear algebra algorithms. Several algorithms for matrix factorization have been implemented and tested. Some of these algorithms are already candidates for the LAPACK library.

Recursion has also been successfully applied to the BLAS (Basic Linear Algebra Subprograms). The ATLAS system (Automatically Tuned Linear Algebra Software) uses a recursive coding of the BLAS.

The Cholesky factorization algorithm for positive definite matrices, LU factorization for general matrices, and LDLT factorization for symmetric indefinite matrices using recursion are formulated in this paper. Performance graphs of our packed Cholesky and LDLT algorithms are presented here.

LAWRA: Linear Algebra with Recursive Algorithms

We present a family of algorithms for computing symmetric rank-revealing VSV decompositions based on triangular factorization of the matrix. The VSV decomposition consists of a middle symmetric matrix that reveals the numerical rank in having three blocks with small norm, plus an orthogonal matrix whose columns span approximations to the numerical range and null space. We show that for semidefinite matrices the VSV decomposition should be computed via the ULV decomposition, while for indefinite matrices it must be computed via a URV-like decomposition that involves hypernormal rotations.

Computing Symmetric Rank-Revealing Decompositions via Triangular Factorization

Stability of the block cyclic reduction

Recursion leads to automatic variable blocking for dense linear‐algebra algorithms. The recursive way of programming algorithms eliminates using BLAS level 2 during the factorization steps. For this and other reasons recursion usually speeds up the algorithms. The Cholesky factorization algorithm for positive definite matrices and LU factorization for general matrices are formulated. Different storage data formats and recursive BLAS are explained in this paper. Performance graphes of packed and recursive Cholesky algorithms are presented.

Lawra – linear algebra with recursive algorithms

Componentwise error analysis for a parallel partitioning algorithm for tridiagonal systems is presented. Bounds on the equivalent perturbations are obtained, depending on three constants. Then bounds on the forward error are presented as well, depending on two types of condition numbers. Estimates on the first constant come directly from the roundoff error analysis of the tridiagonal Gaussian elimination [N. Higham, SIAM J. Matrix Anal. Appl., 11 (1990), pp. 521--530]. In the present paper, the second and third constants are bounded for some special classes of matrices, i.e., diagonally dominant (row or column), symmetric positive definite, M-matrices, and totally nonnegative.
One of the features of the analysis is that the exact forward and backward errors are bounded, not just their first order approximations, with respect to the machine precision. In all the bounds, the linear terms are given separately to show that the terms of higher order are small enough.

Plamen Y. Yalamov

Papers

LAWRA: Linear Algebra with Recursive Algorithms

Computing Symmetric Rank-Revealing Decompositions via Triangular Factorization

Stability of the block cyclic reduction

Lawra – linear algebra with recursive algorithms

On the Stability of a Partitioning Algorithm for Tridiagonal Systems