Parallel reduction of four matrices to condensed form for a generalized matrix eigenvalue algorithm

The VZ algorithm proposed by Charles F. Van Loan (SIMA, 1975) attempts to solve the generalized type of matrix eigenvalue problem A C x = λ B D x , where A , B ∈ R n × m , C , D ∈ R m × n , and m ≥ n , without forming products and inverses. Especially, this algorithm is suitable for solving the generalized singular value problem. Van Loan’s approach first reduces the matrices A , B , C , and D to a condensed form by the finite step initial reduction. The reduction finds orthogonal matrices Q , U , V , and Z , such that QAZ is upper Hessenberg, and QBV , Z T C U , and V T D U are upper triangular. In this initial reduction, A is reduced to upper Hessenberg form, while simultaneously preserving triangularity of other three matrices. This is done by Givens rotations, annihilating one by one element of A , and by generating three more rotations applied to other matrices per each annihilation. Such an algorithm is quite inefficient. In our work, we propose a blocked algorithm for the initial reduction, based on aggregated Givens rotations and matrix–matrix multiplications, which are applied in the outer loop updates. This algorithm has another level of blocking, exploited in the inner loop. Further, we also consider a variant of the algorithm in a hybrid CPU–GPU framework, where compute-intensive outer loop updates are performed on GPU, and can be overlapped with the reduction in the next step performed on CPU. On the other hand, application of a sequence of rotations in the inner loop is parallelized on CPU, with balanced operation count per thread. Since a large number of aggregated rotations are produced in every outer loop step, they are simultaneously accumulated before outer loop updates. These adjustments speed up original initial reduction considerably which is confirmed by numerical experiments, and the efficiency of the whole VZ algorithm is increased.