Purpose
To solve the Sylvester equation -AX + XB = C, where A and B are complex M-by-M and N-by-N matrices, respectively, in Schur form. This routine is intended to be called only by SLICOT Library routine MB03RZ. For efficiency purposes, the computations are aborted when the absolute value of an element of X is greater than a given value PMAX.Specification
SUBROUTINE MB03RW( M, N, PMAX, A, LDA, B, LDB, C, LDC, INFO ) C .. Scalar Arguments .. INTEGER INFO, LDA, LDB, LDC, M, N DOUBLE PRECISION PMAX C .. Array Arguments .. COMPLEX*16 A(LDA,*), B(LDB,*), C(LDC,*)Arguments
Input/Output Parameters
M (input) INTEGER The order of the matrix A and the number of rows of the matrices C and X. M >= 0. N (input) INTEGER The order of the matrix B and the number of columns of the matrices C and X. N >= 0. PMAX (input) DOUBLE PRECISION An upper bound for the absolute value of the elements of X (see METHOD). A (input) COMPLEX*16 array, dimension (LDA,M) The leading M-by-M upper triangular part of this array must contain the matrix A of the Sylvester equation. The elements below the diagonal are not referenced. LDA INTEGER The leading dimension of array A. LDA >= MAX(1,M). B (input) COMPLEX*16 array, dimension (LDB,N) The leading N-by-N upper triangular part of this array must contain the matrix B of the Sylvester equation. The elements below the diagonal are not referenced. LDB INTEGER The leading dimension of array B. LDB >= MAX(1,N). C (input/output) COMPLEX*16 array, dimension (LDC,N) On entry, the leading M-by-N part of this array must contain the matrix C of the Sylvester equation. On exit, if INFO = 0, the leading M-by-N part of this array contains the solution matrix X of the Sylvester equation, and each element of X (see METHOD) has the absolute value less than or equal to PMAX. On exit, if INFO = 1, the solution matrix X has not been computed completely, because an element of X had the absolute value greater than PMAX. Part of the matrix C has possibly been overwritten with the corresponding part of X. LDC INTEGER The leading dimension of array C. LDC >= MAX(1,M).Error Indicator
INFO INTEGER = 0: successful exit; = 1: an element of X had the absolute value greater than the given value PMAX. = 2: A and B have common or very close eigenvalues; perturbed values were used to solve the equation (but the matrices A and B are unchanged). This is a warning.Method
The routine uses an adaptation of the standard method for solving Sylvester equations [1], which controls the magnitude of the individual elements of the computed solution [2]. The equation -AX + XB = C can be rewritten as m l-1 -A X + X B = C + sum A X - sum X B kk kl kl ll kl i=k+1 ki il j=1 kj jl for l = 1:n, and k = m:-1:1, where A , B , C , and X , are the kk ll kl kl elements defined by the partitioning induced by the Schur form of A and B. So, the elements of X are found column by column, starting from the bottom. If any such element has the absolute value greater than the given value PMAX, the calculations are ended.References
[1] Bartels, R.H. and Stewart, G.W. T Solution of the matrix equation A X + XB = C. Comm. A.C.M., 15, pp. 820-826, 1972. [2] Bavely, C. and Stewart, G.W. An Algorithm for Computing Reducing Subspaces by Block Diagonalization. SIAM J. Numer. Anal., 16, pp. 359-367, 1979.Numerical Aspects
2 2 The algorithm requires 0(M N + MN ) operations.Further Comments
Let ( A C ) ( I X ) M = ( ), Y = ( ). ( 0 B ) ( 0 I ) Then -1 ( A 0 ) Y M Y = ( ), ( 0 B ) hence Y is a non-unitary transformation matrix which performs the reduction of M to a block-diagonal form. Bounding a norm of X is equivalent to setting an upper bound to the condition number of the transformation matrix Y.Example
Program Text
NoneProgram Data
NoneProgram Results
None