Bo Kågström

The generalized Schur decomposition of an arbitrary pencil A–λB—robust software with error bounds and applications. Part I: theory and algorithms

Robust software with error bounds for computing the generalized Schur decomposition of an arbitrary matrix pencil A – λB (regular or singular) is presented. The decomposition is a generalization of the Schur canonical form of A – λI to matrix pencils and reveals the Kronecker structure of a singular pencil. Since computing the Kronecker structure of a singular pencil is a potentially ill-posed problem, it is important to be able to compute rigorous and reliable error bounds for the computed features. The error bounds rely on perturbation theory for reducing subspaces and generalized eigenvalues of singular matrix pencils. The first part of this two-part paper presents the theory and algorithms for the decomposition and its error bounds, while the second part describes the computed generalized Schur decomposition and the software, and presents applications and an example of its use.

GEMM-based level 3 BLAS: high-performance model implementations and performance evaluation benchmark

The level 3 Basic Linear Algebra Subprograms (BLAS) are designed to perform various matrix multiply and triangular system solving computations. Due to the complex hardware organization of advanced computer architectures the development of optimal level 3 BLAS code is costly and time consuming. However, it is possible to develop a portable and high-performance level 3 BLAS library mainly relying on a highly optimized GEMM, the routine for the general matrix multiply and add operation. With suitable partitioning, all the other level 3 BLAS can be defined in terms of GEMM and a small amount of level 1 and level 2 computations. Our contribution is twofold. First, the model implementations in Fortran 77 of the GEMM-based level 3 BLAS are structured to reduced effectively data traffic in a memory hierarchy. Second, the GEMM-based level 3 BLAS performance evaluation benchmark is a tool for evaluating and comparing different implementations of the level 3 BLAS with the GEMM-based model implementations.

Recursive blocked algorithms and hybrid data structures for dense matrix library software

Matrix computations are both fundamental and ubiquitous in computational science and its vast application areas. Along with the development of more advanced computer systems with complex memory hierarchies, there is a continuing demand for new algorithms and library software that efficiently utilize and adapt to new architecture features. This article reviews and details some of the recent advances made by applying the paradigm of recursion to dense matrix computations on todays memory-tiered computer systems. Recursion allows for efficient utilization of a memory hierarchy and generalizes existing fixed blocking by introducing automatic variable blocking that has the potential of matching every level of a deep memory hierarchy. Novel recursive blocked algorithms offer new ways to compute factorizations such as Cholesky and QR and to solve matrix equations. In fact, the whole gamut of existing dense linear algebra factorization is beginning to be reexamined in view of the recursive paradigm. Use of recursion has led to using new hybrid data structures and optimized superscalar kernels. The results we survey include new algorithms and library software implementations for level 3 kernels, matrix factorizations, and the solution of general systems of linear equations and several common matrix equations. The software implementations we survey are robust and show impressive performance on todays high performance computing systems.

A Geometric Approach to Perturbation Theory of Matrices and Matrix Pencils. Part I: Versal Deformations

We derive versal deformations of the Kronecker canonical form by deriving the tangent space and orthogonal bases for the normal space to the orbits of strictly equivalent matrix pencils. These deformations reveal the local perturbation theory of matrix pencils related to the Kronecker canonical form. We also obtain a new singular value bound for the distance to the orbits of less generic pencils. The concepts, results, and their derivations are mainly expressed in the language of numerical linear algebra. We conclude with experiments and applications.

Recursive blocked algorithms for solving triangular systems—Part I: one-sided and coupled Sylvester-type matrix equations

Triangular matrix equations appear naturally in estimating the condition numbers of matrix equations and different eigenspace computations, including block-diagonalization of matrices and matrix pairs and computation of functions of matrices. To solve a triangular matrix equation is also a major step in the classical Bartels--Stewart method for solving the standard continuous-time Sylvester equation (AX − XB = C). We present novel recursive blocked algorithms for solving one-sided triangular matrix equations, including the continuous-time Sylvester and Lyapunov equations, and a generalized coupled Sylvester equation. The main parts of the computations are performed as level-3 general matrix multiply and add (GEMM) operations. In contrast to explicit standard blocking techniques, our recursive approach leads to an automatic variable blocking that has the potential of matching the memory hierarchies of todays HPC systems. Different implementation issues are discussed, including when to terminate the recursion, the design of new optimized superscalar kernels for solving leaf-node triangular matrix equations efficiently, and how parallelism is utilized in our implementations. Uniprocessor and SMP parallel performance results of our recursive blocked algorithms and corresponding routines in the state-of-the-art libraries LAPACK and SLICOT are presented. The performance improvements of our recursive algorithms are remarkable, including 10-fold speedups compared to standard algorithms.

Generalized Schur methods with condition estimators for solving the generalized Sylvester equation

Stable algorithms are presented for solving the generalized Sylvester equation. They are based on orthogonal equivalence transformations of the original problem. Perturbation theory and rounding error analysis are included. Condition estimators (dif/sup -1/-estimators) are developed which when substituted into derived error bounds give accuracy estimates of a computed solution. Results from numerical experiments on well-conditioned and ill-conditioned problems are reported. >

An Algorithm for Numerical Computation of the Jordan Normal Form of a Complex Matrix

An algorithm for the computation of the Jordan normal form of a complex square matrix is given. The definition of the Jordan normal form is modified in order to be applicable when working in finiteprecision arithmetic. It is then shown how an accurate and stable algorithm, which computes eigenvalue approximations and chains of principal vectors, can be constructed in this case. The • algorithm is based on a sequence of similarity transformations and successive range-nullspace separations, following a suggestion by Kublanovskaya. It is shown how tolerance parameters in the algorithm should be chosen and how the results of the algorithm should be interpreted and evaluated. This is illustrated by a few numerical examples.

A Geometric Approach to Perturbation Theory of Matrices and Matrix Pencils. Part II: A Stratification-Enhanced Staircase Algorithm

Computing the Jordan form of a matrix or the Kronecker structure of a pencil is a well-known ill-posed problem. We propose that knowledge of the closure relations, i.e., the stratification, of the orbits and bundles of the various forms may be applied in the staircase algorithm. Here we discuss and complete the mathematical theory of these relationships and show how they may be applied to the staircase algorithm. This paper is a continuation of our Part I paper on versal deformations, but it may also be read independently.

The generalized Schur decomposition of an arbitrary pencil A–λB—robust software with error bounds and applications. Part II: software and applications

Robust software with error bounds for computing the generalized Schur decomposition of an arbitrary matrix pencil A – λB (regular or singular) is presented. The decomposition is a generalization of the Schur canonical form of A – λI to matrix pencils and reveals the Kronecker structure of a singular pencil. The second part of this two-part paper describes the computed generalized Schur decomposition in more detail and the software, and presents applications and an example of its use. Background theory and algorithms for the decomposition and its error bounds are presented in Part I of this paper.

Recursive Blocked Data Formats and BLAS's for Dense Linear Algebra Algorithms

Recursive blocked data formats and recursive blocked BLAS’s are introduced and applied to dense linear algebra algorithms that are typified by LAPACK. The new data formats allow for maintaining data locality at every level of the memory hierarchy and hence providing high performance on today’s memory tiered processors. This new data format is hybrid. It contains blocking parameters which are chosen so that the associated submatrices of a block-partitioned A fir into level 1 cache. The recursive part of the data format chooses a linear order of the blocks that maintains a two-dimensional data locality of A in a one-dimensional tiered memory structure. We argue that, out of the NB factorial choices of ordering the NB blocks, our recursive ordering leads to one of the best. This is because our algorithms are also recursive and will do their computations on submatrices that follow the new recursive data structure definition. This is in analogy with the well known principle that the data structure should be matched to the algorithm. Performance results in support for our recursive approach are also presented.