Purpose
To compute orthogonal transformation matrices Q and Z which
reduce the regular pole pencil A-lambda*E of the descriptor system
(A-lambda*E,B,C) to the generalized real Schur form with ordered
generalized eigenvalues. The pair (A,E) is reduced to the form
( A1 * * ) ( E1 * * )
Q'*A*Z = ( 0 A2 * ) , Q'*E*Z = ( 0 E2 * ) , (1)
( 0 0 A3 ) ( 0 0 E2 )
where the subpencils Ak-lambda*Ek, for k = 1, 2, 3, contain the
generalized eigenvalues which belong to certain domains of
interest.
Specification
SUBROUTINE TG01QD( DICO, STDOM, JOBFI, N, M, P, ALPHA, A, LDA,
$ E, LDE, B, LDB, C, LDC, N1, N2, N3, ND, NIBLCK,
$ IBLCK, Q, LDQ, Z, LDZ, ALPHAR, ALPHAI, BETA,
$ TOL, IWORK, DWORK, LDWORK, INFO )
C .. Scalar Arguments ..
CHARACTER DICO, JOBFI, STDOM
INTEGER INFO, LDA, LDB, LDC, LDE, LDQ, LDWORK, LDZ, M, N,
$ N1, N2, N3, ND, NIBLCK, P
DOUBLE PRECISION ALPHA, TOL
C .. Array Arguments ..
INTEGER IBLCK( * ), IWORK(*)
DOUBLE PRECISION A(LDA,*), ALPHAI(*), ALPHAR(*), B(LDB,*),
$ BETA(*), C(LDC,*), DWORK(*), E(LDE,*), Q(LDQ,*),
$ Z(LDZ,*)
Arguments
Mode Parameters
DICO CHARACTER*1
Specifies the type of the descriptor system as follows:
= 'C': continuous-time system;
= 'D': discrete-time system.
STDOM CHARACTER*1
Specifies the type of the domain of interest for the
generalized eigenvalues, as follows:
= 'S': stability type domain (i.e., left part of complex
plane or inside of a circle);
= 'U': instability type domain (i.e., right part of complex
plane or outside of a circle);
= 'N': whole complex domain, excepting infinity.
JOBFI CHARACTER*1
Specifies the type of generalized eigenvalues in the
leading diagonal block(s) as follows:
= 'F': finite generalized eigenvalues are in the
leading diagonal blocks (Af,Ef), and the resulting
transformed pair has the form
( Af * ) ( Ef * )
Q'*A*Z = ( ) , Q'*E*Z = ( ) ;
( 0 Ai ) ( 0 Ei )
= 'I': infinite generalized eigenvalues are in the
leading diagonal blocks (Ai,Ei), and the resulting
transformed pair has the form
( Ai * ) ( Ei * )
Q'*A*Z = ( ) , Q'*E*Z = ( ) .
( 0 Af ) ( 0 Ef )
Input/Output Parameters
N (input) INTEGER
The number of rows of the matrix B, the number of columns
of the matrix C and the order of the square matrices A
and E. N >= 0.
M (input) INTEGER
The number of columns of the matrix B. M >= 0.
P (input) INTEGER
The number of rows of the matrix C. P >= 0.
ALPHA (input) DOUBLE PRECISION
The boundary of the domain of interest for the finite
generalized eigenvalues of the pair (A,E). For a
continuous-time system (DICO = 'C'), ALPHA is the boundary
value for the real parts of the generalized eigenvalues,
while for a discrete-time system (DICO = 'D'), ALPHA >= 0
represents the boundary value for the moduli of the
generalized eigenvalues.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading N-by-N part of this array must
contain the state dynamics matrix A.
On exit, the leading N-by-N part of this array contains
the matrix Q'*A*Z in real Schur form, with the elements
below the first subdiagonal set to zero.
If JOBFI = 'I', the N1-by-N1 pair (A1,E1) contains the
infinite spectrum, the N2-by-N2 pair (A2,E2) contains the
finite spectrum in the domain of interest, and the
N3-by-N3 pair (A3,E3) contains the finite spectrum ouside
of the domain of interest.
If JOBFI = 'F', the N1-by-N1 pair (A1,E1) contains the
finite spectrum in the domain of interest, the N2-by-N2
pair (A2,E2) contains the finite spectrum ouside of the
domain of interest, and the N3-by-N3 pair (A3,E3) contains
the infinite spectrum.
Ai has a block structure as in (2), where A0,0 is ND-by-ND
and Ai,i is IBLCK(i)-by-IBLCK(i), for i = 1, ..., NIBLCK.
The domain of interest for the pair (Af,Ef), containing
the finite generalized eigenvalues, is defined by the
parameters ALPHA, DICO and STDOM as follows:
For DICO = 'C':
Real(eig(Af,Ef)) < ALPHA if STDOM = 'S';
Real(eig(Af,Ef)) > ALPHA if STDOM = 'U'.
For DICO = 'D':
Abs(eig(Af,Ef)) < ALPHA if STDOM = 'S';
Abs(eig(Af,Ef)) > ALPHA if STDOM = 'U'.
LDA INTEGER
The leading dimension of the array A. LDA >= MAX(1,N).
E (input/output) DOUBLE PRECISION array, dimension (LDE,N)
On entry, the leading N-by-N part of this array must
contain the descriptor matrix E.
On exit, the leading N-by-N part of this array contains
the matrix Q'*E*Z in upper triangular form, with the
elements below the diagonal set to zero. Its structure
corresponds to the block structure of the matrix Q'*A*Z
(see description of A).
LDE INTEGER
The leading dimension of the array E. LDE >= MAX(1,N).
B (input/output) DOUBLE PRECISION array, dimension (LDB,M)
On entry, the leading N-by-M part of this array must
contain the input matrix B.
On exit, the leading N-by-M part of this array contains
the transformed input matrix Q'*B.
LDB INTEGER
The leading dimension of the array B. LDB >= MAX(1,N).
C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
On entry, the leading P-by-N part of this array must
contain the output matrix C.
On exit, the leading P-by-N part of this array contains
the transformed output matrix C*Z.
LDC INTEGER
The leading dimension of the array C. LDC >= MAX(1,P).
N1 (output) INTEGER
N2 (output) INTEGER
N3 (output) INTEGER
The number of the generalized eigenvalues of the pairs
(A1,E1), (A2,E2) and (A3,E3), respectively.
ND (output) INTEGER.
The number of non-dynamic infinite eigenvalues of the
pair (A,E). Note: N-ND is the rank of the matrix E.
NIBLCK (output) INTEGER
If ND > 0, the number of infinite blocks minus one.
If ND = 0, then NIBLCK = 0.
IBLCK (output) INTEGER array, dimension (N)
IBLCK(i) contains the dimension of the i-th block in the
staircase form (2), where i = 1,2,...,NIBLCK.
Q (output) DOUBLE PRECISION array, dimension (LDQ,N)
The leading N-by-N part of this array contains the
orthogonal matrix Q, where Q' is the product of orthogonal
transformations applied to A, E, and B on the left.
LDQ INTEGER
The leading dimension of the array Q. LDQ >= MAX(1,N).
Z (output) DOUBLE PRECISION array, dimension (LDZ,N)
The leading N-by-N part of this array contains the
orthogonal matrix Z, which is the product of orthogonal
transformations applied to A, E, and C on the right.
LDZ INTEGER
The leading dimension of the array Z. LDZ >= MAX(1,N).
ALPHAR (output) DOUBLE PRECISION array, dimension (N)
ALPHAI (output) DOUBLE PRECISION array, dimension (N)
BETA (output) DOUBLE PRECISION array, dimension (N)
On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j = 1, ..., N,
are the generalized eigenvalues.
ALPHAR(j) + ALPHAI(j)*i, and BETA(j), j = 1, ..., N, are
the diagonals of the complex Schur form (S,T) that would
result if the 2-by-2 diagonal blocks of the real Schur
form of (A,E) were further reduced to triangular form
using 2-by-2 complex unitary transformations.
If ALPHAI(j) is zero, then the j-th eigenvalue is real;
if positive, then the j-th and (j+1)-st eigenvalues are a
complex conjugate pair, with ALPHAI(j+1) negative.
Tolerances
TOL DOUBLE PRECISION
A tolerance used in rank decisions to determine the
effective rank, which is defined as the order of the
largest leading (or trailing) triangular submatrix in the
QR factorization with column pivoting whose estimated
condition number is less than 1/TOL. If the user sets
TOL <= 0, then an implicitly computed, default tolerance
TOLDEF = N**2*EPS, is used instead, where EPS is the
machine precision (see LAPACK Library routine DLAMCH).
TOL < 1.
Workspace
IWORK INTEGER array, dimension (N)
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, DWORK(1) returns the optimal value
of LDWORK.
LDWORK INTEGER
The length of the array DWORK. LDWORK >= 1, and if N = 0,
LDWORK >= 4*N, if STDOM = 'N';
LDWORK >= 4*N+16, if STDOM = 'S' or 'U'.
For optimum performance LDWORK should be larger.
If LDWORK = -1, then a workspace query is assumed; the
routine only calculates the optimal size of the DWORK
array, returns this value as the first entry of the DWORK
array, and no error message related to LDWORK is issued by
XERBLA.
Error Indicator
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value;
= 1: the pencil A-lambda*E is not regular;
= 2: the QZ algorithm failed to compute all generalized
eigenvalues of the pair (A,E);
= 3: a failure occured during the ordering of the
generalized real Schur form of the pair (A,E).
Method
The separation of the finite and infinite parts is based on the
reduction algorithm of [1].
If JOBFI = 'F', the matrices of the pair (Ai,Ei), containing the
infinite generalized eigenvalues, have the form
( A0,0 A0,k ... A0,1 ) ( 0 E0,k ... E0,1 )
Ai = ( 0 Ak,k ... Ak,1 ) , Ei = ( 0 0 ... Ek,1 ) ; (2)
( : : . : ) ( : : . : )
( 0 0 ... A1,1 ) ( 0 0 ... 0 )
if JOBFI = 'I', the matrices Ai and Ei have the form
( A1,1 ... A1,k A1,0 ) ( 0 ... E1,k E1,0 )
Ai = ( : . : : ) , Ei = ( : . : : ) , (3)
( : ... Ak,k Ak,0 ) ( : ... 0 Ek,0 )
( 0 ... 0 A0,0 ) ( 0 ... 0 0 )
where Ai,i , for i = 0, 1, ..., k, are nonsingular upper
triangular matrices, and A0,0 corresponds to the non-dynamic
infinite modes of the system.
References
[1] Misra, P., Van Dooren, P., and Varga, A.
Computation of structural invariants of generalized
state-space systems.
Automatica, 30, pp. 1921-1936, 1994.
Numerical Aspects
3 The algorithm requires about 25N floating point operations.Further Comments
The number of infinite poles is computed as
NIBLCK
NINFP = Sum IBLCK(i) = N - ND - NF,
i=1
where NF is the number of finite generalized eigenvalues.
The multiplicities of infinite poles can be computed as follows:
there are IBLCK(k)-IBLCK(k+1) infinite poles of multiplicity k,
for k = 1, ..., NIBLCK, where IBLCK(NIBLCK+1) = 0.
Note that each infinite pole of multiplicity k corresponds to an
infinite eigenvalue of multiplicity k+1.
Example
Program Text
* TG01QD EXAMPLE PROGRAM TEXT
* Copyright (c) 2002-2017 NICONET e.V.
*
* .. Parameters ..
DOUBLE PRECISION ZERO
PARAMETER ( ZERO = 0.0D0 )
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER NMAX, MMAX, PMAX
PARAMETER ( NMAX = 20, MMAX = 20, PMAX = 20 )
INTEGER LDA, LDB, LDC, LDE, LDQ, LDZ
PARAMETER ( LDA = NMAX, LDB = NMAX, LDC = PMAX,
$ LDE = NMAX, LDQ = NMAX, LDZ = NMAX )
INTEGER LDWORK
PARAMETER ( LDWORK = 4*NMAX+16 )
* .. Local Scalars ..
CHARACTER*1 DICO, JOBFI, STDOM
INTEGER I, INFO, J, M, N, N1, N2, N3, ND, NIBLCK, P
DOUBLE PRECISION ALPHA, TOL
* .. Local Arrays ..
INTEGER IBLCK(NMAX), IWORK(NMAX)
DOUBLE PRECISION A(LDA,NMAX), ALPHAI(NMAX), ALPHAR(NMAX),
$ B(LDB,MMAX), BETA(NMAX), C(LDC,NMAX),
$ DWORK(LDWORK), E(LDE,NMAX), Q(LDQ,NMAX),
$ Z(LDZ,NMAX)
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL TG01QD
* .. Intrinsic Functions ..
INTRINSIC DCMPLX
* .. Executable Statements ..
*
WRITE ( NOUT, FMT = 99999 )
* Skip the heading in the data file and read the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) N, M, P, DICO, STDOM, JOBFI, ALPHA, TOL
IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99988 ) N
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
READ ( NIN, FMT = * ) ( ( E(I,J), J = 1,N ), I = 1,N )
IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
WRITE ( NOUT, FMT = 99987 ) M
ELSE
READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1,N )
IF ( P.LT.0 .OR. P.GT.PMAX ) THEN
WRITE ( NOUT, FMT = 99986 ) P
ELSE
READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,P )
* Find the reduced descriptor system
* (A-lambda E,B,C).
CALL TG01QD( DICO, STDOM, JOBFI, N, M, P, ALPHA, A, LDA,
$ E, LDE, B, LDB, C, LDC, N1, N2, N3, ND,
$ NIBLCK, IBLCK, Q, LDQ, Z, LDZ, ALPHAR,
$ ALPHAI, BETA, TOL, IWORK, DWORK, LDWORK,
$ INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99994 ) N1, N2, N3
WRITE ( NOUT, FMT = 99983 ) ND
WRITE ( NOUT, FMT = 99989 ) NIBLCK + 1
IF ( NIBLCK.GT.0 ) THEN
WRITE ( NOUT, FMT = 99985 )
$ ( IBLCK(I), I = 1, NIBLCK )
END IF
WRITE ( NOUT, FMT = 99997 )
DO 10 I = 1, N
WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,N )
10 CONTINUE
WRITE ( NOUT, FMT = 99996 )
DO 20 I = 1, N
WRITE ( NOUT, FMT = 99995 ) ( E(I,J), J = 1,N )
20 CONTINUE
WRITE ( NOUT, FMT = 99993 )
DO 30 I = 1, N
WRITE ( NOUT, FMT = 99995 ) ( B(I,J), J = 1,M )
30 CONTINUE
WRITE ( NOUT, FMT = 99992 )
DO 40 I = 1, P
WRITE ( NOUT, FMT = 99995 ) ( C(I,J), J = 1,N )
40 CONTINUE
WRITE ( NOUT, FMT = 99991 )
DO 50 I = 1, N
WRITE ( NOUT, FMT = 99995 ) ( Q(I,J), J = 1,N )
50 CONTINUE
WRITE ( NOUT, FMT = 99990 )
DO 60 I = 1, N
WRITE ( NOUT, FMT = 99995 ) ( Z(I,J), J = 1,N )
60 CONTINUE
WRITE ( NOUT, FMT = 99985 )
DO 70 I = 1, N
IF ( BETA(I).EQ.ZERO .OR. ALPHAI(I).EQ.ZERO ) THEN
WRITE ( NOUT, FMT = 99984 )
$ ALPHAR(I)/BETA(I)
ELSE
WRITE ( NOUT, FMT = 99984 )
$ DCMPLX( ALPHAR(I), ALPHAI(I) )/BETA(I)
END IF
70 CONTINUE
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' TG01QD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from TG01QD = ',I2)
99997 FORMAT (/' The transformed state dynamics matrix Q''*A*Z is ')
99996 FORMAT (/' The transformed descriptor matrix Q''*E*Z is ')
99995 FORMAT (20(1X,F8.4))
99994 FORMAT (' Number of eigenvalues of the three blocks =', 3I5)
99993 FORMAT (/' The transformed input/state matrix Q''*B is ')
99992 FORMAT (/' The transformed state/output matrix C*Z is ')
99991 FORMAT (/' The left transformation matrix Q is ')
99990 FORMAT (/' The right transformation matrix Z is ')
99989 FORMAT ( ' Number of infinite blocks = ',I5)
99988 FORMAT (/' N is out of range.',/' N = ',I5)
99987 FORMAT (/' M is out of range.',/' M = ',I5)
99986 FORMAT (/' P is out of range.',/' P = ',I5)
99985 FORMAT (/' The finite generalized eigenvalues are '/
$ ' real part imag part ')
99984 FORMAT (1X,F9.4,SP,F9.4,S,'i ')
99983 FORMAT ( ' Number of non-dynamic infinite eigenvalues = ',I5)
END
Program Data
TG01QD EXAMPLE PROGRAM DATA
4 2 2 C S F -1.E-7 0.0
-1 0 0 3
0 0 1 2
1 1 0 4
0 0 0 0
1 2 0 0
0 1 0 1
3 9 6 3
0 0 2 0
1 0
0 0
0 1
1 1
-1 0 1 0
0 1 -1 1
Program Results
TG01QD EXAMPLE PROGRAM RESULTS
Number of eigenvalues of the three blocks = 1 2 1
Number of non-dynamic infinite eigenvalues = 1
Number of infinite blocks = 1
The transformed state dynamics matrix Q'*A*Z is
1.6311 2.1641 -0.5848 -3.6517
0.0000 -0.4550 1.0935 1.6851
0.0000 0.0000 0.0000 1.5770
0.0000 0.0000 0.0000 2.2913
The transformed descriptor matrix Q'*E*Z is
-0.4484 9.6340 5.1601 -2.6183
0.0000 -3.3099 -1.6050 -0.2756
0.0000 0.0000 2.3524 -0.6008
0.0000 0.0000 0.0000 0.0000
The transformed input/state matrix Q'*B is
0.0232 -0.9413
0.7251 0.2478
-0.4336 -0.9538
1.1339 0.3780
The transformed state/output matrix C*Z is
0.8621 0.3754 -0.8847 0.5774
0.1511 -1.1192 1.1795 0.5774
The left transformation matrix Q is
0.0232 0.7251 0.3902 0.5669
-0.3369 -0.6425 0.3902 0.5669
-0.9413 0.2478 -0.1301 -0.1890
0.0000 0.0000 -0.8238 0.5669
The right transformation matrix Z is
-0.8621 -0.3754 -0.0843 -0.3299
0.4258 -0.9008 -0.0211 -0.0825
0.0000 0.0000 -0.9689 0.2474
-0.2748 -0.2184 0.2317 0.9073
The finite generalized eigenvalues are
real part imag part
-3.6375
0.1375
0.0000
Infinity