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.

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