Tridiagonal matrix algorithm

Abstract: The Householder reduction of a member of the anti-symmetric Gaussian unitary ensemble gives an anti-symmetric tridiagonal matrix with all independent elements. The random variables permit the introduction of a positive parameter $\beta$, and the eigenvalue probability density function of the corresponding random matrices can be computed explicitly, as can the distribution of $\{q_i\}$, the first components of the eigenvectors. Three proofs are given. One involves an inductive construction based on bordering of a family of random matrices which are shown to have the same distributions as the anti-symmetric tridiagonal matrices. This proof uses the Dixon-Anderson integral from Selberg integral theory. A second proof involves the explicit computation of the Jacobian for the change of variables between real anti-symmetric tridiagonal matrices, its eigenvalues and $\{q_i\}$. The third proof maps matrices from the anti-symmetric Gaussian $\beta$-ensemble to those realizing particular examples of the Laguerre $\beta$-ensemble. In addition to these proofs, we note some simple properties of the shooting eigenvector and associated Prufer phases of the random matrices.

Abstract: Graphical Processing Units (GPUs) potentially promise widespread and inexpensive high performance computation. However, architectural limitations (only some operations and memory access patterns can be performed quickly, partial support for IEEE floating point arithmetic) make it necessary to change existing algorithms to attain high performance and correctness. Here we show how to make the bisection algorithm for eigenvalues of symmetric tridiagonal matrices (sstebz from LAPACK) run both fast and correctly on an ATI Radeon X1900 GPU. Our fastest algorithm takes up to 156! less time than Intel's Math Kernel Library version of sstebz running on the CPU, but does so by doing many redundant floating point operations compared to the CPU version. We use an automatic tuning procedure analogous to ATLAS or PHiPAC to decide the optimal redundancy. Correctness despite partial IEEE floating point semantics required explicitly adding 0 in the inner loop. The problems and solutions discussed here are of interest on other GPU architectures. 1 Motivation and Objectives Modern graphics processors (GPUs) are data parallel architectures that can run general-purpose computations in single precision (so far) at high computational rates. They are capable of achieving 110 GFLOPS in matrix-matrix multiplication [Segal and Peercy 2006] and show 30-40x speedups compared to the recent Intel Xeon processors in computationally intensive applications such as Black-Scholes option pricing [McCool et al. 2006] and gas dynamics solvers [Hagen et al. 2007]. It is tempting to exploit this computational power in solving other common numerical problems. In this work we consider an implementation of another widely used linear algebra routine — the bisection algorithm for finding the eigenvalues of symmetric tridiagonal matrices. A numerically robust, vectorized implementation of this algorithm in single precision is available in LAPACK’s sstebz routine [Anderson et al. 1999]. Our goal is to port the vectorized segments of the code to the GPU. In order to increase the utilization of the parallel resources, we use the Multi-section with Multiple Eigenvalues method used previously by Katagiri et al. [2006]. For the purpose of this study we restrict our attention to finding all eigenvalues of the matrix. The extension to finding a subset of the eigenvalues as done in LAPACK’s sstebz routine, is straightforward. 2 The Bisection Algorithm