Stef Graillat

This article presents an implementation of the power method, enhanced with the ExBLAS library, for accurately computing the dominant eigenvalue of a matrix. The power method is a fundamental algorithm with widespread applications, notably in PageRank for web search and network analysis. Traditional implementations of the power method can suffer from numerical inaccuracies due to the non-associativity in floating-point operations. These inaccuracies are significantly amplified in scenarios with ill-conditioned matrices, such as those with very closely spaced eigenvalues, or in the presence of small perturbations that can significantly affect the computed results. To address these challenges, we leverage ExBLAS to provide bitwise reproducible results, ensuring the reliability and consistency of the computed eigenvalue. This work details the integration of ExBLAS into the power method, analyzes its impact on accuracy, and demonstrates its effectiveness using benchmark matrices from SuiteSparce collection. The reproducible power method offers consistent and reliable results across diverse computing environments, which is critical for applications where eigenvalue calculation is important.

In this article, we derive accurate Horner methods in real and complex floating-point arithmetic. In particular, we show that these methods are as accurate as if computed in k-fold precision, where k ≥ 2, and then rounded into the working precision. When k = 2, our methods are comparable or faster than the existing compensated Horner routines. When compared to multi-precision software, such as MPFR and MPC, our methods are significantly faster for k ≤ 8, that is, up to 489 bits in the significand.

Krylov subspace algorithms are important methods for solving linear systems. In order to efficiently solve large-scale linear systems, parallelism techniques are often applied. However, parallelism often enlarge the non-associativity of floating-point operations, which can lead to non-reproducibility of the computations. This paper compares the performance of the parallel preconditioned BiCGSTAB algorithm implemented with two different libraries (ExBLAS and ReproBLAS) that can ensure the reproducibility of computations. To address the effect of the compiler, we explicitly utilize the FMA instructions. Finally, numerical experiments show that based on two BLAS implementations, the BiCGSTAB algorithms are reproducible. By contrast, the BiCGSTAB algorithm based on ExBLAS is more accurate but more time-consuming than the one based on ReproBLAS.

In this article, we use the complex compensated Horner method to derive a compensated Ehrlich-Aberth method for the accurate computation of all roots, real or complex, of a polynomial. In particular, under suitable conditions, we prove that the limiting accuracy for the compensated Ehrlich-Aberth iterations is as accurate as if computed in twice the working precision and then rounded to the working precision. Moreover, we derive a running error bound for the complex compensated Horner method and use it to form robust stopping criteria for the compensated Ehrlich-Aberth iterations. Finally, extensive numerical experiments illustrate that the backward and forward errors of the root approximations computed via the compensated Ehrlich-Aberth method are similar to those obtained with a quadruple precision implementation of the Ehrlich-Aberth method with a significant speed-up in terms of computation time.

Numerical validation enables one to improve the reliability of numerical computations that rely upon floating-point operations through obtaining trustful results. Discrete Stochastic Arithmetic (DSA) makes it possible to validate the accuracy of floating-point computations using random rounding. However, it may bring a large performance overhead compared with the standard floating-point operations. In this article, we show that with perturbed data it is possible to use standard floating-point arithmetic instead of DSA for the purpose of numerical validation. For instance, for codes including matrix multiplications, we can directly utilize the matrix multiplication routine (GEMM) of level-3 BLAS that is performed with standard floating-point arithmetic. Consequently, we can achieve a significant performance improvement by avoiding the performance overhead of DSA operations as well as by exploiting the speed of highly-optimized BLAS implementations. Finally, we demonstrate the performance gain using Intel MKL routines compared against the DSA version of BLAS routines.