Purpose
To compute equivalence 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 form (if JOB = 'F')
( Af 0 ) ( Ef 0 )
Q*A*Z = ( ) , Q*E*Z = ( ) , (1)
( 0 Ai ) ( 0 Ei )
or to the form (if JOB = 'I')
( Ai 0 ) ( Ei 0 )
Q*A*Z = ( ) , Q*E*Z = ( ) , (2)
( 0 Af ) ( 0 Ef )
where the pair (Af,Ef) is in a generalized real Schur form, with
Ef nonsingular and upper triangular and Af in real Schur form.
The subpencil Af-lambda*Ef contains the finite eigenvalues.
The pair (Ai,Ei) is in a generalized real Schur form with
both Ai and Ei upper triangular. The subpencil Ai-lambda*Ei,
with Ai nonsingular and Ei nilpotent contains the infinite
eigenvalues and is in a block staircase form (see METHOD).
This decomposition corresponds to an additive decomposition of
the transfer-function matrix of the descriptor system as the
sum of a proper term and a polynomial term.
Specification
SUBROUTINE TG01ND( JOB, JOBT, N, M, P, A, LDA, E, LDE, B, LDB,
$ C, LDC, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ,
$ NF, ND, NIBLCK, IBLCK, TOL, IWORK, DWORK,
$ LDWORK, INFO )
C .. Scalar Arguments ..
CHARACTER JOB, JOBT
INTEGER INFO, LDA, LDB, LDC, LDE, LDQ, LDWORK, LDZ, M,
$ N, ND, NF, NIBLCK, P
DOUBLE PRECISION TOL
C .. Array Arguments ..
INTEGER IBLCK( * ), IWORK(*)
DOUBLE PRECISION A(LDA,*), ALPHAR(*), ALPHAI(*), B(LDB,*),
$ BETA(*), C(LDC,*), DWORK(*), E(LDE,*),
$ Q(LDQ,*), Z(LDZ,*)
Arguments
Mode Parameters
JOB CHARACTER*1
= 'F': perform the finite-infinite separation;
= 'I': perform the infinite-finite separation.
JOBT CHARACTER*1
= 'D': compute the direct transformation matrices;
= 'I': compute the inverse transformation matrices
inv(Q) and inv(Z).
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.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading N-by-N part of this array must
contain the N-by-N state matrix A.
On exit, the leading N-by-N part of this array contains
the transformed state matrix Q*A*Z (if JOBT = 'D') or
inv(Q)*A*inv(Z) (if JOBT = 'I') in the form
( Af 0 ) ( Ai 0 )
( ) for JOB = 'F', or ( ) for JOB = 'I',
( 0 Ai ) ( 0 Af )
where Af is an NF-by-NF matrix in real Schur form, and Ai
is an (N-NF)-by-(N-NF) nonsingular and upper triangular
matrix. Ai has a block structure as in (3) or (4), where
A0,0 is ND-by-ND and Ai,i , for i = 1, ..., NIBLCK, is
IBLCK(i)-by-IBLCK(i). (See METHOD.)
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 N-by-N descriptor matrix E.
On exit, the leading N-by-N part of this array contains
the transformed descriptor matrix Q*E*Z (if JOBT = 'D') or
inv(Q)*E*inv(Z) (if JOBT = 'I') in the form
( Ef 0 ) ( Ei 0 )
( ) for JOB = 'F', or ( ) for JOB = 'I',
( 0 Ei ) ( 0 Ef )
where Ef is an NF-by-NF nonsingular and upper triangular
matrix, and Ei is an (N-NF)-by-(N-NF) nilpotent matrix in
an upper triangular block form as in (3) or (4).
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 N-by-M input matrix B.
On exit, the leading N-by-M part of this array contains
the transformed input matrix Q*B (if JOBT = 'D') or
inv(Q)*B (if JOBT = 'I').
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 state/output matrix C.
On exit, the leading P-by-N part of this array contains
the transformed matrix C*Z (if JOBT = 'D') or C*inv(Z)
(if JOBT = 'I').
LDC INTEGER
The leading dimension of the array C. LDC >= MAX(1,P).
ALPHAR (output) DOUBLE PRECISION array, dimension (N)
ALPHAR(1:NF) will be set to the real parts of the diagonal
elements of Af that would result from reducing A and E to
the Schur form, and then further reducing both of them to
triangular form using unitary transformations, subject to
having the diagonal of E positive real. Thus, if Af(j,j)
is in a 1-by-1 block (i.e., Af(j+1,j) = Af(j,j+1) = 0),
then ALPHAR(j) = Af(j,j). Note that the (real or complex)
values (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,NF, are
the finite generalized eigenvalues of the matrix pencil
A - lambda*E.
ALPHAI (output) DOUBLE PRECISION array, dimension (N)
ALPHAI(1:NF) will be set to the imaginary parts of the
diagonal elements of Af that would result from reducing A
and E to Schur form, and then further reducing both of
them to triangular form using unitary transformations,
subject to having the diagonal of E positive real. Thus,
if Af(j,j) is in a 1-by-1 block (see above), then
ALPHAI(j) = 0. Note that the (real or complex) values
(ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,NF, are the
finite generalized eigenvalues of the matrix pencil
A - lambda*E.
BETA (output) DOUBLE PRECISION array, dimension (N)
BETA(1:NF) will be set to the (real) diagonal elements of
Ef that would result from reducing A and E to Schur form,
and then further reducing both of them to triangular form
using unitary transformations, subject to having the
diagonal of E positive real. Thus, if Af(j,j) is in a
1-by-1 block (see above), then BETA(j) = Ef(j,j).
Note that the (real or complex) values
(ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,NF, are the
finite generalized eigenvalues of the matrix pencil
A - lambda*E.
Q (output) DOUBLE PRECISION array, dimension (LDQ,N)
The leading N-by-N part of this array contains the
left transformation matrix Q, if JOBT = 'D', or its
inverse inv(Q), if JOBT = 'I'.
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
right transformation matrix Z, if JOBT = 'D', or its
inverse inv(Z), if JOBT = 'I'.
LDZ INTEGER
The leading dimension of the array Z. LDZ >= MAX(1,N).
NF (output) INTEGER
The order of the reduced matrices Af and Ef; also, the
number of finite generalized eigenvalues of the pencil
A-lambda*E.
ND (output) INTEGER
The number of non-dynamic infinite eigenvalues of the
matrix 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 (3), where i = 1,2,...,NIBLCK.
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+6)
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 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 iteration did not converge;
= 3: (Af,Ef) and (Ai,Ei) have too close generalized
eigenvalues.
Method
For the separation of infinite structure, the reduction algorithm
of [1] is employed. This separation is achieved by computing
orthogonal matrices Q1 and Z1 such that Q1*A*Z1 and Q1*E*Z1
have the form (if JOB = 'F')
( Af Ao ) ( Ef Eo )
Q1*A*Z1 = ( ) , Q1*E*Z1 = ( ) ,
( 0 Ai ) ( 0 Ei )
or to the form (if JOB = 'I')
( Ai Ao ) ( Ei Eo )
Q1*A*Z1 = ( ) , Q1*E*Z1 = ( ) .
( 0 Af ) ( 0 Ef )
If JOB = 'F', the matrices Ai and Ei 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 ) ; (3)
( : : . : ) ( : : . : )
( 0 0 ... A1,1 ) ( 0 0 ... 0 )
if JOB = 'I' the matrices Ai and Ei have the form
( A1,1 ... A1,k A1,0 ) ( 0 ... E1,k E1,0 )
Ai = ( : . : : ) , Ei = ( : . : : ) , (4)
( : ... 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. A0,0 corresponds to the non-dynamic infinite modes of
the system.
In a second step, the transformation matrices Q2 and Z2 are
determined, of the form
( I -X ) ( I Y )
Q2 = ( ) , Z2 = ( )
( 0 I ) ( 0 I )
such that with Q = Q2*Q1 and Z = Z1*Z2, Q*A*Z and Q*E*Z are
block diagonal as in (1) (if JOB = 'F') or in (2) (if JOB = 'I').
X and Y are computed by solving generalized Sylvester equations.
If we partition Q*B and C*Z according to (1) or (2) in the form
( Bf ) and ( Cf Ci ), if JOB = 'F', or ( Bi ) and ( Ci Cf ), if
( Bi ) ( Bf )
JOB = 'I', then (Af-lambda*Ef,Bf,Cf) is the stricly proper part
of the original descriptor system and (Ai-lambda*Ei,Bi,Ci) is its
polynomial part.
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
The algorithm is numerically backward stable and requires 0( N**3 ) floating point operations.Further Comments
The number of infinite poles is computed as
NIBLCK
NINFP = Sum IBLCK(i) = N - ND - NF.
i=1
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
* TG01ND EXAMPLE PROGRAM TEXT
* Copyright (c) 2002-2017 NICONET e.V.
*
* .. Parameters ..
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 )
* .. Local Scalars ..
CHARACTER*1 JOB, JOBT
INTEGER I, INFO, J, M, N, ND, NF, NIBLCK, P
DOUBLE PRECISION TOL
* .. Local Arrays ..
INTEGER IBLCK(NMAX), IWORK(NMAX+6)
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 Subroutines ..
EXTERNAL TG01ND
* .. 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, JOB, JOBT, 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 TG01ND( JOB, JOBT, N, M, P, A, LDA, E, LDE, B, LDB,
$ C, LDC, ALPHAR, ALPHAI, BETA, Q, LDQ, Z,
$ LDZ, NF, ND, NIBLCK, IBLCK, TOL, IWORK,
$ DWORK, LDWORK, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99994 ) NF, 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, NF
WRITE ( NOUT, FMT = 99984 )
$ DCMPLX( ALPHAR(I), ALPHAI(I) )/BETA(I)
70 CONTINUE
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' TG01ND EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from TG01ND = ',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 (' Order of reduced system =', I5/
$ ' Number of non-dynamic infinite eigenvalues =', I5)
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 ')
END
Program Data
TG01ND EXAMPLE PROGRAM DATA
4 2 2 F D 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
TG01ND EXAMPLE PROGRAM RESULTS
Order of reduced system = 3
Number of non-dynamic infinite eigenvalues = 1
Number of infinite blocks = 1
The transformed state dynamics matrix Q*A*Z is
1.2803 -2.3613 -0.9025 0.0000
0.0000 -0.5796 0.8504 0.0000
0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 2.2913
The transformed descriptor matrix Q*E*Z is
9.3142 -4.1463 5.4026 0.0000
0.0000 0.1594 0.1212 0.0000
0.0000 0.0000 2.3524 0.0000
0.0000 0.0000 0.0000 0.0000
The transformed input/state matrix Q*B is
7.7328 1.6760
2.2870 0.4660
-1.2140 -1.2140
1.1339 0.3780
The transformed state/output matrix C*Z is
-0.0469 -0.9391 -0.8847 -6.0622
-1.0697 0.3620 1.1795 -0.0000
The left transformation matrix Q is
3.7620 3.8560 -2.2948 3.9708
1.4909 0.0798 -0.3301 0.7961
-0.0000 -0.0000 0.0000 -1.2140
0.5669 0.5669 -0.1890 0.5669
The right transformation matrix Z is
0.0469 0.9391 -0.0843 6.0622
-0.9962 0.0189 -0.0211 -3.0311
0.0000 -0.0000 -0.9689 -0.0000
-0.0735 0.3432 0.2317 3.0311
The finite generalized eigenvalues are
real part imag part
0.1375 +0.0000i
-3.6375 +0.0000i
0.0000 +0.0000i