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

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,*)

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

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

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.

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

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.

[1] Laub, A.J. Efficient Calculation of Frequency Response Matrices from State-Space Models. ACM TOMS, 12, pp. 26-33, 1986.

3 The algorithm requires 0(N ) operations.

None

**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

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

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)