path: root/modules/cacsd/src/slicot/ib01bd.f

      SUBROUTINE IB01BD( METH, JOB, JOBCK, NOBR, N, M, L, NSMPL, R,
     $                   LDR, A, LDA, C, LDC, B, LDB, D, LDD, Q, LDQ, 
     $                   RY, LDRY, S, LDS, K, LDK, TOL, IWORK, DWORK, 
     $                   LDWORK, BWORK, IWARN, INFO )
C
C     RELEASE 4.0, WGS COPYRIGHT 2000.
C
C     PURPOSE
C
C     To estimate the system matrices A, C, B, and D, the noise
C     covariance matrices Q, Ry, and S, and the Kalman gain matrix K
C     of a linear time-invariant state space model, using the 
C     processed triangular factor R of the concatenated block Hankel
C     matrices, provided by SLICOT Library routine IB01AD.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     METH    CHARACTER*1
C             Specifies the subspace identification method to be used,
C             as follows:
C             = 'M':  MOESP  algorithm with past inputs and outputs;
C             = 'N':  N4SID  algorithm;
C             = 'C':  combined method:  MOESP  algorithm for finding the
C                     matrices A and C, and  N4SID  algorithm for 
C                     finding the matrices B and D.
C
C     JOB     CHARACTER*1
C             Specifies which matrices should be computed, as follows:
C             = 'A':  compute all system matrices, A, B, C, and D;
C             = 'C':  compute the matrices A and C only;
C             = 'B':  compute the matrix B only;
C             = 'D':  compute the matrices B and D only.
C
C     JOBCK   CHARACTER*1
C             Specifies whether or not the covariance matrices and the
C             Kalman gain matrix are to be computed, as follows:
C             = 'C':  the covariance matrices only should be computed;
C             = 'K':  the covariance matrices and the Kalman gain 
C                     matrix should be computed;
C             = 'N':  the covariance matrices and the Kalman gain matrix
C                     should not be computed.
C
C     Input/Output Parameters
C
C     NOBR    (input) INTEGER
C             The number of block rows,  s,  in the input and output
C             Hankel matrices processed by other routines.  NOBR > 1.
C
C     N       (input) INTEGER
C             The order of the system.  NOBR > N > 0.
C
C     M       (input) INTEGER
C             The number of system inputs.  M >= 0.
C
C     L       (input) INTEGER
C             The number of system outputs.  L > 0.
C
C     NSMPL   (input) INTEGER
C             If  JOBCK = 'C' or 'K',  the total number of samples used
C             for calculating the covariance matrices.
C             NSMPL >= 2*(M+L)*NOBR.
C             This parameter is not meaningful if  JOBCK = 'N'.
C
C     R       (input/workspace) DOUBLE PRECISION array, dimension
C             ( LDR,2*(M+L)*NOBR )
C             On entry, the leading  2*(M+L)*NOBR-by-2*(M+L)*NOBR  part
C             of this array must contain the relevant data for the MOESP
C             or N4SID algorithms, as constructed by SLICOT Library
C             routine IB01AD. Let  R_ij,  i,j = 1:4,  be the
C             ij submatrix of  R  (denoted  S  in IB01AD),  partitioned 
C             by  M*NOBR,  L*NOBR,  M*NOBR,  and  L*NOBR  rows and 
C             columns. The submatrix  R_22  contains the matrix of left
C             singular vectors used. Also needed, for  METH = 'N'  or
C             JOBCK <> 'N',  are the submatrices  R_11,  R_14 : R_44,
C             and, for  METH = 'M' or 'C'  and  JOB <> 'C', the 
C             submatrices  R_31  and  R_12,  containing the processed
C             matrices  R_1c  and  R_2c,  respectively, as returned by
C             SLICOT Library routine IB01AD.
C             Moreover, if  METH = 'N'  and  JOB = 'A' or 'C',  the 
C             block-row  R_41 : R_43  must contain the transpose of the
C             block-column  R_14 : R_34  as returned by SLICOT Library
C             routine IB01AD. 
C             The remaining part of  R  is used as workspace.
C             On exit, part of this array is overwritten. Specifically,
C             if  METH = 'M',  R_22  and  R_31  are overwritten if
C                 JOB = 'B' or 'D',  and  R_12,  R_22,  R_14 : R_34,  
C                 and possibly  R_11  are overwritten if  JOBCK <> 'N';
C             if  METH = 'N',  all needed submatrices are overwritten.
C             The details of the contents of  R  need not be known if
C             this routine is called once just after calling the SLICOT
C             Library routine IB01AD.
C
C     LDR     INTEGER
C             The leading dimension of the array  R.
C             LDR >= 2*(M+L)*NOBR.
C
C     A       (input or output) DOUBLE PRECISION array, dimension 
C             (LDA,N)
C             On entry, if  METH = 'N' or 'C'  and  JOB = 'B' or 'D', 
C             the leading N-by-N part of this array must contain the
C             system state matrix.
C             If  METH = 'M'  or  (METH = 'N' or 'C'  and JOB = 'A' 
C             or 'C'),  this array need not be set on input.
C             On exit, if  JOB = 'A' or 'C'  and  INFO = 0,  the
C             leading N-by-N part of this array contains the system 
C             state matrix.
C
C     LDA     INTEGER
C             The leading dimension of the array A.  
C             LDA >= N,  if  JOB = 'A' or 'C',  or  METH = 'N' or 'C'  
C                            and  JOB = 'B' or 'D';
C             LDA >= 1,  otherwise.
C
C     C       (input or output) DOUBLE PRECISION array, dimension 
C             (LDC,N)
C             On entry, if  METH = 'N' or 'C'  and  JOB = 'B' or 'D',
C             the leading L-by-N part of this array must contain the
C             system output matrix.
C             If  METH = 'M'  or  (METH = 'N' or 'C'  and JOB = 'A'
C             or 'C'),  this array need not be set on input.
C             On exit, if  JOB = 'A' or 'C'  and  INFO = 0,  or 
C             INFO = 3  (or  INFO >= 0,  for  METH = 'M'),  the leading 
C             L-by-N part of this array contains the system output
C             matrix.
C
C     LDC     INTEGER
C             The leading dimension of the array C.
C             LDC >= L,  if  JOB = 'A' or 'C',  or  METH = 'N' or 'C'  
C                            and  JOB = 'B' or 'D';
C             LDC >= 1,  otherwise.
C
C     B       (output) DOUBLE PRECISION array, dimension (LDB,M)
C             If  M > 0,  JOB = 'A', 'B', or 'D'  and  INFO = 0,  the 
C             leading N-by-M part of this array contains the system 
C             input matrix. If  M = 0  or  JOB = 'C',  this array is 
C             not referenced.
C
C     LDB     INTEGER
C             The leading dimension of the array B.
C             LDB >= N,  if M > 0 and JOB = 'A', 'B', or 'D';
C             LDB >= 1,  if M = 0 or  JOB = 'C'.
C
C     D       (output) DOUBLE PRECISION array, dimension (LDD,M)
C             If  M > 0,  JOB = 'A' or 'D'  and  INFO = 0,  the leading
C             L-by-M part of this array contains the system input-output
C             matrix. If  M = 0  or  JOB = 'C' or 'B',  this array is
C             not referenced.  
C
C     LDD     INTEGER
C             The leading dimension of the array D.
C             LDD >= L,  if M > 0 and JOB = 'A' or 'D';
C             LDD >= 1,  if M = 0 or  JOB = 'C' or 'B'.
C
C     Q       (output) DOUBLE PRECISION array, dimension (LDQ,N)
C             If  JOBCK = 'C' or 'K',  the leading N-by-N part of this
C             array contains the positive semidefinite state covariance
C             matrix. If  JOBCK = 'K',  this matrix has been used as
C             state weighting matrix for computing the Kalman gain.
C             This parameter is not referenced if JOBCK = 'N'.
C
C     LDQ     INTEGER
C             The leading dimension of the array Q.
C             LDQ >= N,  if JOBCK = 'C' or 'K';
C             LDQ >= 1,  if JOBCK = 'N'.
C
C     RY      (output) DOUBLE PRECISION array, dimension (LDRY,L)
C             If  JOBCK = 'C' or 'K',  the leading L-by-L part of this
C             array contains the positive (semi)definite output
C             covariance matrix. If  JOBCK = 'K',  this matrix has been
C             used as output weighting matrix for computing the Kalman
C             gain.
C             This parameter is not referenced if JOBCK = 'N'.
C
C     LDRY    INTEGER
C             The leading dimension of the array RY.
C             LDRY >= L,  if JOBCK = 'C' or 'K';
C             LDRY >= 1,  if JOBCK = 'N'.
C
C     S       (output) DOUBLE PRECISION array, dimension (LDS,L)
C             If  JOBCK = 'C' or 'K',  the leading N-by-L part of this
C             array contains the state-output cross-covariance matrix.
C             If  JOBCK = 'K',  this matrix has been used as state-
C             output weighting matrix for computing the Kalman gain.
C             This parameter is not referenced if JOBCK = 'N'.
C
C     LDS     INTEGER
C             The leading dimension of the array S.
C             LDS >= N,  if JOBCK = 'C' or 'K';
C             LDS >= 1,  if JOBCK = 'N'.
C
C     K       (output) DOUBLE PRECISION array, dimension ( LDK,L ) 
C             If  JOBCK = 'K',  the leading  N-by-L  part of this array
C             contains the estimated Kalman gain matrix.
C             If  JOBCK = 'C' or 'N',  this array is not referenced.
C
C     LDK     INTEGER
C             The leading dimension of the array  K.
C             LDK >= N,  if JOBCK = 'K';
C             LDK >= 1,  if JOBCK = 'C' or 'N'.
C
C     Tolerances
C
C     TOL     DOUBLE PRECISION
C             The tolerance to be used for estimating the rank of
C             matrices. If the user sets  TOL > 0,  then the given value
C             of  TOL  is used as a lower bound for the reciprocal
C             condition number;  an m-by-n matrix whose estimated
C             condition number is less than  1/TOL  is considered to  
C             be of full rank.  If the user sets  TOL <= 0,  then an 
C             implicitly computed, default tolerance, defined by 
C             TOLDEF = m*n*EPS,  is used instead, where  EPS  is the 
C             relative machine precision (see LAPACK Library routine 
C             DLAMCH).
C
C     Workspace
C
C     IWORK   INTEGER array, dimension (LIWORK)
C             LIWORK >= max(LIW1,LIW2), where
C             LIW1 = N,                     if METH <> 'N' and M = 0
C                                        or JOB = 'C' and JOBCK = 'N';
C             LIW1 = M*NOBR+N,              if METH <> 'N', JOB = 'C', 
C                                           and JOBCK <> 'N';
C             LIW1 = max(L*NOBR,M*NOBR),    if METH = 'M', JOB <> 'C', 
C                                           and JOBCK = 'N';
C             LIW1 = max(L*NOBR,M*NOBR+N),  if METH = 'M', JOB <> 'C',
C                                           and JOBCK = 'C' or 'K';
C             LIW1 = max(M*NOBR+N,M*(N+L)), if METH = 'N', or METH = 'C'
C                                           and JOB  <> 'C';
C             LIW2 = 0,                     if JOBCK <> 'K';
C             LIW2 = N*N,                   if JOBCK =  'K'.
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C             On exit, if  INFO = 0,  DWORK(1) returns the optimal value
C             of LDWORK,  and  DWORK(2),  DWORK(3),  DWORK(4),  and  
C             DWORK(5)  contain the reciprocal condition numbers of the
C             triangular factors of the following matrices (defined in
C             SLICOT Library routine IB01PD and in the lower level
C             routines):
C                GaL  (GaL = Un(1:(s-1)*L,1:n)),
C                R_1c (if  METH = 'M' or 'C'),
C                M    (if  JOBCK = 'C' or 'K'  or  METH = 'N'),  and  
C                Q or T  (see SLICOT Library routine IB01PY or IB01PX),
C             respectively.
C             If  METH = 'N',  DWORK(3)  is set to one without any
C             calculations. Similarly, if  METH = 'M'  and  JOBCK = 'N',
C             DWORK(4)  is set to one. If  M = 0  or  JOB = 'C',
C             DWORK(3)  and  DWORK(5)  are set to one.
C             If  JOBCK = 'K'  and  INFO = 0,  DWORK(6)  to  DWORK(13)
C             contain information about the accuracy of the results when
C             computing the Kalman gain matrix, as follows:
C                DWORK(6)  - reciprocal condition number of the matrix
C                            U11  of the Nth order system of algebraic
C                            equations from which the solution matrix  X
C                            of the Riccati equation is obtained;
C                DWORK(7)  - reciprocal pivot growth factor for the LU
C                            factorization of the matrix  U11;
C                DWORK(8)  - reciprocal condition number of the matrix
C                            As = A - S*inv(Ry)*C,  which is inverted by
C                            the standard Riccati solver;
C                DWORK(9)  - reciprocal pivot growth factor for the LU
C                            factorization of the matrix  As;
C                DWORK(10) - reciprocal condition number of the matrix
C                            Ry;
C                DWORK(11) - reciprocal condition number of the matrix
C                            Ry + C*X*C';
C                DWORK(12) - reciprocal condition number for the Riccati
C                            equation solution;
C                DWORK(13) - forward error bound for the Riccati
C                            equation solution.
C             On exit, if  INFO = -30,  DWORK(1)  returns the minimum
C             value of LDWORK.
C
C     LDWORK  INTEGER
C             The length of the array DWORK.
C             LDWORK >= max( LDW1,LDW2,LDW3 ), where, if METH = 'M',  
C             LDW1 >= max( 2*(L*NOBR-L)*N+2*N, (L*NOBR-L)*N+N*N+7*N ), 
C                     if JOB = 'C' or JOB = 'A' and M = 0;
C             LDW1 >= max( 2*(L*NOBR-L)*N+N*N+7*N,
C                          (L*NOBR-L)*N+N+6*M*NOBR, (L*NOBR-L)*N+N+ 
C                          max( L+M*NOBR, L*NOBR + max( 3*L*NOBR, M ))),
C                     if M > 0 and JOB = 'A', 'B', or 'D'; 
C             LDW2 >= 0,                          if JOBCK = 'N';
C             LDW2 >= L*NOBR*N+
C                     max( (L*NOBR-L)*N+Aw+2*N+max(5*N,(2*M+L)*NOBR+L),
C                          4*(M*NOBR+N), M*NOBR+2*N+L ),
C                                                 if JOBCK = 'C' or 'K',
C             where Aw = N+N*N, if M = 0 or JOB = 'C';
C                   Aw = 0,     otherwise;
C             if METH = 'N',  
C             LDW1 >= L*NOBR*N+max( (L*NOBR-L)*N+2*N+(2*M+L)*NOBR+L,
C                                   2*(L*NOBR-L)*N+N*N+8*N,
C                                   N+4*(M*NOBR+N), M*NOBR+3*N+L );
C             LDW2 >= 0, if M = 0 or JOB = 'C';
C             LDW2 >= L*NOBR*N+M*NOBR*(N+L)*(M*(N+L)+1)+
C                                max( (N+L)**2, 4*M*(N+L)+1 ),
C                     if M > 0 and JOB = 'A', 'B', or 'D'; 
C             and, if METH = 'C', LDW1 as 
C             max( LDW1 for METH = 'M', JOB = 'C', LDW1 for METH = 'N'),
C             and LDW2 for METH = 'N' are used;
C             LDW3 >= 0,                     if JOBCK <> 'K';
C             LDW3 >= max(  4*N*N+2*N*L+L*L+max( 3*L,N*L ),
C                          14*N*N+12*N+5 ),  if JOBCK =  'K'.
C             For good performance,  LDWORK  should be larger.
C
C     BWORK   LOGICAL array, dimension (LBWORK)
C             LBWORK = 2*N, if JOBCK =  'K';
C             LBWORK = 0,   if JOBCK <> 'K'.
C
C     Warning Indicator
C
C     IWARN   INTEGER
C             = 0:  no warning;
C             = 4:  a least squares problem to be solved has a
C                   rank-deficient coefficient matrix;
C             = 5:  the computed covariance matrices are too small.
C                   The problem seems to be a deterministic one; the
C                   gain matrix is set to zero.
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value;
C             = 2:  the singular value decomposition (SVD) algorithm did
C                   not converge;
C             = 3:  a singular upper triangular matrix was found;
C             = 3+i:  if  JOBCK = 'K'  and the associated Riccati
C                   equation could not be solved, where i = 1,...,6;
C                   (see the description of the parameter INFO for the
C                   SLICOT Library routine SB02RD for the meaning of
C                   the i values);
C             = 10: the QR algorithm did not converge.
C
C     METHOD
C
C     In the MOESP approach, the matrices  A  and  C  are first
C     computed from an estimated extended observability matrix [1],
C     and then, the matrices  B  and  D  are obtained by solving an
C     extended linear system in a least squares sense.
C     In the N4SID approach, besides the estimated extended 
C     observability matrix, the solutions of two least squares problems
C     are used to build another least squares problem, whose solution
C     is needed to compute the system matrices  A,  C,  B,  and  D.  The
C     solutions of the two least squares problems are also optionally
C     used by both approaches to find the covariance matrices.
C     The Kalman gain matrix is obtained by solving a discrete-time 
C     algebraic Riccati equation.
C
C     REFERENCES
C
C     [1] Verhaegen M., and Dewilde, P.
C         Subspace Model Identification. Part 1: The output-error 
C         state-space model identification class of algorithms.
C         Int. J. Control, 56, pp. 1187-1210, 1992.
C
C     [2] Van Overschee, P., and De Moor, B.
C         N4SID: Two Subspace Algorithms for the Identification
C         of Combined Deterministic-Stochastic Systems.
C         Automatica, Vol.30, No.1, pp. 75-93, 1994.
C
C     [3] Van Overschee, P.
C         Subspace Identification : Theory - Implementation -
C         Applications.
C         Ph. D. Thesis, Department of Electrical Engineering,
C         Katholieke Universiteit Leuven, Belgium, Feb. 1995.
C
C     [4] Sima, V.
C         Subspace-based Algorithms for Multivariable System 
C         Identification.
C         Studies in Informatics and Control, 5, pp. 335-344, 1996.
C
C     NUMERICAL ASPECTS
C
C     The implemented method consists in numerically stable steps.
C
C     FURTHER COMMENTS
C
C     The covariance matrices are computed using the N4SID approach.
C     Therefore, for efficiency reasons, it is advisable to set 
C     METH = 'N',  if the Kalman gain matrix or covariance matrices 
C     are needed  (JOBCK = 'K', or 'C').  When  JOBCK = 'N',  it could 
C     be more efficient to use the combined method,  METH = 'C'. 
C     Often, this combination will also provide better accuracy than
C     MOESP algorithm.
C     In some applications, it is useful to compute the system matrices
C     using two calls to this routine, the first one with  JOB = 'C', 
C     and the second one with  JOB = 'B' or 'D'.  This is slightly less
C     efficient than using a single call with  JOB = 'A',  because some
C     calculations are repeated. If  METH = 'N',  all the calculations 
C     at the first call are performed again at the second call; 
C     moreover, it is required to save the needed submatrices of  R 
C     before the first call and restore them before the second call.
C     If the covariance matrices and/or the Kalman gain are desired,  
C     JOBCK  should be set to  'C'  or  'K'  at the second call.
C     If  B  and  D  are both needed, they should be computed at once.
C
C     CONTRIBUTOR
C
C     V. Sima, Research Institute for Informatics, Bucharest, Dec. 1999.
C
C     REVISIONS
C
C     March 2000, August 2000.
C
C     KEYWORDS
C
C     Identification methods; least squares solutions; multivariable
C     systems; QR decomposition; singular value decomposition.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO  = 0.0D0, ONE = 1.0D0 )
C     .. Scalar Arguments ..
      DOUBLE PRECISION   TOL
      INTEGER            INFO, IWARN, L, LDA, LDB, LDC, LDD, LDK, LDQ,
     $                   LDR, LDRY, LDS, LDWORK, M, N, NOBR, NSMPL
      CHARACTER          JOB, JOBCK, METH
C     .. Array Arguments ..
      DOUBLE PRECISION   A(LDA, *), B(LDB, *), C(LDC, *), D(LDD, *),
     $                   DWORK(*),  K(LDK, *), Q(LDQ, *), R(LDR, *), 
     $                   RY(LDRY, *), S(LDS, *)
      INTEGER            IWORK( * )
      LOGICAL            BWORK( * )
C     .. Local Scalars ..
      DOUBLE PRECISION   FERR, RCOND, RCONDR, RNORM, SEP
      INTEGER            I, IA, IAW, IC, ID, IERR, IFACT, IG, IK, IO, 
     $                   IQ, IR, IS, IT, IV, IWARNL, IWI, IWR, IX,
     $                   JWORK, LDUNN, LL, LMMNOL, LMNOBR, LNOBR,
     $                   MAXWRK, MINWRK, MNOBR, MNOBRN, N2, NL, NN, NPL,
     $                   NR
      CHARACTER          JOBBD, JOBCOV, JOBCV
      LOGICAL            COMBIN, MOESP, N4SID, WITHAL, WITHB, WITHC,
     $                   WITHCO, WITHD, WITHK
C     .. Local Arrays ..
      DOUBLE PRECISION   RCND(8)
      INTEGER            OUFACT(2)
C     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
C     .. External Subroutines ..
      EXTERNAL           DLACPY, DLASET, IB01PD, MA02AD, SB02MT, SB02ND,
     $                   SB02RD, XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC          INT, MAX
C     .. Executable Statements ..
C
C     Decode the scalar input parameters.
C
      MOESP  = LSAME( METH,  'M' )
      N4SID  = LSAME( METH,  'N' )
      COMBIN = LSAME( METH,  'C' )
      WITHAL = LSAME( JOB,   'A' )
      WITHC  = LSAME( JOB,   'C' ) .OR. WITHAL
      WITHD  = LSAME( JOB,   'D' ) .OR. WITHAL
      WITHB  = LSAME( JOB,   'B' ) .OR. WITHD
      WITHK  = LSAME( JOBCK, 'K' )
      WITHCO = LSAME( JOBCK, 'C' ) .OR. WITHK
      MNOBR  = M*NOBR
      LNOBR  = L*NOBR
      LMNOBR = LNOBR + MNOBR
      MNOBRN = MNOBR + N
      LDUNN  = ( LNOBR - L )*N
      LMMNOL = LNOBR + 2*MNOBR + L
      NR     = LMNOBR + LMNOBR
      NPL    = N + L
      N2     = N + N
      NN     = N*N
      NL     = N*L
      LL     = L*L
      MINWRK = 1
      IWARN  = 0
      INFO   = 0
C
C     Check the scalar input parameters.
C
      IF( .NOT.( MOESP .OR. N4SID .OR. COMBIN ) ) THEN
         INFO = -1
      ELSE IF( .NOT.( WITHB .OR. WITHC ) ) THEN
         INFO = -2
      ELSE IF( .NOT.( WITHCO .OR. LSAME( JOBCK, 'N' ) ) ) THEN
         INFO = -3
      ELSE IF( NOBR.LE.1 ) THEN
         INFO = -4
      ELSE IF( N.LE.0 .OR. N.GE.NOBR ) THEN
         INFO = -5
      ELSE IF( M.LT.0 ) THEN
         INFO = -6
      ELSE IF( L.LE.0 ) THEN
         INFO = -7
      ELSE IF( WITHCO .AND. NSMPL.LT.NR ) THEN
         INFO = -8
      ELSE IF( LDR.LT.NR ) THEN
         INFO = -10
      ELSE IF( LDA.LT.1 .OR. ( ( WITHC .OR. ( WITHB .AND. .NOT.MOESP ) )
     $   .AND. LDA.LT.N ) ) THEN
         INFO = -12
      ELSE IF( LDC.LT.1 .OR. ( ( WITHC .OR. ( WITHB .AND. .NOT.MOESP ) )
     $   .AND. LDC.LT.L ) ) THEN
         INFO = -14
      ELSE IF( LDB.LT.1  .OR. ( WITHB  .AND. LDB.LT.N .AND. M.GT.0 ) )
     $      THEN
         INFO = -16
      ELSE IF( LDD.LT.1  .OR. ( WITHD  .AND. LDD.LT.L .AND. M.GT.0 ) )
     $      THEN
         INFO = -18
      ELSE IF( LDQ.LT.1  .OR. ( WITHCO .AND. LDQ.LT.N ) )  THEN
         INFO = -20
      ELSE IF( LDRY.LT.1 .OR. ( WITHCO .AND. LDRY.LT.L ) ) THEN
         INFO = -22
      ELSE IF( LDS.LT.1  .OR. ( WITHCO .AND. LDS.LT.N ) )  THEN
         INFO = -24
      ELSE IF( LDK.LT.1  .OR. ( WITHK  .AND. LDK.LT.N ) )  THEN
         INFO = -26
      ELSE IF( LDWORK.GE.1 ) THEN
C
C        Compute workspace.
C        (Note: Comments in the code beginning "Workspace:" describe the
C         minimal amount of workspace needed at that point in the code,
C         as well as the preferred amount for good performance.)
C
         IAW    = 0
         MINWRK = LDUNN + 4*N
         IF( .NOT.N4SID ) THEN 
            ID = 0
            IF( WITHC ) THEN 
               MINWRK = MAX( MINWRK, 2*LDUNN + N2, LDUNN + NN + 7*N )
            END IF
         ELSE
            ID = N
         END IF
C
         IF( ( M.GT.0 .AND. WITHB ) .OR. .NOT.MOESP ) THEN 
            MINWRK = MAX( MINWRK, 2*LDUNN + NN + ID + 7*N )
            IF ( MOESP ) 
     $         MINWRK = MAX( MINWRK, LDUNN + N + 6*MNOBR, LDUNN + N + 
     $                       MAX( L + MNOBR, LNOBR + MAX( 3*LNOBR, M ) )
     $                     )
         ELSE 
            IF( .NOT.N4SID ) 
     $         IAW = N + NN
         END IF
C
         IF( .NOT.MOESP .OR. WITHCO ) THEN 
            MINWRK = MAX( MINWRK, LDUNN + IAW + N2 + MAX( 5*N, LMMNOL ),
     $                    ID + 4*MNOBRN, ID + MNOBRN + NPL ) 
            IF( .NOT.MOESP .AND. M.GT.0 .AND. WITHB ) 
     $         MINWRK = MAX( MINWRK, MNOBR*NPL*( M*NPL + 1 ) + 
     $                       MAX( NPL**2, 4*M*NPL + 1 ) )
            MINWRK = LNOBR*N + MINWRK
         END IF
C
         IF( WITHK ) THEN 
            MINWRK = MAX( MINWRK, 4*NN + 2*NL + LL + MAX( 3*L, NL ),
     $                    14*NN + 12*N + 5 ) 
         END IF
C          
         IF ( LDWORK.LT.MINWRK ) THEN
            INFO = -30
            DWORK( 1 ) = MINWRK
         END IF
      END IF
C
C     Return if there are illegal arguments.
C
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'IB01BD', -INFO )
         RETURN
      END IF
C
      IF ( .NOT.WITHK ) THEN
         JOBCV = JOBCK
      ELSE
         JOBCV = 'C'
      END IF
C
      IO = 1
      IF ( .NOT.MOESP .OR. WITHCO ) THEN
         JWORK = IO + LNOBR*N
      ELSE
         JWORK = IO
      END IF
      MAXWRK = MINWRK
C
C     Call the computational routine for estimating system matrices.
C
      IF ( .NOT.COMBIN ) THEN
         CALL IB01PD( METH, JOB, JOBCV, NOBR, N, M, L, NSMPL, R, LDR,
     $                A, LDA, C, LDC, B, LDB, D, LDD, Q, LDQ, RY, LDRY,
     $                S, LDS, DWORK(IO), LNOBR, TOL, IWORK,
     $                DWORK(JWORK), LDWORK-JWORK+1, IWARN, INFO )
C
      ELSE
C
         IF ( WITHC ) THEN
            IF ( WITHAL ) THEN
               JOBCOV = 'N'
            ELSE
               JOBCOV = JOBCV
            END IF
            CALL IB01PD( 'MOESP', 'C and A', JOBCOV, NOBR, N, M, L,
     $                   NSMPL, R, LDR, A, LDA, C, LDC, B, LDB, D, LDD, 
     $                   Q, LDQ, RY, LDRY, S, LDS, DWORK(IO), LNOBR,
     $                   TOL, IWORK, DWORK(JWORK), LDWORK-JWORK+1, 
     $                   IWARNL, INFO )
            IF ( INFO.NE.0 ) 
     $         RETURN 
            IWARN  = MAX( IWARN, IWARNL )  
            MAXWRK = MAX( MAXWRK, INT( DWORK(JWORK) ) + JWORK - 1 )
         END IF
C
         IF ( WITHB ) THEN
            IF ( .NOT.WITHAL ) THEN
               JOBBD = JOB
            ELSE
               JOBBD = 'D'
            END IF
            CALL IB01PD( 'N4SID', JOBBD, JOBCV, NOBR, N, M, L, NSMPL, R,
     $                LDR, A, LDA, C, LDC, B, LDB, D, LDD, Q, LDQ, 
     $                RY, LDRY, S, LDS, DWORK(IO), LNOBR, TOL, IWORK,
     $                DWORK(JWORK), LDWORK-JWORK+1, IWARNL, INFO )
            IWARN  = MAX( IWARN, IWARNL )  
         END IF
      END IF
C
      IF ( INFO.NE.0 ) 
     $   RETURN 
      MAXWRK = MAX( MAXWRK, INT( DWORK(JWORK) ) + JWORK - 1 )
C
      DO 10 I = 1, 4
         RCND(I) = DWORK(JWORK+I)
   10 CONTINUE
C
      IF ( WITHK ) THEN
         IF ( IWARN.EQ.5 ) THEN
C
C           The problem seems to be a deterministic one. Set the Kalman
C           gain to zero, set accuracy parameters and return.
C
            CALL DLASET( 'Full', N, L, ZERO, ZERO, K, LDK )
C
            DO 20 I = 6, 12
               DWORK(I) = ONE
   20       CONTINUE
C
            DWORK(13) = ZERO
         ELSE
C
C           Compute the Kalman gain matrix.
C
C           Convert the optimal problem with coupling weighting terms 
C           to a standard problem.
C           Workspace:  need   4*N*N+2*N*L+L*L+max( 3*L,N*L );
C                       prefer larger.
C
            IX    = 1
            IQ    = IX + NN
            IA    = IQ + NN
            IG    = IA + NN
            IC    = IG + NN
            IR    = IC + NL
            IS    = IR + LL
            JWORK = IS + NL
C
            CALL MA02AD( 'Full',  N, N, A,  LDA,  DWORK(IA), N ) 
            CALL MA02AD( 'Full',  L, N, C,  LDC,  DWORK(IC), N ) 
            CALL DLACPY( 'Upper', N, N, Q,  LDQ,  DWORK(IQ), N ) 
            CALL DLACPY( 'Upper', L, L, RY, LDRY, DWORK(IR), L ) 
            CALL DLACPY( 'Full',  N, L, S,  LDS,  DWORK(IS), N ) 
C
            CALL SB02MT( 'G needed', 'Nonzero S', 'Not factored',
     $                   'Upper', N, L, DWORK(IA), N, DWORK(IC), N, 
     $                   DWORK(IQ), N, DWORK(IR), L, DWORK(IS), N,
     $                   IWORK, IFACT, DWORK(IG), N, IWORK(L+1), 
     $                   DWORK(JWORK), LDWORK-JWORK+1, IERR )
            IF ( IERR.NE.0 ) THEN
               INFO = 3 
               RETURN 
            END IF
            MAXWRK = MAX( MAXWRK, INT( DWORK(JWORK) ) + JWORK - 1 )
            RCONDR = DWORK(JWORK+1)
C
C           Solve the Riccati equation. 
C           Workspace:  need   14*N*N+12*N+5;
C                       prefer larger.
C
            IT    = IC
            IV    = IT  + NN
            IWR   = IV  + NN
            IWI   = IWR + N2
            IS    = IWI + N2
            JWORK = IS  + N2*N2
C
            CALL SB02RD( 'All', 'Discrete', 'Direct', 'NoTranspose',
     $                   'Upper', 'General scaling', 'Unstable first', 
     $                   'Not factored', 'Reduced', N, DWORK(IA), N,
     $                   DWORK(IT), N, DWORK(IV), N, DWORK(IG), N,
     $                   DWORK(IQ), N, DWORK(IX), N, SEP, RCOND, FERR, 
     $                   DWORK(IWR), DWORK(IWI), DWORK(IS), N2, IWORK, 
     $                   DWORK(JWORK), LDWORK-JWORK+1, BWORK, IERR )
C
            IF ( IERR.NE.0 .AND. IERR.LT.7 ) THEN
               INFO = IERR + 3 
               RETURN 
            END IF
            MAXWRK = MAX( MAXWRK, INT( DWORK(JWORK) ) + JWORK - 1 )
C
            DO 30 I = 1, 4
               RCND(I+4) = DWORK(JWORK+I)
   30       CONTINUE
C
C           Compute the gain matrix. 
C           Workspace:  need   2*N*N+2*N*L+L*L+3*L;
C                       prefer larger.
C
            IA    = IX + NN
            IC    = IA + NN
            IR    = IC + NL
            IK    = IR + LL
            JWORK = IK + NL
C
            CALL MA02AD( 'Full',  N, N, A,  LDA,  DWORK(IA), N ) 
            CALL MA02AD( 'Full',  L, N, C,  LDC,  DWORK(IC), N ) 
            CALL DLACPY( 'Upper', L, L, RY, LDRY, DWORK(IR), L ) 
C
            CALL SB02ND( 'Discrete', 'NotFactored', 'Upper',
     $                   'Nonzero S', N, L, 0, DWORK(IA), N, DWORK(IC),
     $                   N, DWORK(IR), L, IWORK, S, LDS, DWORK(IX), N,
     $                   RNORM, DWORK(IK), L, OUFACT, IWORK(L+1),
     $                   DWORK(JWORK), LDWORK-JWORK+1, IERR )   
C
            IF ( IERR.NE.0 ) THEN
               IF ( IERR.LE.L+1 ) THEN
                  INFO = 3 
               ELSE IF ( IERR.EQ.L+2 ) THEN
                  INFO = 10 
               END IF
               RETURN 
            END IF
            MAXWRK = MAX( MAXWRK, INT( DWORK(JWORK) ) + JWORK - 1 )
C
            CALL MA02AD( 'Full', L, N, DWORK(IK), L, K, LDK ) 
C
C           Set the accuracy parameters.
C
            DWORK(11) = DWORK(JWORK+1)
C
            DO 40 I = 6, 9
               DWORK(I) = RCND(I-1)
   40       CONTINUE
C
            DWORK(10) = RCONDR
            DWORK(12) = RCOND
            DWORK(13) = FERR
         END IF
      END IF
C
C     Return optimal workspace in  DWORK(1)  and the remaining
C     reciprocal condition numbers in the next locations.
C
      DWORK(1) = MAXWRK
C
      DO 50 I = 2, 5
         DWORK(I) = RCND(I-1)
   50 CONTINUE
C
      RETURN
C
C *** Last line of IB01BD ***
      END