Umeå University's logo

We present parallel algorithms for triangular periodic Sylvester-type matrix equations, conceptually being the third step of a periodic Bartels-Stewart-like solution method for general periodic Sylvester-type matrix equations based on variants of the periodic Schur decomposition. The presented algorithms are designed and implemented in the framework of the recently developed HPG library SCASY and are based on explicit blocking, 2-dimensional block cyclic data distribution and a wavefront-like traversal of the right hand side matrices. High performance is obtained by rich usage of level 3 BLAS operations. It is also demonstrated how several important key concepts of SCASY regarding communications and the treatment of quasi-triangular coefficient matrices are generalized to the periodic case. Some experimental results from a distributed memory Linux cluster demonstrate are also presented.

This Thesis contains contributions in two different but closely related subfields of Scientific and Parallel Computing which arise in the context of various eigenvalue problems: periodic and parallel eigenvalue reordering and parallel algorithms for Sylvestertype matrix equations with applications in condition estimation.

Many real world phenomena behave periodically, e.g., helicopter rotors, revolving satellites and dynamic systems corresponding to natural processes, like the water flow in a system of connected lakes, and can be described in terms of periodic eigenvalue problems. Typically, eigenvalues and invariant subspaces (or, specifically, eigenvectors) to certain periodic matrix products are of interest and have direct physical interpretations. The eigenvalues of a matrix product can be computed without forming the product explicitly via variants of the periodic Schur decomposition. In the first part of the Thesis, we propose direct methods for eigenvalue reordering in the periodic standard and generalized real Schur forms which extend earlier work on the standard and generalized eigenvalue problems. The core step of the methods consists of solving periodic Sylvester-type equations to high accuracy. Periodic eigenvalue reordering is vital in the computation of periodic eigenspaces corresponding to specified spectra. The proposed direct reordering methods rely on orthogonal transformations and can be generalized to more general periodic matrix products where the factors have varying dimensions and ±1 exponents of arbitrary order.

In the second part, we consider Sylvester-type matrix equations, like the continuoustime Sylvester equation AX −XB =C, where A of size m×m, B of size n×n, and C of size m×n are general matrices with real entries, which have applications in many areas. Examples include eigenvalue problems and condition estimation, and several problems in control system design and analysis. The parallel algorithms presented are based on the well-known Bartels–Stewart’s method and extend earlier work on triangular Sylvester-type matrix equations resulting in a novel software library SCASY. The parallel library provides robust and scalable software for solving 44 sign and transpose variants of eight common Sylvester-type matrix equations. SCASY also includes a parallel condition estimator associated with each matrix equation.

In the last part of the Thesis, we propose parallel variants of the direct eigenvalue reordering method for the standard and generalized real Schur forms. Together with the existing and future parallel implementations of the non-symmetric QR/QZ algorithms and the parallel Sylvester solvers presented in the Thesis, the developed software can be used for parallel computation of invariant and deflating subspaces corresponding to specified spectra and associated reciprocal condition number estimates.

Recently, recursive blocked algorithms for solving triangular one-sided and two-sided Sylvester-type equations were introduced by Jonsson and Kagstrom. This elegant yet simple technique enables an automatic variable blocking that has the potential of matching the memory hierarchies of today's HPC systems. The main parts of the computations are performed as level 3 general matrix multiply and add (GEMM) operations. We extend and apply the recursive blocking technique to solving periodic Sylvester-type matrix equations. Successive recursive splittings are performed on 3-dimensional arrays, where the third dimension represents the periodicity of a matrix equation.

Library software implementing a parallel small-bulge multishift QR algorithm with Aggressive Early Deflation (AED) targeting distributed memory high-performance computing systems is presented. Starting from recent developments of the parallel multishift QR algorithm [Granat et al., SIAM J. Sci. Comput. 32(4), 2010], we describe a number of algorithmic and implementation improvements. These include communication avoiding algorithms via data redistribution and a refined strategy for balancing between multishift QR sweeps and AED. Guidelines concerning several important tunable algorithmic parameters are provided. As a result of these improvements, a computational bottleneck within AED has been removed in the parallel multishift QR algorithm. A performance model is established to explain the scalability behavior of the new parallel multishift QR algorithm. Numerous computational experiments confirm that our new implementation significantly outperforms previous parallel implementations of the QR algorithm.

We continue our presentation of parallel ScaLAPACK-style algorithms for solving Sylvester-type matrix equations. In Part II, we present SCASY (SCAlable SYlvester solvers), a state-of-the-art HPC software library for solving 44 sign and transpose variants of eight common standard and generalized Sylvester-type matrix equations.

A direct method for eigenvalue reordering in a product of a K-periodic matrix sequence in periodic or extended periodic real Schur form is presented and analyzed. Each reordering of two adjacent sequences of diagonal blocks is performed tentatively to guarantee backward stability and involves solving a K-periodic Sylvester equation (PSE) and constructing a K-periodic sequence of orthogonal transformation matrices. An error analysis of the direct reordering method is presented, and results from computational experiments confirm the stability and accuracy of the method for well-conditioned as well as ill-conditioned problems. These include matrix sequences with fixed and time-varying dimensions, and sequences of small and large periodicity.

We discuss parallel algorithms for solving eight common standard and generalized triangular Sylvester-type matrix equation. Our parallel algorithms are based on explicit blocking, 2D block-cyclic data distribution of the matrices and wavefront-like traversal of the right hand side matrices while solving small-sized matrix equations at different nodes and updating the rest of the right hand side using level 3 operations. We apply the triangular solvers in condition estimation, developing parallel sep(-1)-estimators. Some experimental results are presented.

Parallel ScaLAPACK-style algorithms for solving eight common standard and generalized Sylvester-type matrix equations and various sign and transposed variants are presented. All algorithms are blocked variants based on the Bartels--Stewart method and involve four major steps: reduction to triangular form, updating the right-hand side with respect to the reduction, computing the solution to the reduced triangular problem, and transforming the solution back to the original coordinate system. Novel parallel algorithms for solving reduced triangular matrix equations based on wavefront-like traversal of the right-hand side matrices are presented together with a generic scalability analysis. These algorithms are used in condition estimation and new robust parallel sep^− 1-estimators are developed. Experimental results from three parallel platforms, including results from a mixed OpenMP/MPI platform, are presented and analyzed using several performance and accuracy metrics. The analysis includes results regarding general and triangular parallel solvers as well as parallel condition estimators.