1*                  Hermitian matrix.  On entry, Z must contain the
2*                  unitary matrix used to reduce the original matrix
3*                  to tridiagonal form.
4*          = 'I':  Compute eigenvalues and eigenvectors of the
5*                  tridiagonal matrix.  Z is initialized to the identity
6*                  matrix.
7*
8*  N       (input) INTEGER
9*          The order of the matrix.  N >= 0.
10*
11*  D       (input/output) DOUBLE PRECISION array, dimension (N)
12*          On entry, the diagonal elements of the tridiagonal matrix.
13*          On exit, if INFO = 0, the eigenvalues in ascending order.
14*
15*  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
16*          On entry, the (n-1) subdiagonal elements of the tridiagonal
17*          matrix.
18*          On exit, E has been destroyed.
19*
20*  Z       (input/output) COMPLEX*16 array, dimension (LDZ, N)
21*          On entry, if  COMPZ = 'V', then Z contains the unitary
22*          matrix used in the reduction to tridiagonal form.
23*          On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
24*          orthonormal eigenvectors of the original Hermitian matrix,
25*          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
26*          of the symmetric tridiagonal matrix.
27*          If COMPZ = 'N', then Z is not referenced.
28*
29*  LDZ     (input) INTEGER
30*          The leading dimension of the array Z.  LDZ >= 1, and if
31*          eigenvectors are desired, then  LDZ >= max(1,N).
32*
33*  WORK    (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2))
34*          If COMPZ = 'N', then WORK is not referenced.
35*
36*  INFO    (output) INTEGER
37*          = 0:  successful exit
38*          < 0:  if INFO = -i, the i-th argument had an illegal value
39*          > 0:  the algorithm has failed to find all the eigenvalues in
40*                a total of 30*N iterations; if INFO = i, then i
41*                elements of E have not converged to zero; on exit, D
42*                and E contain the elements of a symmetric tridiagonal
43*                matrix which is unitarily similar to the original
44*                matrix.
45*
46*  =====================================================================
47*
48*     .. Parameters ..
49      DOUBLE PRECISION   ZERO, ONE, TWO, THREE
50      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
51     $                   THREE = 3.0D0 )
52      COMPLEX*16         CZERO, CONE
53      PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
54     $                   CONE = ( 1.0D0, 0.0D0 ) )
55      INTEGER            MAXIT, KTOTAL
56      PARAMETER          ( MAXIT = 30 )
57*     ..
58*     .. Local Scalars ..
59      INTEGER            I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND,
60     $                   LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1,
61     $                   NM1, NMAXIT
62      DOUBLE PRECISION   ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2,
63     $                   S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
64*     ..
65*     .. External Functions ..
66      LOGICAL            LSAME
67      DOUBLE PRECISION   DLAMCH, DLANST, DLAPY2
68      EXTERNAL           LSAME, DLAMCH, DLANST, DLAPY2
69*     ..
70*     .. External Subroutines ..
71      EXTERNAL           DLAE2, DLAEV2, DLARTG, DLASCL, DLASRT, XERBLA,
72     $                   ZLASET, ZLASR, ZSWAP
73*     ..
74*     .. Intrinsic Functions ..
75      INTRINSIC          ABS, MAX, SIGN, SQRT
76*     ..
77*     .. Executable Statements ..
78*
79*     Test the input parameters.
80*
81      INFO = 0
82      KTOTAL = 0
83*
84      IF( LSAME( COMPZ, 'N' ) ) THEN
85         ICOMPZ = 0
86      ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
87         ICOMPZ = 1
88      ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
89         ICOMPZ = 2
90      ELSE
91         ICOMPZ = -1
92      END IF
93      IF( ICOMPZ.LT.0 ) THEN
94         INFO = -1
95      ELSE IF( N.LT.0 ) THEN
96         INFO = -2
97      ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
98     $         N ) ) ) THEN
99         INFO = -6
100      END IF
101      IF( INFO.NE.0 ) THEN
102         CALL XERBLA( 'ZSTEQR', -INFO )
103         RETURN
104      END IF
105*
106*     Quick return if possible
107*
108      IF( N.EQ.0 )
109     $   RETURN
110*
111      IF( N.EQ.1 ) THEN
112         IF( ICOMPZ.EQ.2 )
113     $      Z( 1, 1 ) = CONE
114         RETURN
115      END IF
116*
117*     Determine the unit roundoff and over/underflow thresholds.
118*
119      EPS = DLAMCH( 'E' )
120      EPS2 = EPS**2
121      SAFMIN = DLAMCH( 'S' )
122      SAFMAX = ONE / SAFMIN
123      SSFMAX = SQRT( SAFMAX ) / THREE
124      SSFMIN = SQRT( SAFMIN ) / EPS2
125*
126*     Compute the eigenvalues and eigenvectors of the tridiagonal
127*     matrix.
128*
129      IF( ICOMPZ.EQ.2 )
130     $   CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
131*
132      NMAXIT = N*MAXIT
133      JTOT = 0
134*
135*     Determine where the matrix splits and choose QL or QR iteration
136*     for each block, according to whether top or bottom diagonal
137*     element is smaller.
138*
139      L1 = 1
140      NM1 = N - 1
141*
142   10 CONTINUE
143      IF( L1.GT.N )
144     $   GO TO 160
145      IF( L1.GT.1 )
146     $   E( L1-1 ) = ZERO
147      IF( L1.LE.NM1 ) THEN
148         DO 20 M = L1, NM1
149            TST = ABS( E( M ) )
150            IF( TST.EQ.ZERO )
151     $         GO TO 30
152            IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
153     $          1 ) ) ) )*EPS ) THEN
154               E( M ) = ZERO
155               GO TO 30
156            END IF
157   20    CONTINUE
158      END IF
159      M = N
160*
161   30 CONTINUE
162      L = L1
163      LSV = L
164      LEND = M
165      LENDSV = LEND
166      L1 = M + 1
167      IF( LEND.EQ.L )
168     $   GO TO 10
169*
170*     Scale submatrix in rows and columns L to LEND
171*
172      ANORM = DLANST( 'I', LEND-L+1, D( L ), E( L ) )
173      ISCALE = 0
174*      print*,'checking if scaling is needed. L = ', L
175*      print*,'                               LEND = ', LEND
176*      print*,'                               ANORM = ', ANORM
177      IF( ANORM.EQ.ZERO ) THEN
178         print*,'anorm = ', ANORM
179     $   GO TO 10
180      END IF
181      IF( ANORM.GT.SSFMAX ) THEN
182         ISCALE = 1
183         CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
184     $                INFO )
185         CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
186     $                INFO )
187         print*,'iscale = 1. anorm = ', ANORM
188      ELSE IF( ANORM.LT.SSFMIN ) THEN
189         ISCALE = 2
190         CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
191     $                INFO )
192         CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
193     $                INFO )
194         print*,'iscale = 2. anorm = ', ANORM
195      END IF
196*
197*     Choose between QL and QR iteration
198*
199      IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
200         LEND = LSV
201         L = LENDSV
202      END IF
203*
204      IF( .TRUE. ) THEN
205*
206*        QR Iteration
207*
208*        Look for small superdiagonal element.
209*
210   90    CONTINUE
211         IF( L.NE.LEND ) THEN
212            LENDP1 = LEND + 1
213            DO 100 M = L, LENDP1, -1
214               TST = ABS( E( M-1 ) )**2
215               IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+
216     $             SAFMIN )GO TO 110
217  100       CONTINUE
218         END IF
219*
220         M = LEND
221*
222  110    CONTINUE
223         IF( M.GT.LEND )
224     $      E( M-1 ) = ZERO
225         P = D( L )
226         IF( M.EQ.L )
227     $      GO TO 130
228*
229*        If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
230*        to compute its eigensystem.
231*
232         IF( M.EQ.L-1 ) THEN
233            IF( ICOMPZ.GT.0 ) THEN
234               CALL DLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S )
235               WORK( M ) = C
236               WORK( N-1+M ) = S
237               CALL ZLASR( 'R', 'V', 'F', N, 2, WORK( M ),
238     $                     WORK( N-1+M ), Z( 1, L-1 ), LDZ )
239            ELSE
240               CALL DLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 )
241            END IF
242            D( L-1 ) = RT1
243            D( L ) = RT2
244            E( L-1 ) = ZERO
245            L = L - 2
246            IF( L.GE.LEND )
247     $         GO TO 90
248            GO TO 140
249         END IF
250*
251         IF( JTOT.EQ.NMAXIT )
252     $      GO TO 140
253         JTOT = JTOT + 1
254*
255*        Form shift.
256*
257         G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) )
258         R = DLAPY2( G, ONE )
259         G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) )
260*
261         S = ONE
262         C = ONE
263         P = ZERO
264*
265*        Inner loop
266*
267         LM1 = L - 1
268         DO 120 I = M, LM1
269            F = S*E( I )
270            B = C*E( I )
271            CALL DLARTG( G, F, C, S, R )
272            IF( I.NE.M )
273     $         E( I-1 ) = R
274            G = D( I ) - P
275            R = ( D( I+1 )-G )*S + TWO*C*B
276            P = S*R
277            D( I ) = G + P
278            G = C*R - B
279*
280*           If eigenvectors are desired, then save rotations.
281*
282            IF( ICOMPZ.GT.0 ) THEN
283               WORK( I ) = C
284               WORK( N-1+I ) = S
285            END IF
286*
287  120    CONTINUE
288*
289*        If eigenvectors are desired, then apply saved rotations.
290*
291         IF( ICOMPZ.GT.0 ) THEN
292            MM = L - M + 1
293            CALL ZLASR( 'R', 'V', 'F', N, MM, WORK( M ), WORK( N-1+M ),
294     $                  Z( 1, M ), LDZ )
295            KTOTAL = KTOTAL + 1
296         END IF
297*
298         D( L ) = D( L ) - P
299         E( LM1 ) = G
300         GO TO 90
301*
302*        Eigenvalue found.
303*
304  130    CONTINUE
305         D( L ) = P
306*
307         L = L - 1
308         IF( L.GE.LEND )
309     $      GO TO 90
310         GO TO 140
311*
312      END IF
313*
314*     Undo scaling if necessary
315*
316  140 CONTINUE
317      IF( ISCALE.EQ.1 ) THEN
318         CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
319     $                D( LSV ), N, INFO )
320         CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ),
321     $                N, INFO )
322      ELSE IF( ISCALE.EQ.2 ) THEN
323         CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
324     $                D( LSV ), N, INFO )
325         CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ),
326     $                N, INFO )
327      END IF
328*
329*     Check for no convergence to an eigenvalue after a total
330*     of N*MAXIT iterations.
331*
332      IF( JTOT.EQ.NMAXIT ) THEN
333         DO 150 I = 1, N - 1
334            IF( E( I ).NE.ZERO )
335     $         INFO = INFO + 1
336  150    CONTINUE
337         RETURN
338      END IF
339      GO TO 10
340*
341*     Order eigenvalues and eigenvectors.
342*
343  160 CONTINUE
344c      print *, 'K_TOTAL = ', KTOTAL
345      IF( ICOMPZ.EQ.0 ) THEN
346*
347*        Use Quick Sort
348*
349         CALL DLASRT( 'I', N, D, INFO )
350*
351      ELSE
352*
353*        Use Selection Sort to minimize swaps of eigenvectors
354*
355         DO 180 II = 2, N
356            I = II - 1
357            K = I
358            P = D( I )
359            DO 170 J = II, N
360               IF( D( J ).LT.P ) THEN
361                  K = J
362                  P = D( J )
363               END IF
364  170       CONTINUE
365            IF( K.NE.I ) THEN
366               D( K ) = D( I )
367               D( I ) = P
368               CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
369            END IF
370  180    CONTINUE
371      END IF
372      RETURN
373*
374*     End of ZSTEQR
375*
376      END
377