AB09HY

Cholesky factors of the controllability and observability Grammians

Purpose

  To compute the Cholesky factors Su and Ru of the controllability
  Grammian P = Su*Su' and observability Grammian Q = Ru'*Ru,
  respectively, satisfying

         A*P  + P*A' +  scalec^2*B*B'   = 0,       (1)

         A'*Q + Q*A  +  scaleo^2*Cw'*Cw = 0,       (2)

  where
         Cw = Hw - Bw'*X,
         Hw = inv(Dw)*C,
         Bw = (B*D' + P*C')*inv(Dw'),
         D*D' = Dw*Dw' (Dw upper triangular),

  and, with Aw = A - Bw*Hw, X is the stabilizing solution of the
  Riccati equation

         Aw'*X + X*Aw + Hw'*Hw + X*Bw*Bw'*X = 0.   (3)

  The P-by-M matrix D must have full row rank. Matrix A must be
  stable and in a real Schur form.

      SUBROUTINE AB09HY( N, M, P, A, LDA, B, LDB, C, LDC, D, LDD,
     $                   SCALEC, SCALEO, S, LDS, R, LDR, IWORK,
     $                   DWORK, LDWORK, BWORK, INFO )
C     .. Scalar Arguments ..
      INTEGER          INFO, LDA, LDB, LDC, LDD, LDR, LDS, LDWORK, M, N,
     $                 P
      DOUBLE PRECISION SCALEC, SCALEO
C     .. Array Arguments ..
      INTEGER          IWORK(*)
      DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*),
     $                 DWORK(*), R(LDR,*), S(LDS,*)
      LOGICAL          BWORK(*)

Input/Output Parameters

  N       (input) INTEGER
          The order of state-space representation, i.e.,
          the order of the matrix A.  N >= 0.

  M       (input) INTEGER
          The number of system inputs.  M >= 0.

  P       (input) INTEGER
          The number of system outputs.  M >= P >= 0.

  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
          The leading N-by-N part of this array must contain the
          stable state dynamics matrix A in a real Schur canonical
          form.

  LDA     INTEGER
          The leading dimension of 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, corresponding to the Schur matrix A.

  LDB     INTEGER
          The leading dimension of array B.  LDB >= MAX(1,N).

  C       (input) DOUBLE PRECISION array, dimension (LDC,N)
          The leading P-by-N part of this array must contain the
          state/output matrix C, corresponding to the Schur
          matrix A.

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

  D       (input) DOUBLE PRECISION array, dimension (LDD,M)
          The leading P-by-M part of this array must
          contain the full row rank input/output matrix D.

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

  SCALEC  (output) DOUBLE PRECISION
          Scaling factor for the controllability Grammian in (1).

  SCALEO  (output) DOUBLE PRECISION
          Scaling factor for the observability Grammian in (2).

  S       (output) DOUBLE PRECISION array, dimension (LDS,N)
          The leading N-by-N upper triangular part of this array
          contains the Cholesky factor Su of the cotrollability
          Grammian P = Su*Su' satisfying (1).

  LDS     INTEGER
          The leading dimension of array S.  LDS >= MAX(1,N).

  R       (output) DOUBLE PRECISION array, dimension (LDR,N)
          The leading N-by-N upper triangular part of this array
          contains the Cholesky factor Ru of the observability
          Grammian Q = Ru'*Ru satisfying (2).

  LDR     INTEGER
          The leading dimension of array R.  LDR >= MAX(1,N).

  IWORK   INTEGER array, dimension (2*N)

  DWORK   DOUBLE PRECISION array, dimension (LDWORK)
          On exit, if INFO = 0, DWORK(1) returns the optimal value
          of LDWORK and DWORK(2) contains RCOND, the reciprocal
          condition number of the U11 matrix from the expression
          used to compute X = U21*inv(U11). A small value RCOND
          indicates possible ill-conditioning of the Riccati
          equation (3).

  LDWORK  INTEGER
          The length of the array DWORK.
          LDWORK >= MAX( 2, N*(MAX(N,M,P)+5),
                         2*N*P+MAX(P*(M+2),10*N*(N+1) ) ).
          For optimum performance LDWORK should be larger.

  BWORK   LOGICAL array, dimension 2*N

  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value;
          = 1:  the state matrix A is not stable or is not in a
                real Schur form;
          = 2:  the reduction of Hamiltonian matrix to real Schur
                form failed;
          = 3:  the reordering of the real Schur form of the
                Hamiltonian matrix failed;
          = 4:  the Hamiltonian matrix has less than N stable
                eigenvalues;
          = 5:  the coefficient matrix U11 in the linear system
                X*U11 = U21, used to determine X, is singular to
                working precision;
          = 6:  the feedthrough matrix D has not a full row rank P.

  None

Program Text

  None

  None

  None