1      SUBROUTINE HTRID3(NM,N,A,D,E,E2,TAU)
2C***BEGIN PROLOGUE  HTRID3
3C***DATE WRITTEN   760101   (YYMMDD)
4C***REVISION DATE  830518   (YYMMDD)
5C***CATEGORY NO.  D4C1B1
6C***KEYWORDS  EIGENVALUES,EIGENVECTORS,EISPACK
7C***AUTHOR  SMITH, B. T., ET AL.
8C***PURPOSE  Reduces complex Hermitian (packed) matrix to real
9C            symmetric tridiagonal matrix by unitary similarity
10C            transformations.
11C***DESCRIPTION
12C
13C     This subroutine is a translation of a complex analogue of
14C     the ALGOL procedure TRED3, NUM. MATH. 11, 181-195(1968)
15C     by Martin, Reinsch, and Wilkinson.
16C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
17C
18C     This subroutine reduces a COMPLEX HERMITIAN matrix, stored as
19C     a single square array, to a real symmetric tridiagonal matrix
20C     using unitary similarity transformations.
21C
22C     On INPUT
23C
24C        NM must be set to the row dimension of two-dimensional
25C          array parameters as declared in the calling program
26C          dimension statement.
27C
28C        N is the order of the matrix.
29C
30C        A contains the lower triangle of the complex hermitian input
31C          matrix.  The real parts of the matrix elements are stored
32C          in the full lower triangle of A, and the imaginary parts
33C          are stored in the transposed positions of the strict upper
34C          triangle of A.  No storage is required for the zero
35C          imaginary parts of the diagonal elements.
36C
37C     On OUTPUT
38C
39C        A contains information about the unitary transformations
40C          used in the reduction.
41C
42C        D contains the diagonal elements of the the tridiagonal matrix.
43C
44C        E contains the subdiagonal elements of the tridiagonal
45C          matrix in its last N-1 positions.  E(1) is set to zero.
46C
47C        E2 contains the squares of the corresponding elements of E.
48C          E2 may coincide with E if the squares are not needed.
49C
50C        TAU contains further information about the transformations.
51C
52C     Calls PYTHAG(A,B) for sqrt(A**2 + B**2).
53C
54C     Questions and comments should be directed to B. S. Garbow,
55C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
56C     ------------------------------------------------------------------
57C***REFERENCES  B. T. SMITH, J. M. BOYLE, J. J. DONGARRA, B. S. GARBOW,
58C                 Y. IKEBE, V. C. KLEMA, C. B. MOLER, *MATRIX EIGEN-
59C                 SYSTEM ROUTINES - EISPACK GUIDE*, SPRINGER-VERLAG,
60C                 1976.
61C***ROUTINES CALLED  PYTHAG
62C***END PROLOGUE  HTRID3
63C
64      INTEGER I,J,K,L,N,II,NM,JM1,JP1
65      REAL A(NM,N),D(N),E(N),E2(N),TAU(2,N)
66      REAL F,G,H,FI,GI,HH,SI,SCALE
67      REAL PYTHAG
68C
69C***FIRST EXECUTABLE STATEMENT  HTRID3
70      TAU(1,N) = 1.0E0
71      TAU(2,N) = 0.0E0
72C     .......... FOR I=N STEP -1 UNTIL 1 DO -- ..........
73      DO 300 II = 1, N
74         I = N + 1 - II
75         L = I - 1
76         H = 0.0E0
77         SCALE = 0.0E0
78         IF (L .LT. 1) GO TO 130
79C     .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) ..........
80         DO 120 K = 1, L
81  120    SCALE = SCALE + ABS(A(I,K)) + ABS(A(K,I))
82C
83         IF (SCALE .NE. 0.0E0) GO TO 140
84         TAU(1,L) = 1.0E0
85         TAU(2,L) = 0.0E0
86  130    E(I) = 0.0E0
87         E2(I) = 0.0E0
88         GO TO 290
89C
90  140    DO 150 K = 1, L
91            A(I,K) = A(I,K) / SCALE
92            A(K,I) = A(K,I) / SCALE
93            H = H + A(I,K) * A(I,K) + A(K,I) * A(K,I)
94  150    CONTINUE
95C
96         E2(I) = SCALE * SCALE * H
97         G = SQRT(H)
98         E(I) = SCALE * G
99         F = PYTHAG(A(I,L),A(L,I))
100C     .......... FORM NEXT DIAGONAL ELEMENT OF MATRIX T ..........
101         IF (F .EQ. 0.0E0) GO TO 160
102         TAU(1,L) = (A(L,I) * TAU(2,I) - A(I,L) * TAU(1,I)) / F
103         SI = (A(I,L) * TAU(2,I) + A(L,I) * TAU(1,I)) / F
104         H = H + F * G
105         G = 1.0E0 + G / F
106         A(I,L) = G * A(I,L)
107         A(L,I) = G * A(L,I)
108         IF (L .EQ. 1) GO TO 270
109         GO TO 170
110  160    TAU(1,L) = -TAU(1,I)
111         SI = TAU(2,I)
112         A(I,L) = G
113  170    F = 0.0E0
114C
115         DO 240 J = 1, L
116            G = 0.0E0
117            GI = 0.0E0
118            IF (J .EQ. 1) GO TO 190
119            JM1 = J - 1
120C     .......... FORM ELEMENT OF A*U ..........
121            DO 180 K = 1, JM1
122               G = G + A(J,K) * A(I,K) + A(K,J) * A(K,I)
123               GI = GI - A(J,K) * A(K,I) + A(K,J) * A(I,K)
124  180       CONTINUE
125C
126  190       G = G + A(J,J) * A(I,J)
127            GI = GI - A(J,J) * A(J,I)
128            JP1 = J + 1
129            IF (L .LT. JP1) GO TO 220
130C
131            DO 200 K = JP1, L
132               G = G + A(K,J) * A(I,K) - A(J,K) * A(K,I)
133               GI = GI - A(K,J) * A(K,I) - A(J,K) * A(I,K)
134  200       CONTINUE
135C     .......... FORM ELEMENT OF P ..........
136  220       E(J) = G / H
137            TAU(2,J) = GI / H
138            F = F + E(J) * A(I,J) - TAU(2,J) * A(J,I)
139  240    CONTINUE
140C
141         HH = F / (H + H)
142C     .......... FORM REDUCED A ..........
143         DO 260 J = 1, L
144            F = A(I,J)
145            G = E(J) - HH * F
146            E(J) = G
147            FI = -A(J,I)
148            GI = TAU(2,J) - HH * FI
149            TAU(2,J) = -GI
150            A(J,J) = A(J,J) - 2.0E0 * (F * G + FI * GI)
151            IF (J .EQ. 1) GO TO 260
152            JM1 = J - 1
153C
154            DO 250 K = 1, JM1
155               A(J,K) = A(J,K) - F * E(K) - G * A(I,K)
156     1                         + FI * TAU(2,K) + GI * A(K,I)
157               A(K,J) = A(K,J) - F * TAU(2,K) - G * A(K,I)
158     1                         - FI * E(K) - GI * A(I,K)
159  250       CONTINUE
160C
161  260    CONTINUE
162C
163  270    DO 280 K = 1, L
164            A(I,K) = SCALE * A(I,K)
165            A(K,I) = SCALE * A(K,I)
166  280    CONTINUE
167C
168         TAU(2,L) = -SI
169  290    D(I) = A(I,I)
170         A(I,I) = SCALE * SQRT(H)
171  300 CONTINUE
172C
173      RETURN
174      END
175