*  =====================================================================
*
*  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, K, L, LDQ, LDT, N, N1, N2
      REAL           Q( LDQ, * ), T( LDT, * )
      REAL           ZERO, HALF, ONE
      REAL           BETA, CS, EPS, SCALE, SN, T11, T12, T22, T33,
     $               TAU1, TAU2, TEMP, WI1, WI2, WR1, WR2
      REAL           U( 3 ), U2( 3 ), X( 4 )
      PARAMETER      ( ZERO = 0.0E+0, HALF = 0.5E+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 )
*
         T( J1, J1 ) = T22
         T( J2, J2 ) = T11
*
         IF( WANTQ ) THEN
*           Accumulate transformation in the matrix Q.
            CALL SROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
         END IF
*
      ELSE
*
*        Swapping involves at least one 2-by-2 block.
*
*        Solve T11*X - X*T22 = scale*T12 for X.
*
         CALL SLASY2( -1, N1, N2, T( J1, J1 ), LDT,
     $                T( J1, J1+N1 ), LDT, SCALE, X )
*
*        Swap the adjacent diagonal blocks.
*
         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, * )
*
         U( 1 ) = X( 2 )
         U( 2 ) = X( 1 )
         U( 3 ) = SCALE
         CALL SLARFG3( U )
         TAU1 = U( 1 )
         U( 1 ) = U( 3 )
         U( 3 ) = ONE
         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 )
*
         U( 1 ) = -X( 1 )
         U( 2 ) = -X( 2 )
         U( 3 ) = SCALE
         CALL SLARFG3( U )
         TAU1 = U( 1 )
         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 )
*
         U( 1 ) = -X( 1 )
         U( 2 ) = -X( 2 )
         U( 3 ) = SCALE
         CALL SLARFG3( U )
         TAU1 = U( 1 )
         U( 1 ) = ONE
*
         TEMP = -TAU1*( X( 3 )+U( 2 )*X( 4 ) )
         U2( 1 ) = -TEMP*U( 2 ) - X( 4 )
         U2( 2 ) = -TEMP*U( 3 )
         U2( 3 ) = SCALE
         CALL SLARFG3( U2 )
         TAU2 = U2( 1 )
         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 .AND. T(J1,J2).NE.ZERO.AND.T(J2,J1).NE.ZERO ) THEN
         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 )
            IF ( CS.NE.ONE .OR. SN.NE.ZERO ) THEN
               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 )
               IF( WANTQ )
     $           CALL SROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
            END IF
         END IF
*
         IF( N1.EQ.2 ) THEN
*           Standardize new 2-by-2 block T22
            J3 = J1 + N2
            J4 = J3 + 1
            CALL SLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ),
     $                   T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN )
            IF ( CS.NE.ONE .OR. SN.NE.ZERO ) THEN
               IF( J4.LT.N )
     $            CALL SROT( N-J3-1, T( J3, J3+2 ), LDT,
     $                       T( J4, J3+2 ), LDT, CS, SN )
               CALL SROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN )
               IF( WANTQ )
     $           CALL SROT( N, Q( 1, J3 ), 1, Q( 1, J3+1 ), 1, CS, SN )
            END IF
         END IF
*
      END IF
      RETURN
*     End of SLAEXC
      END

      SUBROUTINE SLARFG3( U )
      IMPLICIT NONE
      REAL*4 ONE, BETA, SCL, U( 3 )
      PARAMETER ( ONE = 1.0E+0 )
      BETA = SIGN(SQRT(U(1)*U(1) + U(2)*U(2) + U(3)*U(3)),U(1))
      U(1) = U(1) + BETA
      SCL  = ONE / U(1)
      U(2) = U(2) * SCL
      U(3) = U(3) * SCL
      U(1) = ( U(1) / BETA )
      RETURN
*     End of SLARFG3
      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

* =====================================================================
*
*  Purpose
*  =======
*
*  SLASY2 solves for the N1 by N2 matrix X, 1 <= N1,N2 <= 2, in
*
*         op(TL)*X + ISGN*X*op(TR) = SCALE*B,
*
*  where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or
*  -1.  op(T) = T or T**T, where T**T denotes the transpose of T.
*
*  Arguments
*  =========
*
*  ISGN    (input) INTEGER
*          On entry, ISGN specifies the sign of the equation
*          as described before. ISGN may only be 1 or -1.
*
*  N1      (input) INTEGER
*          On entry, N1 specifies the order of matrix TL.
*          N1 may only be 0, 1 or 2.
*
*  N2      (input) INTEGER
*          On entry, N2 specifies the order of matrix TR.
*          N2 may only be 0, 1 or 2.
*
*  TL      (input) REAL array, dimension (LDTL,2)
*          On entry, TL contains an N1 by N1 matrix.
*
*  LDTL    (input) INTEGER
*          The leading dimension of the matrix TL. LDTL >= max(1,N1).
*
*  TR      (input) REAL array, dimension (LDTR,2)
*          On entry, TR contains an N2 by N2 matrix.
*
*  LDTR    (input) INTEGER
*          The leading dimension of the matrix TR. LDTR >= max(1,N2).
*
*  B       (input) REAL array, dimension (LDB,2)
*          On entry, the N1 by N2 matrix B contains the right-hand
*          side of the equation.
*
*  LDB     (input) INTEGER
*          The leading dimension of the matrix B. LDB >= max(1,N1).
*
*  SCALE   (output) REAL
*          On exit, SCALE contains the scale factor. SCALE is chosen
*          less than or equal to 1 to prevent the solution overflowing.
*
*  X       (output) REAL array, dimension (LDX,2)
*          On exit, X contains the N1 by N2 solution.
*
*  LDX     (input) INTEGER
*          The leading dimension of the matrix X. LDX >= max(1,N1).
*
*  INFO    (output) INTEGER
*          On exit, INFO is set to
*             0: successful exit.
*             1: TL and TR have too close eigenvalues, so TL or
*                TR is perturbed to get a nonsingular equation.
*          NOTE: In the interests of speed, this routine does not
*                check the inputs for errors.
*
* =====================================================================
*
      SUBROUTINE SLASY2( ISGN, N1, N2, T, LDT, B, LDB, SCALE, X )
      IMPLICIT NONE  
      INTEGER          INFO, ISGN, LDB, LDT, N1, N2
      REAL             B( LDB, * ), T( LDT, * ), X( * )
      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, GAM, L21, SGN, SMIN, SMLNUM, TEMP, TII,
     $                 U11, UINV11, U12, U22, SCALE, XMAX
      LOGICAL          BSWPIV( 4 )
      INTEGER          JPIV( 4 ), IPIVOT( 4, 4 )
      REAL             BTMP( 4 ), T16( 4, 4 ), TMP( 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. /
*
*     Quick return if possible
*
      IF( N1.EQ.0 .OR. N2.EQ.0 )
     $   RETURN
*
*     Set constants to control overflow
*
      EPS = 2*EPSILON( ZERO )
      SMLNUM = TINY( ZERO ) / EPS
      SGN = ISGN
      SCALE = ONE
*
      K = N1 + N1 + N2 - 2
      GO TO ( 10, 20, 30, 50 )K
*
*     1 by 1: TL11*X + SGN*X*TR11 = B11
*
   10 CONTINUE
      SMIN = EPS*MAX( ABS( T( 1, 1 ) ), ABS( T( 2, 2 ) ) )
      U11 = T( 1, 1 ) + SGN*T( 2, 2 )
      IF( ABS( U11 ).LE.SMIN )
     $   U11 = SMIN
*
      GAM = ABS( B( 1, 1 ) )
      IF( SMLNUM*GAM.GT.ABS( U11 ) )
     $   SCALE = ONE / GAM
*
      X( 1 ) = ( B( 1, 1 )*SCALE ) / U11
      RETURN
*
*     1 by 2:
*     TL11*[X11 X12] + ISGN*[X11 X12]*op[TR11 TR12]  = [B11 B12]
*                                       [TR21 TR22]
*
   20 CONTINUE
*
      TMP( 1 ) = T( 1, 1 ) + SGN*T( 2, 2 )
      TMP( 4 ) = T( 1, 1 ) + SGN*T( 3, 3 )
      TMP( 2 ) = SGN*T( 2, 3 )
      TMP( 3 ) = SGN*T( 3, 2 )
      BTMP( 1 ) = B( 1, 1 )
      BTMP( 2 ) = B( 1, 2 )
      GO TO 40
*
*     2 by 1:
*          op[TL11 TL12]*[X11] + ISGN* [X11]*TR11  = [B11]
*            [TL21 TL22] [X21]         [X21]         [B21]
*
   30 CONTINUE
      TMP( 1 ) = T( 1, 1 ) + SGN*T( 3, 3 )
      TMP( 4 ) = T( 2, 2 ) + SGN*T( 3, 3 )
      TMP( 2 ) = T( 2, 1 )
      TMP( 3 ) = T( 1, 2 )
      BTMP( 1 ) = B( 1, 1 )
      BTMP( 2 ) = B( 2, 1 )
   40 CONTINUE
*
*     Solve 2 by 2 system using complete pivoting.
*     Set pivots less than SMIN to SMIN.
*
      SMIN = EPS*MAX( ABS(T(1, 1)), ABS(T(2, 2)), ABS(T(3, 3)),
     $                ABS(TMP(2)), ABS(TMP(3)) )
      IPIV = 1
      XMAX = ZERO
      DO 75 IP = 1, 4
         IF( ABS( TMP( IP ) ).GE.XMAX ) THEN
            XMAX = ABS( TMP( IP ) )
            IPIV = IP
         END IF
   75 CONTINUE
      U11 = TMP( IPIV )
      IF( ABS( U11 ).LE.SMIN )
     $   U11 = SMIN
      UINV11 = ONE / U11
      U12 = TMP( IPIVOT( 3, IPIV ) )
      L21 = TMP( IPIVOT( 2, IPIV ) ) * UINV11
      U22 = TMP( IPIVOT( 4, IPIV ) ) - U12*L21
      IF( ABS( U22 ).LE.SMIN )
     $   U22 = SMIN
      XSWAP = 0
      IF ( IPIV.GT.2 )
     $  XSWAP = 1
      IF( BSWPIV( IPIV ) ) THEN
         TEMP = BTMP( 2 )
         BTMP( 2 ) = BTMP( 1 ) - L21*TEMP
         BTMP( 1 ) = TEMP
      ELSE
         BTMP( 2 ) = BTMP( 2 ) - L21*BTMP( 1 )
      END IF
      IF( ( TWO*SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( U22 ) .OR.
     $    ( TWO*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
*
*     2 by 2:
*     op[TL11 TL12]*[X11 X12] +ISGN* [X11 X12]*op[TR11 TR12] = [B11 B12]
*       [TL21 TL22] [X21 X22]        [X21 X22]   [TR21 TR22]   [B21 B22]
*
*     Solve equivalent 4 by 4 system using complete pivoting.
*     Set pivots less than SMIN to SMIN.
*
   50 CONTINUE
      SMIN = MAX( ABS( T( 3, 3 ) ), ABS( T( 3, 4 ) ),
     $       ABS( T( 4, 3 ) ), ABS( T( 4, 4 ) ) )
      SMIN = MAX( SMIN, ABS( T( 1, 1 ) ), ABS( T( 1, 2 ) ),
     $       ABS( T( 2, 1 ) ), ABS( T( 2, 2 ) ) )
      SMIN = EPS*SMIN
      T16( 1, 1 ) = T( 1, 1 ) + SGN*T( 3, 3 )
      T16( 2, 2 ) = T( 2, 2 ) + SGN*T( 3, 3 )
      T16( 3, 3 ) = T( 1, 1 ) + SGN*T( 4, 4 )
      T16( 4, 4 ) = T( 2, 2 ) + SGN*T( 4, 4 )
      T16( 1, 2 ) = T( 1, 2 )
      T16( 2, 1 ) = T( 2, 1 )
      T16( 1, 3 ) = SGN*T( 4, 3 )
      T16( 3, 1 ) = SGN*T( 3, 4 )
      T16( 2, 3 ) = ZERO
      T16( 3, 2 ) = ZERO
*
      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 )
      BTMP( 1 ) = B( 1, 1 )
      BTMP( 2 ) = B( 2, 1 )
      BTMP( 3 ) = B( 1, 2 )
      BTMP( 4 ) = B( 2, 2 )
*
*     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( ( EIGHT*SMLNUM )*ABS( BTMP( 1 ) ).GT.ABS( T16( 1, 1 ) ) .OR.
     $    ( EIGHT*SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( T16( 2, 2 ) ) .OR.
     $    ( EIGHT*SMLNUM )*ABS( BTMP( 3 ) ).GT.ABS( T16( 3, 3 ) ) .OR.
     $    ( EIGHT*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

