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

      SUBROUTINE IB01AD( METH, ALG, JOBD, BATCH, CONCT, CTRL, NOBR, M,
     $                   L, NSMP, U, LDU, Y, LDY, N, R, LDR, SV, RCOND,
     $                   TOL, IWORK, DWORK, LDWORK, IWARN, INFO )
C
C     RELEASE 4.0, WGS COPYRIGHT 2000.
C
C     PURPOSE
C
C     To preprocess the input-output data for estimating the matrices 
C     of a linear time-invariant dynamical system and to find an 
C     estimate of the system order. The input-output data can,
C     optionally, be processed sequentially.
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     ALG     CHARACTER*1
C             Specifies the algorithm for computing the triangular
C             factor R, as follows:
C             = 'C':  Cholesky algorithm applied to the correlation
C                     matrix of the input-output data;
C             = 'F':  Fast QR algorithm;
C             = 'Q':  QR algorithm applied to the concatenated block
C                     Hankel matrices.
C
C     JOBD    CHARACTER*1
C             Specifies whether or not the matrices B and D should later
C             be computed using the MOESP approach, as follows:
C             = 'M':  the matrices B and D should later be computed
C                     using the MOESP approach;
C             = 'N':  the matrices B and D should not be computed using
C                     the MOESP approach.
C             This parameter is not relevant for METH = 'N'.
C
C     BATCH   CHARACTER*1
C             Specifies whether or not sequential data processing is to
C             be used, and, for sequential processing, whether or not
C             the current data block is the first block, an intermediate
C             block, or the last block, as follows:
C             = 'F':  the first block in sequential data processing;
C             = 'I':  an intermediate block in sequential data
C                     processing;
C             = 'L':  the last block in sequential data processing;
C             = 'O':  one block only (non-sequential data processing).
C             NOTE that when  100  cycles of sequential data processing
C                  are completed for  BATCH = 'I',  a warning is
C                  issued, to prevent for an infinite loop.
C
C     CONCT   CHARACTER*1
C             Specifies whether or not the successive data blocks in
C             sequential data processing belong to a single experiment,
C             as follows:
C             = 'C':  the current data block is a continuation of the
C                     previous data block and/or it will be continued
C                     by the next data block;
C             = 'N':  there is no connection between the current data
C                     block and the previous and/or the next ones.
C             This parameter is not used if BATCH = 'O'.
C
C     CTRL    CHARACTER*1
C             Specifies whether or not the user's confirmation of the
C             system order estimate is desired, as follows:
C             = 'C':  user's confirmation;
C             = 'N':  no confirmation.
C             If  CTRL = 'C',  a reverse communication routine,  IB01OY,
C             is indirectly called (by SLICOT Library routine IB01OD), 
C             and, after inspecting the singular values and system order
C             estimate,  n,  the user may accept  n  or set a new value.
C             IB01OY  is not called if CTRL = 'N'.
C
C     Input/Output Parameters
C
C     NOBR    (input) INTEGER
C             The number of block rows,  s,  in the input and output
C             block Hankel matrices to be processed.  NOBR > 0.
C             (In the MOESP theory,  NOBR  should be larger than  n,
C             the estimated dimension of state vector.)
C
C     M       (input) INTEGER
C             The number of system inputs.  M >= 0.
C             When M = 0, no system inputs are processed.
C
C     L       (input) INTEGER
C             The number of system outputs.  L > 0.
C
C     NSMP    (input) INTEGER
C             The number of rows of matrices  U  and  Y  (number of
C             samples,  t). (When sequential data processing is used,
C             NSMP  is the number of samples of the current data
C             block.)
C             NSMP >= 2*(M+L+1)*NOBR - 1,  for non-sequential
C                                          processing;
C             NSMP >= 2*NOBR,  for sequential processing.
C             The total number of samples when calling the routine with
C             BATCH = 'L'  should be at least  2*(M+L+1)*NOBR - 1.
C             The  NSMP  argument may vary from a cycle to another in
C             sequential data processing, but  NOBR, M,  and  L  should
C             be kept constant. For efficiency, it is advisable to use
C             NSMP  as large as possible.
C
C     U       (input) DOUBLE PRECISION array, dimension (LDU,M)
C             The leading NSMP-by-M part of this array must contain the
C             t-by-m input-data sequence matrix  U,
C             U = [u_1 u_2 ... u_m].  Column  j  of  U  contains the
C             NSMP  values of the j-th input component for consecutive
C             time increments.
C             If M = 0, this array is not referenced.
C
C     LDU     INTEGER
C             The leading dimension of the array U.
C             LDU >= NSMP, if M > 0;
C             LDU >= 1,    if M = 0.
C
C     Y       (input) DOUBLE PRECISION array, dimension (LDY,L)
C             The leading NSMP-by-L part of this array must contain the
C             t-by-l output-data sequence matrix  Y,
C             Y = [y_1 y_2 ... y_l].  Column  j  of  Y  contains the
C             NSMP  values of the j-th output component for consecutive
C             time increments.
C
C     LDY     INTEGER
C             The leading dimension of the array Y.  LDY >= NSMP.
C
C     N       (output) INTEGER
C             The estimated order of the system.
C             If  CTRL = 'C',  the estimated order has been reset to a
C             value specified by the user.
C
C     R       (output or input/output) DOUBLE PRECISION array, dimension
C             ( LDR,2*(M+L)*NOBR )
C             On exit, if ALG = 'C' and BATCH = 'F' or 'I', the leading
C             2*(M+L)*NOBR-by-2*(M+L)*NOBR upper triangular part of this
C             array contains the current upper triangular part of the
C             correlation matrix in sequential data processing.
C             If ALG = 'F' and BATCH = 'F' or 'I', the array R is not
C             referenced.
C             On exit, if INFO = 0, ALG = 'Q', and BATCH = 'F' or 'I',
C             the leading 2*(M+L)*NOBR-by-2*(M+L)*NOBR upper triangular
C             part of this array contains the current upper triangular
C             factor R from the QR factorization of the concatenated
C             block Hankel matrices. Denote  R_ij, i,j = 1:4,  the 
C             ij submatrix of  R,  partitioned by M*NOBR,  M*NOBR, 
C             L*NOBR,  and  L*NOBR  rows and columns. 
C             On exit, if INFO = 0 and BATCH = 'L' or 'O', the leading
C             2*(M+L)*NOBR-by-2*(M+L)*NOBR upper triangular part of
C             this array contains the matrix S, the processed upper
C             triangular factor R from the QR factorization of the
C             concatenated block Hankel matrices, as required by other
C             subroutines. Specifically, let  S_ij, i,j = 1:4,  be the
C             ij submatrix of  S,  partitioned by M*NOBR,  L*NOBR, 
C             M*NOBR,  and  L*NOBR  rows and columns. The submatrix 
C             S_22  contains the matrix of left singular vectors needed
C             subsequently. Useful information is stored in  S_11  and 
C             in the block-column  S_14 : S_44.  For METH = 'M' and 
C             JOBD = 'M', the upper triangular part of  S_31  contains
C             the upper triangular factor in the QR factorization of the
C             matrix  R_1c = [ R_12'  R_22'  R_11' ]',  and  S_12 
C             contains the corresponding leading part of the transformed
C             matrix  R_2c = [ R_13'  R_23'  R_14' ]'.  For  METH = 'N',
C             the subarray  S_41 : S_43  contains the transpose of the
C             matrix contained in  S_14 : S_34. 
C             The details of the contents of R need not be known if this
C             routine is followed by SLICOT Library routine IB01BD.
C             On entry, if ALG = 'C', or ALG = 'Q', and BATCH = 'I' or 
C             'L', the leading  2*(M+L)*NOBR-by-2*(M+L)*NOBR  upper
C             triangular part of this array must contain the upper
C             triangular matrix R computed at the previous call of this
C             routine in sequential data processing. The array R need
C             not be set on entry if ALG = 'F' or if BATCH = 'F' or 'O'.
C
C     LDR     INTEGER
C             The leading dimension of the array  R.
C             LDR >= MAX( 2*(M+L)*NOBR, 3*M*NOBR ),
C                                  for METH = 'M' and JOBD = 'M';
C             LDR >= 2*(M+L)*NOBR, for METH = 'M' and JOBD = 'N' or
C                                  for METH = 'N'.
C
C     SV      (output) DOUBLE PRECISION array, dimension ( L*NOBR )
C             The singular values used to estimate the system order.
C
C     Tolerances
C
C     RCOND   DOUBLE PRECISION
C             The tolerance to be used for estimating the rank of
C             matrices. If the user sets  RCOND > 0,  the given value
C             of  RCOND  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/RCOND  is considered to  
C             be of full rank.  If the user sets  RCOND <= 0,  then an 
C             implicitly computed, default tolerance, defined by 
C             RCONDEF = m*n*EPS,  is used instead, where  EPS  is the 
C             relative machine precision (see LAPACK Library routine 
C             DLAMCH).
C             This parameter is not used for  METH = 'M'.
C
C     TOL     DOUBLE PRECISION
C             Absolute tolerance used for determining an estimate of
C             the system order. If  TOL >= 0,  the estimate is
C             indicated by the index of the last singular value greater
C             than or equal to  TOL.  (Singular values less than  TOL
C             are considered as zero.) When  TOL = 0,  an internally
C             computed default value,  TOL = NOBR*EPS*SV(1),  is used,
C             where  SV(1)  is the maximal singular value, and  EPS  is
C             the relative machine precision (see LAPACK Library routine
C             DLAMCH). When  TOL < 0,  the estimate is indicated by the
C             index of the singular value that has the largest
C             logarithmic gap to its successor.
C
C     Workspace
C
C     IWORK   INTEGER array, dimension (LIWORK)
C             LIWORK >= (M+L)*NOBR, if METH = 'N';
C             LIWORK >= M+L, if METH = 'M' and ALG = 'F';
C             LIWORK >= 0,   if METH = 'M' and ALG = 'C' or 'Q'.
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C             On exit, if  INFO = 0,  DWORK(1) returns the optimal value
C             of LDWORK,  and, for  METH = 'N',  and  BATCH = 'L'  or
C             'O',  DWORK(2)  and  DWORK(3)  contain the reciprocal 
C             condition numbers of the triangular factors of the 
C             matrices  U_f  and  r_1  [6].
C             On exit, if  INFO = -23,  DWORK(1)  returns the minimum
C             value of LDWORK.
C             Let
C             k = 0,               if CONCT = 'N' and ALG = 'C' or 'Q';
C             k = 2*NOBR-1,        if CONCT = 'C' and ALG = 'C' or 'Q';
C             k = 2*NOBR*(M+L+1),  if CONCT = 'N' and ALG = 'F';
C             k = 2*NOBR*(M+L+2),  if CONCT = 'C' and ALG = 'F'.
C             The first (M+L)*k elements of  DWORK  should be preserved
C             during successive calls of the routine with  BATCH = 'F'
C             or  'I',  till the final call with  BATCH = 'L'.
C
C     LDWORK  INTEGER
C             The length of the array DWORK.
C             LDWORK >= (4*NOBR-2)*(M+L), if ALG = 'C', BATCH = 'F' or
C                             'I' and CONCT = 'C';
C             LDWORK >= 1, if ALG = 'C', BATCH = 'F' or 'I' and 
C                             CONCT = 'N';
C             LDWORK >= max((4*NOBR-2)*(M+L), 5*L*NOBR), if METH = 'M',
C                             ALG = 'C', BATCH = 'L' and CONCT = 'C';
C             LDWORK >= max((2*M-1)*NOBR, (M+L)*NOBR, 5*L*NOBR), 
C                             if METH = 'M', JOBD = 'M', ALG = 'C',
C                              BATCH = 'O', or 
C                             (BATCH = 'L' and CONCT = 'N');
C             LDWORK >= 5*L*NOBR, if METH = 'M', JOBD = 'N', ALG = 'C',
C                              BATCH = 'O', or 
C                             (BATCH = 'L' and CONCT = 'N');
C             LDWORK >= 5*(M+L)*NOBR, if METH = 'N', ALG = 'C', and
C                             BATCH = 'L' or 'O';
C             LDWORK >= (M+L)*2*NOBR*(M+L+3), if ALG = 'F',
C                             BATCH <> 'O' and CONCT = 'C';
C             LDWORK >= (M+L)*2*NOBR*(M+L+1), if ALG = 'F',
C                             BATCH = 'F', 'I' and CONCT = 'N';
C             LDWORK >= (M+L)*4*NOBR*(M+L+1)+(M+L)*2*NOBR, if ALG = 'F',
C                             BATCH = 'L' and CONCT = 'N', or
C                             BATCH = 'O';
C             LDWORK >= 4*(M+L)*NOBR, if ALG = 'Q', BATCH = 'F', and 
C                             LDR >= NS = NSMP - 2*NOBR + 1;
C             LDWORK >= max(4*(M+L)*NOBR, 5*L*NOBR), if METH = 'M',
C                             ALG = 'Q', BATCH = 'O', and LDR >= NS;
C             LDWORK >= 5*(M+L)*NOBR, if METH = 'N', ALG = 'Q',
C                             BATCH = 'O', and LDR >= NS;
C             LDWORK >= 6*(M+L)*NOBR, if ALG = 'Q', (BATCH = 'F' or 'O',
C                             and LDR < NS), or (BATCH = 'I' or 
C                             'L' and CONCT = 'N');
C             LDWORK >= 4*(NOBR+1)*(M+L)*NOBR, if ALG = 'Q', BATCH = 'I'
C                             or 'L' and CONCT = 'C'.
C             The workspace used for ALG = 'Q' is 
C                       LDRWRK*2*(M+L)*NOBR + 4*(M+L)*NOBR, 
C             where LDRWRK = LDWORK/(2*(M+L)*NOBR) - 2; recommended 
C             value LDRWRK = NS, assuming a large enough cache size.
C             For good performance,  LDWORK  should be larger.
C
C     Warning Indicator
C
C     IWARN   INTEGER
C             = 0:  no warning;
C             = 1:  the number of 100 cycles in sequential data
C                   processing has been exhausted without signaling
C                   that the last block of data was get; the cycle
C                   counter was reinitialized;
C             = 2:  a fast algorithm was requested (ALG = 'C' or 'F'),
C                   but it failed, and the QR algorithm was then used
C                   (non-sequential data processing);
C             = 3:  all singular values were exactly zero, hence  N = 0
C                   (both input and output were identically zero);
C             = 4:  the least squares problems with coefficient matrix
C                   U_f,  used for computing the weighted oblique
C                   projection (for METH = 'N'), have a rank-deficient
C                   coefficient matrix;
C             = 5:  the least squares problem with coefficient matrix
C                   r_1  [6], used for computing the weighted oblique
C                   projection (for METH = 'N'), has a rank-deficient
C                   coefficient matrix.
C             NOTE: the values 4 and 5 of IWARN have no significance
C                   for the identification problem. 
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             = 1:  a fast algorithm was requested (ALG = 'C', or 'F')
C                   in sequential data processing, but it failed; the
C                   routine can be repeatedly called again using the
C                   standard QR algorithm;
C             = 2:  the singular value decomposition (SVD) algorithm did
C                   not converge.
C
C     METHOD
C
C     The procedure consists in three main steps, the first step being
C     performed by one of the three algorithms included.
C
C     1.a) For non-sequential data processing using QR algorithm, a
C     t x 2(m+l)s  matrix H is constructed, where  
C
C          H = [ Uf'         Up'      Y'      ],  for METH = 'M',
C                  s+1,2s,t    1,s,t   1,2s,t
C
C          H = [ U'       Y'      ],              for METH = 'N',
C                 1,2s,t   1,2s,t
C
C     and  Up     , Uf        , U      , and  Y        are block Hankel
C            1,s,t    s+1,2s,t   1,2s,t        1,2s,t
C     matrices defined in terms of the input and output data [3].
C     A QR factorization is used to compress the data.
C     The fast QR algorithm uses a QR factorization which exploits 
C     the block-Hankel structure. Actually, the Cholesky factor of H'*H
C     is computed.
C
C     1.b) For sequential data processing using QR algorithm, the QR
C     decomposition is done sequentially, by updating the upper
C     triangular factor  R.  This is also performed internally if the
C     workspace is not large enough to accommodate an entire batch.
C
C     1.c) For non-sequential or sequential data processing using
C     Cholesky algorithm, the correlation matrix of input-output data is
C     computed (sequentially, if requested), taking advantage of the
C     block Hankel structure [7].  Then, the Cholesky factor of the
C     correlation matrix is found, if possible.
C
C     2) A singular value decomposition (SVD) of a certain matrix is
C     then computed, which reveals the order  n  of the system as the
C     number of "non-zero" singular values. For the MOESP approach, this
C     matrix is  [ R_24'  R_34' ]' := R(ms+1:(2m+l)s,(2m+l)s+1:2(m+l)s),
C     where  R  is the upper triangular factor  R  constructed by SLICOT
C     Library routine  IB01MD.  For the N4SID approach, a weighted
C     oblique projection is computed from the upper triangular factor  R
C     and its SVD is then found.
C
C     3) The singular values are compared to the given, or default TOL,
C     and the estimated order  n  is returned, possibly after user's
C     confirmation.
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] Verhaegen M.
C         Subspace Model Identification. Part 3: Analysis of the
C         ordinary output-error state-space model identification
C         algorithm.
C         Int. J. Control, 58, pp. 555-586, 1993.
C
C     [3] Verhaegen M.
C         Identification of the deterministic part of MIMO state space
C         models given in innovations form from input-output data.
C         Automatica, Vol.30, No.1, pp.61-74, 1994.
C
C     [4] Van Overschee, P., and De Moor, B.
C         N4SID: Subspace Algorithms for the Identification of
C         Combined Deterministic-Stochastic Systems.
C         Automatica, Vol.30, No.1, pp. 75-93, 1994.
C
C     [5] Peternell, K., Scherrer, W. and Deistler, M.
C         Statistical Analysis of Novel Subspace Identification Methods.
C         Signal Processing, 52, pp. 161-177, 1996.
C
C     [6] Sima, V.
C         Subspace-based Algorithms for Multivariable System 
C         Identification.
C         Studies in Informatics and Control, 5, pp. 335-344, 1996.
C
C     [7] Sima, V.
C         Cholesky or QR Factorization for Data Compression in 
C         Subspace-based Identification ?
C         Proceedings of the Second NICONET Workshop on ``Numerical 
C         Control Software: SLICOT, a Useful Tool in Industry'', 
C         December 3, 1999, INRIA Rocquencourt, France, pp. 75-80, 1999.
C
C     NUMERICAL ASPECTS
C
C     The implemented method is numerically stable (when QR algorithm is
C     used), reliable and efficient. The fast Cholesky or QR algorithms
C     are more efficient, but the accuracy could diminish by forming the
C     correlation matrix.
C     The most time-consuming computational step is step 1:
C                                        2 
C     The QR algorithm needs 0(t(2(m+l)s) ) floating point operations.
C                                           2              3
C     The Cholesky algorithm needs 0(2t(m+l) s)+0((2(m+l)s) ) floating 
C     point operations.
C                                          2           3 2 
C     The fast QR algorithm needs 0(2t(m+l) s)+0(4(m+l) s ) floating 
C     point operations.
C                                                3
C     Step 2 of the algorithm requires 0(((m+l)s) ) floating point 
C     operations.
C
C     FURTHER COMMENTS
C
C     For ALG = 'Q', BATCH = 'O' and LDR < NS, or BATCH <> 'O', the  
C     calculations could be rather inefficient if only minimal workspace
C     (see argument LDWORK) is provided. It is advisable to provide as
C     much workspace as possible. Almost optimal efficiency can be 
C     obtained for  LDWORK = (NS+2)*(2*(M+L)*NOBR),  assuming that the 
C     cache size is large enough to accommodate R, U, Y, and DWORK.
C
C     CONTRIBUTOR
C
C     V. Sima, Katholieke Universiteit Leuven, Feb. 2000.
C
C     REVISIONS
C
C     August 2000.
C
C     KEYWORDS
C
C     Cholesky decomposition, Hankel matrix, identification methods, 
C     multivariable systems, QR decomposition, singular value 
C     decomposition.
C
C     ******************************************************************
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION   RCOND, TOL
      INTEGER            INFO, IWARN, L, LDR, LDU, LDWORK, LDY, M, N,
     $                   NOBR, NSMP
      CHARACTER          ALG, BATCH, CONCT, CTRL, JOBD, METH
C     .. Array Arguments ..
      INTEGER            IWORK(*)
      DOUBLE PRECISION   DWORK(*), R(LDR, *), SV(*), U(LDU, *), 
     $                   Y(LDY, *)
C     .. Local Scalars ..
      INTEGER            IWARNL, LMNOBR, LNOBR, MAXWRK, MINWRK, MNOBR,
     $                   NOBR21, NR, NS, NSMPSM
      LOGICAL            CHALG, CONNEC, CONTRL, FIRST, FQRALG, INTERM,
     $                   JOBDM, LAST, MOESP, N4SID, ONEBCH, QRALG
C     .. External Functions ..
      LOGICAL            LSAME 
      EXTERNAL           LSAME
C     .. External Subroutines ..
      EXTERNAL           IB01MD, IB01ND, IB01OD, XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC          MAX
C     .. Save Statement ..
C        MAXWRK  is used to store the optimal workspace.
C        NSMPSM  is used to sum up the  NSMP  values for  BATCH <> 'O'.
      SAVE               MAXWRK, NSMPSM
C     ..
C     .. Executable Statements ..
C
C     Decode the scalar input parameters.
C
      MOESP  = LSAME( METH,  'M' )
      N4SID  = LSAME( METH,  'N' )
      FQRALG = LSAME( ALG,   'F' )
      QRALG  = LSAME( ALG,   'Q' )
      CHALG  = LSAME( ALG,   'C' )
      JOBDM  = LSAME( JOBD,  'M' )
      ONEBCH = LSAME( BATCH, 'O' )
      FIRST  = LSAME( BATCH, 'F' ) .OR. ONEBCH
      INTERM = LSAME( BATCH, 'I' )
      LAST   = LSAME( BATCH, 'L' ) .OR. ONEBCH
      CONTRL = LSAME( CTRL,  'C' )
C
      IF( .NOT.ONEBCH ) THEN
         CONNEC = LSAME( CONCT, 'C' )
      ELSE
         CONNEC = .FALSE.
      END IF
C
      MNOBR  = M*NOBR
      LNOBR  = L*NOBR
      LMNOBR = LNOBR  + MNOBR
      NR     = LMNOBR + LMNOBR
      NOBR21 = 2*NOBR - 1
      IWARN  = 0
      INFO   = 0
      IF( FIRST ) THEN
         MAXWRK = 1
         NSMPSM = 0
      END IF
      NSMPSM = NSMPSM + NSMP
C
C     Check the scalar input parameters.
C
      IF( .NOT.( MOESP .OR. N4SID ) ) THEN
         INFO = -1
      ELSE IF( .NOT.( FQRALG .OR. QRALG .OR. CHALG ) ) THEN
         INFO = -2
      ELSE IF( MOESP .AND. .NOT.( JOBDM .OR. LSAME( JOBD, 'N' ) ) ) THEN
         INFO = -3
      ELSE IF( .NOT.( FIRST .OR. INTERM .OR. LAST ) ) THEN
         INFO = -4
      ELSE IF( .NOT. ONEBCH ) THEN
         IF( .NOT.( CONNEC .OR. LSAME( CONCT, 'N' ) ) )
     $      INFO = -5
      END IF
      IF( INFO.EQ.0 ) THEN
         IF( .NOT.( CONTRL .OR. LSAME( CTRL, 'N' ) ) ) THEN
            INFO = -6
         ELSE IF( NOBR.LE.0 ) THEN
            INFO = -7
         ELSE IF( M.LT.0 ) THEN
            INFO = -8
         ELSE IF( L.LE.0 ) THEN
            INFO = -9
         ELSE IF( NSMP.LT.2*NOBR .OR.
     $            ( LAST .AND. NSMPSM.LT.NR+NOBR21 ) ) THEN
            INFO = -10
         ELSE IF( LDU.LT.1 .OR. ( M.GT.0 .AND. LDU.LT.NSMP ) ) THEN
            INFO = -12
         ELSE IF( LDY.LT.NSMP ) THEN
            INFO = -14
         ELSE IF( LDR.LT.NR .OR. ( MOESP .AND. JOBDM .AND. 
     $            LDR.LT.3*MNOBR ) ) THEN
            INFO = -17
         ELSE
C
C           Compute workspace.
C           (Note: Comments in the code beginning "Workspace:" describe
C           the minimal amount of workspace needed at that point in the
C           code, as well as the preferred amount for good performance.)
C
            NS = NSMP - NOBR21
            IF ( CHALG ) THEN
               IF ( .NOT.LAST ) THEN
                  IF ( CONNEC ) THEN
                     MINWRK = 2*( NR - M - L )
                  ELSE
                     MINWRK = 1
                  END IF
               ELSE IF ( MOESP ) THEN
                  IF ( CONNEC .AND. .NOT.ONEBCH ) THEN
                     MINWRK = MAX( 2*( NR - M - L ), 5*LNOBR )
                  ELSE
                     MINWRK = 5*LNOBR
                     IF ( JOBDM ) 
     $                  MINWRK = MAX( 2*MNOBR - NOBR, LMNOBR, MINWRK )
                  END IF
               ELSE
                  MINWRK = 5*LMNOBR
               END IF
            ELSE IF ( FQRALG ) THEN
               IF ( .NOT.ONEBCH .AND. CONNEC ) THEN
                  MINWRK = NR*( M + L + 3 )
               ELSE IF ( FIRST .OR. INTERM ) THEN
                  MINWRK = NR*( M + L + 1 )
               ELSE 
                  MINWRK = 2*NR*( M + L + 1 ) + NR
               END IF
            ELSE 
               MINWRK = 2*NR
               IF ( ONEBCH .AND. LDR.GE.NS ) THEN
                  IF ( MOESP ) THEN
                     MINWRK = MAX( MINWRK, 5*LNOBR )
                  ELSE
                     MINWRK = 5*LMNOBR
                  END IF
               END IF
               IF ( FIRST ) THEN
                  IF ( LDR.LT.NS ) THEN
                     MINWRK = MINWRK + NR
                  END IF
               ELSE
                  IF ( CONNEC ) THEN
                     MINWRK = MINWRK*( NOBR + 1 )
                  ELSE
                     MINWRK = MINWRK + NR
                  END IF
               END IF
            END IF
C
            MAXWRK = MINWRK
C
            IF( LDWORK.LT.MINWRK ) THEN
               INFO = -23
               DWORK( 1 ) = MINWRK
            END IF
         END IF
      END IF       
C
C     Return if there are illegal arguments.
C
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'IB01AD', -INFO )
         RETURN
      END IF
C
C     Compress the input-output data.
C     Workspace: need   c*(M+L)*NOBR, where c is a constant depending
C                       on the algorithm and the options used
C                       (see SLICOT Library routine IB01MD);
C                prefer larger.
C
      CALL IB01MD( METH, ALG, BATCH, CONCT, NOBR, M, L, NSMP, U, LDU, Y,
     $             LDY, R, LDR, IWORK, DWORK, LDWORK, IWARN, INFO )
C
      IF ( INFO.EQ.1 ) THEN
C
C        Error return: A fast algorithm was requested (ALG = 'C', 'F')
C        in sequential data processing, but it failed. 
C
         RETURN
      END IF
C
      MAXWRK = MAX( MAXWRK, INT( DWORK( 1 ) ) )
C
      IF ( .NOT.LAST ) THEN
C
C        Return to get new data.
C
         RETURN
      END IF
C
C     Find the singular value decomposition (SVD) giving the system
C     order, and perform related preliminary calculations needed for
C     computing the system matrices.
C     Workspace: need   max( (2*M-1)*NOBR, (M+L)*NOBR, 5*L*NOBR ),
C                                          if METH = 'M';
C                            5*(M+L)*NOBR, if METH = 'N';
C                prefer larger.
C
      CALL IB01ND( METH, JOBD, NOBR, M, L, R, LDR, SV, RCOND, IWORK,
     $             DWORK, LDWORK, IWARNL, INFO )
      IWARN = MAX( IWARN, IWARNL )
C
      IF ( INFO.EQ.2 ) THEN
C
C        Error return: the singular value decomposition (SVD) algorithm
C        did not converge.
C
         RETURN
      END IF
C
C     Estimate the system order.
C
      CALL IB01OD( CTRL, NOBR, L, SV, N, TOL, IWARNL, INFO )
      IWARN = MAX( IWARN, IWARNL )
C
C     Return optimal workspace in  DWORK(1).
C
      DWORK( 1 ) = MAX( MAXWRK,  INT( DWORK( 1 ) ) )
      RETURN
C
C *** Last line of IB01AD ***
      END