*  =====================================================================
*
*  Purpose
*  =======
*
*  SLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in
*  an upper quasi-triangular matrix T by an orthogonal similarity
*  transformation.
*
*  T must be in Schur canonical form, that is, block upper triangular
*  with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block
*  has its diagonal elemnts equal and its off-diagonal elements of
*  opposite sign.
*
*  Arguments
*  =========
*
*  WANTQ   (input) LOGICAL
*          = .TRUE. : accumulate the transformation in the matrix Q;
*          = .FALSE.: do not accumulate the transformation.
*
*  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
*          canonical form.
*          On exit, the updated matrix T, again in Schur canonical form.
*
*  LDT     (input)  INTEGER
*          The leading dimension of the array T. LDT >= max(1,N).
*
*  Q       (input/output) REAL array, dimension (LDQ,N)
*          On entry, if WANTQ is .TRUE., the orthogonal matrix Q.
*          On exit, if WANTQ is .TRUE., the updated matrix Q.
*          If WANTQ is .FALSE., Q is not referenced.
*
*  LDQ     (input) INTEGER
*          The leading dimension of the array Q.
*          LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N.
*
*  J1      (input) INTEGER
*          The index of the first row of the first block T11.
*
*  N1      (input) INTEGER
*          The order of the first block T11. N1 = 0, 1 or 2.
*
*  N2      (input) INTEGER
*          The order of the second block T22. N2 = 0, 1 or 2.
*
*  INFO    (output) INTEGER
*          = 0: successful exit
*          = 1: the transformed matrix T would be too far from Schur
*               form; the blocks are not swapped and T and Q are
*               unchanged.
*
*  =====================================================================
*
      SUBROUTINE SLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, INFO )
      IMPLICIT NONE
      LOGICAL        WANTQ
      INTEGER        INFO, J, J1, J2, J3, J4, J5, K, LDQ, LDT, N, N1, N2
      REAL           Q( LDQ, * ), T( LDT, * )
      REAL           ZERO, ONE
      REAL           BETA, EPS, SCALE, SCALE2, CS, SN, CTEMP, STEMP,
     $               T11, T12, T22, T33, SCL, TAU1, TAU2, TEMP,
     $               WI1, WI2, WR1, WR2
      REAL           B0, B1, SMIN, TMP( 4 )
      REAL           T16( 4, 4 ) 
      REAL           U( 3 ), U2( 3 ), X( 4 )
      PARAMETER      ( ZERO = 0.0E+0, ONE = 1.0E+0 )
*
      INFO = 0
*
*     Quick return if possible
*
      IF( N.EQ.0 .OR. N1.EQ.0 .OR. N2.EQ.0 )
     $   RETURN
      IF( J1+N1.GT.N )
     $   RETURN
*
      J2 = J1 + 1
      J3 = J1 + 2
      J4 = J1 + 3
*
      IF( N1.EQ.1 .AND. N2.EQ.1 ) THEN
*
*        Swap two 1-by-1 blocks.
*
         T11 = T( J1, J1 )
         T12 = T( J1, J2 )
         T22 = T( J2, J2 )
*
*        Determine the transformation to perform the interchange.
*
         CS = T12
         SN = T22-T11
         TEMP = ONE / SQRT( CS*CS+SN*SN )
         IF( ABS( CS ).GT.ABS( SN ) .AND. CS.LT.ZERO ) THEN
            TEMP = -TEMP
         END IF
         CS = CS * TEMP
         SN = SN * TEMP
*
*        Apply transformation to the matrix T.
*
*         IF( J3.LE.N )
*     $      CALL SROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT,
*     $                 CS, SN )
*         CALL SROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
         DO J = J3, N
            CTEMP = CS*T( J1, J ) + SN*T( J2, J )
            STEMP = CS*T( J2, J ) - SN*T( J1, J )
            T( J1, J ) = CTEMP
            T( J2, J ) = STEMP
         END DO
         DO J = 1, J1-1
            CTEMP = CS*T( J, J1 ) + SN*T( J, J2 )
            STEMP = CS*T( J, J2 ) - SN*T( J, J1 )
            T( J, J1 ) = CTEMP
            T( J, J2 ) = STEMP
         END DO
*
         T( J1, J1 ) = T22
         T( J2, J2 ) = T11
*
         IF( WANTQ ) THEN
*           Accumulate transformation in the matrix Q.
*DEC$ NOVECTOR
            DO J = 1, N
               CTEMP = CS*Q( J, J1 ) + SN*Q( J, J2 )
               STEMP = CS*Q( J, J2 ) - SN*Q( J, J1 )
               Q( J, J1 ) = CTEMP
               Q( J, J2 ) = STEMP
            END DO
         END IF
*
      ELSE
*
*        Swapping involves at least one 2-by-2 block.
*
*        Solve T11*X - X*T22 = scale*T12 for X.
*
         EPS = 2*EPSILON(ZERO)
         K = N1 + N1 + N2 - 3
         IF ( K.NE.3 ) THEN
            IF ( N1.EQ.1 ) THEN
                TMP( 1 ) =  T( J1, J1 ) - T( J2, J2 )
                TMP( 4 ) =  T( J1, J1 ) - T( J3, J3 )
                TMP( 2 ) = -T( J2, J3 )
                TMP( 3 ) = -T( J3, J2 )
                B0 = T( J1, J2 )
                B1 = T( J1, J3 )
            ELSE
                TMP( 1 ) = T( J1, J1 ) - T( J3, J3 )
                TMP( 4 ) = T( J2, J2 ) - T( J3, J3 )
                TMP( 2 ) = T( J2, J1 )
                TMP( 3 ) = T( J1, J2 )
                B0 = T( J1, J3 )
                B1 = T( J2, J3 )
            END IF
            SMIN = MAX( ABS(T(J1,J1)), ABS(T(J2,J2)), ABS(T(J3,J3)),
     $                  ABS(TMP(2)), ABS(TMP(3)) )
            SMIN = EPS*SMIN
            CALL SLASY12( TMP, B0, B1, SMIN, SCALE, X )
         ELSE
            T16( 1, 1 ) = T( J1, J1 ) - T( J3, J3 )
            T16( 2, 2 ) = T( J2, J2 ) - T( J3, J3 )
            T16( 3, 3 ) = T( J1, J1 ) - T( J4, J4 )
            T16( 4, 4 ) = T( J2, J2 ) - T( J4, J4 )
            T16( 1, 2 ) =  T( J1, J2 )
            T16( 2, 1 ) =  T( J2, J1 )
            T16( 1, 3 ) = -T( J4, J3 )
            T16( 3, 1 ) = -T( J3, J4 )
            T16( 2, 3 ) = ZERO
            T16( 3, 2 ) = ZERO
            TMP( 1 ) = T( J1, J3 )
            TMP( 2 ) = T( J2, J3 )
            TMP( 3 ) = T( J1, J4 )
            TMP( 4 ) = T( J2, J4 )
            SMIN = MAX( ABS(T(J1,J1)), ABS(T(J2,J2)), ABS(T(J3,J3)),
     $                  ABS(T(J4,J4)) )
            SMIN = MAX( SMIN, ABS( T(J1,J2) ), ABS( T(J2,J1) ),
     $                        ABS( T(J3,J4) ), ABS( T(J4,J3) ) )
            SMIN = EPS*SMIN
            CALL SLASY22( T16, 4, TMP, SMIN, SCALE, X )
         END IF
*
*        Swap the adjacent diagonal blocks.
*
         SCALE2 = SCALE*SCALE
         BETA = SIGN(SQRT(X(1)*X(1) + X(2)*X(2) + SCALE2),-X(3-N1))
         K = N1 + N1 + N2 - 3
         GO TO ( 10, 20, 30 )K
*
   10    CONTINUE
*
*        N1 = 1, N2 = 2: generate elementary reflector H so that:
*
*        ( scale, X11, X12 ) H = ( 0, 0, * )
*
         TAU1 = X(2) - BETA
         SCL  = ONE / TAU1
         U( 3 ) = ONE
         U( 2 ) =  X(1) * SCL
         U( 1 ) = SCALE * SCL
         TAU1 = -( TAU1 / BETA )
         T11 = T( J1, J1 )
*
*        Apply transformation to the entire matrix T.
*
         DO 15 J = J1, N
            TEMP = TAU1*(U(1)*T(J1, J) + U(2)*T(J2, J) + T(J3, J))
            T( J1, J ) = T( J1, J ) - TEMP*U(1)
            T( J2, J ) = T( J2, J ) - TEMP*U(2)
            T( J3, J ) = T( J3, J ) - TEMP
   15    CONTINUE
         CALL SLARFX3R( J2, U, TAU1, T(  1, J1 ), LDT )
*
         T( J3, J1 ) = ZERO
         T( J3, J2 ) = ZERO
         T( J3, J3 ) = T11
         GO TO 40
*
   20    CONTINUE
*
*        N1 = 2, N2 = 1: generate elementary reflector H so that:
*
*        H (  -X11 ) = ( * )
*          (  -X21 ) = ( 0 )
*          ( scale ) = ( 0 )
*
         TAU1 = -X(1) + BETA
         SCL  = ONE / TAU1
         U(2) = -X(2) * SCL
         U(3) = SCALE * SCL
         TAU1 = TAU1 / BETA
         U(1) = ONE
         T33 = T( J3, J3 )
*
*        Apply transformation to the entire matrix T.
*
         CALL SLARFX3R( J3, U, TAU1, T( 1, J1 ), LDT )
         DO 25 J = J2, N
            TEMP = TAU1*(U(1)*T(J1, J) + U(2)*T(J2, J) + U(3)*T(J3, J))
            T( J1, J ) = T( J1, J ) - TEMP*U(1)
            T( J2, J ) = T( J2, J ) - TEMP*U(2)
            T( J3, J ) = T( J3, J ) - TEMP*U(3)
   25    CONTINUE
*
         T( J1, J1 ) = T33
         T( J2, J1 ) = ZERO
         T( J3, J1 ) = ZERO
         GO TO 40
*
   30    CONTINUE
*
*        N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so
*        that:
*
*        H(2) H(1) (  -X11  -X12 ) = (  *  * )
*                  (  -X21  -X22 )   (  0  * )
*                  ( scale    0  )   (  0  0 )
*                  (    0  scale )   (  0  0 )
*
         TAU1 = -X(1) + BETA
         SCL  = ONE / TAU1
         U(2) = -X(2) * SCL
         U(3) = SCALE * SCL
         TAU1 = TAU1 / BETA
         U(1) = ONE
*
         TEMP = -TAU1*( X( 3 )+U( 2 )*X( 4 ) )
         U2( 1 ) = -TEMP*U( 2 ) - X( 4 )
         U2( 2 ) = -TEMP*U( 3 )
         BETA = SIGN(SQRT(U2(1)*U2(1) + U2(2)*U2(2) + SCALE2), U2(1))
         TAU2 = U2(1) + BETA
         SCL  = ONE / TAU2
         U2(2) = U2(2) * SCL
         U2(3) = SCALE * SCL
         TAU2 = TAU2 / BETA
         U2(1) = ONE
*
*        Apply transformation to the entire matrix T.
*
         DO 35 J = J1, N
            TEMP = TAU1*(U(1)*T(J1, J) + U(2)*T(J2, J) + U(3)*T(J3, J))
            T( J1, J ) = T( J1, J ) - TEMP*U(1)
            T( J2, J ) = T( J2, J ) - TEMP*U(2)
            T( J3, J ) = T( J3, J ) - TEMP*U(3)
   35    CONTINUE
         CALL SLARFX3R( J4, U, TAU1, T( 1, J1 ), LDT )
         DO 45 J = J1, N
            TEMP = TAU2*(U2(1)*T(J2, J)+U2(2)*T(J3, J)+U2(3)*T(J4, J))
            T( J2, J ) = T( J2, J ) - TEMP*U2(1)
            T( J3, J ) = T( J3, J ) - TEMP*U2(2)
            T( J4, J ) = T( J4, J ) - TEMP*U2(3)
   45    CONTINUE
         CALL SLARFX3R( J4, U2, TAU2, T( 1, J2 ), LDT )
*
         T( J3, J1 ) = ZERO
         T( J3, J2 ) = ZERO
         T( J4, J1 ) = ZERO
         T( J4, J2 ) = ZERO
*
   40    CONTINUE
*
         IF( WANTQ ) THEN
*           Accumulate transformation in the matrix Q.
            CALL SLARFX3R( N, U, TAU1, Q( 1, J1 ), LDQ )
            IF( K.EQ.3 )
     $         CALL SLARFX3R( N, U2, TAU2, Q( 1, J2 ), LDQ )
         END IF
*
         IF( N2.EQ.2 ) THEN
*           Standardize new 2-by-2 block T11
            CALL SLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ),
     $                   T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN )
*            CALL SROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT,
*     $                 CS, SN )
*            CALL SROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
            DO J = J3, N
               CTEMP = CS*T( J1, J ) + SN*T( J2, J )
               STEMP = CS*T( J2, J ) - SN*T( J1, J )
               T( J1, J ) = CTEMP
               T( J2, J ) = STEMP
            END DO
            DO J = 1, J1-1
               CTEMP = CS*T( J, J1 ) + SN*T( J, J2 )
               STEMP = CS*T( J, J2 ) - SN*T( J, J1 )
               T( J, J1 ) = CTEMP
               T( J, J2 ) = STEMP
            END DO
            IF( WANTQ )
     $         CALL SROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
         END IF
*
         IF( N1.EQ.2 ) THEN
*           Standardize new 2-by-2 block T22
            J3 = J1 + N2
            J4 = J3 + 1
            J5 = J3 + 2
            CALL SLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ),
     $                   T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN )
*            IF( J4.LT.N )
*     $         CALL SROT( N-J3-1, T( J3, J5 ), LDT, T( J4, J5 ),
*     $                    LDT, CS, SN )
*            CALL SROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN )
            DO J = J5, N
               CTEMP = CS*T( J3, J ) + SN*T( J4, J )
               STEMP = CS*T( J4, J ) - SN*T( J3, J )
               T( J3, J ) = CTEMP
               T( J4, J ) = STEMP
            END DO
            DO J = 1, J3-1
               CTEMP = CS*T( J, J3 ) + SN*T( J, J4 )
               STEMP = CS*T( J, J4 ) - SN*T( J, J3 )
               T( J, J3 ) = CTEMP
               T( J, J4 ) = STEMP
            END DO
            IF( WANTQ )
     $         CALL SROT( N, Q( 1, J3 ), 1, Q( 1, J3+1 ), 1, CS, SN )
         END IF
*
      END IF
      RETURN
*     End of SLAEXC
      END

*  =====================================================================
*
*  Purpose
*  =======
*
*  SLARFX applies a real elementary reflector H to a real m by n
*  matrix C, from either the left or the right. H is represented in the
*  form
*
*        H = I - tau * v * v**T
*
*  where tau is a real scalar and v is a real vector.
*
*  If tau = 0, then H is taken to be the unit matrix
*
*  Arguments
*  =========
*
*  SIDE    (input) CHARACTER*1
*          = 'L': form  H * C
*          = 'R': form  C * H
*
*  N       (input) INTEGER
*          side = 'R': The number of columns of the matrix C.
*          side = 'L': The number of rows of the matrix C.
*
*  V       (input) REAL array, dimension (3)
*          The vector v in the representation of H.
*
*  TAU     (input) REAL
*          The value tau in the representation of H.
*
*  C       (input/output) REAL array, dimension (LDC,N)
*          On entry, the m by n matrix C.
*          On exit, C is overwritten by the matrix H * C if SIDE = 'L',
*          or C * H if SIDE = 'R'.
*
*  LDC     (input) INTEGER
*          The leading dimension of the array C. LDA >= (1,M).
*
*  =====================================================================
*
      SUBROUTINE SLARFX3R( N, V, TAU, A, LDA )
      IMPLICIT NONE
      INTEGER    LDA, J, N
      REAL       A( LDA, * ), SUM, TAU, V( 3 ), V1, V2, V3, ZERO
      PARAMETER  ( ZERO = 0.0E+0 )
*
      IF( TAU.EQ.ZERO )
     $   RETURN
*     Form  C * H, where H has order 3.
      V1 = V( 1 )
      V2 = V( 2 )
      V3 = V( 3 )
*DEC$ NOVECTOR
      DO 260 J = 1, N
         SUM = TAU*(V1*A( J, 1 ) + V2*A( J, 2 ) + V3*A( J, 3 ))
         A( J, 1 ) = A( J, 1 ) - SUM*V1
         A( J, 2 ) = A( J, 2 ) - SUM*V2
         A( J, 3 ) = A( J, 3 ) - SUM*V3
  260 CONTINUE
      GO TO 410
  410 RETURN
*     End of SLARFX3
      END

      SUBROUTINE SLASY12( TMP, B0, B1, SMIN, SCALE, X )
      IMPLICIT NONE
      REAL               B0, B1, SCALE, TMP( * ), X( * )
      REAL               ZERO, ONE, HALF
      REAL               BTMP( 2 )
      INTEGER            I, II, IP, IPIV, J, K, XSWAP
      REAL               EPS, SMIN, SMLNUM, TEMP, TII,
     $                   U11, UINV11, L21, U12, U22, XMAX
      LOGICAL            BSWAP, BSWPIV( 4 )
      INTEGER            IPIVOT( 4, 4 )
      DATA               IPIVOT / 1,2,3,4, 2,1,4,3, 3,4,1,2, 4,3,2,1 /
      DATA               BSWPIV / .FALSE., .TRUE., .FALSE., .TRUE. /
      PARAMETER          ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0 )
*
*     Set constants to control overflow
*
      SCALE = ONE
      EPS = 2*EPSILON(ZERO)
      SMLNUM = TINY(ZERO) / EPS
      SMLNUM = SMLNUM + SMLNUM
*     Solve 2 by 2 system using complete pivoting.
*     Set pivots less than SMIN to SMIN.
*      SMIN = EPS*MAX( ABS(TL(1, 1)), ABS(TL(2, 2)), ABS(TL(3, 3)),
*     $                ABS(TMP(2)), ABS(TMP(3)) )
      BTMP( 1 ) = B0
      BTMP( 2 ) = B1
      IPIV = 1
      XMAX = ABS(TMP(1))
      DO IP = 2, 4
         IF (ABS(TMP(IP)).GT.XMAX) THEN
            XMAX = ABS(TMP(IP))
            IPIV = IP
         END IF
      END DO
      U11 = TMP( IPIV )
      IF( ABS( U11 ).LE.SMIN ) THEN
         U11 = SMIN
      END IF
      UINV11 = ONE / U11
      U12 = TMP( IPIVOT( 3, IPIV ) )
      L21 = TMP( IPIVOT( 2, IPIV ) ) * UINV11
      U22 = TMP( IPIVOT( 4, IPIV ) ) - U12*L21
      XSWAP = 0
      IF ( IPIV.GT.2 )
     $  XSWAP = 1
      IF( ABS( U22 ).LE.SMIN ) THEN
         U22 = SMIN
      END IF
      BSWAP = BSWPIV( IPIV )
      IF( BSWAP ) THEN
         TEMP = BTMP( 2 )
         BTMP( 2 ) = BTMP( 1 ) - L21*TEMP
         BTMP( 1 ) = TEMP
      ELSE
         BTMP( 2 ) = BTMP( 2 ) - L21*BTMP( 1 )
      END IF
      SCALE = ONE
      IF( ( SMLNUM*ABS( BTMP( 2 ) ).GT.ABS( U22 ) ) .OR.
     $    ( SMLNUM*ABS( BTMP( 1 ) ).GT.ABS( U11 ) ) ) THEN
         SCALE = HALF / MAX( ABS( BTMP( 1 ) ), ABS( BTMP( 2 ) ) )
         BTMP( 1 ) = BTMP( 1 )*SCALE
         BTMP( 2 ) = BTMP( 2 )*SCALE
      END IF
      BTMP( 2 ) = BTMP( 2 ) / U22
      BTMP( 1 ) = BTMP( 1 ) * UINV11 - ( U12 * UINV11 )*BTMP( 2 )
      X( 1 ) = BTMP( 1+XSWAP )
      X( 2 ) = BTMP( 2-XSWAP )
      RETURN
*     End of SLASY2
      END

      SUBROUTINE SLASY22( T16, LDT, BTMP, SMIN, SCALE, X )
      IMPLICIT NONE  
      INTEGER          LDT
      REAL             BTMP( 4 ), T16( LDT, * ), X( * )
      REAL             SCALE, SMIN
      REAL             ZERO, ONE, TWO, HALF, EIGHT, EIGHTH
      PARAMETER      ( ZERO = 0.0E+0,ONE = 1.0E+0,TWO = 2.0E+0 )
      PARAMETER      ( HALF = 0.5E+0,EIGHT = 8.0E+0,EIGHTH = 0.125E+0 )
      INTEGER          I, II, IP, IPIV, IPSV, J, JP, JPSV, K, XSWAP
      REAL             EPS, SMLNUM, TEMP, TII, XMAX
      INTEGER          JPIV( 4 )
*
*     Set constants to control overflow
*
      SCALE = ONE
      EPS = 2*EPSILON(ZERO)
      SMLNUM = EIGHT*TINY(ZERO) / EPS
*
*     Copy some values
*
      T16( 1, 4 ) = ZERO
      T16( 2, 4 ) = T16( 1, 3 )
      T16( 3, 4 ) = T16( 1, 2 )
      T16( 4, 1 ) = ZERO
      T16( 4, 2 ) = T16( 3, 1 )
      T16( 4, 3 ) = T16( 2, 1 )
*
*     Perform elimination
*
      DO 100 I = 1, 3
         XMAX = ZERO
         DO 70 IP = I, 4
            DO 60 JP = I, 4
               IF( ABS( T16( IP, JP ) ).GE.XMAX ) THEN
                  XMAX = ABS( T16( IP, JP ) )
                  IPSV = IP
                  JPSV = JP
               END IF
   60       CONTINUE
   70    CONTINUE
         IF( IPSV.NE.I ) THEN
*            CALL SSWAP( 4, T16( IPSV, 1 ), 4, T16( I, 1 ), 4 )
            DO II = 1,4
               TEMP = T16(IPSV,II)
               T16(IPSV,II) = T16(I,II)
               T16(I,II) = TEMP
            END DO
            TEMP = BTMP( I )
            BTMP( I ) = BTMP( IPSV )
            BTMP( IPSV ) = TEMP
         END IF
*         IF( JPSV.NE.I )
*     $      CALL SSWAP( 4, T16( 1, JPSV ), 1, T16( 1, I ), 1 )
         IF( JPSV.NE.I ) THEN
*            CALL SSWAP( 4, T16( 1, JPSV ), 1, T16( 1, I ), 1 )
            DO II = 1,4
               TEMP = T16(II,JPSV)
               T16(II,JPSV) = T16(II,I)
               T16(II,I) = TEMP
            END DO
         END IF
         JPIV( I ) = JPSV
         IF( ABS( T16( I, I ) ).LT.SMIN )
     $      T16( I, I ) = SMIN
         TII = ONE / T16( I, I )
         DO 90 J = I + 1, 4
            T16( J, I ) = T16( J, I ) * TII
            BTMP( J ) = BTMP( J ) - T16( J, I )*BTMP( I )
            DO 80 K = I + 1, 4
               T16( J, K ) = T16( J, K ) - T16( J, I )*T16( I, K )
   80       CONTINUE
   90    CONTINUE
  100 CONTINUE
      IF( ABS( T16( 4, 4 ) ).LT.SMIN )
     $   T16( 4, 4 ) = SMIN
      IF( ( SMLNUM )*ABS( BTMP( 1 ) ).GT.ABS( T16( 1, 1 ) ) .OR.
     $    ( SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( T16( 2, 2 ) ) .OR.
     $    ( SMLNUM )*ABS( BTMP( 3 ) ).GT.ABS( T16( 3, 3 ) ) .OR.
     $    ( SMLNUM )*ABS( BTMP( 4 ) ).GT.ABS( T16( 4, 4 ) ) ) THEN
         SCALE = EIGHTH / MAX( ABS( BTMP( 1 ) ), ABS( BTMP( 2 ) ),
     $                         ABS( BTMP( 3 ) ), ABS( BTMP( 4 ) ) )
         BTMP( 1 ) = BTMP( 1 )*SCALE
         BTMP( 2 ) = BTMP( 2 )*SCALE
         BTMP( 3 ) = BTMP( 3 )*SCALE
         BTMP( 4 ) = BTMP( 4 )*SCALE
      END IF
      DO 120 I = 1, 4
         K = 5 - I
         TEMP = ONE / T16( K, K )
         X( K ) = BTMP( K )*TEMP
         DO 110 J = K + 1, 4
            X( K ) = X( K ) - ( TEMP*T16( K, J ) )*X( J )
  110    CONTINUE
  120 CONTINUE
      DO 130 I = 1, 3
         IF( JPIV( 4-I ).NE.4-I ) THEN
            TEMP = X( 4-I )
            X( 4-I ) = X( JPIV( 4-I ) )
            X( JPIV( 4-I ) ) = TEMP
         END IF
  130 CONTINUE
      RETURN
*     End of SLASY2
      END

