1      SUBROUTINE TRED3(N,NV,A,D,E,E2)
2C***BEGIN PROLOGUE  TRED3
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  Reduce real symmetric matrix stored in packed form to
9C            symmetric tridiagonal matrix using orthogonal
10C            transformations.
11C***DESCRIPTION
12C
13C     This subroutine is a translation of the ALGOL procedure TRED3,
14C     NUM. MATH. 11, 181-195(1968) by Martin, Reinsch, and Wilkinson.
15C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
16C
17C     This subroutine reduces a REAL SYMMETRIC matrix, stored as
18C     a one-dimensional array, to a symmetric tridiagonal matrix
19C     using orthogonal similarity transformations.
20C
21C     On Input
22C
23C        n is the order of the matrix.
24C
25C        NV must be set to the dimension of the array parameter A
26C          as declared in the calling program dimension statement.
27C
28C        A contains the lower triangle of the real symmetric
29C          input matrix, stored row-wise as a one-dimensional
30C          array, in its first N*(N+1)/2 positions.
31C
32C     On Output
33C
34C        A contains information about the orthogonal
35C          transformations used in the reduction.
36C
37C        D contains the diagonal elements of the tridiagonal matrix.
38C
39C        E contains the subdiagonal elements of the tridiagonal
40C          matrix in its last N-1 positions.  E(1) is set to zero.
41C
42C        E2 contains the squares of the corresponding elements of E.
43C          E2 may coincide with E if the squares are not needed.
44C
45C     Questions and comments should be directed to B. S. Garbow,
46C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
47C     ------------------------------------------------------------------
48C***REFERENCES  B. T. SMITH, J. M. BOYLE, J. J. DONGARRA, B. S. GARBOW,
49C                 Y. IKEBE, V. C. KLEMA, C. B. MOLER, *MATRIX EIGEN-
50C                 SYSTEM ROUTINES - EISPACK GUIDE*, SPRINGER-VERLAG,
51C                 1976.
52C***ROUTINES CALLED  (NONE)
53C***END PROLOGUE  TRED3
54C
55      INTEGER I,J,K,L,N,II,IZ,JK,NV
56      REAL A(NV),D(N),E(N),E2(N)
57      REAL F,G,H,HH,SCALE
58C
59C     .......... FOR I=N STEP -1 UNTIL 1 DO -- ..........
60C***FIRST EXECUTABLE STATEMENT  TRED3
61      DO  300 II = 1, N
62         I = N + 1 - II
63         L = I - 1
64         IZ = (I * L) / 2
65         H = 0.0E0
66         SCALE = 0.0E0
67         IF (L .LT. 1) GO TO 130
68C     .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) ..........
69         DO 120 K = 1, L
70            IZ = IZ + 1
71            D(K) = A(IZ)
72            SCALE = SCALE + ABS(D(K))
73  120    CONTINUE
74C
75         IF (SCALE .NE. 0.0E0) GO TO 140
76  130    E(I) = 0.0E0
77         E2(I) = 0.0E0
78         GO TO 290
79C
80  140    DO 150 K = 1, L
81            D(K) = D(K) / SCALE
82            H = H + D(K) * D(K)
83  150    CONTINUE
84C
85         E2(I) = SCALE * SCALE * H
86         F = D(L)
87         G = -SIGN(SQRT(H),F)
88         E(I) = SCALE * G
89         H = H - F * G
90         D(L) = F - G
91         A(IZ) = SCALE * D(L)
92         IF (L .EQ. 1) GO TO 290
93         F = 0.0E0
94C
95         DO 240 J = 1, L
96            G = 0.0E0
97            JK = (J * (J-1)) / 2
98C     .......... FORM ELEMENT OF A*U ..........
99            DO 180 K = 1, L
100               JK = JK + 1
101               IF (K .GT. J) JK = JK + K - 2
102               G = G + A(JK) * D(K)
103  180       CONTINUE
104C     .......... FORM ELEMENT OF P ..........
105            E(J) = G / H
106            F = F + E(J) * D(J)
107  240    CONTINUE
108C
109         HH = F / (H + H)
110         JK = 0
111C     .......... FORM REDUCED A ..........
112         DO 260 J = 1, L
113            F = D(J)
114            G = E(J) - HH * F
115            E(J) = G
116C
117            DO 260 K = 1, J
118               JK = JK + 1
119               A(JK) = A(JK) - F * E(K) - G * D(K)
120  260    CONTINUE
121C
122  290    D(I) = A(IZ+1)
123         A(IZ+1) = SCALE * SQRT(H)
124  300 CONTINUE
125C
126      RETURN
127      END
128