<HTML>
<HEAD><TITLE>IB01CD - SLICOT Library Routine Documentation</TITLE>
</HEAD>
<BODY>
<H2><A Name="IB01CD">IB01CD</A></H2>
<H3>
Estimating the initial state and system matrices B and D using A, B, and input-output data (driver)
</H3>
<A HREF ="#Specification"><B>[Specification]</B></A>
<A HREF ="#Arguments"><B>[Arguments]</B></A>
<A HREF ="#Method"><B>[Method]</B></A>
<A HREF ="#References"><B>[References]</B></A>
<A HREF ="#Comments"><B>[Comments]</B></A>
<A HREF ="#Example"><B>[Example]</B></A>
<P>
<B><FONT SIZE="+1">Purpose</FONT></B>
<PRE>
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.
</PRE>
<A name="Specification"><B><FONT SIZE="+1">Specification</FONT></B></A>
<PRE>
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(*)
</PRE>
<A name="Arguments"><B><FONT SIZE="+1">Arguments</FONT></B></A>
<P>
<B>Mode Parameters</B>
<PRE>
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).
</PRE>
<B>Input/Output Parameters</B>
<PRE>
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'.
</PRE>
<B>Tolerances</B>
<PRE>
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.
</PRE>
<B>Workspace</B>
<PRE>
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.
</PRE>
<B>Warning Indicator</B>
<PRE>
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.
</PRE>
<B>Error Indicator</B>
<PRE>
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.
</PRE>
<A name="Method"><B><FONT SIZE="+1">Method</FONT></B></A>
<PRE>
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.
</PRE>
<A name="References"><B><FONT SIZE="+1">References</FONT></B></A>
<PRE>
[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.
</PRE>
<A name="Numerical Aspects"><B><FONT SIZE="+1">Numerical Aspects</FONT></B></A>
<PRE>
The implemented method is numerically stable.
</PRE>
<A name="Comments"><B><FONT SIZE="+1">Further Comments</FONT></B></A>
<PRE>
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.
</PRE>
<A name="Example"><B><FONT SIZE="+1">Example</FONT></B></A>
<P>
<B>Program Text</B>
<PRE>
* IB01CD EXAMPLE PROGRAM TEXT
*
* .. 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
</PRE>
<B>Program Data</B>
<PRE>
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
</PRE>
<B>Program Results</B>
<PRE>
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
</PRE>
<HR>
<p>
<A HREF=..\libindex.html><B>Return to index</B></A></BODY>
</HTML>