Efficient parallel reduction to bidiagonal form

Numerical Algorithms, 1995

Journal of Computational and Applied Mathematics, 1989

Solving dense symmetric eigenvalue problems and computing singular value decompositions continue to be two of the most dominating tasks in numerous scientific applications. With the advent of multiprocessor computer systems, the design of efficient parallel algorithms to determine these solutions becomes of paramount importance. In this paper, we discuss two fundamental approaches for determining the eigenvalues and eigenvectors of dense symmetric matrices on machines such as the Alliant FX/8 and CRAY X-MP. One approach capitalizes upon the inherent parallelism offered by Jacobi methods, while the other relies upon an efficient reduction to tridiagonal form via Householder's transformations followed by a multisectioning technique to obtain the eigenvalues and eigenvectors of the corresponding symmetric tridiagonal matrix. For the singular value decomposition, we discuss an efficient method for rectangular matrices in which the number of rows is substantially larger or smaller than the number of columns. This scheme performs an initial orthogonal factorization using block Householder transformation, followed by a parallel one-sided Jacobi method to obtain the singular values and singular vectors of the resulting upper-triangular matrix. Exceptional performance for this SVD scheme is demonstrated for tall matrices of full or deficient rank having clustered or multiple singular values. A hybrid method that combines one-and two-sided Jacobi schemes is also discussed. Performance results for each of the above algorithms on the Alliant FX/8 and CRAY X-MP computer systems will be presented with particular emphasis given to speedups obtained over such classical EISPACK algorithms .

2021

This paper discusses the parallelisation (using the programming language Chapel) of one of the most often used algorithms for the Singular V alue Decomposition of a matrix. The mathematical alterations and programming constructs used to achie ve that parallelisation are examined at length. Performance and validation tests showing the increased parallel performance that occurs when run on a multi-core computer architecture are documented. Chapel’ s serial performance is compared against that of Fortran that implements the original algorithm. Only parallelisation in a Symmetric Multi-Processor environment is explored.

ACM Transactions on Mathematical Software, 2020

The Singular Value Decomposition (SVD) is widely used in numerical analysis and scientific computing applications, including dimensionality reduction, data compression and clustering, and computation of pseudoinverses. In many cases, a crucial part of the SVD of a general matrix is to find the SVD of an associated bidiagonal matrix. This article discusses an algorithm to compute the SVD of a bidiagonal matrix through the eigenpairs of an associated symmetric tridiagonal matrix. The algorithm enables the computation of only a subset of singular values and corresponding vectors, with potential performance gains. The article focuses on a sequential version of the algorithm, and discusses special cases and implementation details. The implementation, called BDSVDX, has been included in the LAPACK library. We use a large set of bidiagonal matrices to assess the accuracy of the implementation, both in single and double precision, as well as to identify potential shortcomings. The results show that BDSVDX can be up to three orders of magnitude faster than existing algorithms, which are limited to the computation of a full SVD. We also show comparisons of an implementation that uses BDSVDX as a building block for the computation of the SVD of general matrices. CCS Concepts: • Mathematics of computing → Solvers; Computations on matrices;

1993

In this paper a new implementation of a divide and conquer algorithm will be considered. This algorithm, in contrast to the LAPACK algorithm, uses a di erent formulation of the update problem, and extended precision in order to maintain accuracy and orthogonality. O u r I n tel Paragon implementation shows, in contrast to the Hypercube implementation by Ipsen and Jessup 14], that good speedups can be obtained from a distributed memory parallel version of the divide and conquer eigenvalue algorithm .

SIAM Journal on Matrix Analysis and Applications, 2007

Two new algorithms for one-sided bidiagonalization are presented. The first is a block version which improves execution time by improving cache utilization from the use of BLAS 2.5 operations and more BLAS 3 operations. The second is adapted to parallel computation. When incorporated into singular value decomposition software, the second algorithm is faster than the corresponding ScaLAPACK routine in most cases. An error analysis is presented for the first algorithm. Numerical results and timings are presented for both algorithms.

SIAM Journal on Scientific Computing, 1999

A new algorithm for the orthogonal reduction of a symmetric matrix to tridiagonal form is developed and analysed. It uses a Cholesky factorization of the original matrix and the rotations are applied to the factors. The idea is similar to the one used for the one-sided Jacobi algorithms B. Zhou and R. Brent, A Parallel Ordering Algorithm for E cient One-Sided Jacobi SVD Computations, Proc. Sixth IASTED-ISMM International Conference on Parallel and Distributed Computing and Systems, pp. 369 372, 1994.. The algorithm uses little communication, accesses data with stride one and is to a large extent independent of data distribution. It has been implemented on the Fujitsu VPP 500. The algorithm is designed to be the rst step of an eigensolver so the procedure for accumulating transforms for eventual calculation of eigenvectors is given.

Procedia Computer Science

A challenging class of problems arising in many GPU applications, called batched problems, involves linear algebra operations on many small-sized matrices. We designed batched BLAS (Basic Linear Algebra Subroutines) routines, and in particular the Level-2 BLAS GEMV and the Level-3 BLAS GEMM routines, to solve them. We proposed device functions and big-tile settings in our batched BLAS design. We adopted auto-tuning to optimize different instances of GEMV routines. We illustrated our batched BLAS approach to optimize batched bi-diagonalization progressively on a K40c GPU. The optimization techniques in this paper are applicable to the other two-sided factorizations as well.

SIAM Journal on Matrix Analysis and Applications, 1997

In this paper we present a new parallel algorithm for computing the LL T decomposition of real symmetric positive de nite tridiagonal matrices. The algorithm consists of a preprocessing and a factoring stage. In the preprocessing stage it determines a rank-(p?1) correction to the original matrix (p = number of processors) by precomputing selected components x k of the L factor, k = 1; : : : ; p ? 1. In the factoring stage it performs independent factorizations of p matrices of order n=p. The algorithm is especially suited for machines with both vector and processor parallelism, as con rmed by the experiments carried out on a Connection Machine CM5 with 32 nodes. Letx k andx 0 k denote the components computed in the preprocessing stage and the corresponding values (re)computed in the factorization stage, respectively. Assuming that jx k =x 0 k j is small, k = 1; : : : ; p ? 1, we are able to prove that the algorithm is stable in the backward sense. The above assumption is justi ed both experimentally and theoretically. In fact we found experimentally that jx k =x 0 k j is small even for illconditioned matrices, and we have proved by an a priori analysis that Research produced with the help of the National Science Foundation, Infrastructure Grant number CDA-8722788, and the Cooperation Agreement CNR-MOSA.

Linear Algebra and its Applications, 2005

A new bidiagonal reduction method is proposed for X ∈ R m×n. For m n, it decomposes X into the product X = UBV T where U ∈ R m×n has orthonormal columns, V ∈ R n×n is orthogonal, and B ∈ R n×n is upper bidiagonal. The matrix V is computed as a product of Householder transformations. The matrices U and B are constructed using a recurrence. If U is desired from the computation, the new procedure requires fewer operations than the Golub-Kahan procedure [SIAM J. Num. Anal. Ser. B 2 (1965) 205] and similar procedures. In floating point arithmetic, the columns of U may be far from orthonormal, but that departure from orthonormality is structured. The application of any backward stable singular value decomposition procedure to B recovers the left singular vectors associated with the leading (largest) singular values of X to near orthogonality. The singular values of B are those of X perturbed by no more than f (m, n)ε M X F where f (m, n) is a modestly growing function and ε M is the machine unit. Under certain assumptions, relative error bounds on the singular values are possible.