Contents
dtgsyl - solve the generalized Sylvester equation
SUBROUTINE DTGSYL(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD,
E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO)
CHARACTER * 1 TRANS
INTEGER IJOB, M, N, LDA, LDB, LDC, LDD, LDE, LDF, LWORK,
INFO
INTEGER IWORK(*)
DOUBLE PRECISION SCALE, DIF
DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*),
E(LDE,*), F(LDF,*), WORK(*)
SUBROUTINE DTGSYL_64(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO)
CHARACTER * 1 TRANS
INTEGER*8 IJOB, M, N, LDA, LDB, LDC, LDD, LDE, LDF, LWORK,
INFO
INTEGER*8 IWORK(*)
DOUBLE PRECISION SCALE, DIF
DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*),
E(LDE,*), F(LDF,*), WORK(*)
F95 INTERFACE
SUBROUTINE TGSYL(TRANS, IJOB, [M], [N], A, [LDA], B, [LDB], C, [LDC],
D, [LDD], E, [LDE], F, [LDF], SCALE, DIF, [WORK], [LWORK], [IWORK],
[INFO])
CHARACTER(LEN=1) :: TRANS
INTEGER :: IJOB, M, N, LDA, LDB, LDC, LDD, LDE, LDF, LWORK,
INFO
INTEGER, DIMENSION(:) :: IWORK
REAL(8) :: SCALE, DIF
REAL(8), DIMENSION(:) :: WORK
REAL(8), DIMENSION(:,:) :: A, B, C, D, E, F
SUBROUTINE TGSYL_64(TRANS, IJOB, [M], [N], A, [LDA], B, [LDB], C,
[LDC], D, [LDD], E, [LDE], F, [LDF], SCALE, DIF, [WORK], [LWORK],
[IWORK], [INFO])
CHARACTER(LEN=1) :: TRANS
INTEGER(8) :: IJOB, M, N, LDA, LDB, LDC, LDD, LDE, LDF,
LWORK, INFO
INTEGER(8), DIMENSION(:) :: IWORK
REAL(8) :: SCALE, DIF
REAL(8), DIMENSION(:) :: WORK
REAL(8), DIMENSION(:,:) :: A, B, C, D, E, F
C INTERFACE
#include <sunperf.h>
void dtgsyl(char trans, int ijob, int m, int n, double *a,
int lda, double *b, int ldb, double *c, int ldc,
double *d, int ldd, double *e, int lde, double *f,
int ldf, double *scale, double *dif, int *info);
void dtgsyl_64(char trans, long ijob, long m, long n, double
*a, long lda, double *b, long ldb, double *c, long
ldc, double *d, long ldd, double *e, long lde,
double *f, long ldf, double *scale, double *dif,
long *info);
dtgsyl solves the generalized Sylvester equation:
A * R - L * B = scale * C (1)
D * R - L * E = scale * F
where R and L are unknown m-by-n matrices, (A, D), (B, E)
and (C, F) are given matrix pairs of size m-by-m, n-by-n and
m-by-n, respectively, with real entries. (A, D) and (B, E)
must be in generalized (real) Schur canonical form, i.e. A,
B are upper quasi triangular and D, E are upper triangular.
The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an
output scaling factor chosen to avoid overflow.
In matrix notation (1) is equivalent to solve Zx = scale b,
where Z is defined as
Z = [ kron(In, A) -kron(B', Im) ] (2)
[ kron(In, D) -kron(E', Im) ].
Here Ik is the identity matrix of size k and X' is the tran-
spose of X. kron(X, Y) is the Kronecker product between the
matrices X and Y.
If TRANS = 'T', DTGSYL solves the transposed system Z'*y =
scale*b, which is equivalent to solve for R and L in
A' * R + D' * L = scale * C (3)
R * B' + L * E' = scale * (-F)
This case (TRANS = 'T') is used to compute an one-norm-based
estimate of Dif[(A,D), (B,E)], the separation between the
matrix pairs (A,D) and (B,E), using SLACON.
If IJOB >= 1, DTGSYL computes a Frobenius norm-based esti-
mate of Dif[(A,D),(B,E)]. That is, the reciprocal of a lower
bound on the reciprocal of the smallest singular value of Z.
See [1-2] for more information.
This is a level 3 BLAS algorithm.
TRANS (input)
= 'N', solve the generalized Sylvester equation
(1). = 'T', solve the 'transposed' system (3).
IJOB (input)
Specifies what kind of functionality to be per-
formed. =0: solve (1) only.
=1: The functionality of 0 and 3.
=2: The functionality of 0 and 4.
=3: Only an estimate of Dif[(A,D), (B,E)] is com-
puted. (look ahead strategy IJOB = 1 is used).
=4: Only an estimate of Dif[(A,D), (B,E)] is com-
puted. ( SGECON on sub-systems is used ). Not
referenced if TRANS = 'T'.
M (input) The order of the matrices A and D, and the row
dimension of the matrices C, F, R and L.
N (input) The order of the matrices B and E, and the column
dimension of the matrices C, F, R and L.
A (input) The upper quasi triangular matrix A.
LDA (input)
The leading dimension of the array A. LDA >=
max(1, M).
B (input) The upper quasi triangular matrix B.
LDB (input)
The leading dimension of the array B. LDB >=
max(1, N).
C (input/output)
On entry, C contains the right-hand-side of the
first matrix equation in (1) or (3). On exit, if
IJOB = 0, 1 or 2, C has been overwritten by the
solution R. If IJOB = 3 or 4 and TRANS = 'N', C
holds R, the solution achieved during the computa-
tion of the Dif-estimate.
LDC (input)
The leading dimension of the array C. LDC >=
max(1, M).
D (input) The upper triangular matrix D.
LDD (input)
The leading dimension of the array D. LDD >=
max(1, M).
E (input) The upper triangular matrix E.
LDE (input)
The leading dimension of the array E. LDE >=
max(1, N).
F (input/output)
On entry, F contains the right-hand-side of the
second matrix equation in (1) or (3). On exit, if
IJOB = 0, 1 or 2, F has been overwritten by the
solution L. If IJOB = 3 or 4 and TRANS = 'N', F
holds L, the solution achieved during the computa-
tion of the Dif-estimate.
LDF (input)
The leading dimension of the array F. LDF >=
max(1, M).
SCALE (output)
On exit SCALE is the reciprocal of a lower bound
of the reciprocal of the Dif-function, i.e. SCALE
is an upper bound of Dif[(A,D), (B,E)] =
sigma_min(Z), where Z as in (2). If IJOB = 0 or
TRANS = 'T', SCALE is not touched.
DIF (output)
On exit SCALE is the reciprocal of a lower bound
of the reciprocal of the Dif-function, i.e. SCALE
is an upper bound of Dif[(A,D), (B,E)] =
sigma_min(Z), where Z as in (2). If IJOB = 0 or
TRANS = 'T', SCALE is not touched.
WORK (workspace)
If IJOB = 0, WORK is not referenced. Otherwise,
on exit, if INFO = 0, WORK(1) returns the optimal
LWORK.
LWORK (input)
The dimension of the array WORK. LWORK > = 1. If
IJOB = 1 or 2 and TRANS = 'N', LWORK >= 2*M*N.
If LWORK = -1, then a workspace query is assumed;
the routine only calculates the optimal size of
the WORK array, returns this value as the first
entry of the WORK array, and no error message
related to LWORK is issued by XERBLA.
IWORK (workspace)
dimension(M+N+2)
INFO (output)
=0: successful exit
<0: If INFO = -i, the i-th argument had an illegal
value.
>0: (A, D) and (B, E) have common or close eigen-
values.
Based on contributions by
Bo Kagstrom and Peter Poromaa, Department of Computing
Science,
Umea University, S-901 87 Umea, Sweden.
[1] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and
Software
for Solving the Generalized Sylvester Equation and
Estimating the
Separation between Regular Matrix Pairs, Report UMINF -
93.23,
Department of Computing Science, Umea University, S-901
87 Umea,
Sweden, December 1993, Revised April 1994, Also as
LAPACK Working
Note 75. To appear in ACM Trans. on Math. Software, Vol
22,
No 1, 1996.
[2] B. Kagstrom, A Perturbation Analysis of the Generalized
Sylvester
Equation (AR - LB, DR - LE ) = (C, F), SIAM J. Matrix
Anal.
Appl., 15(4):1045-1060, 1994
[3] B. Kagstrom and L. Westin, Generalized Schur Methods
with
Condition Estimators for Solving the Generalized Sylves-
ter
Equation, IEEE Transactions on Automatic Control, Vol.
34, No. 7,
July 1989, pp 745-751.