Each of the generalized Lyapunov matrix equations can be seen simply as a system of linear equations in the elements of the unknown matrix X. For numerical solutions, however, it is usually kept in its matrix form. The continuous time equation, A’XE + E’XA = Y, for example, can be reduced to a quasi-triangular form by solving the generalized eigenvalue problem, Au = λ Eu and then using the matrix of eigenvectors. The A and E are real square matrices, and Y is real symmetric. In the reduced form A and E are triangular except for 2-by-2 blocks on the diagonal of E representing pairs of complex eigenvalues. This form allows for a variety of solution methods.

The algorithms discussed in this paper are derived from the work of Bartels and Stewart [1] with generalization by Penzl [2] and have a common initial phase. This phase begins by applying the above-mentioned reduction. Then, given a nominal block size, p>2, E is partitioned more or less uniformly so that each diagonal block is of size p, p + 1, or p - 1 with the size adjusted if necessary so that each 2-by-2 block is contained in a larger one.

The first phase now applies the block structure to the entire reduced matrix equation and then solves for the blocks of X by a forward block substitution scheme. This scheme requires, at each step, the solution of a Sylvester matrix equation of order p more or less.

The authors have implemented this scheme in three ways using calls to level-3 basic linear algebra subprograms (BLAS), and using different algorithms for solving the Sylvester equations. They used an Intel Fortran 90 compiler linked with an Intel LAPACK library for testing. The test system is described in detail with numerical results and times. Tests were run for a variety of problem sizes and block sizes as well as single and multithread calls to the BLAS. They also ran tests showing the speed-up to be gained by controlling the memory alignment of the arrays. The tests used randomly generated problems as well as a set problems with known solutions, all beginning in quasi-triangular form.

The algorithms, including two of the three methods for the Sylvester equations, are given in detailed pseudocode. There is no Fortran code included in the paper, but the pseudocode together with the level-3 BLAS should translate into a clear and simple Fortran 90 program.