| # SCASY - Parallel ScaLAPACK-style High Performance Library for Sylvester-Type Matrix Equations Robert Granat and Bo Kågström
### Introduction Matrix equations arise in a number of applications in science and engineering, for example in eigenvalue problems, condition estimation and control system design and analysis. SCASY is a HPC parallel ScaLAPACK-style library for solving 44 sign and transpose variants of eight common general Sylvester-type matrix equations. SCASY implements Bartels-Stewart's-like methods for the considered matrix equations which include four major steps: Reduction to triangular forms Update of the right hand side matrix or matrices with respect to the reduction(s) in step 1 Computing the solution of the triangular matrix equation Transforming the solution to the triangular problem back to the original coordinate system The reductions in step 1 and the transformations in step 2 and 4 are performed by orthogonal transformations only. Step 1 is carried out by using existing software in ScaLAPACK or new preliminary developments of fast and reliable algorithms for unsymmetric eigenvalue problems (parallel variants of the unsymmetric QR and QZ algorithms). Steps 2 and 4 are parallel matrix multiplications, which are implemented using PBLAS. In the 3rd step, specially designed parallel wave-front algorithms are employed that compute the solution to the triangular problems in a fixed number of operations; the major contributions from SCASY belong to this kind of algorithms and software. ### Sylvester-Type Matrix Equations The following general Sylvester-type matrix equations are considered in SCASY.
Name | Matrix equation | Acronym
Standard Sylvester (CT) | op(A)X ± Xop(B) = C | SYCT
Standard Lyapunov (CT) | op(A)X + Xop(AT) = C | LYCT
Generalized coupled Sylvester | ( op(A)X ± Yop(B) = C, op(D)X ± Yop(E) = F ) | GCSY
Standard Sylvester (DT) | op(A)Xop(B) ± X = C | SYDT
Standard Lyapunov (DT) | op(A)Xop(AT) - X = C | LYDT
Generalized Sylvester | op(A)Xop(B) ± op(C)Xop(D) = E | GSYL
Generalized Lyapunov (CT) | op(A)Xop(ET) + op(E)Xop(AT) = C | GLYCT
Generalized Lyapunov (DT) | op(A)Xop(AT) - op(E)Xop(ET) = C | GLYDT
Table 1. The Sylvester-type matrix equations considered in SCASY For example, SCASY solves for the matrix X in eight sign and transpose variants of the general standard continuous-time Sylvester equation op(A)X ± Xop(B) = scale C, where 'op()' denotes a possible transpose and scale is a scalar to prevent overflow in the right hand side. Moreover, for all Lyapunov equations, the right hand side and the solution matrices may be symmetric. We differentiate between one-sided and two-sided matrix equations, i.e, matrix equations where the solution is multiplied by another matrix from one or two sides. We also differentiate between standard and generalized matrix equations. In particular, in Table 1, the matrix equations SYCT, LYCT and GCSY are one-sided, while the remaining five are two-sided matrix equations. In addition, the matrix equations GCSY, GSYL, GLYCT and GLYDT are generalized variants of the standard matrix equations SYCT, LYCT, SYDT and LYDT. With the solver routines in SCASY, it is possible to work with submatrices of distributed global matrices, as well. Below we denote the M-by-N submatrix of the global distributed matrix A starting at row and column indicies (IA,JA) by sub(A), where sub(A)=A(IA:IA+M-1,JA:JA+N-1). ### Solver routines The following routines are available in SCASY:
Name | Matrix equation
PGESYCTD | SYCT
PGELYCTD | LYCT
PGESYDTD | SYDT
PGELYDTD | LYDT
PGEGCSYD | GCSY
PGEGSYLD | GSYL
PGEGLYCTD | GLYCT
PGEGLYDTD | GLYDT
Table 2. The solver routines in SCASY In general, the solver PGE[ACRO]D solves for the matrix equation with the corresponding acronym ACRO. The routines calls eight corresponding triangular solvers PTR[ACRO]D which implements the parallel wave-front algorithms decribed above. However, it is not recommended that the user calls these triangular solves directly, since the major part of the error checking is performed by the general solver. The reduction and transformation steps can also be omitted by setting the in parameters to the general routine properly. Below follows the subroutine headings for each routine. SUBROUTINE PGESYCTD( JOB, ASCHUR, BSCHUR, TRANSA, TRANSB, ISGN, COMM, M, N, A, IA, JA, DESCA, B, IB, JB, DESCB, C, IC , JC, DESCC, MBNB2, DWORK, LDWORK, IWORK, LIWORK, NOEXSY, SCALE, INFO ) Solves the continuous-time Sylvester equation (SYCT) where sub(A) is M-by-M, sub(B) is N-by-N and sub(C)/X is M-by-N. SUBROUTINE PGELYCTD( JOB, SYMM, TRANSA, ASCHUR, M, A, IA, JA, DESCA, C, IC, JC, DESCC, MB2, DWORK, LDWORK, IWORK, LIWORK, NOEXSY, SCALE, INFO ) Solves the continuous-time Lyapunov equation (LYCT) where sub(A) is M-by-M and sub(C)/X is M-by-M. SUBROUTINE PGESYDTD( JOB, ASCHUR, BSCHUR, TRANSA, TRANSB, ISGN, COMM, M, N, A, IA, JA, DESCA, B, IB, JB, DESCB, C, IC, JC, DESCC, MBNB2, DWORK, LDWORK, IWORK, LIWORK, NOEXSY, SCALE, INFO ) Solves the discrete-time Sylvester equation (SYDT) where sub(A) is M-by-M, sub(B) is N-by-N and sub(C)/X is M-by-N. SUBROUTINE PGELYDTD( JOB, SYMM, TRANSA, ASCHUR, M, A, IA, JC, DESCA, C, IC, JC, DESCC, MB2, DWORK, LDWORK, IWORK, LIWORK, NOEXSY, SCALE, INFO ) Solves the discrete-time Lyapunov equation (LYDT) where sub(A) is M-by-M and sub(C)/X is M-by-M. SUBROUTINE PGEGCSYD( JOB, TRANZ, ADSCHR, BESCHR, TRANAD, TRANBE, ISGN, COMM, M, N, A, IA, JA, DESCA, B, IB, JB, DESCB, C, IC, JC, DESCC, D, ID, JD, DESCD, E, IE, JE, DESCE, F, IF, JF, DESCF, MBNB2, DWORK, LDWORK, IWORK, LIWORK, NOEXSY, SCALE, INFO ) Solves the generalized coupled Sylvester equation (GCSY) where sub(A) and sub(D) are M-by-M, sub(B) and sub(E) are N-by-N and sub(C)/X and sub(F)/Y are M-by-N. SUBROUTINE PGEGSYLD( JOB, ACSCHR, BDSCHR, TRANAC, TRANBD, ISGN, COMM, M, N, A, IA, JA, DESCA, B, IB, JB, DESCB, C, IC, JC, DESCC, D, ID, JD, DESCD, E, IE, JE, DESCE, MBNB2, DWORK, LDWORK, IWORK, LIWORK, NOEXSY, SCALE, INFO ) Solves the generalized Sylvester equation (GSYL) where sub(A) and sub(C) are M-by-M, sub(B) and sub(D) are N-by-N and sub(E)/X is M-by-N. SUBROUTINE PGEGLYCTD( JOB, SYMM, OP, AESCHR, M, A, IA, JA, DESCA, E, IE, JE, DESCE, C, IC, JC, DESCC, MB2, DWORK, LDWORK, IWORK, LIWORK, NOEXSY, SCALE, INFO ) Solves the continuous-time generalized Lyapunov equation (GLYCT) where sub(A) and sub(E) are M-by-M and sub(C)/X is M-by-M. SUBROUTINE PGEGLYDTD( JOB, SYMM, OP, AESCHR, M, A, IA, JA, DESCA, E, IE, JE, DESCE, C, IC, JC, DESCC, MB2, DWORK, LDWORK, IWORK, LIWORK, NOEXSY, SCALE, INFO ) Solves the discrete-time generalized Lyapunov equation (GLYDT) where sub(A) and sub(E) are M-by-M and sub(C)/X is M-by-M. Each parameter is decribed in the following table.
Mode parameters | Type | Meaning | Remark
JOB | CHARACTER*1 | Specifies solving mode ('S') or reduction mode ('R'). |
SYMM | CHARACTER*1 | Specifies if the right hand side is symmetric ('S') or not ('N'). | 1.
OP | CHARACTER*1 | Specifies the transpose variant of the corresponding Lyapunov equation ('N' or 'T'). | 1.
TRANZ | CHARACTER*1 | Specifies if the transposed Kronecker product representation of the equation should be solved ('T') or not ('N'). Only applies to GCSY. | 2.
_SCHUR | CHARACTER*1 | Specifies if the matrix _ is in real Schur form ('S') or not ('N'). |
__SCHR | CHARACTER*1 | Specifies if the matrix pair (_,_) is in generalized Schur form ('S') or not ('N'). |
TRANS_ | CHARACTER*1 | Specifies the transpose variant of the matrix _ ('N' or 'T'). |
TRAN__ | CHARACTER*1 | Specifies the transpose variant of the matrix pair (_,_) ('N' or 'T'). |
ISGN | INTEGER | Specifies the sign variant of the equation (-1 or 1). |
Input/output parameters |  |  |
COMM | CHARACTER*1 | Sets communication pattern in PTR[ACRO]D. | 3.
M,N | INTEGER | (Sub)matrix dimensions (see above). |
A,B,C,D,E,F | DOUBLE PRECISION(*) | One-dimensional arrays corresponding to the local parts of the globally distributed matrices involved in the actual equation. On output, each right hand side matrix or matrix pair is returned in real or generalized Schur form. On output and in solving mode, the right hand side matrix or matrix pair is overwritten with the solution. In reduction mode, the right hand side is not referenced. |
IA,IB,IC,ID,IE,IF | INTEGER | Row starting indices for submatrices to operate on. | 4.
JA,JB,JC,JD,JE,JF | INTEGER | Column starting indices for submatrices to operate on. | 4.
DESC_ | INTEGER(*) | ScaLAPACK matrix descriptor arrays. |
MBNB2, MB2 | INTEGER(2), INTEGER | Internal block factors used in pipelining. | 5.
Workspace |  |  |
DWORK | DOUBLE PRECISION(*) | Double precision workspace. |
LDWORK | INTEGER | Length of DWORK. |
IWORK | INTEGER(*) | Integer workspace. |
LIWORK | INTEGER | Length of IWORK. |
Output information |  |  |
NOEXSY | INTEGER | Counts the number of extended subsystems solved in PTR[ACRO]D. | 6.
SCALE | DOUBLE PRECISION*1 | Scaling factor, 0 < SCALE <= 1.0 |
Error handling |  |  |
INFO | INTEGER | Return error information to the calling routine, including overflow warning |
Table 3. Parameters in the solver routines Remarks: This argument is only used for the LYCT, LYDT, GLYCT and GLYDT equations. This argument is only used in condition estimation of GCSY. Input/output parameters COMM chooses between two communication schemes in the triangular solver PTR[ACRO]D, "on demand" communication or "matrix shifts" communication. "On demand" is allowed for all combination of mode parameters, whereas "matrix shifts" is only allowed for a subset the combinations. COMM is also output from the routines, showing which scheme that was actually used in the triangular solver. See the references for details. For general problems (non-triangular), SCASY can not work on submatrices – this limitation comes from missing functionality in the utilized ScaLAPACK routine PDLAHQR. This might be changed in a future release. Currently not used for LYDT, GLYCT, GLYDT. The triangular solvers handle splitted 2-by-2 blocks on the diagonals of the Schur forms via an implicit redistribution [1,2] of the involved matrices. Some subsystems will then be extended one row and column and some other subsystems will be diminished. ### Condition estimators SCASY provides condition estimators for the triangular counterparts of the Sylvester-type matrix equations in Table 1. The reduction to (quasi-)triangular form of the actual equation can be conducted by invoking the standard solver routine with the inargument 'JOB' to reduction mode. All condition estimators build on the ScaLAPACK condition estimator PDLACON. In analogy with the naming scheme for the solver routines, the condition estimators go under the names P[ACRO]CON, where ACRO is replaced by the corresponding matrix equation acronym. ### Matrix generators SCASY includes two generator routines for generating matrices and matrix pairs with specified standard and generalized eigenvalues. These can be utilized to construct test examples with varying conditioning in the sense of sep[ACRO]. ### New versions of the reduction routines Future releases of SCASY will include new complete implementations of the unsymmetric multishift QR and QZ algorithms with advanced deflations techniques which are one to two orders of magnitude faster than the existing ScaLAPACK implementation. However, preliminary implementations of the new QR algorithm can be tested in the reduction step to triangular form of the standard equations SYCT, LYCT, SYDT and LYDT, see above, by a specific pre-processing options (see the README file and the Users' Guide in the software package for details). The routines connected to the QZ algorithm are in the current release empty wrappers which perform no computations, but specify the assumed interface of the future production code. ### Prerequisites To use SCASY on your system you need the following. A Fortran 90/95 compiler. We have tested the software with Portland Groups and the Pathscale Compiler Suites, as well as the NAG Fortran Compiler, IBMs XLFortran Compiler, and gfortran. You should have no problem with compiling the software as long as you follow the instructions in the README file and use the proper compilation flags. A good implementation of the BLAS, e.g., ATLAS or the GOTO-BLAS. The LAPACK library. The RECSY library (which also depends on a few routines from SLICOT). The ScaLAPACK library, including PBLAS and BLACS. The following is a subroutine callgraph of the library. ### Download The latest version of the complete library (SCASY 1.0beta) can be downloaded as a tar.gz-file here and as a zip-file right here. Notice that this version is a beta version - for bug reports and related issues, see below. ### Installation instruction s There is currently no installation routine. However, the built library should be moved to a proper directory, e.g., /usr/local/lib. After unpacking, the library is built as follows. Move into the SCASY directory. Then, edit the Make_include file to fit for your system regarding compiler, linker, flags, library paths, pre-processing options etc. (Make_include templates are provided for the Portland Group and the Pathscale Compiler Suites, the GNU Fortran (gfortran) compiler and the NAG Fortran Compiler.) Then try host:~> make all to build the library with a test program. To build only the library try host:~> make libscasy ### Revision history 2007-09-10: First release, 0.10beta. 2009-08-31: Second release, 1.0beta. Please send comments, suggestions and bug reports to granat@cs.umu.se and bokg@cs.umu.se. ### References Robert Granat, A Parallel ScaLAPACK-style Sylvester Solver, Master's Thesis, Report UMNAD 435/03, Dept. Computing Science, Umeå University, Sweden, January, 2003. Robert Granat, Bo Kågström and Peter Poromaa, Parallell ScaLAPACK-style Algorithms for Solving Continous-Time Sylvester Equations. In H. Kosch et al (Eds), Euro-Par 2003 Parallel Processing. Lecture Notes in Computer Science, Springer Verlag, Vol. 2790, pp. 800-809, 2003. Robert Granat and Bo Kågström. Evaluating Parallel Algorithms for Solving Sylvester-Type Matrix Equations: Direct Transformation-Based versus Iterative Matrix-Sign-Function-Based Methods. To appear in PARA'04 State-of-the-Art in Scientific Computing Conference Proceedings, LCNS, Springer Verlag, 2004. Robert Granat, Isak Jonsson and Bo Kågström, Combining Explicit and Recursive Blocking for Solving Triangular Sylvester-Type Matrix Equations on Distributed Memory Platforms. In M. Danelutto, D. Laforenza, M. Vanneschi (Eds), Euro-Par 2004. Lecture Notes in Computer Science, Springer Verlag, Vol. 3149, pp. 742-750, 2004. Robert Granat and Bo Kågström, Parallel Solvers for Sylvester-Type Matrix Equations with Applications in Condition Estimation, Part I: Theory and Algorithms, ACM TOMS, submitted June 2007, revised January and August 2009. Robert Granat and Bo Kågström, ALGORITHM XXX: The SCASY Library -- Parallel Solvers for Sylvester-Type Matrix Equations with Applications in Condition Estimation, Part II, ACM TOMS, submitted June 2007, revised January and August 2009. Robert Granat and Bo Kågström, SUG - SCASY Users' Guide, Tech Report UMINF-09.10, Dept. Computing Science, Umeå University, Sweden, August 2009. ### Acknowledgements This research is financially supported by the Swedish Research Council under grant VR 621-2001-3284 and the Swedish Foundation for Strategic Research under grant A3 02:128 .The research is conducted using the resources of the High Performance Computing Center North (HPC 2N). ### Terms of Usage Authors and main contributors are Robert Granat and Bo Kågström, Department of Computing Science and HPC 2N, Umeå University, SE-901 87 Umeå, Sweden. Direct and indirect co-contributors are Björn Adlerborn, Isak Jonsson and Daniel Kressner (see the references above and the software documentation for details). The library is freely available for academic (non-commercial) use, and is provided on an "as is" basis. Any use of the SCASY library should be acknowledged by citing the corresponding papers [5-6] above and by refering to this homepage. This homepage was last updated 090819 by Robert Granat.