SUBROUTINE BSTREXC( N, T, LDT, IFST, ILST, NITRAF, ITRAF,
     $                    NDTRAF, DTRAF, WORK, INFO )
      IMPLICIT NONE
*
*
*  -- LAPACK routine (version 3.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     March 31, 1993
*
*     .. Scalar Arguments ..
      INTEGER            IFST, ILST, INFO, LDT, N, NDTRAF, NITRAF
*     ..
*     .. Array Arguments ..
      INTEGER            ITRAF( * )
      REAL               DTRAF( * ), T( LDT, * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  BSTREXC reorders the real Schur factorization of a real matrix
*  A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
*  moved to row ILST.
*
*  The real Schur form T is reordered by an orthogonal similarity
*  transformation Z**T*T*Z. In contrast to the LAPACK routine DTREXC,
*  the orthogonal matrix Z is not explicitly constructed but
*  represented by paramaters contained in the arrays ITRAF and DTRAF,
*  see further details.
*
*  T must be in Schur canonical form (as returned by DHSEQR), that is,
*  block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
*  2-by-2 diagonal block has its diagonal elements equal and its
*  off-diagonal elements of opposite sign.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The order of the matrix T. N >= 0.
*
*  T       (input/output) REAL             array, dimension (LDT,N)
*          On entry, the upper quasi-triangular matrix T, in Schur
*          Schur canonical form.
*          On exit, the reordered upper quasi-triangular matrix, again
*          in Schur canonical form.
*
*  LDT     (input) INTEGER
*          The leading dimension of the array T. LDT >= max(1,N).
*
*  IFST    (input/output) INTEGER
*  ILST    (input/output) INTEGER
*          Specify the reordering of the diagonal blocks of T.
*          The block with row index IFST is moved to row ILST, by a
*          sequence of transpositions between adjacent blocks.
*          On exit, if IFST pointed on entry to the second row of a
*          2-by-2 block, it is changed to point to the first row; ILST
*          always points to the first row of the block in its final
*          position (which may differ from its input value by +1 or -1).
*          1 <= IFST <= N; 1 <= ILST <= N.
*
*  NITRAF  (input/output) INTEGER
*          On entry, length of the array ITRAF.
*          As a minimum requirement, NITRAF >= max(1,|ILST-IFST|).
*          If there are 2-by-2 blocks in T then NITRAF must be larger;
*          a safe choice is NITRAF >= max(1,2*|ILST-IFST|).
*          On exit, actual length of the array ITRAF.
*
*  ITRAF   (output) INTEGER array, length NITRAF
*          List of parameters for representing the transformation
*          matrix Z, see further details.
*
*  NDTRAF  (input/output) INTEGER
*          On entry, length of the array DTRAF.
*          As a minimum requirement, NDTRAF >= max(1,2*|ILST-IFST|).
*          If there are 2-by-2 blocks in T then NDTRAF must be larger;
*          a safe choice is NDTRAF >= max(1,5*|ILST-IFST|).
*          On exit, actual length of the array DTRAF.
*
*  DTRAF   (output) REAL             array, length NDTRAF
*          List of parameters for representing the transformation
*          matrix Z, see further details.
*
*  WORK    (workspace) REAL             array, dimension (N)
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
*          = 1:  two adjacent blocks were too close to swap (the problem
*                is very ill-conditioned); T may have been partially
*                reordered, and ILST points to the first row of the
*                current position of the block being moved.
*          = 2:  the 2 by 2 block to be reordered split into two 1 by 1
*                blocks and the second block failed to swap with an
*                adjacent block. ILST points to the first row of the
*                current position of the whole block being moved.
*
*  Further Details
*  ===============
*
*  The orthogonal transformation matrix Z is a product of NITRAF
*  elementary orthogonal transformations. The parameters defining these
*  transformations are stored in the arrays ITRAF and DTRAF as follows:
*
*  Consider the i-th transformation acting on rows/columns POS,
*  POS+1, ... If this transformation is
*
*     (1) a Givens rotation with cosine C and sine S then
*
*           ITRAF(i) = POS,
*           DTRAF(j) = C,    DTRAF(j+1) = S;
*
*     (2) a Householder reflector H = I - tau * v * v' with
*         v = [ 1; v2; v3 ] then
*
*           ITRAF(i) = N + POS,
*           DTRAF(j) = tau,  DTRAF(j+1) = v2,  DTRAF(j+2) = v3;
*
*     (3) a Householder reflector H = I - tau * v * v' with
*         v = [ v1; v2; 1 ] then
*
*           ITRAF(i) = 2*N + POS,
*           DTRAF(j) = v1,  DTRAF(j+1) = v2,  DTRAF(j+2) = tau;
*
*  Note that the parameters in DTRAF are stored consecutively.
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO
      PARAMETER          ( ZERO = 0.0E+0 )
      INTEGER            DLNGTH(3), ILNGTH(3)
*     ..
*     .. Local Scalars ..
      INTEGER            CDTRAF, CITRAF, LDTRAF, LITRAF, HERE, I, NBF,
     $                   NBL, NBNEXT
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     ..
*     .. External Subroutines ..
      EXTERNAL           BSLAEXC, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX
*     .. Data Statements ..
c      DATA ( ILNGTH(I), I = 1, 3 ) / 1, 2, 4 /
c      DATA ( DLNGTH(I), I = 1, 3 ) / 2, 5, 10 /
      DATA ILNGTH(1)/1/, ILNGTH(2)/2/, ILNGTH(3)/4/
      DATA DLNGTH(1)/2/, DLNGTH(2)/5/, DLNGTH(3)/10/
*     ..
*     .. Executable Statements ..
*
*     Decode and test the input arguments.
*
      INFO = 0
      IF( N.LT.0 ) THEN
         INFO = -1
      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
         INFO = -3
      ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
         INFO = -4
      ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
         INFO = -5
      ELSE IF ( NITRAF.LT.MAX( 1, ABS( ILST-IFST ) ) ) THEN
         INFO = -6
      ELSE IF ( NDTRAF.LT.MAX( 1, 2*ABS( ILST-IFST ) ) ) THEN
         INFO = -8
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'DTREXC', -INFO )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( N.LE.1 )
     $   RETURN
      CITRAF = 1
      CDTRAF = 1
*
*     Determine the first row of specified block
*     and find out it is 1 by 1 or 2 by 2.
*
      IF( IFST.GT.1 ) THEN
         IF( T( IFST, IFST-1 ).NE.ZERO )
     $      IFST = IFST - 1
      END IF
      NBF = 1
      IF( IFST.LT.N ) THEN
         IF( T( IFST+1, IFST ).NE.ZERO )
     $      NBF = 2
      END IF
*
*     Determine the first row of the final block
*     and find out it is 1 by 1 or 2 by 2.
*
      IF( ILST.GT.1 ) THEN
         IF( T( ILST, ILST-1 ).NE.ZERO )
     $      ILST = ILST - 1
      END IF
      NBL = 1
      IF( ILST.LT.N ) THEN
         IF( T( ILST+1, ILST ).NE.ZERO )
     $      NBL = 2
      END IF
*
      IF( IFST.EQ.ILST )
     $   RETURN
*
      IF( IFST.LT.ILST ) THEN
*
*        Update ILST
*
         IF( NBF.EQ.2 .AND. NBL.EQ.1 )
     $      ILST = ILST - 1
         IF( NBF.EQ.1 .AND. NBL.EQ.2 )
     $      ILST = ILST + 1
*
         HERE = IFST
*
   10    CONTINUE
*
*        Swap block with next one below
*
         IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
*
*           Current block either 1 by 1 or 2 by 2
*
            NBNEXT = 1
            IF( HERE+NBF+1.LE.N ) THEN
               IF( T( HERE+NBF+1, HERE+NBF ).NE.ZERO )
     $            NBNEXT = 2
            END IF
*
            LITRAF = ILNGTH(NBF+NBNEXT-1)
            LDTRAF = DLNGTH(NBF+NBNEXT-1)
            IF( CITRAF+LITRAF-1.GT.NITRAF ) THEN
               INFO = -6
               CALL XERBLA( 'BSTREXC', -INFO )
               RETURN
            END IF
            IF( CDTRAF+LDTRAF-1.GT.NDTRAF ) THEN
               INFO = -8
               CALL XERBLA( 'BSTREXC', -INFO )
               RETURN
            END IF
            CALL BSLAEXC( N, T, LDT, HERE, NBF, NBNEXT, ITRAF(CITRAF),
     $                    DTRAF(CDTRAF), WORK, INFO )
            IF( INFO.NE.0 ) THEN
               ILST = HERE
               NITRAF = CITRAF - 1
               NDTRAF = CDTRAF - 1
               RETURN
            END IF
            CITRAF = CITRAF + LITRAF
            CDTRAF = CDTRAF + LDTRAF
            HERE = HERE + NBNEXT
*
*           Test if 2 by 2 block breaks into two 1 by 1 blocks
*
            IF( NBF.EQ.2 ) THEN
               IF( T( HERE+1, HERE ).EQ.ZERO )
     $            NBF = 3
            END IF
*
         ELSE
*
*           Current block consists of two 1 by 1 blocks each of which
*           must be swapped individually
*
            NBNEXT = 1
            IF( HERE+3.LE.N ) THEN
               IF( T( HERE+3, HERE+2 ).NE.ZERO )
     $            NBNEXT = 2
            END IF
            LITRAF = ILNGTH(NBNEXT)
            LDTRAF = DLNGTH(NBNEXT)
            IF( CITRAF+LITRAF-1.GT.NITRAF ) THEN
               INFO = -6
               CALL XERBLA( 'BSTREXC', -INFO )
               RETURN
            END IF
            IF( CDTRAF+LDTRAF-1.GT.NDTRAF ) THEN
               INFO = -8
               CALL XERBLA( 'BSTREXC', -INFO )
               RETURN
            END IF
            CALL BSLAEXC( N, T, LDT, HERE+1, 1, NBNEXT, ITRAF(CITRAF),
     $                    DTRAF(CDTRAF), WORK, INFO )
            IF( INFO.NE.0 ) THEN
               ILST = HERE
               NITRAF = CITRAF - 1
               NDTRAF = CDTRAF - 1
               RETURN
            END IF
            CITRAF = CITRAF + LITRAF
            CDTRAF = CDTRAF + LDTRAF
*
            IF( NBNEXT.EQ.1 ) THEN
*
*              Swap two 1 by 1 blocks, no problems possible
*
               LITRAF = ILNGTH(1)
               LDTRAF = DLNGTH(1)
               IF( CITRAF+LITRAF-1.GT.NITRAF ) THEN
                  INFO = -6
                  CALL XERBLA( 'BSTREXC', -INFO )
                  RETURN
               END IF
               IF( CDTRAF+LDTRAF-1.GT.NDTRAF ) THEN
                  INFO = -8
                  CALL XERBLA( 'BSTREXC', -INFO )
                  RETURN
               END IF
               CALL BSLAEXC( N, T, LDT, HERE, 1, NBNEXT, ITRAF(CITRAF),
     $                       DTRAF(CDTRAF), WORK, INFO )
               CITRAF = CITRAF + LITRAF
               CDTRAF = CDTRAF + LDTRAF
               HERE = HERE + 1
            ELSE
*
*              Recompute NBNEXT in case 2 by 2 split
*
               IF( T( HERE+2, HERE+1 ).EQ.ZERO )
     $            NBNEXT = 1
               IF( NBNEXT.EQ.2 ) THEN
*
*                 2 by 2 Block did not split
*
                  LITRAF = ILNGTH(2)
                  LDTRAF = DLNGTH(2)
                  IF( CITRAF+LITRAF-1.GT.NITRAF ) THEN
                     INFO = -6
                     CALL XERBLA( 'BSTREXC', -INFO )
                     RETURN
                  END IF
                  IF( CDTRAF+LDTRAF-1.GT.NDTRAF ) THEN
                     INFO = -8
                     CALL XERBLA( 'BSTREXC', -INFO )
                     RETURN
                  END IF
                  CALL BSLAEXC( N, T, LDT, HERE, 1, NBNEXT,
     $                          ITRAF(CITRAF), DTRAF(CDTRAF), WORK,
     $                          INFO )
                  IF( INFO.NE.0 ) THEN
                     INFO = 2
                     ILST = HERE
                     NITRAF = CITRAF - 1
                     NDTRAF = CDTRAF - 1
                     RETURN
                  END IF
                  CITRAF = CITRAF + LITRAF
                  CDTRAF = CDTRAF + LDTRAF
                  HERE = HERE + 2
               ELSE
*
*                 2 by 2 Block did split
*
                  LITRAF = ILNGTH(1)
                  LDTRAF = DLNGTH(1)
                  IF( CITRAF+2*LITRAF-1.GT.NITRAF ) THEN
                     INFO = -6
                     CALL XERBLA( 'BSTREXC', -INFO )
                     RETURN
                  END IF
                  IF( CDTRAF+2*LDTRAF-1.GT.NDTRAF ) THEN
                     INFO = -8
                     CALL XERBLA( 'BSTREXC', -INFO )
                     RETURN
                  END IF
                  CALL BSLAEXC( N, T, LDT, HERE, 1, 1, ITRAF(CITRAF),
     $                          DTRAF(CDTRAF), WORK, INFO )
                  CITRAF = CITRAF + LITRAF
                  CDTRAF = CDTRAF + LDTRAF
                  CALL BSLAEXC( N, T, LDT, HERE+1, 1, 1, ITRAF(CITRAF),
     $                          DTRAF(CDTRAF), WORK, INFO )
                  CITRAF = CITRAF + LITRAF
                  CDTRAF = CDTRAF + LDTRAF
                  HERE = HERE + 2
               END IF
            END IF
         END IF
         IF( HERE.LT.ILST )
     $      GO TO 10
*
      ELSE
*
         HERE = IFST
   20    CONTINUE
*
*        Swap block with next one above
*
         IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
*
*           Current block either 1 by 1 or 2 by 2
*
            NBNEXT = 1
            IF( HERE.GE.3 ) THEN
               IF( T( HERE-1, HERE-2 ).NE.ZERO )
     $            NBNEXT = 2
            END IF
*
            LITRAF = ILNGTH(NBF+NBNEXT-1)
            LDTRAF = DLNGTH(NBF+NBNEXT-1)
            IF( CITRAF+LITRAF-1.GT.NITRAF ) THEN
               INFO = -6
               CALL XERBLA( 'BSTREXC', -INFO )
               RETURN
            END IF
            IF( CDTRAF+LDTRAF-1.GT.NDTRAF ) THEN
               INFO = -8
               CALL XERBLA( 'BSTREXC', -INFO )
               RETURN
            END IF
            CALL BSLAEXC( N, T, LDT, HERE-NBNEXT, NBNEXT, NBF,
     $                    ITRAF(CITRAF), DTRAF(CDTRAF), WORK, INFO )
            IF( INFO.NE.0 ) THEN
               ILST = HERE
               NITRAF = CITRAF - 1
               NDTRAF = CDTRAF - 1
               RETURN
            END IF
            CITRAF = CITRAF + LITRAF
            CDTRAF = CDTRAF + LDTRAF
            HERE = HERE - NBNEXT
*
*           Test if 2 by 2 block breaks into two 1 by 1 blocks
*
            IF( NBF.EQ.2 ) THEN
               IF( T( HERE+1, HERE ).EQ.ZERO )
     $            NBF = 3
            END IF
*
         ELSE
*
*           Current block consists of two 1 by 1 blocks each of which
*           must be swapped individually
*
            NBNEXT = 1
            IF( HERE.GE.3 ) THEN
               IF( T( HERE-1, HERE-2 ).NE.ZERO )
     $            NBNEXT = 2
            END IF
            LITRAF = ILNGTH(NBNEXT)
            LDTRAF = DLNGTH(NBNEXT)
            IF( CITRAF+LITRAF-1.GT.NITRAF ) THEN
               INFO = -6
               CALL XERBLA( 'BSTREXC', -INFO )
               RETURN
            END IF
            IF( CDTRAF+LDTRAF-1.GT.NDTRAF ) THEN
               INFO = -8
               CALL XERBLA( 'BSTREXC', -INFO )
               RETURN
            END IF
            CALL BSLAEXC( N, T, LDT, HERE-NBNEXT, NBNEXT, 1,
     $                    ITRAF(CITRAF), DTRAF(CDTRAF), WORK, INFO )
            IF( INFO.NE.0 ) THEN
               ILST = HERE
               NITRAF = CITRAF - 1
               NDTRAF = CDTRAF - 1
               RETURN
            END IF
            CITRAF = CITRAF + LITRAF
            CDTRAF = CDTRAF + LDTRAF
*
            IF( NBNEXT.EQ.1 ) THEN
*
*              Swap two 1 by 1 blocks, no problems possible
*
               LITRAF = ILNGTH(1)
               LDTRAF = DLNGTH(1)
               IF( CITRAF+LITRAF-1.GT.NITRAF ) THEN
                  INFO = -6
                  CALL XERBLA( 'BSTREXC', -INFO )
                  RETURN
               END IF
               IF( CDTRAF+LDTRAF-1.GT.NDTRAF ) THEN
                  INFO = -8
                  CALL XERBLA( 'BSTREXC', -INFO )
                  RETURN
               END IF
               CALL BSLAEXC( N, T, LDT, HERE, NBNEXT, 1, ITRAF(CITRAF),
     $                       DTRAF(CDTRAF), WORK, INFO )
               CITRAF = CITRAF + LITRAF
               CDTRAF = CDTRAF + LDTRAF
               HERE = HERE - 1
            ELSE
*
*              Recompute NBNEXT in case 2 by 2 split
*
               IF( T( HERE, HERE-1 ).EQ.ZERO )
     $            NBNEXT = 1
               IF( NBNEXT.EQ.2 ) THEN
*
*                 2 by 2 Block did not split
*
                  LITRAF = ILNGTH(2)
                  LDTRAF = DLNGTH(2)
                  IF( CITRAF+LITRAF-1.GT.NITRAF ) THEN
                     INFO = -6
                     CALL XERBLA( 'BSTREXC', -INFO )
                     RETURN
                  END IF
                  IF( CDTRAF+LDTRAF-1.GT.NDTRAF ) THEN
                     INFO = -8
                     CALL XERBLA( 'BSTREXC', -INFO )
                     RETURN
                  END IF
                  CALL BSLAEXC( N, T, LDT, HERE-1, 2, 1, ITRAF(CITRAF),
     $                          DTRAF(CDTRAF), WORK, INFO )
                  IF( INFO.NE.0 ) THEN
                     INFO = 2
                     ILST = HERE
                     NITRAF = CITRAF - 1
                     NDTRAF = CDTRAF - 1
                     RETURN
                  END IF
                  CITRAF = CITRAF + LITRAF
                  CDTRAF = CDTRAF + LDTRAF
                  HERE = HERE - 2
               ELSE
*
*                 2 by 2 Block did split
*
                  LITRAF = ILNGTH(1)
                  LDTRAF = DLNGTH(1)
                  IF( CITRAF+2*LITRAF-1.GT.NITRAF ) THEN
                     INFO = -6
                     CALL XERBLA( 'BSTREXC', -INFO )
                     RETURN
                  END IF
                  IF( CDTRAF+2*LDTRAF-1.GT.NDTRAF ) THEN
                     INFO = -8
                     CALL XERBLA( 'BSTREXC', -INFO )
                     RETURN
                  END IF
                  CALL BSLAEXC( N, T, LDT, HERE, 1, 1, ITRAF(CITRAF),
     $                          DTRAF(CDTRAF), WORK, INFO )
                  CITRAF = CITRAF + LITRAF
                  CDTRAF = CDTRAF + LDTRAF
                  CALL BSLAEXC( N, T, LDT, HERE-1, 1, 1, ITRAF(CITRAF),
     $                          DTRAF(CDTRAF), WORK, INFO )
                  CITRAF = CITRAF + LITRAF
                  CDTRAF = CDTRAF + LDTRAF
                  HERE = HERE - 2
               END IF
            END IF
         END IF
         IF( HERE.GT.ILST )
     $      GO TO 20
      END IF
      ILST = HERE
      NITRAF = CITRAF - 1
      NDTRAF = CDTRAF - 1
*
      RETURN
*
*     End of BSTREXC
*
      END