An efficient parallel version of the householder-QL matrix diagonalisation algorithm

Abstract

In this paper we report an effective parallelisation of the Householder routine for the reduction of a real symmetric matrix to tri-diagonal form and the QL algorithm for the diagonalisation of the resulting matrix. The Householder algorithm scales like $\text{αN}{}^3\text{/P+βN}{}^2\text{log}{}_2\text{(P)}$ and the QL algorithm like γN²+δN³/P as the number of processors P is increased for fixed problem size. The constant parameters α,β,γ and δ are obtained empirically. When the eigenvalues only are required the Householder method scales as above while the QL algorithm remains sequential. The code is implemented in c in conjunction with the message passing interface (MPI) libraries and verified on a sixteen node IBM SP2 and for real matrices that occur in the simulation of properties of crystaline materials.

Introduction

The parallelisation of methods for finding the eigenvalues and eigenvectors for real symetric matrices has received little attention over the years, for two reasons. Firstly for very large systems quite often only a few eigenvalues and eigenvectors are needed. Secondly, the parallelisation of the best sequential techniques for matrix diagonalisation has proved difficult and ineffective.

In the simulation of the dynamics of large structures, only the lowest few eigenvalues that correspond to the most likely occuring low frequency modes of vibration are needed. Similarly, in solid state physics, states at or near the ground state are those most likely to be occupied, at least at low temperatures, and so make the largest contribution to the material properties. For these reasons, in subjects like structural mechanics, techniques have concentrated in reducing the degrees of freedom in the system, from perhaps millions to a few hundred, enabling classical sequential techniques to be used. In Physics and Chemistry, and where the original matrix is regularly structured and sparse, iterative techniques like the Lanczos [9] or Davidson [2] methods, that recover the lowest few eigenvalues and their corresponding eigenvectors can be used effectively.

There are however many problems for which it is desirable to have all of the eigenvalues and eigenvectors of systems. Physical simulations of quantum mechanical systems that have a few particles present at finite temperature cannot be properly analysed by methods for large numbers of particles like statistical mechanics and are sufficiently far from their ground states as to be influenced by all available states. Similarly in the design of mechanical structures acoustic resonances might occur over a broad bandwidth. In such cases all eigenvalues and eigenvectors of large symmetric matrices are needed, for which no algorithm of complexity less than O(N³), where N is the dimension of the matrix to be diagonalised, exists, regardless of the sparsity of the matrix.

The most efficient algorithm for the diagonalisation of large symmetric matrices is the Householder tri-diagonalisation procedure followed by the QL diagonalisation algorithm. The operation count for the combined Householder method is approximately N³ whether or not eigenvectors are required and the operation count for the QL method is roughly 30N² when the eigenvectors are not required and 3N³ when they are. These counts presume that the average number of iterations is about 1.7 [10], an assumption that is confirmed for the test matrices that we use. The only viable alternative is the Jacobi method which has about 20N³ operations when eigenvectors are not required and about 40N³ when the are assuming less that 10 Jacobi cycles for convergence [10]. The Jacobi method might readily parallelise [11] but requires more than 20 processors before becoming competitive with the sequential Householder-QL algorithm. Previous attempts at parallelising the Householder and QL algorithms have resulted in inefficient solutions [4] or have altered the algorithms in some way 12, 3, 1, 8