Purpose
To find the complex frequency response matrix (transfer matrix)
G(freq) of the state-space representation (A,B,C) given by
-1
G(freq) = C * ((freq*I - A) ) * B
where A, B and C are real N-by-N, N-by-M and P-by-N matrices
respectively and freq is a complex scalar.
Specification
SUBROUTINE TB05AD( BALEIG, INITA, N, M, P, FREQ, A, LDA, B, LDB,
$ C, LDC, RCOND, G, LDG, EVRE, EVIM, HINVB,
$ LDHINV, IWORK, DWORK, LDWORK, ZWORK, LZWORK,
$ INFO )
C .. Scalar Arguments ..
CHARACTER BALEIG, INITA
INTEGER INFO, LDA, LDB, LDC, LDG, LDHINV, LDWORK,
$ LZWORK, M, N, P
DOUBLE PRECISION RCOND
COMPLEX*16 FREQ
C .. Array Arguments ..
INTEGER IWORK(*)
DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*), EVIM(*),
$ EVRE(*)
COMPLEX*16 ZWORK(*), G(LDG,*), HINVB(LDHINV,*)
Arguments
Mode Parameters
BALEIG CHARACTER*1
Determines whether the user wishes to balance matrix A
and/or compute its eigenvalues and/or estimate the
condition number of the problem as follows:
= 'N': The matrix A should not be balanced and neither
the eigenvalues of A nor the condition number
estimate of the problem are to be calculated;
= 'C': The matrix A should not be balanced and only an
estimate of the condition number of the problem
is to be calculated;
= 'B' or 'E' and INITA = 'G': The matrix A is to be
balanced and its eigenvalues calculated;
= 'A' and INITA = 'G': The matrix A is to be balanced,
and its eigenvalues and an estimate of the
condition number of the problem are to be
calculated.
INITA CHARACTER*1
Specifies whether or not the matrix A is already in upper
Hessenberg form as follows:
= 'G': The matrix A is a general matrix;
= 'H': The matrix A is in upper Hessenberg form and
neither balancing nor the eigenvalues of A are
required.
INITA must be set to 'G' for the first call to the
routine, unless the matrix A is already in upper
Hessenberg form and neither balancing nor the eigenvalues
of A are required. Thereafter, it must be set to 'H' for
all subsequent calls.
Input/Output Parameters
N (input) INTEGER
The number of states, i.e. the order of the state
transition matrix A. N >= 0.
M (input) INTEGER
The number of inputs, i.e. the number of columns in the
matrix B. M >= 0.
P (input) INTEGER
The number of outputs, i.e. the number of rows in the
matrix C. P >= 0.
FREQ (input) COMPLEX*16
The frequency freq at which the frequency response matrix
(transfer matrix) is to be evaluated.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading N-by-N part of this array must
contain the state transition matrix A.
If INITA = 'G', then, on exit, the leading N-by-N part of
this array contains an upper Hessenberg matrix similar to
(via an orthogonal matrix consisting of a sequence of
Householder transformations) the original state transition
matrix A.
LDA INTEGER
The leading dimension of array A. LDA >= 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/state matrix B.
If INITA = 'G', then, on exit, the leading N-by-M part of
this array contains the product of the transpose of the
orthogonal transformation matrix used to reduce A to upper
Hessenberg form and the original input/state matrix B.
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.
If INITA = 'G', then, on exit, the leading P-by-N part of
this array contains the product of the original output/
state matrix C and the orthogonal transformation matrix
used to reduce A to upper Hessenberg form.
LDC INTEGER
The leading dimension of array C. LDC >= MAX(1,P).
RCOND (output) DOUBLE PRECISION
If BALEIG = 'C' or BALEIG = 'A', then RCOND contains an
estimate of the reciprocal of the condition number of
matrix H with respect to inversion (see METHOD).
G (output) COMPLEX*16 array, dimension (LDG,M)
The leading P-by-M part of this array contains the
frequency response matrix G(freq).
LDG INTEGER
The leading dimension of array G. LDG >= MAX(1,P).
EVRE, (output) DOUBLE PRECISION arrays, dimension (N)
EVIM If INITA = 'G' and BALEIG = 'B' or 'E' or BALEIG = 'A',
then these arrays contain the real and imaginary parts,
respectively, of the eigenvalues of the matrix A.
Otherwise, these arrays are not referenced.
HINVB (output) COMPLEX*16 array, dimension (LDHINV,M)
The leading N-by-M part of this array contains the
-1
product H B.
LDHINV INTEGER
The leading dimension of array HINVB. LDHINV >= MAX(1,N).
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 >= MAX(1, N - 1 + MAX(N,M,P)),
if INITA = 'G' and BALEIG = 'N', or 'B', or 'E';
LDWORK >= MAX(1, N + MAX(N,M-1,P-1)),
if INITA = 'G' and BALEIG = 'C', or 'A';
LDWORK >= MAX(1, 2*N),
if INITA = 'H' and BALEIG = 'C', or 'A';
LDWORK >= 1, otherwise.
For optimum performance when INITA = 'G' LDWORK should be
larger.
ZWORK COMPLEX*16 array, dimension (LZWORK)
LZWORK INTEGER
The length of the array ZWORK.
LZWORK >= MAX(1,N*N+2*N), if BALEIG = 'C', or 'A';
LZWORK >= MAX(1,N*N), otherwise.
Error Indicator
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value;
= 1: if more than 30*N iterations are required to
isolate all the eigenvalues of the matrix A; the
computations are continued;
= 2: if either FREQ is too near to an eigenvalue of the
matrix A, or RCOND is less than EPS, where EPS is
the machine precision (see LAPACK Library routine
DLAMCH).
Method
The matrix A is first balanced (if BALEIG = 'B' or 'E', or
BALEIG = 'A') and then reduced to upper Hessenberg form; the same
transformations are applied to the matrix B and the matrix C.
The complex Hessenberg matrix H = (freq*I - A) is then used
-1
to solve for C * H * B.
Depending on the input values of parameters BALEIG and INITA,
the eigenvalues of matrix A and the condition number of
matrix H with respect to inversion are also calculated.
References
[1] Laub, A.J.
Efficient Calculation of Frequency Response Matrices from
State-Space Models.
ACM TOMS, 12, pp. 26-33, 1986.
Numerical Aspects
3 The algorithm requires 0(N ) operations.Further Comments
NoneExample
Program Text
* TB05AD 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, LDG, LDHINV
PARAMETER ( LDA = NMAX, LDB = NMAX, LDC = PMAX, LDG = PMAX,
$ LDHINV = NMAX )
INTEGER LIWORK
PARAMETER ( LIWORK = NMAX )
INTEGER LDWORK
PARAMETER ( LDWORK = 2*NMAX )
INTEGER LZWORK
PARAMETER ( LZWORK = NMAX*( NMAX+2 ) )
* .. Local Scalars ..
COMPLEX*16 FREQ
DOUBLE PRECISION RCOND
INTEGER I, INFO, J, M, N, P
CHARACTER*1 BALEIG, INITA
LOGICAL LBALBA, LBALEA, LBALEB, LBALEC, LINITA
* .. Local Arrays ..
COMPLEX*16 G(LDG,MMAX), HINVB(LDHINV,MMAX), ZWORK(LZWORK)
DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,NMAX),
$ DWORK(LDWORK), EVIM(NMAX), EVRE(NMAX)
INTEGER IWORK(LIWORK)
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL TB05AD
* .. 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, FREQ, INITA, BALEIG
LBALEC = LSAME( BALEIG, 'C' )
LBALEB = LSAME( BALEIG, 'B' ) .OR. LSAME( BALEIG, 'E' )
LBALEA = LSAME( BALEIG, 'A' )
LBALBA = LBALEB.OR.LBALEA
LINITA = LSAME( INITA, 'G' )
IF ( N.LE.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99992 ) N
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
IF ( M.LE.0 .OR. M.GT.MMAX ) THEN
WRITE ( NOUT, FMT = 99991 ) M
ELSE
READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1,N )
IF ( P.LE.0 .OR. P.GT.PMAX ) THEN
WRITE ( NOUT, FMT = 99990 ) P
ELSE
READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,P )
* Find the frequency response matrix of the ssr (A,B,C).
CALL TB05AD( BALEIG, INITA, N, M, P, FREQ, A, LDA, B,
$ LDB, C, LDC, RCOND, G, LDG, EVRE, EVIM,
$ HINVB, LDHINV, IWORK, DWORK, LDWORK, ZWORK,
$ LZWORK, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
IF ( ( LBALEC ) .OR. ( LBALEA ) ) WRITE ( NOUT,
$ FMT = 99997 ) RCOND
IF ( ( LINITA ) .AND. ( LBALBA ) )
$ WRITE ( NOUT, FMT = 99996 )
$ ( EVRE(I), EVIM(I), I = 1,N )
WRITE ( NOUT, FMT = 99995 )
DO 20 I = 1, P
WRITE ( NOUT, FMT = 99994 ) ( G(I,J), J = 1,M )
20 CONTINUE
WRITE ( NOUT, FMT = 99993 )
DO 40 I = 1, N
WRITE ( NOUT, FMT = 99994 ) ( HINVB(I,J), J = 1,M )
40 CONTINUE
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' TB05AD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from TB05AD = ',I2)
99997 FORMAT (' RCOND = ',F4.2)
99996 FORMAT (/' Eigenvalues of the state transmission matrix A are ',
$ /(1X,2F7.2,'*j'))
99995 FORMAT (/' The frequency response matrix G(freq) is ')
99994 FORMAT (20(' (',F5.2,',',F5.2,') ',:))
99993 FORMAT (/' H(inverse)*B is ')
99992 FORMAT (/' N is out of range.',/' N = ',I5)
99991 FORMAT (/' M is out of range.',/' M = ',I5)
99990 FORMAT (/' P is out of range.',/' P = ',I5)
END
Program Data
TB05AD EXAMPLE PROGRAM DATA 3 1 2 (0.0,0.5) G A 1.0 2.0 0.0 4.0 -1.0 0.0 0.0 0.0 1.0 1.0 0.0 1.0 1.0 0.0 -1.0 0.0 0.0 1.0Program Results
TB05AD EXAMPLE PROGRAM RESULTS
RCOND = 0.22
Eigenvalues of the state transmission matrix A are
3.00 0.00*j
-3.00 0.00*j
1.00 0.00*j
The frequency response matrix G(freq) is
( 0.69, 0.35)
(-0.80,-0.40)
H(inverse)*B is
(-0.11,-0.05)
(-0.43, 0.00)
(-0.80,-0.40)