Benefits of IEEE‐754 Features in Modern Symmetric Tridiagonal Eigensolvers

Abstract

Abstract. Bisection is one of the most common methods used to compute the eigenvalues of symmetric tridiagonal matrices. Bisection relies on the Sturm count: for a given shift σ, the number of negative pivots in the factorization T − σI = LDL T equals the number of eigenvalues of T that are smaller than σ. In IEEE-754 arithmetic, the value ∞ permits the computation to continue past a zero pivot, producing a correct Sturm count when T is unreduced. Demmel and Li showed in the 90s that using ∞ rather than testing for zero pivots within the loop could improve performance significantly on certain architectures.When eigenvalues are to be computed to high relative accuracy, it is often preferable to work with LDL T factorizations instead of the original tridiagonal T . One important example is the MRRR algorithm. When bisection is applied to the factored matrix, the Sturm count is computed from LDL T which makes differential stationary and progressive qds algorithms the methods of choice.While it seems trivial to replace T by LDL T , in reality these algorithms are more complicated: in IEEE-754 arithmetic, a zero pivot produces an overflow followed by an invalid exception (NaN) that renders the Sturm count incorrect.We present alternative, safe formulations that are guaranteed to produce the correct result. Benchmarking these algorithms on a variety of platforms shows that the original formulation without tests is always faster provided no exception occurs. The transforms see speed-ups of up to 2.6× over the careful formulations.Tests on industrial matrices show that encountering exceptions in practice is rare. This leads to the following design: First, compute the Sturm count by the fast but unsafe algorithm. Then, if an exception occurred, recompute the count by a safe, slower alternative.The new Sturm count algorithms improve the speed of bisection by up to 2× on our test matrices. Furthermore, unlike the traditional tiny-pivot substitution, proper use of IEEE-754 features provides a careful formulation that imposes no input range restrictions.

High-Performance Solvers for Dense Hermitian Eigenproblems

“…As a second step, bisection is used to compute approximations to each eigenvalueλ j . At the cost of O(n) flops, bisection guarantees that each eigenvalue is computed with high relative accuracy [31].…”

Proposed Consistent Exception Handling for the BLAS and LAPACK

“…SSTEMR computes selected eigenvalues, and optionally eigenvectors, of a symmetric tridiagonal matrix T. One internal operation is counting the number of eigenvalues of T that are < s, for various values of s. Letting D be the array of diagonal entries of T, and E be the array of offdiagonal entries, the inner loop (which appears in SLANEG) that does the counting looks roughly like this (an analysis is provided elsewhere [33]):…”

“…As described in the symmetric tridiagonal eigensolvers report [33], it is possible for a tiny DPLUS to cause T to overflow to Inf, which makes the next DPLUS equal to Inf, which makes the next T = Inf/Inf = NaN, which then continues to propagate. Checking for this rare event in the inner loop would be expensive, so SLANEG only checks for T being a NaN every 128 iterations, yielding significant speedups in the most common cases.…”

Algorithm 880

“…(1) to experiment with and select algorithmic variants as in Marques et al [2006], (2) to experiment with parameter and threshold choices such as those to be set in LAPACK's ilaenv, (3) to carry out large scale performance comparisons such as Marques et al [2006]; Demmel et al [2006], and (4) to validate and benchmark new algorithms such as Matsekh [2005].…”