Purpose
To find a reduced descriptor representation (Ar-lambda*Er,Br,Cr) without non-dynamic modes for a descriptor representation (A-lambda*E,B,C). Optionally, the reduced descriptor system can be put into a standard form with the leading diagonal block of Er identity.Specification
SUBROUTINE TG01GD( JOBS, L, N, M, P, A, LDA, E, LDE, B, LDB,
$ C, LDC, D, LDD, LR, NR, RANKE, INFRED, TOL,
$ IWORK, DWORK, LDWORK, INFO )
C .. Scalar Arguments ..
CHARACTER JOBS
INTEGER INFO, INFRED, L, LDA, LDB, LDC, LDD, LDE,
$ LDWORK, LR, M, N, NR, P, RANKE
DOUBLE PRECISION TOL
C .. Array Arguments ..
INTEGER IWORK(*)
DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*),
$ DWORK(*), E(LDE,*)
Arguments
Mode Parameters
JOBS CHARACTER*1
Indicates whether the user wishes to transform the leading
diagonal block of Er to an identity matrix, as follows:
= 'S': make Er with leading diagonal identity;
= 'D': keep Er unreduced or upper triangular.
Input/Output Parameters
L (input) INTEGER
The number of rows of the matrices A, E, and B;
also the number of differential equations. L >= 0.
N (input) INTEGER
The number of columns of the matrices A, E, and C;
also the dimension of descriptor state vector. N >= 0.
M (input) INTEGER
The number of columns of the matrix B;
also the dimension of the input vector. M >= 0.
P (input) INTEGER
The number of rows of the matrix C.
also the dimension of the output vector. P >= 0.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading L-by-N part of this array must
contain the state dynamics matrix A.
On exit, if NR < N, the leading LR-by-NR part of this
array contains the reduced order state matrix Ar of a
descriptor realization without non-dynamic modes.
Array A contains the original state dynamics matrix if
INFRED < 0.
LDA INTEGER
The leading dimension of the array A. LDA >= MAX(1,L).
E (input/output) DOUBLE PRECISION array, dimension (LDE,N)
On entry, the leading L-by-N part of this array must
contain the descriptor matrix E.
On exit, if INFRED >= 0, the leading LR-by-NR part of this
array contains the reduced order descriptor matrix Er of a
descriptor realization without non-dynamic modes.
In this case, only the leading RANKE-by-RANKE submatrix
of Er is nonzero and this submatrix is nonsingular and
upper triangular. Array E contains the original descriptor
matrix if INFRED < 0. If JOBS = 'S', then the leading
RANKE-by-RANKE submatrix results in an identity matrix.
LDE INTEGER
The leading dimension of the array E. LDE >= MAX(1,L).
B (input/output) DOUBLE PRECISION array, dimension (LDB,M)
On entry, the leading L-by-M part of this array must
contain the input matrix B.
On exit, the leading LR-by-M part of this array contains
the reduced order input matrix Br of a descriptor
realization without non-dynamic modes. Array B contains
the original input matrix if INFRED < 0.
LDB INTEGER
The leading dimension of the array B. LDB >= MAX(1,L).
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-NR part of this array contains
the reduced order output matrix Cr of a descriptor
realization without non-dynamic modes. Array C contains
the original output matrix if INFRED < 0.
LDC INTEGER
The leading dimension of the array C. LDC >= MAX(1,P).
D (input/output) DOUBLE PRECISION array, dimension (LDD,M)
On entry, the leading P-by-M part of this array must
contain the original feedthrough matrix D.
On exit, the leading P-by-M part of this array contains
the feedthrough matrix Dr of a reduced descriptor
realization without non-dynamic modes.
LDD INTEGER
The leading dimension of the array D. LDD >= MAX(1,P).
LR (output) INTEGER
The number of reduced differential equations.
NR (output) INTEGER
The dimension of the reduced descriptor state vector.
RANKE (output) INTEGER
The estimated rank of the matrix E.
INFRED (output) INTEGER
This parameter contains information on performed reduction
and on structure of resulting system matrices, as follows:
INFRED >= 0 the reduced system is in an SVD-like
coordinate form with Er upper triangular;
INFRED is the achieved order reduction.
INFRED < 0 no reduction achieved and the original
system has been restored.
Tolerances
TOL DOUBLE PRECISION
The tolerance to be used in rank determinations when
transforming (A-lambda*E). 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 = L*N*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, if MIN(L,N) = 0; otherwise,
LDWORK >= MAX( N+P, MIN(L,N)+MAX(3*N-1,M,L) ).
If LDWORK >= 2*L*N+L*M+N*P+
MAX( 1, N+P, MIN(L,N)+MAX(3*N-1,M,L) ) then
the original matrices are restored if no order reduction
is possible. This is achieved by saving system matrices
before reduction and restoring them if no order reduction
took place.
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. The optimal size does not necessarily include the
space needed for saving the original system matrices.
Error Indicator
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value.
Method
The subroutine elliminates the non-dynamics modes in two steps:
Step 1: Reduce the system to the SVD-like coordinate form
(Q'*A*Z-lambda*Q'*E*Z, Q'*B, C*Z) , where
( A11 A12 A13 ) ( E11 0 0 ) ( B1 )
Q'*A*Z = ( A21 A22 0 ), Q'*E*Z = ( 0 0 0 ), Q'*B = ( B2 ),
( A31 0 0 ) ( 0 0 0 ) ( B3 )
C*Z = ( C1 C2 C3 ),
where E11 and A22 are upper triangular invertible matrices.
Step 2: Compute the reduced system as (Ar-lambda*Er,Br,Cr,Dr),
where
( A11 - A12*inv(A22)*A21, A13 ) ( E11 0 )
Ar = ( ), Er = ( ),
( A31 0 ) ( 0 0 )
( B1 - A12*inv(A22)*B2 )
Br = ( ), Cr = ( C1 - C2*inv(A22)*A21, C3 ),
( B3 )
Dr = D - C2*inv(A22)*B2.
Step 3: If desired (JOBS = 'S'), reduce the descriptor system to
the standard form
Ar <- diag(inv(Er),I)*Ar; Br <- diag(inv(Er),I)*Br;
Er = diag(I,0).
If L = N and LR = NR = RANKE, then if Step 3 is performed,
the resulting system is a standard state space system.
Numerical Aspects
If L = N, the algorithm requires 0( N**3 ) floating point operations.Further Comments
NoneExample
Program Text
* TG01GD 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, LDD, LDE
PARAMETER ( LDA = LMAX, LDB = LMAX, LDC = PMAX,
$ LDD = PMAX, LDE = LMAX )
INTEGER LDWORK
PARAMETER ( LDWORK = MIN( LMAX, NMAX ) +
$ MAX( 3*NMAX - 1, MMAX, LMAX ) +
$ 2*LMAX*NMAX + LMAX*MMAX + PMAX*NMAX )
* .. Local Scalars ..
CHARACTER*1 JOBS
INTEGER I, INFO, INFRED, J, L, LR, M, N, NR, P, RANKE
DOUBLE PRECISION TOL
* .. Local Arrays ..
INTEGER IWORK(NMAX)
DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,NMAX),
$ D(LDD,MMAX), DWORK(LDWORK), E(LDE,NMAX)
* .. External Subroutines ..
EXTERNAL TG01GD
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN
* .. Executable Statements ..
*
WRITE ( NOUT, FMT = 99999 )
* Skip the heading in the data file and read the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) L, N, M, P, JOBS, TOL
IF ( L.LT.0 .OR. L.GT.LMAX ) THEN
WRITE ( NOUT, FMT = 99989 ) L
ELSE
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,L )
READ ( NIN, FMT = * ) ( ( E(I,J), J = 1,N ), I = 1,L )
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,L )
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 )
READ ( NIN, FMT = * ) ( ( D(I,J), J = 1,M ), I = 1,P )
* Find the reduced descriptor system
* (A-lambda E,B,C,D).
CALL TG01GD( JOBS, L, N, M, P, A, LDA, E, LDE, B, LDB,
$ C, LDC, D, LDD, LR, NR, RANKE, INFRED,
$ TOL, IWORK, DWORK, LDWORK, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99994 ) RANKE
WRITE ( NOUT, FMT = 99997 )
DO 10 I = 1, LR
WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,NR )
10 CONTINUE
WRITE ( NOUT, FMT = 99996 )
DO 20 I = 1, LR
WRITE ( NOUT, FMT = 99995 ) ( E(I,J), J = 1,NR )
20 CONTINUE
WRITE ( NOUT, FMT = 99993 )
DO 30 I = 1, LR
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,NR )
40 CONTINUE
WRITE ( NOUT, FMT = 99991 )
DO 50 I = 1, P
WRITE ( NOUT, FMT = 99995 ) ( D(I,J), J = 1,M )
50 CONTINUE
END IF
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' TG01GD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from TG01GD = ',I2)
99997 FORMAT (/' The reduced state dynamics matrix is ')
99996 FORMAT (/' The reduced descriptor matrix is ')
99995 FORMAT (20(1X,F8.4))
99994 FORMAT (' Rank of matrix E =', I5)
99993 FORMAT (/' The reduced input/state matrix is ')
99992 FORMAT (/' The reduced state/output matrix is ')
99991 FORMAT (/' The transformed feedthrough matrix 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)
END
Program Data
TG01GD EXAMPLE PROGRAM DATA
4 4 2 2 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
1 0
1 1
Program Results
TG01GD EXAMPLE PROGRAM RESULTS Rank of matrix E = 3 The reduced state dynamics matrix is 2.5102 -3.8550 -11.4533 -0.0697 0.0212 0.7015 0.3798 -0.1156 -3.8250 The reduced descriptor matrix is 10.1587 5.8230 1.3021 0.0000 -2.4684 -0.1896 0.0000 0.0000 1.0338 The reduced input/state matrix is 7.7100 1.6714 0.7678 1.1070 2.5428 0.6935 The reduced state/output matrix is 0.5477 -2.5000 -6.2610 -1.0954 1.0000 -0.8944 The transformed feedthrough matrix is 4.0000 1.0000 1.0000 1.0000