TB05AD

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

Return to index