### Frequency response matrix of a given state-space representation (A,B,C)

[Specification] [Arguments] [Method] [References] [Comments] [Example]

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
```   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.

```
```  None
```
Example

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 ..
*     .. 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.0
```
Program 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)
```