IB01CD

Estimating the initial state and system matrices B and D using A, B, and input-output data (driver)

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

Purpose

  To estimate the initial state and, optionally, the system matrices
  B  and  D  of a linear time-invariant (LTI) discrete-time system,
  given the system matrices  (A,B,C,D),  or (when  B  and  D  are
  estimated) only the matrix pair  (A,C),  and the input and output
  trajectories of the system. The model structure is :

        x(k+1) = Ax(k) + Bu(k),   k >= 0,
        y(k)   = Cx(k) + Du(k),

  where  x(k)  is the  n-dimensional state vector (at time k),
         u(k)  is the  m-dimensional input vector,
         y(k)  is the  l-dimensional output vector,
  and  A, B, C, and D  are real matrices of appropriate dimensions.
  The input-output data can internally be processed sequentially.

Specification
      SUBROUTINE IB01CD( JOBX0, COMUSE, JOB, N, M, L, NSMP, A, LDA, B,
     $                   LDB, C, LDC, D, LDD, U, LDU, Y, LDY, X0, V,
     $                   LDV, TOL, IWORK, DWORK, LDWORK, IWARN, INFO )
C     .. Scalar Arguments ..
      DOUBLE PRECISION   TOL
      INTEGER            INFO, IWARN, L, LDA, LDB, LDC, LDD, LDU, LDV,
     $                   LDWORK, LDY, M, N, NSMP
      CHARACTER          COMUSE, JOB, JOBX0
C     .. Array Arguments ..
      DOUBLE PRECISION   A(LDA, *), B(LDB, *), C(LDC, *), D(LDD, *),
     $                   DWORK(*),  U(LDU, *), V(LDV, *), X0(*),
     $                   Y(LDY, *)
      INTEGER            IWORK(*)

Arguments

Mode Parameters

  JOBX0   CHARACTER*1
          Specifies whether or not the initial state should be
          computed, as follows:
          = 'X':  compute the initial state x(0);
          = 'N':  do not compute the initial state (possibly,
                  because x(0) is known to be zero).

  COMUSE  CHARACTER*1
          Specifies whether the system matrices B and D should be
          computed or used, as follows:
          = 'C':  compute the system matrices B and D, as specified
                  by JOB;
          = 'U':  use the system matrices B and D, as specified by
                  JOB;
          = 'N':  do not compute/use the matrices B and D.
          If  JOBX0 = 'N'  and  COMUSE <> 'N',  then  x(0)  is set
          to zero.
          If  JOBX0 = 'N'  and  COMUSE =  'N',  then  x(0)  is
          neither computed nor set to zero.

  JOB     CHARACTER*1
          If  COMUSE = 'C'  or  'U',  specifies which of the system
          matrices  B and D  should be computed or used, as follows:
          = 'B':  compute/use the matrix B only (D is known to be
                  zero);
          = 'D':  compute/use the matrices B and D.
          The value of  JOB  is irrelevant if  COMUSE = 'N'  or if
          JOBX0 = 'N'  and  COMUSE = 'U'.
          The combinations of options, the data used, and the
          returned results, are given in the table below, where
          '*'  denotes an irrelevant value.

           JOBX0   COMUSE    JOB     Data used    Returned results
          ----------------------------------------------------------
             X       C        B       A,C,u,y          x,B
             X       C        D       A,C,u,y          x,B,D
             N       C        B       A,C,u,y          x=0,B
             N       C        D       A,C,u,y          x=0,B,D
          ----------------------------------------------------------
             X       U        B      A,B,C,u,y            x
             X       U        D      A,B,C,D,u,y          x
             N       U        *          -               x=0
          ----------------------------------------------------------
             X       N        *        A,C,y              x
             N       N        *          -                -
          ----------------------------------------------------------

          For  JOBX0 = 'N'  and  COMUSE = 'N',  the routine just
          sets  DWORK(1)  to 2 and  DWORK(2)  to 1, and returns
          (see the parameter DWORK).

Input/Output Parameters
  N       (input) INTEGER
          The order of the system.  N >= 0.

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

  L       (input) INTEGER
          The number of system outputs.  L > 0.

  NSMP    (input) INTEGER
          The number of rows of matrices  U  and  Y  (number of
          samples,  t).
          NSMP >= 0,            if  JOBX0 = 'N'  and  COMUSE <> 'C';
          NSMP >= N,            if  JOBX0 = 'X'  and  COMUSE <> 'C';
          NSMP >= N*M + a + e,  if  COMUSE = 'C',
          where   a = 0,  if  JOBX0 = 'N';
                  a = N,  if  JOBX0 = 'X';
                  e = 0,  if  JOBX0 = 'X'  and  JOB = 'B';
                  e = 1,  if  JOBX0 = 'N'  and  JOB = 'B';
                  e = M,  if  JOB   = 'D'.

  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
          If  JOBX0 = 'X'  or  COMUSE = 'C',  the leading N-by-N
          part of this array must contain the system state matrix A.
          If  N = 0,  or  JOBX0 = 'N'  and  COMUSE <> 'C',  this
          array is not referenced.

  LDA     INTEGER
          The leading dimension of the array A.
          LDA >= MAX(1,N),  if  JOBX0 = 'X'  or   COMUSE =  'C';
          LDA >= 1,         if  JOBX0 = 'N'  and  COMUSE <> 'C'.

  B       (input or output) DOUBLE PRECISION array, dimension
          (LDB,M)
          If  JOBX0 = 'X'  and  COMUSE = 'U',  B  is an input
          parameter and, on entry, the leading N-by-M part of this
          array must contain the system input matrix  B.
          If  COMUSE = 'C',  B  is an output parameter and, on exit,
          if  INFO = 0,  the leading N-by-M part of this array
          contains the estimated system input matrix  B.
          If  min(N,M) = 0,  or  JOBX0 = 'N'  and  COMUSE = 'U',
          or  COMUSE = 'N',  this array is not referenced.

  LDB     INTEGER
          The leading dimension of the array B.
          LDB >= MAX(1,N),  if  M > 0,  COMUSE = 'U',  JOBX0 = 'X',
                            or  M > 0,  COMUSE = 'C';
          LDB >= 1,         if  min(N,M) = 0,  or  COMUSE = 'N',
                            or  JOBX0  = 'N'  and  COMUSE = 'U'.

  C       (input) DOUBLE PRECISION array, dimension (LDC,N)
          If  JOBX0 = 'X'  or  COMUSE = 'C',  the leading L-by-N
          part of this array must contain the system output
          matrix  C.
          If  N = 0,  or  JOBX0 = 'N'  and  COMUSE <> 'C',  this
          array is not referenced.

  LDC     INTEGER
          The leading dimension of the array C.
          LDC >= L,  if  N > 0, and  JOBX0 = 'X'  or  COMUSE = 'C';
          LDC >= 1,  if  N = 0, or  JOBX0 = 'N'  and  COMUSE <> 'C'.

  D       (input or output) DOUBLE PRECISION array, dimension
          (LDD,M)
          If  JOBX0 = 'X',  COMUSE = 'U',  and  JOB = 'D',  D  is an
          input parameter and, on entry, the leading L-by-M part of
          this array must contain the system input-output matrix  D.
          If  COMUSE = 'C'  and  JOB = 'D',  D  is an output
          parameter and, on exit, if  INFO = 0,  the leading
          L-by-M part of this array contains the estimated system
          input-output matrix  D.
          If  M = 0,  or  JOBX0 = 'N'  and  COMUSE = 'U',  or
          COMUSE = 'N',  or  JOB = 'B',  this array is not
          referenced.

  LDD     INTEGER
          The leading dimension of the array D.
          LDD >= L,  if  M > 0,   JOBX0 = 'X',  COMUSE = 'U',  and
                                                JOB = 'D',  or
                     if  M > 0,  COMUSE = 'C',  and  JOB = 'D';
          LDD >= 1,  if  M = 0,  or  JOBX0 = 'N'  and  COMUSE = 'U',
                     or  COMUSE = 'N',  or  JOB = 'B'.

  U       (input or input/output) DOUBLE PRECISION array, dimension
          (LDU,M)
          On entry, if  COMUSE = 'C',  or  JOBX0 = 'X'  and
          COMUSE = 'U',  the leading NSMP-by-M part of this array
          must contain the t-by-m input-data sequence matrix  U,
          U = [u_1 u_2 ... u_m].  Column  j  of  U  contains the
          NSMP  values of the j-th input component for consecutive
          time increments.
          On exit, if  COMUSE = 'C'  and  JOB = 'D',  the leading
          NSMP-by-M part of this array contains details of the
          QR factorization of the t-by-m matrix  U,  possibly
          computed sequentially (see METHOD).
          If  COMUSE = 'C'  and  JOB = 'B',  or  COMUSE = 'U',  this
          array is unchanged on exit.
          If  M = 0,  or  JOBX0 = 'N'  and  COMUSE = 'U',  or
          COMUSE = 'N',  this array is not referenced.

  LDU     INTEGER
          The leading dimension of the array U.
          LDU >= MAX(1,NSMP),  if  M > 0    and  COMUSE = 'C'  or
                               JOBX0 = 'X'  and  COMUSE = 'U;
          LDU >= 1,            if  M = 0,   or   COMUSE = 'N',  or
                               JOBX0 = 'N'  and  COMUSE = 'U'.

  Y       (input) DOUBLE PRECISION array, dimension (LDY,L)
          On entry, if  JOBX0 = 'X'  or  COMUSE = 'C',  the leading
          NSMP-by-L part of this array must contain the t-by-l
          output-data sequence matrix  Y,  Y = [y_1 y_2 ... y_l].
          Column  j  of  Y  contains the  NSMP  values of the j-th
          output component for consecutive time increments.
          If  JOBX0 = 'N'  and  COMUSE <> 'C',  this array is not
          referenced.

  LDY     INTEGER
          The leading dimension of the array Y.
          LDY >= MAX(1,NSMP),  if  JOBX0 = 'X'  or   COMUSE = 'C;
          LDY >= 1,            if  JOBX0 = 'N'  and  COMUSE <> 'C'.

  X0      (output) DOUBLE PRECISION array, dimension (N)
          If  INFO = 0  and  JOBX0 = 'X',  this array contains the
          estimated initial state of the system,  x(0).
          If  JOBX0 = 'N'  and  COMUSE = 'C',  this array is used as
          workspace and finally it is set to zero.
          If  JOBX0 = 'N'  and  COMUSE = 'U',  then  x(0)  is set to
          zero without any calculations.
          If  JOBX0 = 'N'  and  COMUSE = 'N',  this array is not
          referenced.

  V       (output) DOUBLE PRECISION array, dimension (LDV,N)
          On exit, if  INFO = 0  or 2,  JOBX0 = 'X'  or
          COMUSE = 'C',  the leading N-by-N part of this array
          contains the orthogonal matrix V of a real Schur
          factorization of the matrix  A.
          If  JOBX0 = 'N'  and  COMUSE <> 'C',  this array is not
          referenced.

  LDV     INTEGER
          The leading dimension of the array V.
          LDV >= MAX(1,N),  if  JOBX0 = 'X'  or   COMUSE =  'C;
          LDV >= 1,         if  JOBX0 = 'N'  and  COMUSE <> 'C'.

Tolerances
  TOL     DOUBLE PRECISION
          The tolerance to be used for estimating the rank of
          matrices. If the user sets  TOL > 0,  then the given value
          of  TOL  is used as a lower bound for the reciprocal
          condition number;  a matrix whose estimated condition
          number is less than  1/TOL  is considered to be of full
          rank.  If the user sets  TOL <= 0,  then  EPS  is used
          instead, where  EPS  is the relative machine precision
          (see LAPACK Library routine DLAMCH).  TOL <= 1.

Workspace
  IWORK   INTEGER array, dimension (LIWORK), where
          LIWORK >= 0,          if  JOBX0 = 'N'  and  COMUSE <> 'C';
          LIWORK >= N,          if  JOBX0 = 'X'  and  COMUSE <> 'C';
          LIWORK >= N*M + a,        if COMUSE = 'C' and JOB = 'B',
          LIWORK >= max(N*M + a,M), if COMUSE = 'C' and JOB = 'D',
          with  a = 0,  if  JOBX0 = 'N';
                a = N,  if  JOBX0 = 'X'.

  DWORK   DOUBLE PRECISION array, dimension (LDWORK)
          On exit, if  INFO = 0,  DWORK(1) returns the optimal value
          of LDWORK;  DWORK(2)  contains the reciprocal condition
          number of the triangular factor of the QR factorization of
          the matrix  W2,  if  COMUSE = 'C',  or of the matrix
          Gamma,  if  COMUSE = 'U'  (see METHOD); if  JOBX0 = 'N'
          and  COMUSE <> 'C',  DWORK(2)  is set to one;
          if  COMUSE = 'C',  M > 0,  and  JOB = 'D',   DWORK(3)
          contains the reciprocal condition number of the triangular
          factor of the QR factorization of  U;  denoting
             g = 2,  if  JOBX0  = 'X'  and  COMUSE <> 'C'  or
                         COMUSE = 'C'  and  M = 0  or   JOB = 'B',
             g = 3,  if  COMUSE = 'C'  and  M > 0  and  JOB = 'D',
          then  DWORK(i), i = g+1:g+N*N,
                DWORK(j), j = g+1+N*N:g+N*N+L*N,  and
                DWORK(k), k = g+1+N*N+L*N:g+N*N+L*N+N*M,
          contain the transformed system matrices  At, Ct, and Bt,
          respectively, corresponding to the real Schur form of the
          given system state matrix  A,  i.e.,
             At = V'*A*V,  Bt = V'*B,  Ct = C*V.
          The matrices  At, Ct, Bt  are not computed if  JOBX0 = 'N'
          and  COMUSE <> 'C'.
          On exit, if  INFO = -26,  DWORK(1)  returns the minimum
          value of LDWORK.

  LDWORK  INTEGER
          The length of the array DWORK.
          LDWORK >= 2,  if  JOBX0 = 'N'  and  COMUSE <> 'C',  or
                        if  max( N, M ) = 0.
          Otherwise,
          LDWORK >= LDW1 + N*( N + M + L ) +
                           max( 5*N, LDW1, min( LDW2, LDW3 ) ),
          where, if  COMUSE = 'C',  then
          LDW1 = 2,          if  M = 0  or   JOB = 'B',
          LDW1 = 3,          if  M > 0  and  JOB = 'D',
          LDWa = t*L*(r + 1) + max( N + max( d, f ), 6*r ),
          LDW2 = LDWa,       if  M = 0  or  JOB = 'B',
          LDW2 = max( LDWa, t*L*(r + 1) + 2*M*M + 6*M ),
                             if  M > 0  and JOB = 'D',
          LDWb = (b + r)*(r + 1) +
                  max( q*(r + 1) + N*N*M + c + max( d, f ), 6*r ),
          LDW3 = LDWb,       if  M = 0  or  JOB = 'B',
          LDW3 = max( LDWb, (b + r)*(r + 1) + 2*M*M + 6*M ),
                             if  M > 0  and JOB = 'D',
             r = N*M + a,
             a = 0,                  if  JOBX0 = 'N',
             a = N,                  if  JOBX0 = 'X';
             b = 0,                  if  JOB   = 'B',
             b = L*M,                if  JOB   = 'D';
             c = 0,                  if  JOBX0 = 'N',
             c = L*N,                if  JOBX0 = 'X';
             d = 0,                  if  JOBX0 = 'N',
             d = 2*N*N + N,          if  JOBX0 = 'X';
             f = 2*r,                if  JOB   = 'B'   or  M = 0,
             f = M + max( 2*r, M ),  if  JOB   = 'D'  and  M > 0;
             q = b + r*L;
          and, if  JOBX0 = 'X'  and  COMUSE <> 'C',  then
          LDW1 = 2,
          LDW2 = t*L*(N + 1) + 2*N + max( 2*N*N, 4*N ),
          LDW3 = N*(N + 1) + 2*N + max( q*(N + 1) + 2*N*N + L*N,
                                        4*N ),
             q = N*L.
          For good performance,  LDWORK  should be larger.
          If  LDWORK >= LDW2,  or if  COMUSE = 'C'  and
              LDWORK >= t*L*(r + 1) + (b + r)*(r + 1) + N*N*M + c +
                        max( d, f ),
          then standard QR factorizations of the matrices  U  and/or
          W2,  if  COMUSE = 'C',  or of the matrix  Gamma,  if
          JOBX0 = 'X'  and  COMUSE <> 'C'  (see METHOD), are used.
          Otherwise, the QR factorizations are computed sequentially
          by performing  NCYCLE  cycles, each cycle (except possibly
          the last one) processing  s < t  samples, where  s  is
          chosen by equating  LDWORK  to the first term of  LDWb,
          if  COMUSE = 'C',  or of  LDW3,  if  COMUSE <> 'C',  for
          q  replaced by  s*L.  (s  is larger than or equal to the
          minimum value of  NSMP.)  The computational effort may
          increase and the accuracy may slightly decrease with the
          decrease of  s.  Recommended value is  LDWORK = LDW2,
          assuming a large enough cache size, to also accommodate
          A,  (B,)  C,  (D,)  U,  and  Y.

Warning Indicator
  IWARN   INTEGER
          = 0:  no warning;
          = 4:  the least squares problem to be solved has a
                rank-deficient coefficient matrix;
          = 6:  the matrix  A  is unstable;  the estimated  x(0)
                and/or  B and D  could be inaccurate.
          NOTE: the value 4 of  IWARN  has no significance for the
                identification problem.

Error Indicator
  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value;
          = 1:  if the QR algorithm failed to compute all the
                eigenvalues of the matrix A (see LAPACK Library
                routine DGEES); the locations  DWORK(i),  for
                i = g+1:g+N*N,  contain the partially converged
                Schur form;
          = 2:  the singular value decomposition (SVD) algorithm did
                not converge.

Method
  Matrix  A  is initially reduced to a real Schur form, A = V*At*V',
  and the given system matrices are transformed accordingly. For the
  reduced system, an extension and refinement of the method in [1,2]
  is used. Specifically, for  JOBX0 = 'X',  COMUSE = 'C',  and
  JOB = 'D',  denoting

        X = [ vec(D')' vec(B)' x0' ]',

  where  vec(M)  is the vector obtained by stacking the columns of
  the matrix  M,  then  X  is the least squares solution of the
  system  S*X = vec(Y),  with the matrix  S = [ diag(U)  W ],
  defined by

        ( U         |     | ... |     |     | ... |     |         )
        (   U       |  11 | ... |  n1 |  12 | ... |  nm |         )
    S = (     :     | y   | ... | y   | y   | ... | y   | P*Gamma ),
        (       :   |     | ... |     |     | ... |     |         )
        (         U |     | ... |     |     | ... |     |         )
                                                                  ij
  diag(U)  having  L  block rows and columns.  In this formula,  y
  are the outputs of the system for zero initial state computed
  using the following model, for j = 1:m, and for i = 1:n,
         ij          ij                    ij
        x  (k+1) = Ax  (k) + e_i u_j(k),  x  (0) = 0,

         ij          ij
        y  (k)   = Cx  (k),

  where  e_i  is the i-th n-dimensional unit vector,  Gamma  is
  given by

             (     C     )
             (    C*A    )
     Gamma = (   C*A^2   ),
             (     :     )
             ( C*A^(t-1) )

  and  P  is a permutation matrix that groups together the rows of
  Gamma  depending on the same row of  C,  namely
  [ c_j;  c_j*A;  c_j*A^2; ...  c_j*A^(t-1) ],  for j = 1:L.
  The first block column,  diag(U),  is not explicitly constructed,
  but its structure is exploited. The last block column is evaluated
  using powers of A with exponents 2^k. No interchanges are applied.
  A special QR decomposition of the matrix  S  is computed. Let
  U = q*[ r' 0 ]'  be the QR decomposition of  U,  if  M > 0,  where
  r  is  M-by-M.   Then,  diag(q')  is applied to  W  and  vec(Y).
  The block-rows of  S  and  vec(Y)  are implicitly permuted so that
  matrix  S  becomes

     ( diag(r)  W1 )
     (    0     W2 ),

  where  W1  has L*M rows. Then, the QR decomposition of  W2 is
  computed (sequentially, if  M > 0) and used to obtain  B  and  x0.
  The intermediate results and the QR decomposition of  U  are
  needed to find  D.  If a triangular factor is too ill conditioned,
  then singular value decomposition (SVD) is employed. SVD is not
  generally needed if the input sequence is sufficiently
  persistently exciting and  NSMP  is large enough.
  If the matrix  W  cannot be stored in the workspace (i.e.,
  LDWORK < LDW2),  the QR decompositions of  W2  and  U  are
  computed sequentially.
  For  JOBX0 = 'N'  and  COMUSE = 'C',  or  JOB = 'B',  a simpler
  problem is solved efficiently.

  For  JOBX0 = 'X'  and  COMUSE <> 'C',  a simpler method is used.
  Specifically, the output y0(k) of the system for zero initial
  state is computed for k = 0, 1, ...,  t-1 using the given model.
  Then the following least squares problem is solved for x(0)

                      (   y(0) - y0(0)   )
                      (   y(1) - y0(1)   )
     Gamma * x(0)  =  (        :         ).
                      (        :         )
                      ( y(t-1) - y0(t-1) )

  The coefficient matrix  Gamma  is evaluated using powers of A with
  exponents 2^k. The QR decomposition of this matrix is computed.
  If its triangular factor  R  is too ill conditioned, then singular
  value decomposition of  R  is used.
  If the coefficient matrix cannot be stored in the workspace (i.e.,
  LDWORK < LDW2),  the QR decomposition is computed sequentially.

References
  [1] Verhaegen M., and Varga, A.
      Some Experience with the MOESP Class of Subspace Model
      Identification Methods in Identifying the BO105 Helicopter.
      Report TR R165-94, DLR Oberpfaffenhofen, 1994.

  [2] Sima, V., and Varga, A.
      RASP-IDENT : Subspace Model Identification Programs.
      Deutsche Forschungsanstalt fur Luft- und Raumfahrt e. V.,
      Report TR R888-94, DLR Oberpfaffenhofen, Oct. 1994.

Numerical Aspects
  The implemented method is numerically stable.

Further Comments
  The algorithm for computing the system matrices  B  and  D  is
  less efficient than the MOESP or N4SID algorithms implemented in
  SLICOT Library routines IB01BD/IB01PD, because a large least
  squares problem has to be solved, but the accuracy is better, as
  the computed matrices  B  and  D  are fitted to the input and
  output trajectories. However, if matrix  A  is unstable, the
  computed matrices  B  and  D  could be inaccurate.

Example

Program Text

*     IB01CD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          LDA, LDB, LDC, LDD, LDR, LDU, LDV, LDW1, LDW2,
     $                 LDW4, LDW5, LDWORK, LDY, LIWORK, LMAX, MMAX,
     $                 NMAX, NOBRMX, NSMPMX
      PARAMETER        ( LMAX = 5, MMAX = 5, NOBRMX = 20, NSMPMX = 2000,
     $                   NMAX = NOBRMX - 1, LDA = NMAX, LDB = NMAX,
     $                   LDC  = LMAX, LDD = LMAX, LDV = NMAX,
     $                   LDR  = MAX( 2*( MMAX + LMAX )*NOBRMX,
     $                               3*MMAX*NOBRMX ), LDU = NSMPMX,
     $                   LDW1 = MAX( LMAX*( NOBRMX - 1 )*NMAX + NMAX +
     $                               MAX( 6*MMAX, 4*LMAX )*NOBRMX,
     $                               LMAX*NOBRMX*NMAX +
     $                               MAX( LMAX*( NOBRMX - 1 )*NMAX +
     $                                    3*NMAX + LMAX +
     $                                    ( 2*MMAX + LMAX )*NOBRMX,
     $                                    2*LMAX*( NOBRMX - 1 )*NMAX +
     $                                    NMAX*NMAX + 8*NMAX,
     $                                    NMAX +
     $                                    4*( MMAX*NOBRMX + NMAX ) ) ),
     $                   LDW2 = LMAX*NOBRMX*NMAX +
     $                          MMAX*NOBRMX*( NMAX + LMAX )*
     $                          ( MMAX*( NMAX + LMAX ) + 1 ) +
     $                          MAX( ( NMAX + LMAX )**2,
     $                          4*MMAX*( NMAX + LMAX ) + 1 ),
     $                   LDW4 = NSMPMX*LMAX*NMAX*( MMAX + 1 ) +
     $                          MAX( NMAX +
     $                               MAX( 2*NMAX*NMAX + NMAX,
     $                                    MMAX +
     $                                    MAX( 2*NMAX*( MMAX + 1 ),
     $                                         MMAX ),
     $                                    6*NMAX*( MMAX + 1 ) ),
     $                               2*MMAX*MMAX*NMAX + 6*MMAX ),
     $                   LDW5 = ( LMAX*MMAX + NMAX*( MMAX + 1 ) )*
     $                          NMAX*( MMAX + 1 ) +
     $                          MAX( ( LMAX*MMAX +
     $                               LMAX*NMAX*( MMAX + 1 ) )*
     $                               NMAX*( MMAX + 1 ) +
     $                               NMAX*NMAX*MMAX + LMAX*NMAX +
     $                               MAX( 2*NMAX*NMAX + NMAX,
     $                                    MMAX +
     $                                    MAX( 2*NMAX*( MMAX + 1 ),
     $                                         MMAX ),
     $                                    6*NMAX*( MMAX + 1 ) ),
     $                               2*MMAX*MMAX*NMAX + 6*MMAX ),
     $                   LDWORK = MAX( 6*( MMAX + LMAX )*NOBRMX,
     $                                 ( MMAX + LMAX )*( 4*NOBRMX*
     $                                 ( MMAX + LMAX + 2 ) - 2 ),
     $                                 ( MMAX + LMAX )*4*NOBRMX*
     $                                 ( NOBRMX + 1 ), LDW1, LDW2,
     $                                 3 + ( NMAX + MMAX + LMAX )*NMAX +
     $                                 MAX( 5*NMAX, 3,
     $                                      MIN( LDW4, LDW5 ) ) ),
     $                   LDY = NSMPMX,
     $                   LIWORK = MAX( ( MMAX + LMAX )*NOBRMX,
     $                                 MMAX*NOBRMX + NMAX,
     $                                 MMAX*( NMAX + LMAX ),
     $                                 NMAX*MMAX + NMAX, MMAX )
     $                 )
*     .. Local Scalars ..
      LOGICAL          NGIVEN
      CHARACTER        ALG, BATCH, COMUSE, CONCT, CTRL, JOB, JOBBD,
     $                 JOBCK, JOBD, JOBDA, JOBX0, METH, METHA
      INTEGER          I, ICYCLE, II, INFO, IWARN, J, L, M, N, NCYCLE,
     $                 NGIV, NOBR, NSAMPL, NSMP
      DOUBLE PRECISION RCOND, TOL
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA, NMAX), B(LDB, MMAX), C(LDC, NMAX),
     $                 D(LDD, MMAX), DUM(1), DWORK(LDWORK),
     $                 R(LDR, 2*(MMAX+LMAX)*NOBRMX),
     $                 SV(LMAX*NOBRMX), U(LDU, MMAX), V(LDV, NMAX),
     $                 X0(NMAX), Y(LDY, LMAX)
      INTEGER          IWORK(LIWORK)
      LOGICAL          BWORK(1)
*     .. External Functions ..
      LOGICAL          LSAME
      EXTERNAL         LSAME
*     .. External Subroutines ..
      EXTERNAL         IB01AD, IB01BD, IB01CD
*     .. Intrinsic Functions ..
      INTRINSIC        MAX, MIN
*     .. Executable Statements ..
*
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
*     If the value of N is positive, it will be taken as system order.
      READ ( NIN, FMT = '()' )
      READ ( NIN, FMT = * ) NOBR, N, M, L, NSMP, RCOND, TOL
      READ ( NIN, FMT = * ) METH, ALG, JOBD, BATCH, CONCT, CTRL, JOB,
     $                      COMUSE, JOBX0
      IF ( LSAME( BATCH, 'F' ) ) THEN
         READ ( NIN, FMT = * ) NCYCLE
      ELSE
         NCYCLE = 1
      END IF
      NSAMPL = NCYCLE*NSMP
*
      NGIVEN = N.GT.0
      IF( NGIVEN )
     $   NGIV = N
      IF ( NOBR.LE.0 .OR. NOBR.GT.NOBRMX ) THEN
         WRITE ( NOUT, FMT = 99997 ) NOBR
      ELSE IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
         WRITE ( NOUT, FMT = 99996 ) M
      ELSE IF ( L.LE.0 .OR. L.GT.LMAX ) THEN
         WRITE ( NOUT, FMT = 99995 ) L
      ELSE IF ( NSMP.LT.0 .OR. NSMP.GT.NSMPMX .OR.
     $        ( NSMP.LT.2*( M + L + 1 )*NOBR - 1 .AND.
     $          LSAME( BATCH, 'O' ) ) .OR.
     $        ( NSAMPL.LT.2*( M + L + 1 )*NOBR - 1 .AND.
     $          LSAME( BATCH, 'L' ) ) .OR.
     $          NSMP.LT.2*NOBR .AND. ( LSAME( BATCH, 'F' ) .OR.
     $                                 LSAME( BATCH, 'I' ) ) ) THEN
         WRITE ( NOUT, FMT = 99994 ) NSMP
      ELSE IF ( NCYCLE.LE.0 .OR. NSAMPL.GT.NSMPMX ) THEN
         WRITE ( NOUT, FMT = 99993 ) NCYCLE
      ELSE IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
         WRITE ( NOUT, FMT = 99983 ) N
      ELSE
*        Read the matrices U and Y from the input file.
         IF ( M.GT.0 )
     $      READ ( NIN, FMT = * )
     $                         ( ( U(I,J), J = 1, M ), I = 1, NSAMPL )
         READ ( NIN, FMT = * ) ( ( Y(I,J), J = 1, L ), I = 1, NSAMPL )
*        Force some options, depending on the specifications.
         IF ( LSAME( METH, 'C' ) ) THEN
            METHA = 'M'
            JOBDA = 'N'
         ELSE
            METHA = METH
            JOBDA = JOBD
         END IF
*        The covariances and Kalman gain matrix are not computed.
         JOBCK = 'N'
         IF ( LSAME( JOB, 'A' ) .OR. LSAME( JOB, 'C' ) ) THEN
            JOBBD = 'D'
         ELSE
            JOBBD = JOB
         END IF
         IF ( LSAME( COMUSE, 'C' ) ) THEN
            JOB = 'C'
         ELSE IF ( LSAME( COMUSE, 'U' ) ) THEN
            JOB = 'A'
         END IF
*        Compute the  R  factor from a QR (or Cholesky) factorization
*        of the Hankel-like matrix (or correlation matrix).
         DO 10 ICYCLE = 1, NCYCLE
            II = ( ICYCLE - 1 )*NSMP + 1
            IF ( NCYCLE.GT.1 ) THEN
               IF ( ICYCLE.GT.1 )      BATCH = 'I'
               IF ( ICYCLE.EQ.NCYCLE ) BATCH = 'L'
            END IF
            CALL IB01AD( METHA, ALG, JOBDA, BATCH, CONCT, CTRL, NOBR, M,
     $                   L, NSMP, U(II,1), LDU, Y(II,1), LDY, N, R, LDR,
     $                   SV, RCOND, TOL, IWORK, DWORK, LDWORK, IWARN,
     $                   INFO )
   10    CONTINUE
         IF ( INFO.NE.0 ) THEN
            WRITE ( NOUT, FMT = 99998 ) INFO
         ELSE
            IF ( IWARN.NE.0 )
     $         WRITE ( NOUT, FMT = 99990 ) IWARN
            IF( NGIVEN )
     $         N = NGIV
*           Compute the system matrices and x0.
            CALL IB01BD( METH, JOB, JOBCK, NOBR, N, M, L, NSMP, R,
     $                   LDR, A, LDA, C, LDC, B, LDB, D, LDD, DUM, 1,
     $                   DUM, 1, DUM, 1, DUM, 1, RCOND, IWORK, DWORK,
     $                   LDWORK, BWORK, IWARN, INFO )
            IF ( INFO.NE.0 ) THEN
               WRITE ( NOUT, FMT = 99982 ) INFO
            ELSE
               IF ( IWARN.NE.0 )
     $            WRITE ( NOUT, FMT = 99981 ) IWARN
               CALL IB01CD( JOBX0, COMUSE, JOBBD, N, M, L, NSMP, A, LDA,
     $                      B, LDB, C, LDC, D, LDD, U, LDU, Y, LDY, X0,
     $                      V, LDV, RCOND, IWORK, DWORK, LDWORK, IWARN,
     $                      INFO )
               IF ( INFO.NE.0 ) THEN
                  WRITE ( NOUT, FMT = 99992 ) INFO
               ELSE
                  IF ( IWARN.NE.0 )
     $               WRITE ( NOUT, FMT = 99991 ) IWARN
                  IF ( LSAME( JOB, 'A' ) .OR. LSAME( JOB, 'C' ) ) THEN
                     WRITE ( NOUT, FMT = 99989 )
                     DO 20 I = 1, N
                        WRITE ( NOUT, FMT = 99988 ) ( A(I,J), J = 1,N )
   20                CONTINUE
                     WRITE ( NOUT, FMT = 99987 )
                     DO 30 I = 1, L
                        WRITE ( NOUT, FMT = 99988 ) ( C(I,J), J = 1,N )
   30                CONTINUE
                  END IF
                  IF ( LSAME( COMUSE, 'C' ) ) THEN
                     WRITE ( NOUT, FMT = 99986 )
                     DO 40 I = 1, N
                        WRITE ( NOUT, FMT = 99988 ) ( B(I,J), J = 1,M )
   40                CONTINUE
                     IF ( LSAME( JOBBD, 'D' ) ) THEN
                        WRITE ( NOUT, FMT = 99985 )
                        DO 50 I = 1, L
                           WRITE ( NOUT, FMT = 99988 )
     $                           ( D(I,J), J = 1,M )
   50                   CONTINUE
                     END IF
                  END IF
                  IF ( LSAME( JOBX0, 'X' ) ) THEN
                     WRITE ( NOUT, FMT = 99984 )
                     WRITE ( NOUT, FMT = 99988 ) ( X0(I), I = 1,N )
                  END IF
               END IF
            END IF
         END IF
      END IF
      STOP
99999 FORMAT ( ' IB01CD EXAMPLE PROGRAM RESULTS', /1X)
99998 FORMAT ( ' INFO on exit from IB01AD = ',I2)
99997 FORMAT (/' NOBR is out of range.',/' NOBR = ', I5)
99996 FORMAT (/' M is out of range.',/' M = ', I5)
99995 FORMAT (/' L is out of range.',/' L = ', I5)
99994 FORMAT (/' NSMP is out of range.',/' NSMP = ', I5)
99993 FORMAT (/' NCYCLE is out of range.',/' NCYCLE = ', I5)
99992 FORMAT ( ' INFO on exit from IB01CD = ',I2)
99991 FORMAT ( ' IWARN on exit from IB01CD = ',I2)
99990 FORMAT ( ' IWARN on exit from IB01AD = ',I2)
99989 FORMAT (/' The system state matrix A is ')
99988 FORMAT (20(1X,F8.4))
99987 FORMAT (/' The system output matrix C is ')
99986 FORMAT (/' The system input matrix B is ')
99985 FORMAT (/' The system input-output matrix D is ')
99984 FORMAT (/' The initial state vector x0 is ')
99983 FORMAT (/' N is out of range.',/' N = ', I5)
99982 FORMAT ( ' INFO on exit from IB01BD = ',I2)
99981 FORMAT ( ' IWARN on exit from IB01BD = ',I2)
      END
Program Data
 IB01CD EXAMPLE PROGRAM DATA
  15     0     1     1  1000    0.0   -1.0
   C     C     N     O     N     N     A     C     X
   6.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   6.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   3.41
   3.41
   3.41
   3.41
   6.41
   3.41
   3.41
   3.41
   6.41
   6.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   6.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   6.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   3.41
   6.41
   3.41
   3.41
   3.41
   6.41
   3.41
   3.41
   3.41
   6.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   6.41
   6.41
   3.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   3.41
   3.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   3.41
   3.41
   3.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   6.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   3.41
   6.41
   6.41
   6.41
   6.41
   6.41
   6.41
   4.766099
   4.763659
   4.839359
   5.002979
   5.017629
   5.056699
   5.154379
   5.361949
   5.425439
   5.569519
   5.681849
   5.742899
   5.803949
   5.918729
   5.821049
   5.447419
   5.061589
   4.629349
   4.267939
   4.011519
   3.850349
   3.711159
   3.569519
   3.518239
   3.652549
   3.818609
   3.862559
   4.011519
   4.353409
   4.705049
   5.083559
   5.344859
   5.274039
   5.127519
   4.761219
   4.451089
   4.221539
   4.045709
   3.874769
   3.730689
   3.662319
   3.576849
   3.542659
   3.479169
   3.454749
   3.359509
   3.298459
   3.225199
   3.200779
   3.225199
   3.227639
   3.274039
   3.457189
   3.867449
   4.321659
   4.492599
   4.431549
   4.243519
   4.050599
   3.857679
   3.730689
   3.791739
   3.921169
   3.955359
   3.847909
   3.725809
   3.611039
   3.716039
   4.092109
   4.480389
   4.814939
   5.054259
   5.303339
   5.486489
   5.672089
   5.779529
   5.799069
   5.664759
   5.291129
   4.880879
   4.558529
   4.184909
   3.889419
   3.708719
   3.623249
   3.569519
   3.718479
   4.033499
   4.412009
   4.629349
   4.558529
   4.394919
   4.180019
   4.197119
   4.431549
   4.714819
   4.961459
   5.300899
   5.567079
   5.681849
   5.545099
   5.188569
   4.883319
   4.600049
   4.270379
   4.038389
   3.838139
   3.711159
   3.591499
   3.535329
   3.486489
   3.476729
   3.425439
   3.381489
   3.369279
   3.364389
   3.347299
   3.381489
   3.420559
   3.413229
   3.452309
   3.635459
   4.038389
   4.375379
   4.727029
   5.056699
   5.298459
   5.532889
   5.466959
   5.195899
   4.885759
   4.763659
   4.875989
   5.042049
   5.283809
   5.491379
   5.596379
   5.672089
   5.772209
   5.830819
   5.933379
   5.899189
   5.935819
   5.894309
   5.918729
   5.994429
   5.957799
   6.031059
   6.062809
   6.040829
   6.096999
   6.123859
   6.162929
   6.040829
   5.845469
   5.772209
   5.799069
   5.923609
   5.928499
   6.001759
   6.001759
   6.060369
   5.882099
   5.510909
   5.322879
   5.371719
   5.454749
   5.437649
   5.159269
   4.902859
   4.587839
   4.502369
   4.595159
   4.824709
   5.064029
   5.271599
   5.466959
   5.615919
   5.528009
   5.254499
   4.883319
   4.517019
   4.197119
   4.001759
   3.806399
   3.904079
   3.923609
   3.869889
   3.806399
   3.720929
   3.818609
   4.140949
   4.529229
   4.805179
   5.086009
   5.339969
   5.532889
   5.576849
   5.667199
   5.791739
   5.850349
   5.923609
   5.921169
   5.977339
   5.740459
   5.388809
   5.000539
   4.849129
   4.944369
   5.173919
   5.369279
   5.447419
   5.603709
   5.730689
   5.850349
   5.979779
   5.991989
   6.084789
   5.940709
   5.803949
   5.791739
   5.603709
   5.264269
   4.946809
   4.619579
   4.514579
   4.433989
   4.285029
   4.121419
   3.945589
   3.984659
   4.219099
   4.546319
   4.873549
   5.154379
   5.388809
   5.613479
   5.835699
   5.884539
   5.955359
   5.762439
   5.459629
   5.061589
   4.707499
   4.458409
   4.267939
   4.053039
   3.943149
   3.825929
   3.967569
   4.280149
   4.480389
   4.492599
   4.390039
   4.197119
   4.111649
   3.982219
   3.867449
   3.767319
   3.872329
   4.236189
   4.663539
   4.971229
   5.066469
   4.902859
   4.675749
   4.392479
   4.099439
   4.114089
   4.326539
   4.643999
   4.971229
   5.159269
   5.388809
   5.576849
   5.652549
   5.803949
   5.913839
   5.886979
   5.799069
   5.730689
   5.762439
   5.813719
   5.821049
   5.928499
   6.013969
   5.764879
   5.413229
   5.098219
   4.678189
   4.372939
   4.392479
   4.590279
   4.919949
   5.017629
   4.858899
   4.675749
   4.619579
   4.834479
   5.090889
   5.376599
   5.681849
   5.823489
   5.952919
   6.062809
   6.089669
   6.075019
   6.026179
   5.994429
   6.077459
   5.857679
   5.701389
   5.730689
   5.784419
   5.823489
   5.894309
   5.762439
   5.415679
   4.961459
   4.595159
   4.331429
   4.297239
   4.582949
   4.861339
   5.173919
   5.166589
   4.919949
   4.607369
   4.370499
   4.182469
   4.038389
   4.145839
   4.431549
   4.556089
   4.480389
   4.375379
   4.370499
   4.558529
   4.858899
   4.895529
   4.741679
   4.744129
   4.875989
   5.105539
   5.239849
   5.518239
   5.652549
   5.723369
   5.855239
   5.962679
   5.984659
   5.984659
   6.055479
   6.062809
   6.055479
   6.070129
   5.784419
   5.440099
   5.056699
   4.941929
   5.010299
   5.134849
   5.313109
   5.479169
   5.623249
   5.562199
   5.330209
   5.010299
   4.665979
   4.414459
   4.201999
   4.048159
   4.079899
   4.189789
   4.131179
   4.004199
   3.916289
   3.960239
   4.199559
   4.624469
   4.883319
   5.137289
   5.379049
   5.623249
   5.762439
   5.833259
   5.686739
   5.366839
   5.225199
   5.239849
   5.354629
   5.508469
   5.596379
   5.752669
   5.874769
   5.906519
   5.894309
   5.742899
   5.447419
   5.024959
   4.883319
   4.885759
   4.893089
   4.714819
   4.451089
   4.233749
   4.043269
   3.864999
   3.757559
   3.669639
   3.593939
   3.547539
   3.506029
   3.454749
   3.398579
   3.361949
   3.339969
   3.374159
   3.520679
   3.713599
   3.757559
   3.779529
   3.696509
   3.777089
   3.886979
   3.904079
   3.850349
   3.965129
   4.282589
   4.521899
   4.714819
   4.971229
   5.220319
   5.532889
   5.652549
   5.781979
   5.955359
   6.035939
   6.118969
   6.133629
   6.153159
   6.192229
   6.143389
   6.167809
   5.991989
   5.652549
   5.459629
   5.437649
   5.339969
   5.098219
   4.785639
   4.492599
   4.236189
   4.067689
   3.933379
   3.823489
   3.730689
   3.611039
   3.564639
   3.549989
   3.557309
   3.513359
   3.515799
   3.694059
   4.072579
   4.480389
   4.705049
   4.612259
   4.385149
   4.201999
   4.026179
   3.904079
   3.774649
   3.691619
   3.845469
   4.201999
   4.585399
   4.902859
   5.256949
   5.510909
   5.640339
   5.843029
   5.974889
   5.935819
   5.821049
   5.528009
   5.171479
   4.810059
   4.453529
   4.380269
   4.565859
   4.805179
   5.125079
   5.354629
   5.589059
   5.764879
   5.923609
   5.940709
   5.857679
   5.694059
   5.486489
   5.149499
   4.844249
   4.541439
   4.267939
   4.060369
   3.960239
   3.789299
   3.642779
   3.525569
   3.498699
   3.454749
   3.408349
   3.379049
   3.376599
   3.361949
   3.359509
   3.369279
   3.398579
   3.579289
   3.948029
   4.412009
   4.585399
   4.514579
   4.343639
   4.155599
   3.984659
   4.043269
   4.307009
   4.421779
   4.353409
   4.223979
   4.053039
   3.940709
   3.838139
   3.730689
   3.652549
   3.611039
   3.564639
   3.496259
   3.462069
   3.454749
   3.425439
   3.379049
   3.432769
   3.623249
   3.974889
   4.380269
   4.714819
   5.073799
   5.369279
   5.603709
   5.745349
   5.652549
   5.401019
   5.015189
   4.709939
   4.416899
   4.236189
   4.236189
   4.248399
   4.221539
   4.297239
   4.590279
   4.893089
   5.134849
   5.427889
   5.379049
   5.364389
   5.452309
   5.567079
   5.672089
   5.769769
   5.830819
   5.923609
   5.965129
   6.057919
   6.050599
   6.072579
   6.111649
   6.070129
   5.896749
   5.755109
   5.718479
   5.821049
   6.001759
   6.001759
   5.901629
   5.557309
   5.173919
   4.800289
   4.431549
   4.194679
   4.006639
   3.850349
   3.747789
   3.642779
   3.591499
   3.569519
   3.528009
   3.537779
   3.554869
   3.493819
   3.447419
   3.440099
   3.408349
   3.410789
   3.452309
   3.681849
   4.060369
   4.441319
   4.854019
   5.154379
   5.425439
   5.596379
   5.586619
   5.354629
   5.027399
   4.863779
   4.761219
   4.570739
   4.368059
   4.397359
   4.573189
   4.841809
   5.203219
   5.452309
   5.652549
   5.855239
   5.906519
   5.952919
   5.828369
   5.791739
   5.799069
   5.813719
   5.877209
   5.955359
   5.781979
   5.518239
   5.127519
   4.763659
   4.492599
   4.233749
   4.011519
   3.855239
   3.691619
   3.635459
   3.818609
   4.155599
   4.590279
   4.988329
   5.076239
   4.907739
   4.648889
   4.377829
   4.216649
   4.287469
   4.590279
   4.846689
   5.139729
   5.388809
   5.689179
   5.884539
   6.043269
   6.170259
   6.211769
   6.250839
   6.209329
   6.013969
   5.701389
   5.469399
   5.479169
   5.557309
   5.728249
   5.882099
   5.984659
   5.901629
   5.581729
   5.371719
   5.418119
   5.510909
   5.667199
   5.791739
   5.698949
   5.484049
   5.154379
   4.980999
   5.061589
   5.195899
   5.359509
   5.615919
   5.762439
   5.857679
   5.948029
   5.835699
   5.706269
   5.498699
   5.188569
   5.117749
   5.191009
   5.315549
   5.532889
   5.444979
   5.396139
   5.274039
   5.027399
   4.744129
   4.668419
   4.651329
   4.514579
   4.267939
   4.260609
   4.263049
   4.189789
   4.277699
   4.600049
   4.932159
   5.283809
   5.528009
   5.740459
   5.874769
   5.955359
   5.991989
   5.845469
   5.528009
   5.061589
   4.734359
   4.534109
   4.534109
   4.697729
   4.744129
   4.619579
   4.643999
   4.832039
   5.132399
   5.410789
   5.625689
   5.603709
   5.315549
   4.961459
   4.619579
   4.358289
   4.155599
   4.033499
   3.886979
   3.772209
   3.640339
   3.532889
   3.435209
   3.427889
   3.422999
   3.398579
   3.603709
   4.023729
   4.451089
   4.792969
   4.902859
   4.780759
   4.590279
   4.336309
   4.145839
   4.216649
   4.433989
   4.714819
   5.098219
   5.359509
   5.569519
   5.772209
   5.921169
   6.055479
   5.962679
   5.642779
   5.435209
   5.388809
   5.537779
   5.681849
   5.701389
   5.615919
   5.667199
   5.740459
   5.803949
   5.882099
   5.950469
   6.072579
   6.148279
   6.116529
   6.177579
   6.201999
   6.206889
   5.991989
   5.564639
   5.178799
   4.998089
   5.051819
   5.232529
   5.484049
   5.686739
   5.899189
   5.869889
   5.977339
   6.053039
   6.079899
   6.128739
   6.079899
   6.167809
   6.194679
   6.236189
   6.053039
   5.652549
   5.274039
   4.858899
   4.534109
   4.455969
   4.619579
   4.866229
   5.117749
   5.166589
   5.056699
   5.002979
   5.098219
   5.325319
   5.567079
   5.466959
   5.252059
   4.946809
   4.880879
   4.980999
   5.225199
   5.459629
   5.723369
   5.791739
   5.906519
   5.991989
   5.835699
   5.528009
   5.142169
   4.775869
   4.490159
   4.236189
   4.023729
   3.886979
   3.752669
   3.681849
   3.806399
   4.145839
   4.600049
   5.002979
   5.303339
   5.552429
   5.615919
   5.523119
   5.611039
   5.713599
   5.845469
   5.899189
   5.994429
   6.092109
   6.092109
   6.143389
   6.153159
   6.233749
   6.187349
   6.013969
   5.835699
   5.774649
   5.686739
   5.537779
   5.327759
   5.054259
   4.700169
   4.394919
   4.180019
   4.043269
   3.877209
   3.752669
   3.728249
   3.869889
   4.206889
   4.355849
   4.426669
   4.453529
   4.521899
   4.392479
   4.155599
   3.965129
   3.877209
   3.970009
   4.258169
   4.421779
   4.336309
   4.299679
   4.392479
   4.675749
   4.761219
   4.658659
   4.490159
   4.307009
   4.126299
   3.972449
   4.077459
   4.372939
   4.741679
   5.088449
   5.186129
   5.037169
   4.785639
   4.563419
   4.534109
   4.705049
   4.741679
   4.648889
   4.431549
   4.238629
   4.065249
   3.943149
   3.811279
   3.691619
   3.652549
   3.825929
   4.223979
   4.424219
   4.429109
   4.319219
   4.138509
   3.965129
   3.886979
   3.801509
   3.701389
   3.640339
   3.767319
   4.150719
   4.648889
   4.990769
   5.088449
   5.022509
   4.783199
   4.685519
   4.665979
   4.707499
   4.912619
   5.195899
   5.415679
   5.623249
   5.740459
   5.899189
   5.928499
   6.050599
   6.153159
   5.965129
   5.586619
   5.381489
   5.371719
   5.486489
   5.567079
   5.821049
   5.913839
   5.994429
   6.011519
   5.999309
   6.018849
   5.821049
   5.728249
   5.740459
   5.764879
   5.882099
   5.926049
   5.750229
   5.415679
   4.995649
   4.861339
   4.902859
   5.103099
   5.364389
   5.596379
   5.752669
   5.845469
   5.928499
   6.006639
   5.840579
   5.518239
   5.173919
   4.739239
   4.458409
   4.426669
   4.602489
   4.822269
   5.183689
   5.430329
   5.652549
   5.821049
   5.706269
   5.369279
   5.027399
   4.705049
   4.414459
   4.145839
   3.965129
   4.033499
   4.372939
   4.683079
Program Results
 IB01CD EXAMPLE PROGRAM RESULTS


 The system state matrix A is 
   0.8924   0.3887   0.1285   0.1716
  -0.0837   0.6186  -0.6273  -0.4582
   0.0052   0.1307   0.6685  -0.6755
   0.0055   0.0734  -0.2148   0.4788

 The system output matrix C is 
  -0.4442   0.6663   0.3961   0.4102

 The system input matrix B is 
  -0.2150
  -0.1962
   0.0511
   0.0373

 The system input-output matrix D is 
  -0.0018

 The initial state vector x0 is 
 -11.4329  -0.6767   0.0472   0.3600

Return to index