Purpose
To compute orthogonal transformation matrices Q and Z which
reduce the N-th order descriptor system (A-lambda*E,B,C)
to the form
( Ac * ) ( Ec * ) ( Bc )
Q'*A*Z = ( ) , Q'*E*Z = ( ) , Q'*B = ( ) ,
( 0 Anc ) ( 0 Enc ) ( 0 )
C*Z = ( Cc Cnc ) ,
where the NCONT-th order descriptor system (Ac-lambda*Ec,Bc,Cc)
is a finite and/or infinite controllable. The pencil
Anc - lambda*Enc is regular of order N-NCONT and contains the
uncontrollable finite and/or infinite eigenvalues of the pencil
A-lambda*E.
For JOBCON = 'C' or 'I', the pencil ( Bc Ec-lambda*Ac ) has full
row rank NCONT for all finite lambda and is in a staircase form
with
_ _ _ _
( E1,0 E1,1 ... E1,k-1 E1,k )
( _ _ _ )
( Bc Ec ) = ( 0 E2,1 ... E2,k-1 E2,k ) , (1)
( ... _ _ )
( 0 0 ... Ek,k-1 Ek,k )
_ _ _
( A1,1 ... A1,k-1 A1,k )
( _ _ )
Ac = ( 0 ... A2,k-1 A2,k ) , (2)
( ... _ )
( 0 ... 0 Ak,k )
_
where Ei,i-1 is an rtau(i)-by-rtau(i-1) full row rank matrix
_
(with rtau(0) = M) and Ai,i is an rtau(i)-by-rtau(i)
upper triangular matrix.
For JOBCON = 'F', the pencil ( Bc Ac-lambda*Ec ) has full
row rank NCONT for all finite lambda and is in a staircase form
with
_ _ _ _
( A1,0 A1,1 ... A1,k-1 A1,k )
( _ _ _ )
( Bc Ac ) = ( 0 A2,1 ... A2,k-1 A2,k ) , (3)
( ... _ _ )
( 0 0 ... Ak,k-1 Ak,k )
_ _ _
( E1,1 ... E1,k-1 E1,k )
( _ _ )
Ec = ( 0 ... E2,k-1 E2,k ) , (4)
( ... _ )
( 0 ... 0 Ek,k )
_
where Ai,i-1 is an rtau(i)-by-rtau(i-1) full row rank matrix
_
(with rtau(0) = M) and Ei,i is an rtau(i)-by-rtau(i)
upper triangular matrix.
For JOBCON = 'C', the (N-NCONT)-by-(N-NCONT) regular pencil
Anc - lambda*Enc has the form
( Ainc - lambda*Einc * )
Anc - lambda*Enc = ( ) ,
( 0 Afnc - lambda*Efnc )
where:
1) the NIUCON-by-NIUCON regular pencil Ainc - lambda*Einc,
with Ainc upper triangular and nonsingular, contains the
uncontrollable infinite eigenvalues of A - lambda*E;
2) the (N-NCONT-NIUCON)-by-(N-NCONT-NIUCON) regular pencil
Afnc - lambda*Efnc, with Efnc upper triangular and
nonsingular, contains the uncontrollable finite
eigenvalues of A - lambda*E.
Note: The significance of the two diagonal blocks can be
interchanged by calling the routine with the
arguments A and E interchanged. In this case,
Ainc - lambda*Einc contains the uncontrollable zero
eigenvalues of A - lambda*E, while Afnc - lambda*Efnc
contains the uncontrollable nonzero finite and infinite
eigenvalues of A - lambda*E.
For JOBCON = 'F', the pencil Anc - lambda*Enc has the form
Anc - lambda*Enc = Afnc - lambda*Efnc ,
where the regular pencil Afnc - lambda*Efnc, with Efnc
upper triangular and nonsingular, contains the uncontrollable
finite eigenvalues of A - lambda*E.
For JOBCON = 'I', the pencil Anc - lambda*Enc has the form
Anc - lambda*Enc = Ainc - lambda*Einc ,
where the regular pencil Ainc - lambda*Einc, with Ainc
upper triangular and nonsingular, contains the uncontrollable
nonzero finite and infinite eigenvalues of A - lambda*E.
The left and/or right orthogonal transformations Q and Z
performed to reduce the system matrices can be optionally
accumulated.
The reduced order descriptor system (Ac-lambda*Ec,Bc,Cc) has
the same transfer-function matrix as the original system
(A-lambda*E,B,C).
Specification
SUBROUTINE TG01HD( JOBCON, COMPQ, COMPZ, N, M, P, A, LDA, E, LDE,
$ B, LDB, C, LDC, Q, LDQ, Z, LDZ, NCONT, NIUCON,
$ NRBLCK, RTAU, TOL, IWORK, DWORK, INFO )
C .. Scalar Arguments ..
CHARACTER COMPQ, COMPZ, JOBCON
INTEGER INFO, LDA, LDB, LDC, LDE, LDQ, LDZ,
$ M, N, NCONT, NIUCON, NRBLCK, P
DOUBLE PRECISION TOL
C .. Array Arguments ..
INTEGER IWORK( * ), RTAU( * )
DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
$ DWORK( * ), E( LDE, * ), Q( LDQ, * ),
$ Z( LDZ, * )
Arguments
Mode Parameters
JOBCON CHARACTER*1
= 'C': separate both finite and infinite uncontrollable
eigenvalues;
= 'F': separate only finite uncontrollable eigenvalues:
= 'I': separate only nonzero finite and infinite
uncontrollable eigenvalues.
COMPQ CHARACTER*1
= 'N': do not compute Q;
= 'I': Q is initialized to the unit matrix, and the
orthogonal matrix Q is returned;
= 'U': Q must contain an orthogonal matrix Q1 on entry,
and the product Q1*Q is returned.
COMPZ CHARACTER*1
= 'N': do not compute Z;
= 'I': Z is initialized to the unit matrix, and the
orthogonal matrix Z is returned;
= 'U': Z must contain an orthogonal matrix Z1 on entry,
and the product Z1*Z is returned.
Input/Output Parameters
N (input) INTEGER
The dimension of the descriptor state vector; also the
order of square matrices A and E, the number of rows of
matrix B, and the number of columns of matrix C. N >= 0.
M (input) INTEGER
The dimension of descriptor system input vector; also the
number of columns of matrix B. M >= 0.
P (input) INTEGER
The dimension of descriptor system output vector; also the
number of rows of 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,
( Ac * )
Q'*A*Z = ( ) ,
( 0 Anc )
where Ac is NCONT-by-NCONT and Anc is
(N-NCONT)-by-(N-NCONT).
If JOBCON = 'F', the matrix ( Bc Ac ) is in the
controllability staircase form (3).
If JOBCON = 'C' or 'I', the submatrix Ac is upper
triangular.
If JOBCON = 'C', the Anc matrix has the form
( Ainc * )
Anc = ( ) ,
( 0 Afnc )
where the NIUCON-by-NIUCON matrix Ainc is nonsingular and
upper triangular.
If JOBCON = 'I', Anc is nonsingular and upper triangular.
LDA INTEGER
The leading dimension of 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,
( Ec * )
Q'*E*Z = ( ) ,
( 0 Enc )
where Ec is NCONT-by-NCONT and Enc is
(N-NCONT)-by-(N-NCONT).
If JOBCON = 'C' or 'I', the matrix ( Bc Ec ) is in the
controllability staircase form (1).
If JOBCON = 'F', the submatrix Ec is upper triangular.
If JOBCON = 'C', the Enc matrix has the form
( Einc * )
Enc = ( ) ,
( 0 Efnc )
where the NIUCON-by-NIUCON matrix Einc is nilpotent
and the (N-NCONT-NIUCON)-by-(N-NCONT-NIUCON) matrix Efnc
is nonsingular and upper triangular.
If JOBCON = 'F', Enc is nonsingular and upper triangular.
LDE INTEGER
The leading dimension of 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
( Bc )
Q'*B = ( ) ,
( 0 )
where Bc is NCONT-by-M.
For JOBCON = 'C' or 'I', the matrix ( Bc Ec ) is in the
controllability staircase form (1).
For JOBCON = 'F', the matrix ( Bc Ac ) is in the
controllability staircase form (3).
LDB INTEGER
The leading dimension of 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.
LDC INTEGER
The leading dimension of array C. LDC >= MAX(1,P).
Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
If COMPQ = 'N': Q is not referenced.
If COMPQ = 'I': on entry, Q need not be set;
on exit, the leading N-by-N part of this
array contains the orthogonal matrix Q,
where Q' is the product of transformations
which are applied to A, E, and B on
the left.
If COMPQ = 'U': on entry, the leading N-by-N part of this
array must contain an orthogonal matrix
Qc;
on exit, the leading N-by-N part of this
array contains the orthogonal matrix
Qc*Q.
LDQ INTEGER
The leading dimension of array Q.
LDQ >= 1, if COMPQ = 'N';
LDQ >= MAX(1,N), if COMPQ = 'U' or 'I'.
Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
If COMPZ = 'N': Z is not referenced.
If COMPZ = 'I': on entry, Z need not be set;
on exit, the leading N-by-N part of this
array contains the orthogonal matrix Z,
which is the product of transformations
applied to A, E, and C on the right.
If COMPZ = 'U': on entry, the leading N-by-N part of this
array must contain an orthogonal matrix
Zc;
on exit, the leading N-by-N part of this
array contains the orthogonal matrix
Zc*Z.
LDZ INTEGER
The leading dimension of array Z.
LDZ >= 1, if COMPZ = 'N';
LDZ >= MAX(1,N), if COMPZ = 'U' or 'I'.
NCONT (output) INTEGER
The order of the reduced matrices Ac and Ec, and the
number of rows of reduced matrix Bc; also the order of
the controllable part of the pair (A-lambda*E,B).
NIUCON (output) INTEGER
For JOBCON = 'C', the order of the reduced matrices
Ainc and Einc; also the number of uncontrollable
infinite eigenvalues of the pencil A - lambda*E.
For JOBCON = 'F' or 'I', NIUCON has no significance
and is set to zero.
NRBLCK (output) INTEGER
For JOBCON = 'C' or 'I', the number k, of full row rank
_
blocks Ei,i in the staircase form of the pencil
(Bc Ec-lambda*Ac) (see (1) and (2)).
For JOBCON = 'F', the number k, of full row rank blocks
_
Ai,i in the staircase form of the pencil (Bc Ac-lambda*Ec)
(see (3) and (4)).
RTAU (output) INTEGER array, dimension (N)
RTAU(i), for i = 1, ..., NRBLCK, is the row dimension of
_ _
the full row rank block Ei,i-1 or Ai,i-1 in the staircase
form (1) or (3) for JOBCON = 'C' or 'I', or
for JOBCON = 'F', respectively.
Tolerances
TOL DOUBLE PRECISION
The tolerance to be used in rank determinations when
transforming (A-lambda*E, B). If the user sets TOL > 0,
then the given value of TOL is used as a lower bound for
reciprocal condition numbers in rank determinations; a
(sub)matrix whose estimated condition number is less than
1/TOL is considered to be of full rank. If the user sets
TOL <= 0, then an implicitly computed, default tolerance,
defined by TOLDEF = N*N*EPS, is used instead, where EPS
is the machine precision (see LAPACK Library routine
DLAMCH). TOL < 1.
Workspace
IWORK INTEGER array, dimension (M) DWORK DOUBLE PRECISION array, dimension (MAX(N,2*M))Error Indicator
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value.
Method
The subroutine is based on the reduction algorithms of [1].References
[1] A. Varga
Computation of Irreducible Generalized State-Space
Realizations.
Kybernetika, vol. 26, pp. 89-106, 1990.
Numerical Aspects
The algorithm is numerically backward stable and requires 0( N**3 ) floating point operations.Further Comments
If the system matrices A, E and B are badly scaled, it is generally recommendable to scale them with the SLICOT routine TG01AD, before calling TG01HD.Example
Program Text
* TG01HD EXAMPLE PROGRAM TEXT
* Copyright (c) 2002-2017 NICONET e.V.
*
* .. Parameters ..
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER LMAX, NMAX, MMAX, PMAX
PARAMETER ( LMAX = 20, NMAX = 20, MMAX = 20, PMAX = 20 )
INTEGER LDA, LDB, LDC, LDE, LDQ, LDZ
PARAMETER ( LDA = LMAX, LDB = LMAX, LDC = PMAX,
$ LDE = LMAX, LDQ = LMAX, LDZ = NMAX )
INTEGER LDWORK
PARAMETER ( LDWORK = MAX( 1, NMAX, 2*MMAX ) )
* .. Local Scalars ..
CHARACTER*1 COMPQ, COMPZ, JOBCO
INTEGER I, INFO, J, M, N, NCONT, NIUCON, NRBLCK, P
DOUBLE PRECISION TOL
* .. Local Arrays ..
INTEGER IWORK(MMAX), RTAU(NMAX)
DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,NMAX),
$ DWORK(LDWORK), E(LDE,NMAX), Q(LDQ,LMAX),
$ Z(LDZ,NMAX)
* .. External Subroutines ..
EXTERNAL TG01HD
* .. Intrinsic Functions ..
INTRINSIC MAX
* .. 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, TOL, JOBCO
COMPQ = 'I'
COMPZ = 'I'
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 transformed descriptor system (A-lambda E,B,C).
CALL TG01HD( JOBCO, COMPQ, COMPZ, N, M, P, A, LDA,
$ E, LDE, B, LDB, C, LDC, Q, LDQ, Z, LDZ,
$ NCONT, NIUCON, NRBLCK, RTAU, TOL, IWORK,
$ DWORK, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99994 ) NCONT, NIUCON
WRITE ( NOUT, FMT = 99985 )
WRITE ( NOUT, FMT = 99984 ) ( RTAU(I), I = 1,NRBLCK )
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
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' TG01HD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from TG01HD = ',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 (' Dimension of controllable part =', I5/
$ ' Number of uncontrollable 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 (/' L is out of range.',/' L = ',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 staircase form row dimensions are ' )
99984 FORMAT (10I5)
END
Program Data
TG01HD EXAMPLE PROGRAM DATA
7 3 2 0.0 C
2 0 2 0 -1 3 1
0 1 0 0 1 0 0
0 0 0 1 0 0 1
0 0 2 0 -1 3 1
0 0 0 1 0 0 1
0 1 0 0 1 0 0
0 0 0 1 0 0 1
0 0 1 0 0 0 0
0 0 0 0 0 1 0
0 0 0 0 0 0 1
0 0 0 0 0 0 1
0 0 0 1 0 0 0
0 0 1 0 -1 0 0
1 3 0 2 0 0 0
2 1 0
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
1 2 3
1 0 0 1 0 0 1
0 -1 1 0 -1 1 0
Program Results
TG01HD EXAMPLE PROGRAM RESULTS
Dimension of controllable part = 3
Number of uncontrollable infinite eigenvalues = 1
The staircase form row dimensions are
2 1
The transformed state dynamics matrix Q'*A*Z is
0.0000 0.0000 0.0000 0.0000 -1.2627 0.4334 0.4666
0.0000 2.0000 0.0000 -3.7417 -0.8520 0.2924 -0.4342
0.0000 0.0000 1.7862 0.3780 -0.2651 -0.7723 0.0000
0.0000 0.0000 0.0000 3.7417 0.8520 -0.2924 0.4342
0.0000 0.0000 0.0000 0.0000 -1.5540 0.5334 0.5742
0.0000 0.0000 0.0000 0.0000 -0.6533 0.2242 0.2414
0.0000 0.0000 0.0000 0.0000 -0.5892 0.2022 0.2177
The transformed descriptor matrix Q'*E*Z is
-1.8325 1.0000 2.3752 0.0000 -0.8214 0.2819 1.8016
0.4887 0.0000 0.3770 -0.5345 0.1874 0.5461 0.0000
-0.1728 0.0000 -0.1333 -1.1339 0.1325 0.3861 0.0000
0.0000 0.0000 0.0000 0.0000 0.8520 -0.2924 0.4342
0.0000 0.0000 0.0000 0.0000 -1.0260 -0.1496 0.0000
0.0000 0.0000 0.0000 0.0000 0.0000 1.1937 0.0000
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000
The transformed input/state matrix Q'*B is
1.0000 2.0000 3.0000
2.0000 1.0000 0.0000
0.0000 0.0000 0.0000
0.0000 0.0000 0.0000
0.0000 0.0000 0.0000
0.0000 0.0000 0.0000
0.0000 0.0000 0.0000
The transformed state/output matrix C*Z is
0.0000 1.0000 0.0000 0.0000 -1.2627 0.4334 0.4666
0.3665 0.0000 -0.9803 -1.6036 0.1874 0.5461 0.0000
The left transformation matrix Q is
0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.7071 0.0000 0.2740 -0.6519 0.0000
0.0000 0.0000 0.0000 0.0000 0.8304 0.3491 -0.4342
0.0000 0.0000 0.0000 -1.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000 0.4003 0.1683 0.9008
0.0000 0.0000 0.7071 0.0000 -0.2740 0.6519 0.0000
1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
The right transformation matrix Z is
0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000
-0.6108 0.0000 0.7917 0.0000 0.0000 0.0000 0.0000
0.4887 0.0000 0.3770 -0.5345 0.1874 0.5461 0.0000
0.0000 0.0000 0.0000 0.0000 -0.4107 0.1410 0.9008
0.6108 0.0000 0.4713 0.2673 -0.1874 -0.5461 0.0000
-0.1222 0.0000 -0.0943 -0.8018 -0.1874 -0.5461 0.0000
0.0000 0.0000 0.0000 0.0000 -0.8520 0.2924 -0.4342