```      SUBROUTINE SLASQ2( N, Z, INFO )
*
*  -- LAPACK routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     Modified to call SLAZQ3 in place of SLASQ3, 13 Feb 03, SJH.
*
*     .. Scalar Arguments ..
INTEGER            INFO, N
*     ..
*     .. Array Arguments ..
REAL               Z( * )
*     ..
*
*  Purpose
*  =======
*
*  SLASQ2 computes all the eigenvalues of the symmetric positive
*  definite tridiagonal matrix associated with the qd array Z to high
*  relative accuracy are computed to high relative accuracy, in the
*  absence of denormalization, underflow and overflow.
*
*  To see the relation of Z to the tridiagonal matrix, let L be a
*  unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
*  let U be an upper bidiagonal matrix with 1's above and diagonal
*  Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
*  symmetric tridiagonal to which it is similar.
*
*  Note : SLASQ2 defines a logical variable, IEEE, which is true
*  on machines which follow ieee-754 floating-point standard in their
*  handling of infinities and NaNs, and false otherwise. This variable
*  is passed to SLAZQ3.
*
*  Arguments
*  =========
*
*  N     (input) INTEGER
*        The number of rows and columns in the matrix. N >= 0.
*
*  Z     (workspace) REAL array, dimension (4*N)
*        On entry Z holds the qd array. On exit, entries 1 to N hold
*        the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
*        trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
*        N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
*        holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
*        shifts that failed.
*
*  INFO  (output) INTEGER
*        = 0: successful exit
*        < 0: if the i-th argument is a scalar and had an illegal
*             value, then INFO = -i, if the i-th argument is an
*             array and the j-entry had an illegal value, then
*             INFO = -(i*100+j)
*        > 0: the algorithm failed
*              = 1, a split was marked by a positive value in E
*              = 2, current block of Z not diagonalized after 30*N
*                   iterations (in inner while loop)
*              = 3, termination criterion of outer while loop not met
*                   (program created more than N unreduced blocks)
*
*  Further Details
*  ===============
*  Local Variables: I0:N0 defines a current unreduced segment of Z.
*  The shifts are accumulated in SIGMA. Iteration count is in ITER.
*  Ping-pong is controlled by PP (alternates between 0 and 1).
*
*  =====================================================================
*
*     .. Parameters ..
REAL               CBIAS
PARAMETER          ( CBIAS = 1.50E0 )
REAL               ZERO, HALF, ONE, TWO, FOUR, HUNDRD
PARAMETER          ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0,
\$                     TWO = 2.0E0, FOUR = 4.0E0, HUNDRD = 100.0E0 )
*     ..
*     .. Local Scalars ..
LOGICAL            IEEE
INTEGER            I0, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, K,
\$                   N0, NBIG, NDIV, NFAIL, PP, SPLT, TTYPE
REAL               D, DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, E,
\$                   EMAX, EMIN, EPS, OLDEMN, QMAX, QMIN, S, SAFMIN,
\$                   SIGMA, T, TAU, TEMP, TOL, TOL2, TRACE, ZMAX
*     ..
*     .. External Subroutines ..
EXTERNAL           SLAZQ3, SLASRT, XERBLA
*     ..
*     .. External Functions ..
INTEGER            ILAENV
REAL               SLAMCH
EXTERNAL           ILAENV, SLAMCH
*     ..
*     .. Intrinsic Functions ..
INTRINSIC          ABS, MAX, MIN, REAL, SQRT
*     ..
*     .. Executable Statements ..
*
*     Test the input arguments.
*     (in case SLASQ2 is not called by SLASQ1)
*
INFO = 0
EPS = SLAMCH( 'Precision' )
SAFMIN = SLAMCH( 'Safe minimum' )
TOL = EPS*HUNDRD
TOL2 = TOL**2
*
IF( N.LT.0 ) THEN
INFO = -1
CALL XERBLA( 'SLASQ2', 1 )
RETURN
ELSE IF( N.EQ.0 ) THEN
RETURN
ELSE IF( N.EQ.1 ) THEN
*
*        1-by-1 case.
*
IF( Z( 1 ).LT.ZERO ) THEN
INFO = -201
CALL XERBLA( 'SLASQ2', 2 )
END IF
RETURN
ELSE IF( N.EQ.2 ) THEN
*
*        2-by-2 case.
*
IF( Z( 2 ).LT.ZERO .OR. Z( 3 ).LT.ZERO ) THEN
INFO = -2
CALL XERBLA( 'SLASQ2', 2 )
RETURN
ELSE IF( Z( 3 ).GT.Z( 1 ) ) THEN
D = Z( 3 )
Z( 3 ) = Z( 1 )
Z( 1 ) = D
END IF
Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 )
IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN
T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) )
S = Z( 3 )*( Z( 2 ) / T )
IF( S.LE.T ) THEN
S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) )
ELSE
S = Z( 3 )*( Z( 2 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
END IF
T = Z( 1 ) + ( S+Z( 2 ) )
Z( 3 ) = Z( 3 )*( Z( 1 ) / T )
Z( 1 ) = T
END IF
Z( 2 ) = Z( 3 )
Z( 6 ) = Z( 2 ) + Z( 1 )
RETURN
END IF
*
*     Check for negative data and compute sums of q's and e's.
*
Z( 2*N ) = ZERO
EMIN = Z( 2 )
QMAX = ZERO
ZMAX = ZERO
D = ZERO
E = ZERO
*
DO 10 K = 1, 2*( N-1 ), 2
IF( Z( K ).LT.ZERO ) THEN
INFO = -( 200+K )
CALL XERBLA( 'SLASQ2', 2 )
RETURN
ELSE IF( Z( K+1 ).LT.ZERO ) THEN
INFO = -( 200+K+1 )
CALL XERBLA( 'SLASQ2', 2 )
RETURN
END IF
D = D + Z( K )
E = E + Z( K+1 )
QMAX = MAX( QMAX, Z( K ) )
EMIN = MIN( EMIN, Z( K+1 ) )
ZMAX = MAX( QMAX, ZMAX, Z( K+1 ) )
10 CONTINUE
IF( Z( 2*N-1 ).LT.ZERO ) THEN
INFO = -( 200+2*N-1 )
CALL XERBLA( 'SLASQ2', 2 )
RETURN
END IF
D = D + Z( 2*N-1 )
QMAX = MAX( QMAX, Z( 2*N-1 ) )
ZMAX = MAX( QMAX, ZMAX )
*
*     Check for diagonality.
*
IF( E.EQ.ZERO ) THEN
DO 20 K = 2, N
Z( K ) = Z( 2*K-1 )
20    CONTINUE
CALL SLASRT( 'D', N, Z, IINFO )
Z( 2*N-1 ) = D
RETURN
END IF
*
TRACE = D + E
*
*     Check for zero data.
*
IF( TRACE.EQ.ZERO ) THEN
Z( 2*N-1 ) = ZERO
RETURN
END IF
*
*     Check whether the machine is IEEE conformable.
*
IEEE = ILAENV( 10, 'SLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND.
\$       ILAENV( 11, 'SLASQ2', 'N', 1, 2, 3, 4 ).EQ.1
*
*     Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
*
DO 30 K = 2*N, 2, -2
Z( 2*K ) = ZERO
Z( 2*K-1 ) = Z( K )
Z( 2*K-2 ) = ZERO
Z( 2*K-3 ) = Z( K-1 )
30 CONTINUE
*
I0 = 1
N0 = N
*
*     Reverse the qd-array, if warranted.
*
IF( CBIAS*Z( 4*I0-3 ).LT.Z( 4*N0-3 ) ) THEN
IPN4 = 4*( I0+N0 )
DO 40 I4 = 4*I0, 2*( I0+N0-1 ), 4
TEMP = Z( I4-3 )
Z( I4-3 ) = Z( IPN4-I4-3 )
Z( IPN4-I4-3 ) = TEMP
TEMP = Z( I4-1 )
Z( I4-1 ) = Z( IPN4-I4-5 )
Z( IPN4-I4-5 ) = TEMP
40    CONTINUE
END IF
*
*     Initial split checking via dqd and Li's test.
*
PP = 0
*
DO 80 K = 1, 2
*
D = Z( 4*N0+PP-3 )
DO 50 I4 = 4*( N0-1 ) + PP, 4*I0 + PP, -4
IF( Z( I4-1 ).LE.TOL2*D ) THEN
Z( I4-1 ) = -ZERO
D = Z( I4-3 )
ELSE
D = Z( I4-3 )*( D / ( D+Z( I4-1 ) ) )
END IF
50    CONTINUE
*
*        dqd maps Z to ZZ plus Li's test.
*
EMIN = Z( 4*I0+PP+1 )
D = Z( 4*I0+PP-3 )
DO 60 I4 = 4*I0 + PP, 4*( N0-1 ) + PP, 4
Z( I4-2*PP-2 ) = D + Z( I4-1 )
IF( Z( I4-1 ).LE.TOL2*D ) THEN
Z( I4-1 ) = -ZERO
Z( I4-2*PP-2 ) = D
Z( I4-2*PP ) = ZERO
D = Z( I4+1 )
ELSE IF( SAFMIN*Z( I4+1 ).LT.Z( I4-2*PP-2 ) .AND.
\$               SAFMIN*Z( I4-2*PP-2 ).LT.Z( I4+1 ) ) THEN
TEMP = Z( I4+1 ) / Z( I4-2*PP-2 )
Z( I4-2*PP ) = Z( I4-1 )*TEMP
D = D*TEMP
ELSE
Z( I4-2*PP ) = Z( I4+1 )*( Z( I4-1 ) / Z( I4-2*PP-2 ) )
D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) )
END IF
EMIN = MIN( EMIN, Z( I4-2*PP ) )
60    CONTINUE
Z( 4*N0-PP-2 ) = D
*
*        Now find qmax.
*
QMAX = Z( 4*I0-PP-2 )
DO 70 I4 = 4*I0 - PP + 2, 4*N0 - PP - 2, 4
QMAX = MAX( QMAX, Z( I4 ) )
70    CONTINUE
*
*        Prepare for the next iteration on K.
*
PP = 1 - PP
80 CONTINUE
*
*     Initialise variables to pass to SLAZQ3
*
TTYPE = 0
DMIN1 = ZERO
DMIN2 = ZERO
DN    = ZERO
DN1   = ZERO
DN2   = ZERO
TAU   = ZERO
*
ITER = 2
NFAIL = 0
NDIV = 2*( N0-I0 )
*
DO 140 IWHILA = 1, N + 1
IF( N0.LT.1 )
\$      GO TO 150
*
*        While array unfinished do
*
*        E(N0) holds the value of SIGMA when submatrix in I0:N0
*        splits from the rest of the array, but is negated.
*
DESIG = ZERO
IF( N0.EQ.N ) THEN
SIGMA = ZERO
ELSE
SIGMA = -Z( 4*N0-1 )
END IF
IF( SIGMA.LT.ZERO ) THEN
INFO = 1
RETURN
END IF
*
*        Find last unreduced submatrix's top index I0, find QMAX and
*        EMIN. Find Gershgorin-type bound if Q's much greater than E's.
*
EMAX = ZERO
IF( N0.GT.I0 ) THEN
EMIN = ABS( Z( 4*N0-5 ) )
ELSE
EMIN = ZERO
END IF
QMIN = Z( 4*N0-3 )
QMAX = QMIN
DO 90 I4 = 4*N0, 8, -4
IF( Z( I4-5 ).LE.ZERO )
\$         GO TO 100
IF( QMIN.GE.FOUR*EMAX ) THEN
QMIN = MIN( QMIN, Z( I4-3 ) )
EMAX = MAX( EMAX, Z( I4-5 ) )
END IF
QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) )
EMIN = MIN( EMIN, Z( I4-5 ) )
90    CONTINUE
I4 = 4
*
100    CONTINUE
I0 = I4 / 4
*
*        Store EMIN for passing to SLAZQ3.
*
Z( 4*N0-1 ) = EMIN
*
*        Put -(initial shift) into DMIN.
*
DMIN = -MAX( ZERO, QMIN-TWO*SQRT( QMIN )*SQRT( EMAX ) )
*
*        Now I0:N0 is unreduced. PP = 0 for ping, PP = 1 for pong.
*
PP = 0
*
NBIG = 30*( N0-I0+1 )
DO 120 IWHILB = 1, NBIG
IF( I0.GT.N0 )
\$         GO TO 130
*
*           While submatrix unfinished take a good dqds step.
*
CALL SLAZQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL,
\$                   ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1,
\$                   DN2, TAU )
*
PP = 1 - PP
*
*           When EMIN is very small check for splits.
*
IF( PP.EQ.0 .AND. N0-I0.GE.3 ) THEN
IF( Z( 4*N0 ).LE.TOL2*QMAX .OR.
\$             Z( 4*N0-1 ).LE.TOL2*SIGMA ) THEN
SPLT = I0 - 1
QMAX = Z( 4*I0-3 )
EMIN = Z( 4*I0-1 )
OLDEMN = Z( 4*I0 )
DO 110 I4 = 4*I0, 4*( N0-3 ), 4
IF( Z( I4 ).LE.TOL2*Z( I4-3 ) .OR.
\$                   Z( I4-1 ).LE.TOL2*SIGMA ) THEN
Z( I4-1 ) = -SIGMA
SPLT = I4 / 4
QMAX = ZERO
EMIN = Z( I4+3 )
OLDEMN = Z( I4+4 )
ELSE
QMAX = MAX( QMAX, Z( I4+1 ) )
EMIN = MIN( EMIN, Z( I4-1 ) )
OLDEMN = MIN( OLDEMN, Z( I4 ) )
END IF
110             CONTINUE
Z( 4*N0-1 ) = EMIN
Z( 4*N0 ) = OLDEMN
I0 = SPLT + 1
END IF
END IF
*
120    CONTINUE
*
INFO = 2
RETURN
*
*        end IWHILB
*
130    CONTINUE
*
140 CONTINUE
*
INFO = 3
RETURN
*
*     end IWHILA
*
150 CONTINUE
*
*     Move q's to the front.
*
DO 160 K = 2, N
Z( K ) = Z( 4*K-3 )
160 CONTINUE
*
*     Sort and compute sum of eigenvalues.
*
CALL SLASRT( 'D', N, Z, IINFO )
*
E = ZERO
DO 170 K = N, 1, -1
E = E + Z( K )
170 CONTINUE
*
*     Store trace, sum(eigenvalues) and information on performance.
*
Z( 2*N+1 ) = TRACE
Z( 2*N+2 ) = E
Z( 2*N+3 ) = REAL( ITER )
Z( 2*N+4 ) = REAL( NDIV ) / REAL( N**2 )
Z( 2*N+5 ) = HUNDRD*NFAIL / REAL( ITER )
RETURN
*
*     End of SLASQ2
*
END

```