AB09JV

State-space representation of a projection of a left weighted transfer-function matrix

Purpose

  To construct a state-space representation (A,BS,CS,DS) of the
  projection of V*G or conj(V)*G containing the poles of G, from the
  state-space representations (A,B,C,D) and (AV-lambda*EV,BV,CV,DV),
  of the transfer-function matrices G and V, respectively.
  G is assumed to be a stable transfer-function matrix and
  the state matrix A must be in a real Schur form.
  When computing the stable projection of V*G, it is assumed
  that G and V have completely distinct poles.
  When computing the stable projection of conj(V)*G, it is assumed
  that G and conj(V) have completely distinct poles.

  Note: For a transfer-function matrix G, conj(G) denotes the
  conjugate of G given by G'(-s) for a continuous-time system or
  G'(1/z) for a discrete-time system.

      SUBROUTINE AB09JV( JOB, DICO, JOBEV, STBCHK, N, M, P, NV, PV,
     $                   A, LDA, B, LDB, C, LDC, D, LDD, AV, LDAV,
     $                   EV, LDEV, BV, LDBV, CV, LDCV, DV, LDDV, IWORK,
     $                   DWORK, LDWORK, INFO )
C     .. Scalar Arguments ..
      CHARACTER         DICO, JOB, JOBEV, STBCHK
      INTEGER           INFO, LDA, LDAV, LDB, LDBV, LDC, LDCV,
     $                  LDD, LDDV, LDEV, LDWORK, M, N, NV, P, PV
C     .. Array Arguments ..
      INTEGER           IWORK(*)
      DOUBLE PRECISION  A(LDA,*), AV(LDAV,*), B(LDB,*), BV(LDBV,*),
     $                  C(LDC,*), CV(LDCV,*), D(LDD,*), DV(LDDV,*),
     $                  DWORK(*), EV(LDEV,*)

Mode Parameters

  JOB     CHARACTER*1
          Specifies the projection to be computed as follows:
          = 'V':  compute the projection of V*G containing
                  the poles of G;
          = 'C':  compute the projection of conj(V)*G containing
                  the poles of G.

  DICO    CHARACTER*1
          Specifies the type of the systems as follows:
          = 'C':  G and V are continuous-time systems;
          = 'D':  G and V are discrete-time systems.

  JOBEV   CHARACTER*1
          Specifies whether EV is a general square or an identity
          matrix as follows:
          = 'G':  EV is a general square matrix;
          = 'I':  EV is the identity matrix.

  STBCHK  CHARACTER*1
          Specifies whether stability/antistability of V is to be
          checked as follows:
          = 'C':  check stability if JOB = 'C' or antistability if
                  JOB = 'V';
          = 'N':  do not check stability or antistability.

  N       (input) INTEGER
          The dimension of the state vector of the system with
          the transfer-function matrix G.  N >= 0.

  M       (input) INTEGER
          The dimension of the input vector of the system with
          the transfer-function matrix G.  M >= 0.

  P       (input) INTEGER
          The dimension of the output vector of the system with the
          transfer-function matrix G, and also the dimension of
          the input vector if JOB = 'V', or of the output vector
          if JOB = 'C', of the system with the transfer-function
          matrix V.  P >= 0.

  NV      (input) INTEGER
          The dimension of the state vector of the system with
          the transfer-function matrix V.  NV >= 0.

  PV      (input) INTEGER
          The dimension of the output vector, if JOB = 'V', or
          of the input vector, if JOB = 'C', of the system with
          the transfer-function matrix V.  PV >= 0.

  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
          The leading N-by-N part of this array must contain the
          state matrix A of the system with the transfer-function
          matrix G in a real Schur form.

  LDA     INTEGER
          The leading dimension of the array A.  LDA >= MAX(1,N).

  B       (input) DOUBLE PRECISION array, dimension (LDB,M)
          The leading N-by-M part of this array must contain
          the input/state matrix B of the system with the
          transfer-function matrix G. The matrix BS is equal to B.

  LDB     INTEGER
          The leading dimension of the 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 output matrix C of the system with the
          transfer-function matrix G.
          On exit, if INFO = 0, the leading PV-by-N part of this
          array contains the output matrix CS of the projection of
          V*G, if JOB = 'V', or of conj(V)*G, if JOB = 'C'.

  LDC     INTEGER
          The leading dimension of the array C.  LDC >= MAX(1,P,PV).

  D       (input/output) DOUBLE PRECISION array, dimension (LDD,M)
          On entry, the leading P-by-M part of this array must
          contain the feedthrough matrix D of the system with the
          transfer-function matrix G.
          On exit, if INFO = 0, the leading PV-by-M part of
          this array contains the feedthrough matrix DS of the
          projection of V*G, if JOB = 'V', or of conj(V)*G,
          if JOB = 'C'.

  LDD     INTEGER
          The leading dimension of the array D.  LDD >= MAX(1,P,PV).

  AV      (input/output) DOUBLE PRECISION array, dimension (LDAV,NV)
          On entry, the leading NV-by-NV part of this array must
          contain the state matrix AV of the system with the
          transfer-function matrix V.
          On exit, if INFO = 0, the leading NV-by-NV part of this
          array contains a condensed matrix as follows:
          if JOBEV = 'I', it contains the real Schur form of AV;
          if JOBEV = 'G' and JOB = 'V', it contains a quasi-upper
          triangular matrix representing the real Schur matrix
          in the real generalized Schur form of the pair (AV,EV);
          if JOBEV = 'G', JOB = 'C' and DICO = 'C', it contains a
          quasi-upper triangular matrix corresponding to the
          generalized real Schur form of the pair (AV',EV');
          if JOBEV = 'G', JOB = 'C' and DICO = 'D', it contains an
          upper triangular matrix corresponding to the generalized
          real Schur form of the pair (EV',AV').

  LDAV    INTEGER
          The leading dimension of the array AV.  LDAV >= MAX(1,NV).

  EV      (input/output) DOUBLE PRECISION array, dimension (LDEV,NV)
          On entry, if JOBEV = 'G', the leading NV-by-NV part of
          this array must contain the descriptor matrix EV of the
          system with the transfer-function matrix V.
          If JOBEV = 'I', EV is assumed to be an identity matrix
          and is not referenced.
          On exit, if INFO = 0 and JOBEV = 'G', the leading NV-by-NV
          part of this array contains a condensed matrix as follows:
          if JOB = 'V', it contains an upper triangular matrix
          corresponding to the real generalized Schur form of the
          pair (AV,EV);
          if JOB = 'C' and DICO = 'C', it contains an upper
          triangular matrix corresponding to the generalized real
          Schur form of the pair (AV',EV');
          if JOB = 'C' and DICO = 'D', it contains a quasi-upper
          triangular matrix corresponding to the generalized
          real Schur form of the pair (EV',AV').

  LDEV    INTEGER
          The leading dimension of the array EV.
          LDEV >= MAX(1,NV), if JOBEV = 'G';
          LDEV >= 1,         if JOBEV = 'I'.

  BV      (input/output) DOUBLE PRECISION array,
          dimension (LDBV,MBV), where MBV = P, if JOB = 'V', and
          MBV = PV, if JOB = 'C'.
          On entry, the leading NV-by-MBV part of this array must
          contain the input matrix BV of the system with the
          transfer-function matrix V.
          On exit, if INFO = 0, the leading NV-by-MBV part of this
          array contains Q'*BV, where Q is the orthogonal matrix
          that reduces AV to the real Schur form or the left
          orthogonal matrix used to reduce the pair (AV,EV),
          (AV',EV') or (EV',AV') to the generalized real Schur form.

  LDBV    INTEGER
          The leading dimension of the array BV.  LDBV >= MAX(1,NV).

  CV      (input/output) DOUBLE PRECISION array, dimension (LDCV,NV)
          On entry, the leading PCV-by-NV part of this array must
          contain the output matrix CV of the system with the
          transfer-function matrix V, where PCV = PV, if JOB = 'V',
          or PCV = P, if JOB = 'C'.
          On exit, if INFO = 0, the leading PCV-by-NV part of this
          array contains CV*Q, where Q is the orthogonal matrix that
          reduces AV to the real Schur form, or CV*Z, where Z is the
          right orthogonal matrix used to reduce the pair (AV,EV),
          (AV',EV') or (EV',AV') to the generalized real Schur form.

  LDCV    INTEGER
          The leading dimension of the array CV.
          LDCV >= MAX(1,PV) if JOB = 'V';
          LDCV >= MAX(1,P)  if JOB = 'C'.

  DV      (input) DOUBLE PRECISION array,
          dimension (LDDV,MBV), where MBV = P, if JOB = 'V', and
          MBV = PV, if JOB = 'C'.
          The leading PCV-by-MBV part of this array must contain
          the feedthrough matrix DV of the system with the
          transfer-function matrix V, where PCV = PV, if JOB = 'V',
          or PCV = P, if JOB = 'C'.

  LDDV    INTEGER
          The leading dimension of the array DV.
          LDDV >= MAX(1,PV) if JOB = 'V';
          LDDV >= MAX(1,P)  if JOB = 'C'.

  IWORK   INTEGER array, dimension (LIWORK)
          LIWORK =   0,    if JOBEV = 'I';
          LIWORK = NV+N+6, if JOBEV = 'G'.

  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 >= LW1, if JOBEV = 'I',
          LDWORK >= LW2, if JOBEV = 'G', where
            LW1 = MAX( 1, NV*(NV+5), NV*N + MAX( a, PV*N, PV*M ) )
                  a = 0,    if DICO = 'C' or  JOB = 'V',
                  a = 2*NV, if DICO = 'D' and JOB = 'C';
            LW2 = MAX( 2*NV*NV + MAX( 11*NV+16, P*NV, PV*NV ),
                       NV*N + MAX( NV*N+N*N, PV*N, PV*M ) ).
          For good performance, LDWORK should be larger.

  INFO    INTEGER
          =  0:  successful exit;
          <  0:  if INFO = -i, the i-th argument had an illegal
                 value;
          =  1:  the reduction of the pair (AV,EV) to the real
                 generalized Schur form failed (JOBEV = 'G'),
                 or the reduction of the matrix AV to the real
                 Schur form failed (JOBEV = 'I);
          =  2:  the solution of the Sylvester equation failed
                 because the matrix A and the pencil AV-lambda*EV
                 have common eigenvalues (if JOB = 'V'), or the
                 pencil -AV-lambda*EV and A have common eigenvalues
                 (if JOB = 'C' and DICO = 'C'), or the pencil
                 AV-lambda*EV has an eigenvalue which is the
                 reciprocal of one of eigenvalues of A
                 (if JOB = 'C' and DICO = 'D');
          =  3:  the solution of the Sylvester equation failed
                 because the matrices A and AV have common
                 eigenvalues (if JOB = 'V'), or the matrices A
                 and -AV have common eigenvalues (if JOB = 'C' and
                 DICO = 'C'), or the matrix A has an eigenvalue
                 which is the reciprocal of one of eigenvalues of AV
                 (if JOB = 'C' and DICO = 'D');
          =  4:  JOB = 'V' and the pair (AV,EV) has not completely
                 unstable generalized eigenvalues, or JOB = 'C' and
                 the pair (AV,EV) has not completely stable
                 generalized eigenvalues.

  If JOB = 'V', the matrices of the stable projection of V*G are
  computed as

    BS = B,  CS = CV*X + DV*C,  DS = DV*D,

  where X satisfies the generalized Sylvester equation

    AV*X - EV*X*A + BV*C = 0.

  If JOB = 'C', the matrices of the stable projection of conj(V)*G
  are computed using the following formulas:

  - for a continuous-time system, the matrices BS, CS and DS of
    the stable projection are computed as

      BS = B,  CS = BV'*X + DV'*C,  DS = DV'*D,

    where X satisfies the generalized Sylvester equation

      AV'*X + EV'*X*A + CV'*C = 0.

  - for a discrete-time system, the matrices BS, CS and DS of
    the stable projection are computed as

      BS = B,  CS = BV'*X*A + DV'*C,  DS = DV'*D + BV'*X*B,

    where X satisfies the generalized Sylvester equation

      EV'*X - AV'*X*A = CV'*C.

  [1] Varga, A.
      Efficient and numerically reliable implementation of the
      frequency-weighted Hankel-norm approximation model reduction
      approach.
      Proc. 2001 ECC, Porto, Portugal, 2001.

  [2] Zhou, K.
      Frequency-weighted H-infinity norm and optimal Hankel norm
      model reduction.
      IEEE Trans. Autom. Control, vol. 40, pp. 1687-1699, 1995.

  The implemented methods rely on numerically stable algorithms.

  None

Program Text

  None

  None

  None