SUBROUTINE IB01MD( METH, ALG, BATCH, CONCT, NOBR, M, L, NSMP, U,
     $                   LDU, Y, LDY, R, LDR, IWORK, DWORK, LDWORK,
     $                   IWARN, INFO )
C
C     SLICOT RELEASE 5.0.
C
C     Copyright (c) 2002-2009 NICONET e.V.
C
C     This program is free software: you can redistribute it and/or
C     modify it under the terms of the GNU General Public License as
C     published by the Free Software Foundation, either version 2 of
C     the License, or (at your option) any later version.
C
C     This program is distributed in the hope that it will be useful,
C     but WITHOUT ANY WARRANTY; without even the implied warranty of
C     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C     GNU General Public License for more details.
C
C     You should have received a copy of the GNU General Public License
C     along with this program.  If not, see
C     <http://www.gnu.org/licenses/>.
C
C     PURPOSE
C
C     To construct an upper triangular factor  R  of the concatenated
C     block Hankel matrices using input-output data.  The input-output
C     data can, 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     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     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     R       (output or input/output) DOUBLE PRECISION array, dimension
C             ( LDR,2*(M+L)*NOBR )
C             On exit, if INFO = 0 and ALG = 'Q', or (ALG = 'C' or 'F',
C             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 (current) upper triangular factor
C             R from the QR factorization of the concatenated block
C             Hankel matrices. The diagonal elements of R are positive
C             when the Cholesky algorithm was successfully used.
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 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 >= 2*(M+L)*NOBR.
C
C     Workspace
C
C     IWORK   INTEGER array, dimension (LIWORK)
C             LIWORK >= M+L, if ALG = 'F';
C             LIWORK >= 0,   if ALG = 'C' or 'Q'.
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C             On exit, if  INFO = 0,  DWORK(1)  returns the optimal
C             value of LDWORK.
C             On exit, if  INFO = -17,  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 <> 'O' and
C                                     CONCT = 'C';
C             LDWORK >= 1,            if ALG = 'C', BATCH = 'O' or
C                                     CONCT = 'N';
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' or 'O',
C                                     and LDR >= NS = NSMP - 2*NOBR + 1;
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
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
C     METHOD
C
C     1) 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     2) 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     3) 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     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                                        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
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, Research Institute for Informatics, Bucharest, Aug. 1999.
C
C     REVISIONS
C
C     Feb. 2000, Aug. 2000, Feb. 2004.
C
C     KEYWORDS
C
C     Cholesky decomposition, Hankel matrix, identification methods,
C     multivariable systems, QR decomposition.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
      INTEGER            MAXCYC
      PARAMETER          ( MAXCYC = 100 )
C     .. Scalar Arguments ..
      INTEGER            INFO, IWARN, L, LDR, LDU, LDWORK, LDY, M, NOBR,
     $                   NSMP
      CHARACTER          ALG, BATCH, CONCT, METH
C     .. Array Arguments ..
      INTEGER            IWORK(*)
      DOUBLE PRECISION   DWORK(*), R(LDR, *), U(LDU, *), Y(LDY, *)
C     .. Local Scalars ..
      DOUBLE PRECISION   UPD, TEMP
      INTEGER            I, ICOL, ICYCLE, ID, IERR, II, INICYC, INIT,
     $                   INITI, INU, INY, IREV, ISHFT2, ISHFTU, ISHFTY,
     $                   ITAU, J, JD, JWORK, LDRWMX, LDRWRK, LLDRW,
     $                   LMNOBR, LNOBR, MAXWRK, MINWRK, MLDRW, MMNOBR,
     $                   MNOBR, NCYCLE, NICYCL, NOBR2, NOBR21, NOBRM1,
     $                   NR, NS, NSF, NSL, NSLAST, NSMPSM
      LOGICAL            CHALG, CONNEC, FIRST, FQRALG, INTERM, LAST,
     $                   LINR, MOESP, N4SID, ONEBCH, QRALG
C     .. Local Arrays ..
      DOUBLE PRECISION   DUM( 1 )
C     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ILAENV
      EXTERNAL           ILAENV, LSAME
C     .. External Subroutines ..
      EXTERNAL           DAXPY, DCOPY, DGEMM, DGEQRF, DGER, DLACPY,
     $                   DLASET, DPOTRF, DSWAP, DSYRK, IB01MY, MB04OD,
     $                   XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC          INT, MAX, MIN
C     .. Save Statement ..
C        ICYCLE  is used to count the cycles for  BATCH = 'I'. It is
C                reinitialized at each MAXCYC cycles.
C        MAXWRK  is used to store the optimal workspace.
C        NSMPSM  is used to sum up the  NSMP  values for  BATCH <> 'O'.
      SAVE               ICYCLE, 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' )
      ONEBCH = LSAME( BATCH, 'O' )
      FIRST  = LSAME( BATCH, 'F' ) .OR. ONEBCH
      INTERM = LSAME( BATCH, 'I' )
      LAST   = LSAME( BATCH, 'L' ) .OR. ONEBCH
      IF( .NOT.ONEBCH ) THEN
         CONNEC = LSAME( CONCT, 'C' )
      ELSE
         CONNEC = .FALSE.
      END IF
C
      MNOBR  = M*NOBR
      LNOBR  = L*NOBR
      LMNOBR = LNOBR + MNOBR
      MMNOBR = MNOBR + MNOBR
      NOBRM1 = NOBR - 1
      NOBR21 = NOBR + NOBRM1
      NOBR2  = NOBR21 + 1
      IWARN  = 0
      INFO   = 0
      IERR   = 0
      IF( FIRST ) THEN
         ICYCLE = 1
         MAXWRK = 1
         NSMPSM = 0
      END IF
      NSMPSM = NSMPSM + NSMP
      NR = LMNOBR + LMNOBR
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( .NOT.( FIRST .OR. INTERM .OR. LAST ) ) THEN
         INFO = -3
      ELSE IF( .NOT. ONEBCH ) THEN
         IF( .NOT.( CONNEC .OR. LSAME( CONCT, 'N' ) ) )
     $      INFO = -4
      END IF
      IF( INFO.EQ.0 ) THEN
         IF( NOBR.LE.0 ) THEN
            INFO = -5
         ELSE IF( M.LT.0 ) THEN
            INFO = -6
         ELSE IF( L.LE.0 ) THEN
            INFO = -7
         ELSE IF( NSMP.LT.NOBR2 .OR.
     $          ( LAST .AND. NSMPSM.LT.NR+NOBR21 ) ) THEN
            INFO = -8
         ELSE IF( LDU.LT.1 .OR. ( M.GT.0 .AND. LDU.LT.NSMP ) ) THEN
            INFO = -10
         ELSE IF( LDY.LT.NSMP ) THEN
            INFO = -12
         ELSE IF( LDR.LT.NR ) THEN
            INFO = -14
         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           NB refers to the optimal block size for the immediately
C           following subroutine, as returned by ILAENV.)
C
            NS = NSMP - NOBR21
            IF ( CHALG ) THEN
               IF ( .NOT.ONEBCH .AND. CONNEC ) THEN
                  MINWRK = 2*( NR - M - L )
               ELSE
                  MINWRK = 1
               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
               MAXWRK = NR + NR*ILAENV( 1, 'DGEQRF', ' ', NS, NR, -1,
     $                                  -1 )
               IF ( FIRST ) THEN
                  IF ( LDR.LT.NS ) THEN
                     MINWRK = MINWRK + NR
                     MAXWRK = NS*NR + MAXWRK
                  END IF
               ELSE
                  IF ( CONNEC ) THEN
                     MINWRK = MINWRK*( NOBR + 1 )
                  ELSE
                     MINWRK = MINWRK + NR
                  END IF
                  MAXWRK = NS*NR + MAXWRK
               END IF
            END IF
            MAXWRK = MAX( MINWRK, MAXWRK )
C
            IF( LDWORK.LT.MINWRK ) THEN
               INFO = -17
               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( 'IB01MD', -INFO )
         RETURN
      END IF
C
      IF ( CHALG ) THEN
C
C        Compute the  R  factor from a Cholesky factorization of the
C        input-output data correlation matrix.
C
C        Set the parameters for constructing the correlations of the
C        current block.
C
         LDRWRK = 2*NOBR2 - 2
         IF( FIRST ) THEN
            UPD = ZERO
         ELSE
            UPD = ONE
         END IF
C
         IF( .NOT.FIRST .AND. CONNEC ) THEN
C
C           Restore the saved (M+L)*(2*NOBR-1) "connection" elements of
C           U  and  Y  into their appropriate position in sequential
C           processing. The process is performed column-wise, in
C           reverse order, first for  Y  and then for  U.
C           Workspace: need   (4*NOBR-2)*(M+L).
C
            IREV =     NR - M - L   - NOBR21 + 1
            ICOL = 2*( NR - M - L ) - LDRWRK + 1
C
            DO 10 J = 2, M + L
               DO 5 I = NOBR21 - 1, 0, -1
                  DWORK(ICOL+I) = DWORK(IREV+I)
    5          CONTINUE
               IREV = IREV - NOBR21
               ICOL = ICOL - LDRWRK
   10       CONTINUE
C
            IF ( M.GT.0 )
     $         CALL DLACPY( 'Full', NOBR21, M, U, LDU, DWORK(NOBR2),
     $                      LDRWRK )
            CALL DLACPY( 'Full', NOBR21, L, Y, LDY,
     $                   DWORK(LDRWRK*M+NOBR2), LDRWRK )
         END IF
C
         IF ( M.GT.0 ) THEN
C
C           Let  Guu(i,j) = Guu0(i,j) + u_i*u_j' + u_(i+1)*u_(j+1)' +
C                                 ... + u_(i+NS-1)*u_(j+NS-1)',
C           where  u_i'  is the i-th row of  U,  j = 1 : 2s,  i = 1 : j,
C           NS = NSMP - 2s + 1,  and  Guu0(i,j)  is a zero matrix for
C           BATCH = 'O' or 'F', and it is the matrix Guu(i,j) computed
C           till the current block for BATCH = 'I' or 'L'. The matrix
C           Guu(i,j)  is  m-by-m,  and  Guu(j,j)  is symmetric. The
C           upper triangle of the U-U correlations,  Guu,  is computed
C           (or updated) column-wise in the array  R,  that is, in the
C           order  Guu(1,1),  Guu(1,2),  Guu(2,2),  ...,  Guu(2s,2s).
C           Only the submatrices of the first block-row are fully
C           computed (or updated). The remaining ones are determined
C           exploiting the block-Hankel structure, using the updating
C           formula
C
C           Guu(i+1,j+1) = Guu0(i+1,j+1) - Guu0(i,j) + Guu(i,j) +
C                                 u_(i+NS)*u_(j+NS)' - u_i*u_j'.
C
            IF( .NOT.FIRST ) THEN
C
C              Subtract the contribution of the previous block of data
C              in sequential processing. The columns must be processed
C              in backward order.
C
               DO 20 I = NOBR21*M, 1, -1
                  CALL DAXPY( I, -ONE, R(1,I), 1, R(M+1,M+I), 1 )
   20          CONTINUE
C
            END IF
C
C           Compute/update  Guu(1,1).
C
            IF( .NOT.FIRST .AND. CONNEC )
     $         CALL DSYRK( 'Upper', 'Transpose', M, NOBR21, ONE, DWORK,
     $                     LDRWRK, UPD, R, LDR )
            CALL DSYRK( 'Upper', 'Transpose', M, NS, ONE, U, LDU, UPD,
     $                  R, LDR )
C
            JD = 1
C
            IF( FIRST .OR. .NOT.CONNEC ) THEN
C
               DO 70 J = 2, NOBR2
                  JD = JD + M
                  ID = M  + 1
C
C                 Compute/update  Guu(1,j).
C
                  CALL DGEMM( 'Transpose', 'NoTranspose', M, M, NS, ONE,
     $                        U, LDU, U(J,1), LDU, UPD, R(1,JD), LDR )
C
C                 Compute/update  Guu(2:j,j), exploiting the
C                 block-Hankel structure.
C
                  IF( FIRST ) THEN
C
                     DO 30 I = JD - M, JD - 1
                        CALL DCOPY( I, R(1,I), 1, R(M+1,M+I), 1 )
   30                CONTINUE
C
                  ELSE
C
                     DO 40 I = JD - M, JD - 1
                        CALL DAXPY( I, ONE, R(1,I), 1, R(M+1,M+I), 1 )
   40                CONTINUE
C
                  END IF
C
                  DO 50 I = 2, J - 1
                     CALL DGER( M, M, ONE, U(NS+I-1,1), LDU,
     $                          U(NS+J-1,1), LDU, R(ID,JD), LDR )
                     CALL DGER( M, M, -ONE, U(I-1,1), LDU, U(J-1,1),
     $                          LDU, R(ID,JD), LDR )
                     ID = ID + M
   50             CONTINUE
C
                  DO 60 I = 1, M
                     CALL DAXPY( I, U(NS+J-1,I), U(NS+J-1,1), LDU,
     $                           R(JD,JD+I-1), 1 )
                     CALL DAXPY( I, -U(J-1,I), U(J-1,1), LDU,
     $                           R(JD,JD+I-1), 1 )
   60             CONTINUE
C
   70          CONTINUE
C
            ELSE
C
               DO 120 J = 2, NOBR2
                  JD = JD + M
                  ID = M  + 1
C
C                 Compute/update  Guu(1,j)  for sequential processing
C                 with connected blocks.
C
                  CALL DGEMM( 'Transpose', 'NoTranspose', M, M, NOBR21,
     $                        ONE, DWORK, LDRWRK, DWORK(J), LDRWRK, UPD,
     $                        R(1,JD), LDR )
                  CALL DGEMM( 'Transpose', 'NoTranspose', M, M, NS, ONE,
     $                        U, LDU, U(J,1), LDU, ONE, R(1,JD), LDR )
C
C                 Compute/update  Guu(2:j,j)  for sequential processing
C                 with connected blocks, exploiting the block-Hankel
C                 structure.
C
                  IF( FIRST ) THEN
C
                     DO 80 I = JD - M, JD - 1
                        CALL DCOPY( I, R(1,I), 1, R(M+1,M+I), 1 )
   80                CONTINUE
C
                  ELSE
C
                     DO 90 I = JD - M, JD - 1
                        CALL DAXPY( I, ONE, R(1,I), 1, R(M+1,M+I), 1 )
   90                CONTINUE
C
                  END IF
C
                  DO 100 I = 2, J - 1
                     CALL DGER( M, M, ONE, U(NS+I-1,1), LDU,
     $                          U(NS+J-1,1), LDU, R(ID,JD), LDR )
                     CALL DGER( M, M, -ONE, DWORK(I-1), LDRWRK,
     $                          DWORK(J-1), LDRWRK, R(ID,JD), LDR )
                     ID = ID + M
  100             CONTINUE
C
                  DO 110 I = 1, M
                     CALL DAXPY( I, U(NS+J-1,I), U(NS+J-1,1), LDU,
     $                           R(JD,JD+I-1), 1 )
                     CALL DAXPY( I, -DWORK((I-1)*LDRWRK+J-1),
     $                           DWORK(J-1), LDRWRK, R(JD,JD+I-1), 1 )
  110             CONTINUE
C
  120          CONTINUE
C
            END IF
C
            IF ( LAST .AND. MOESP ) THEN
C
C              Interchange past and future parts for MOESP algorithm.
C              (Only the upper triangular parts are interchanged, and
C              the (1,2) part is transposed in-situ.)
C
               TEMP = R(1,1)
               R(1,1) = R(MNOBR+1,MNOBR+1)
               R(MNOBR+1,MNOBR+1) = TEMP
C
               DO 130 J = 2, MNOBR
                  CALL DSWAP( J, R(1,J), 1, R(MNOBR+1,MNOBR+J), 1 )
                  CALL DSWAP( J-1, R(1,MNOBR+J), 1, R(J,MNOBR+1), LDR )
  130          CONTINUE
C
            END IF
C
C           Let  Guy(i,j) = Guy0(i,j) + u_i*y_j' + u_(i+1)*y_(j+1)' +
C                                 ... + u_(i+NS-1)*y_(j+NS-1)',
C           where  u_i'  is the i-th row of  U,  y_j'  is the j-th row
C           of  Y,  j = 1 : 2s,  i = 1 : 2s,  NS = NSMP - 2s + 1,  and
C           Guy0(i,j)  is a zero matrix for  BATCH = 'O' or 'F', and it
C           is the matrix Guy(i,j) computed till the current block for
C           BATCH = 'I' or 'L'.  Guy(i,j) is m-by-L. The U-Y
C           correlations,  Guy,  are computed (or updated) column-wise
C           in the array  R. Only the submatrices of the first block-
C           column and block-row are fully computed (or updated). The
C           remaining ones are determined exploiting the block-Hankel
C           structure, using the updating formula
C
C           Guy(i+1,j+1) = Guy0(i+1,j+1) - Guy0(i,j) + Guy(i,j) +
C                                 u_(i+NS)*y(j+NS)' - u_i*y_j'.
C
            II = MMNOBR - M
            IF( .NOT.FIRST ) THEN
C
C              Subtract the contribution of the previous block of data
C              in sequential processing. The columns must be processed
C              in backward order.
C
               DO 140 I = NR - L, MMNOBR + 1, -1
                  CALL DAXPY( II, -ONE, R(1,I), 1, R(M+1,L+I), 1 )
  140          CONTINUE
C
            END IF
C
C           Compute/update the first block-column of  Guy,  Guy(i,1).
C
            IF( FIRST .OR. .NOT.CONNEC ) THEN
C
               DO 150 I = 1, NOBR2
                  CALL DGEMM( 'Transpose', 'NoTranspose', M, L, NS, ONE,
     $                        U(I,1), LDU, Y, LDY, UPD,
     $                        R((I-1)*M+1,MMNOBR+1), LDR )
  150          CONTINUE
C
            ELSE
C
               DO 160 I = 1, NOBR2
                  CALL DGEMM( 'Transpose', 'NoTranspose', M, L, NOBR21,
     $                        ONE, DWORK(I), LDRWRK, DWORK(LDRWRK*M+1),
     $                        LDRWRK, UPD, R((I-1)*M+1,MMNOBR+1), LDR )
                  CALL DGEMM( 'Transpose', 'NoTranspose', M, L, NS, ONE,
     $                        U(I,1), LDU, Y, LDY, ONE,
     $                        R((I-1)*M+1,MMNOBR+1), LDR )
  160          CONTINUE
C
            END IF
C
            JD = MMNOBR + 1
C
            IF( FIRST .OR. .NOT.CONNEC ) THEN
C
               DO 200 J = 2, NOBR2
                  JD = JD + L
                  ID = M  + 1
C
C                 Compute/update  Guy(1,j).
C
                  CALL DGEMM( 'Transpose', 'NoTranspose', M, L, NS, ONE,
     $                        U, LDU, Y(J,1), LDY, UPD, R(1,JD), LDR )
C
C                 Compute/update  Guy(2:2*s,j), exploiting the
C                 block-Hankel structure.
C
                  IF( FIRST ) THEN
C
                     DO 170 I = JD - L, JD - 1
                        CALL DCOPY( II, R(1,I), 1, R(M+1,L+I), 1 )
  170                CONTINUE
C
                  ELSE
C
                     DO 180 I = JD - L, JD - 1
                        CALL DAXPY( II, ONE, R(1,I), 1, R(M+1,L+I), 1 )
  180                CONTINUE
C
                  END IF
C
                  DO 190 I = 2, NOBR2
                     CALL DGER( M, L, ONE, U(NS+I-1,1), LDU,
     $                          Y(NS+J-1,1), LDY, R(ID,JD), LDR )
                     CALL DGER( M, L, -ONE, U(I-1,1), LDU, Y(J-1,1),
     $                          LDY, R(ID,JD), LDR )
                     ID = ID + M
  190             CONTINUE
C
  200          CONTINUE
C
            ELSE
C
               DO 240 J = 2, NOBR2
                  JD = JD + L
                  ID = M  + 1
C
C                 Compute/update  Guy(1,j)  for sequential processing
C                 with connected blocks.
C
                  CALL DGEMM( 'Transpose', 'NoTranspose', M, L, NOBR21,
     $                        ONE, DWORK, LDRWRK, DWORK(LDRWRK*M+J),
     $                        LDRWRK, UPD, R(1,JD), LDR )
                  CALL DGEMM( 'Transpose', 'NoTranspose', M, L, NS, ONE,
     $                        U, LDU, Y(J,1), LDY, ONE, R(1,JD), LDR )
C
C                 Compute/update  Guy(2:2*s,j)  for sequential
C                 processing with connected blocks, exploiting the
C                 block-Hankel structure.
C
                  IF( FIRST ) THEN
C
                     DO 210 I = JD - L, JD - 1
                        CALL DCOPY( II, R(1,I), 1, R(M+1,L+I), 1 )
  210                CONTINUE
C
                  ELSE
C
                     DO 220 I = JD - L, JD - 1
                        CALL DAXPY( II, ONE, R(1,I), 1, R(M+1,L+I), 1 )
  220                CONTINUE
C
                  END IF
C
                  DO 230 I = 2, NOBR2
                     CALL DGER( M, L, ONE, U(NS+I-1,1), LDU,
     $                          Y(NS+J-1,1), LDY, R(ID,JD), LDR )
                     CALL DGER( M, L, -ONE, DWORK(I-1), LDRWRK,
     $                          DWORK(LDRWRK*M+J-1), LDRWRK, R(ID,JD),
     $                          LDR )
                     ID = ID + M
  230             CONTINUE
C
  240          CONTINUE
C
            END IF
C
            IF ( LAST .AND. MOESP ) THEN
C
C              Interchange past and future parts of U-Y correlations
C              for MOESP algorithm.
C
               DO 250 J = MMNOBR + 1, NR
                  CALL DSWAP( MNOBR, R(1,J), 1, R(MNOBR+1,J), 1 )
  250          CONTINUE
C
            END IF
         END IF
C
C        Let  Gyy(i,j) = Gyy0(i,j) + y_i*y_i' + y_(i+1)*y_(i+1)' + ... +
C                                    y_(i+NS-1)*y_(i+NS-1)',
C        where  y_i'  is the i-th row of  Y,  j = 1 : 2s,  i = 1 : j,
C        NS = NSMP - 2s + 1,  and  Gyy0(i,j)  is a zero matrix for
C        BATCH = 'O' or 'F', and it is the matrix Gyy(i,j) computed till
C        the current block for BATCH = 'I' or 'L'.  Gyy(i,j) is L-by-L,
C        and  Gyy(j,j)  is symmetric. The upper triangle of the Y-Y
C        correlations,  Gyy,  is computed (or updated) column-wise in
C        the corresponding part of the array  R,  that is, in the order
C        Gyy(1,1),  Gyy(1,2),  Gyy(2,2),  ...,  Gyy(2s,2s).  Only the
C        submatrices of the first block-row are fully computed (or
C        updated). The remaining ones are determined exploiting the
C        block-Hankel structure, using the updating formula
C
C        Gyy(i+1,j+1) = Gyy0(i+1,j+1) - Gyy0(i,j) + Gyy(i,j) +
C                              y_(i+NS)*y_(j+NS)' - y_i*y_j'.
C
         JD = MMNOBR + 1
C
         IF( .NOT.FIRST ) THEN
C
C           Subtract the contribution of the previous block of data
C           in sequential processing. The columns must be processed in
C           backward order.
C
            DO 260 I = NR - L, MMNOBR + 1, -1
               CALL DAXPY( I-MMNOBR, -ONE, R(JD,I), 1, R(JD+L,L+I), 1 )
  260       CONTINUE
C
         END IF
C
C        Compute/update  Gyy(1,1).
C
         IF( .NOT.FIRST .AND. CONNEC )
     $      CALL DSYRK( 'Upper', 'Transpose', L, NOBR21, ONE,
     $                  DWORK(LDRWRK*M+1), LDRWRK, UPD, R(JD,JD), LDR )
         CALL DSYRK( 'Upper', 'Transpose', L, NS, ONE, Y, LDY, UPD,
     $               R(JD,JD), LDR )
C
         IF( FIRST .OR. .NOT.CONNEC ) THEN
C
            DO 310 J = 2, NOBR2
               JD = JD + L
               ID = MMNOBR + L + 1
C
C              Compute/update  Gyy(1,j).
C
               CALL DGEMM( 'Transpose', 'NoTranspose', L, L, NS, ONE, Y,
     $                     LDY, Y(J,1), LDY, UPD, R(MMNOBR+1,JD), LDR )
C
C              Compute/update  Gyy(2:j,j), exploiting the block-Hankel
C              structure.
C
               IF( FIRST ) THEN
C
                  DO 270 I = JD - L, JD - 1
                     CALL DCOPY( I-MMNOBR, R(MMNOBR+1,I), 1,
     $                           R(MMNOBR+L+1,L+I), 1 )
  270             CONTINUE
C
               ELSE
C
                  DO 280 I = JD - L, JD - 1
                     CALL DAXPY( I-MMNOBR, ONE, R(MMNOBR+1,I), 1,
     $                           R(MMNOBR+L+1,L+I), 1 )
  280             CONTINUE
C
               END IF
C
               DO 290 I = 2, J - 1
                  CALL DGER( L, L, ONE, Y(NS+I-1,1), LDY, Y(NS+J-1,1),
     $                       LDY, R(ID,JD), LDR )
                  CALL DGER( L, L, -ONE, Y(I-1,1), LDY, Y(J-1,1), LDY,
     $                       R(ID,JD), LDR )
                  ID = ID + L
  290          CONTINUE
C
               DO 300 I = 1, L
                  CALL DAXPY( I, Y(NS+J-1,I), Y(NS+J-1,1), LDY,
     $                        R(JD,JD+I-1), 1 )
                  CALL DAXPY( I, -Y(J-1,I), Y(J-1,1), LDY, R(JD,JD+I-1),
     $                        1 )
  300          CONTINUE
C
  310       CONTINUE
C
         ELSE
C
            DO 360 J = 2, NOBR2
               JD = JD + L
               ID = MMNOBR + L + 1
C
C              Compute/update  Gyy(1,j)  for sequential processing with
C              connected blocks.
C
               CALL DGEMM( 'Transpose', 'NoTranspose', L, L, NOBR21,
     $                     ONE, DWORK(LDRWRK*M+1), LDRWRK,
     $                     DWORK(LDRWRK*M+J), LDRWRK, UPD,
     $                     R(MMNOBR+1,JD), LDR )
               CALL DGEMM( 'Transpose', 'NoTranspose', L, L, NS, ONE, Y,
     $                     LDY, Y(J,1), LDY, ONE, R(MMNOBR+1,JD), LDR )
C
C              Compute/update  Gyy(2:j,j)  for sequential processing
C              with connected blocks, exploiting the block-Hankel
C              structure.
C
               IF( FIRST ) THEN
C
                  DO 320 I = JD - L, JD - 1
                     CALL DCOPY( I-MMNOBR, R(MMNOBR+1,I), 1,
     $                           R(MMNOBR+L+1,L+I), 1 )
  320             CONTINUE
C
               ELSE
C
                  DO 330 I = JD - L, JD - 1
                     CALL DAXPY( I-MMNOBR, ONE, R(MMNOBR+1,I), 1,
     $                           R(MMNOBR+L+1,L+I), 1 )
  330             CONTINUE
C
               END IF
C
               DO 340 I = 2, J - 1
                  CALL DGER( L, L, ONE, Y(NS+I-1,1), LDY, Y(NS+J-1,1),
     $                       LDY, R(ID,JD), LDR )
                  CALL DGER( L, L, -ONE, DWORK(LDRWRK*M+I-1), LDRWRK,
     $                       DWORK(LDRWRK*M+J-1), LDRWRK, R(ID,JD),
     $                       LDR )
                  ID = ID + L
  340          CONTINUE
C
               DO 350 I = 1, L
                  CALL DAXPY( I, Y(NS+J-1,I), Y(NS+J-1,1), LDY,
     $                        R(JD,JD+I-1), 1 )
                  CALL DAXPY( I, -DWORK(LDRWRK*(M+I-1)+J-1),
     $                        DWORK(LDRWRK*M+J-1), LDRWRK, R(JD,JD+I-1),
     $                        1 )
  350          CONTINUE
C
  360       CONTINUE
C
         END IF
C
         IF ( .NOT.LAST ) THEN
            IF ( CONNEC ) THEN
C
C              For sequential processing with connected data blocks,
C              save the remaining ("connection") elements of  U  and  Y
C              in the first  (M+L)*(2*NOBR-1)  locations of  DWORK.
C
               IF ( M.GT.0 )
     $            CALL DLACPY( 'Full', NOBR21, M, U(NS+1,1), LDU, DWORK,
     $                         NOBR21 )
               CALL DLACPY( 'Full', NOBR21, L, Y(NS+1,1), LDY,
     $                      DWORK(MMNOBR-M+1), NOBR21 )
            END IF
C
C           Return to get new data.
C
            ICYCLE = ICYCLE + 1
            IF ( ICYCLE.GT.MAXCYC )
     $         IWARN = 1
            RETURN
C
         ELSE
C
C           Try to compute the Cholesky factor of the correlation
C           matrix.
C
            CALL DPOTRF( 'Upper', NR, R, LDR, IERR )
            GO TO 370
         END IF
      ELSE IF ( FQRALG ) THEN
C
C        Compute the  R  factor from a fast QR factorization of the
C        input-output data correlation matrix.
C
         CALL IB01MY( METH, BATCH, CONCT, NOBR, M, L, NSMP, U, LDU,
     $                Y, LDY, R, LDR, IWORK, DWORK, LDWORK, IWARN,
     $                IERR )
         IF( .NOT.LAST )
     $      RETURN
         MAXWRK = MAX( MAXWRK, INT( DWORK(1) ) )
      END IF
C
  370 CONTINUE
C
      IF( IERR.NE.0 ) THEN
C
C        Error return from a fast factorization algorithm of the
C        input-output data correlation matrix.
C
         IF( ONEBCH ) THEN
            QRALG  = .TRUE.
            IWARN  = 2
            MINWRK = 2*NR
            MAXWRK = NR + NR*ILAENV( 1, 'DGEQRF', ' ', NS, NR, -1,
     $                               -1 )
            IF ( LDR.LT.NS ) THEN
               MINWRK = MINWRK + NR
               MAXWRK = NS*NR + MAXWRK
            END IF
            MAXWRK = MAX( MINWRK, MAXWRK )
C
            IF( LDWORK.LT.MINWRK ) THEN
               INFO = -17
C
C              Return: Not enough workspace.
C
               DWORK( 1 ) = MINWRK
               CALL XERBLA( 'IB01MD', -INFO )
               RETURN
            END IF
         ELSE
            INFO = 1
            RETURN
         END IF
      END IF
C
      IF ( QRALG ) THEN
C
C        Compute the  R  factor from a QR factorization of the matrix  H
C        of concatenated block Hankel matrices.
C
C        Construct the matrix  H.
C
C        Set the parameters for constructing the current segment of the
C        Hankel matrix, taking the available memory space into account.
C        INITI+1 points to the beginning rows of  U  and  Y  from which
C                data are taken when NCYCLE > 1 inner cycles are needed,
C                or for sequential processing with connected blocks.
C        LDRWMX is the number of rows that can fit in the working space.
C        LDRWRK is the actual number of rows processed in this space.
C        NSLAST is the number of samples to be processed at the last
C               inner cycle.
C
         INITI  = 0
         LDRWMX = LDWORK / NR - 2
         NCYCLE = 1
         NSLAST = NSMP
         LINR   = .FALSE.
         IF ( FIRST ) THEN
            LINR   = LDR.GE.NS
            LDRWRK = NS
         ELSE IF ( CONNEC ) THEN
            LDRWRK = NSMP
         ELSE
            LDRWRK = NS
         END IF
         INICYC = 1
C
         IF ( .NOT.LINR ) THEN
            IF ( LDRWMX.LT.LDRWRK ) THEN
C
C              Not enough working space for doing a single inner cycle.
C              NCYCLE inner cycles are to be performed for the current
C              data block using the working space.
C
               NCYCLE = LDRWRK / LDRWMX
               NSLAST = MOD( LDRWRK, LDRWMX )
               IF ( NSLAST.NE.0 ) THEN
                  NCYCLE = NCYCLE + 1
               ELSE
                  NSLAST = LDRWMX
               END IF
               LDRWRK = LDRWMX
               NS = LDRWRK
               IF ( FIRST ) INICYC = 2
            END IF
            MLDRW = M*LDRWRK
            LLDRW = L*LDRWRK
            INU = MLDRW*NOBR  + 1
            INY = MLDRW*NOBR2 + 1
         END IF
C
C        Process the data given at the current call.
C
         IF ( .NOT.FIRST .AND. CONNEC ) THEN
C
C           Restore the saved (M+L)*(2*NOBR-1) "connection" elements of
C           U  and  Y  into their appropriate position in sequential
C           processing. The process is performed column-wise, in
C           reverse order, first for  Y  and then for  U.
C
            IREV = NR - M - L - NOBR21 + 1
            ICOL = INY + LLDRW - LDRWRK
C
            DO 380 J = 1, L
               DO 375 I = NOBR21 - 1, 0, -1
                  DWORK(ICOL+I) = DWORK(IREV+I)
  375          CONTINUE
               IREV = IREV - NOBR21
               ICOL = ICOL - LDRWRK
  380       CONTINUE
C
            IF( MOESP ) THEN
               ICOL = INU + MLDRW - LDRWRK
            ELSE
               ICOL = MLDRW - LDRWRK + 1
            END IF
C
            DO 390 J = 1, M
               DO 385 I = NOBR21 - 1, 0, -1
                  DWORK(ICOL+I) = DWORK(IREV+I)
  385          CONTINUE
               IREV = IREV - NOBR21
               ICOL = ICOL - LDRWRK
  390       CONTINUE
C
            IF( MOESP )
     $         CALL DLACPY( 'Full', NOBRM1, M, DWORK(INU+NOBR), LDRWRK,
     $                      DWORK, LDRWRK )
         END IF
C
C        Data compression using QR factorization.
C
         IF ( FIRST ) THEN
C
C           Non-sequential data processing or first block in
C           sequential data processing:
C           Use the general QR factorization algorithm.
C
            IF ( LINR ) THEN
C
C              Put the input-output data in the array  R.
C
               IF( M.GT.0 ) THEN
                  IF( MOESP ) THEN
C
                     DO 400 I = 1, NOBR
                        CALL DLACPY( 'Full', NS, M, U(NOBR+I,1), LDU,
     $                               R(1,M*(I-1)+1), LDR )
  400                CONTINUE
C
                     DO 410 I = 1, NOBR
                        CALL DLACPY( 'Full', NS, M, U(I,1), LDU,
     $                               R(1,MNOBR+M*(I-1)+1), LDR )
  410                CONTINUE
C
                  ELSE
C
                     DO 420 I = 1, NOBR2
                        CALL DLACPY( 'Full', NS, M, U(I,1), LDU,
     $                               R(1,M*(I-1)+1), LDR )
  420                CONTINUE
C
                  END IF
               END IF
C
               DO 430 I = 1, NOBR2
                  CALL DLACPY( 'Full', NS, L, Y(I,1), LDY,
     $                         R(1,MMNOBR+L*(I-1)+1), LDR )
  430          CONTINUE
C
C              Workspace: need   4*(M+L)*NOBR,
C                         prefer 2*(M+L)*NOBR+2*(M+L)*NOBR*NB.
C
               ITAU  = 1
               JWORK = ITAU + NR
               CALL DGEQRF( NS, NR, R, LDR, DWORK(ITAU), DWORK(JWORK),
     $                      LDWORK-JWORK+1, IERR )
            ELSE
C
C              Put the input-output data in the array  DWORK.
C
               IF( M.GT.0 ) THEN
                  ISHFTU = 1
                  IF( MOESP ) THEN
                     ISHFT2 = INU
C
                     DO 440 I = 1, NOBR
                        CALL DLACPY( 'Full', NS, M, U(NOBR+I,1), LDU,
     $                               DWORK(ISHFTU), LDRWRK )
                        ISHFTU = ISHFTU + MLDRW
  440                CONTINUE
C
                     DO 450 I = 1, NOBR
                        CALL DLACPY( 'Full', NS, M, U(I,1), LDU,
     $                               DWORK(ISHFT2), LDRWRK )
                        ISHFT2 = ISHFT2 + MLDRW
  450                CONTINUE
C
                  ELSE
C
                     DO 460 I = 1, NOBR2
                        CALL DLACPY( 'Full', NS, M, U(I,1), LDU,
     $                               DWORK(ISHFTU), LDRWRK )
                        ISHFTU = ISHFTU + MLDRW
  460                CONTINUE
C
                  END IF
               END IF
C
               ISHFTY = INY
C
               DO 470 I = 1, NOBR2
                  CALL DLACPY( 'Full', NS, L, Y(I,1), LDY,
     $                         DWORK(ISHFTY), LDRWRK )
                  ISHFTY = ISHFTY + LLDRW
  470          CONTINUE
C
C              Workspace: need   2*(M+L)*NOBR + 4*(M+L)*NOBR,
C                         prefer NS*2*(M+L)*NOBR + 2*(M+L)*NOBR
C                                                + 2*(M+L)*NOBR*NB,
C                         used LDRWRK*2*(M+L)*NOBR + 4*(M+L)*NOBR,
C                         where  NS = NSMP - 2*NOBR + 1,
C                            LDRWRK = min(NS, LDWORK/(2*(M+L)*NOBR)-2).
C
               ITAU  = LDRWRK*NR + 1
               JWORK = ITAU + NR
               CALL DGEQRF( NS, NR, DWORK, LDRWRK, DWORK(ITAU),
     $                      DWORK(JWORK), LDWORK-JWORK+1, IERR )
               CALL DLACPY( 'Upper ', MIN(NS,NR), NR, DWORK, LDRWRK, R,
     $                      LDR )
            END IF
C
            IF ( NS.LT.NR )
     $         CALL DLASET( 'Upper ', NR - NS, NR - NS, ZERO, ZERO,
     $                      R(NS+1,NS+1), LDR )
            INITI = INITI + NS
         END IF
C
         IF ( NCYCLE.GT.1 .OR. .NOT.FIRST ) THEN
C
C           Remaining segments of the first data block or
C           remaining segments/blocks in sequential data processing:
C           Use a structure-exploiting QR factorization algorithm.
C
            NSL = LDRWRK
            IF ( .NOT.CONNEC ) NSL = NS
            ITAU  = LDRWRK*NR + 1
            JWORK = ITAU + NR
C
            DO 560 NICYCL = INICYC, NCYCLE
C
C              INIT  denotes the beginning row where new data are put.
C
               IF ( CONNEC .AND. NICYCL.EQ.1 ) THEN
                  INIT = NOBR2
               ELSE
                  INIT = 1
               END IF
               IF ( NCYCLE.GT.1 .AND. NICYCL.EQ.NCYCLE ) THEN
C
C                 Last samples in the last data segment of a block.
C
                  NS  = NSLAST
                  NSL = NSLAST
               END IF
C
C              Put the input-output data in the array  DWORK.
C
               NSF = NS
               IF ( INIT.GT.1 .AND. NCYCLE.GT.1 ) NSF = NSF - NOBR21
               IF ( M.GT.0 ) THEN
                  ISHFTU = INIT
C
                  IF( MOESP ) THEN
                     ISHFT2 = INIT + INU - 1
C
                     DO 480 I = 1, NOBR
                        CALL DLACPY( 'Full', NSF, M, U(INITI+NOBR+I,1),
     $                               LDU, DWORK(ISHFTU), LDRWRK )
                        ISHFTU = ISHFTU + MLDRW
  480                CONTINUE
C
                     DO 490 I = 1, NOBR
                        CALL DLACPY( 'Full', NSF, M, U(INITI+I,1), LDU,
     $                               DWORK(ISHFT2), LDRWRK )
                        ISHFT2 = ISHFT2 + MLDRW
  490                CONTINUE
C
                  ELSE
C
                     DO 500 I = 1, NOBR2
                        CALL DLACPY( 'Full', NSF, M, U(INITI+I,1), LDU,
     $                               DWORK(ISHFTU), LDRWRK )
                        ISHFTU = ISHFTU + MLDRW
  500                CONTINUE
C
                  END IF
               END IF
C
               ISHFTY = INIT + INY - 1
C
               DO 510 I = 1, NOBR2
                  CALL DLACPY( 'Full', NSF, L, Y(INITI+I,1), LDY,
     $                         DWORK(ISHFTY), LDRWRK )
                  ISHFTY = ISHFTY + LLDRW
  510          CONTINUE
C
               IF ( INIT.GT.1 ) THEN
C
C                 Prepare the connection to the previous block of data
C                 in sequential processing.
C
                  IF( MOESP .AND. M.GT.0 )
     $               CALL DLACPY( 'Full', NOBR, M, U, LDU, DWORK(NOBR),
     $                            LDRWRK )
C
C                 Shift the elements from the connection to the previous
C                 block of data in sequential processing.
C
                  IF ( M.GT.0 ) THEN
                     ISHFTU = MLDRW + 1
C
                     IF( MOESP ) THEN
                        ISHFT2 = MLDRW + INU
C
                        DO 520 I = 1, NOBRM1
                           CALL DLACPY( 'Full', NOBR21, M,
     $                                  DWORK(ISHFTU-MLDRW+1), LDRWRK,
     $                                  DWORK(ISHFTU), LDRWRK )
                           ISHFTU = ISHFTU + MLDRW
  520                  CONTINUE
C
                        DO 530 I = 1, NOBRM1
                           CALL DLACPY( 'Full', NOBR21, M,
     $                                  DWORK(ISHFT2-MLDRW+1), LDRWRK,
     $                                  DWORK(ISHFT2), LDRWRK )
                          ISHFT2 = ISHFT2 + MLDRW
  530                  CONTINUE
C
                     ELSE
C
                        DO 540 I = 1, NOBR21
                           CALL DLACPY( 'Full', NOBR21, M,
     $                                  DWORK(ISHFTU-MLDRW+1), LDRWRK,
     $                                  DWORK(ISHFTU), LDRWRK )
                           ISHFTU = ISHFTU + MLDRW
  540                   CONTINUE
C
                     END IF
                  END IF
C
                  ISHFTY = LLDRW + INY
C
                  DO 550 I = 1, NOBR21
                     CALL DLACPY( 'Full', NOBR21, L,
     $                            DWORK(ISHFTY-LLDRW+1), LDRWRK,
     $                            DWORK(ISHFTY), LDRWRK )
                     ISHFTY = ISHFTY + LLDRW
  550             CONTINUE
C
               END IF
C
C              Workspace: need LDRWRK*2*(M+L)*NOBR + 4*(M+L)*NOBR.
C
               CALL MB04OD( 'Full', NR, 0, NSL, R, LDR, DWORK, LDRWRK,
     $                      DUM, NR, DUM, NR, DWORK(ITAU), DWORK(JWORK)
     $                    )
               INITI = INITI + NSF
  560       CONTINUE
C
         END IF
C
         IF ( .NOT.LAST ) THEN
            IF ( CONNEC ) THEN
C
C              For sequential processing with connected data blocks,
C              save the remaining ("connection") elements of  U  and  Y
C              in the first  (M+L)*(2*NOBR-1)  locations of  DWORK.
C
               IF ( M.GT.0 )
     $            CALL DLACPY( 'Full', NOBR21, M, U(INITI+1,1), LDU,
     $                         DWORK, NOBR21 )
               CALL DLACPY( 'Full', NOBR21, L, Y(INITI+1,1), LDY,
     $                      DWORK(MMNOBR-M+1), NOBR21 )
            END IF
C
C           Return to get new data.
C
            ICYCLE = ICYCLE + 1
            IF ( ICYCLE.LE.MAXCYC )
     $         RETURN
            IWARN  = 1
            ICYCLE = 1
C
         END IF
C
      END IF
C
C     Return optimal workspace in  DWORK(1).
C
      DWORK( 1 ) = MAXWRK
      IF ( LAST ) THEN
         ICYCLE = 1
         MAXWRK = 1
         NSMPSM = 0
      END IF
      RETURN
C
C *** Last line of IB01MD ***
      END