1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
|
*DECK TQLRAT
SUBROUTINE TQLRAT (N, D, E2, IERR)
C***BEGIN PROLOGUE TQLRAT
C***PURPOSE Compute the eigenvalues of symmetric tridiagonal matrix
C using a rational variant of the QL method.
C***LIBRARY SLATEC (EISPACK)
C***CATEGORY D4A5, D4C2A
C***TYPE SINGLE PRECISION (TQLRAT-S)
C***KEYWORDS EIGENVALUES OF A SYMMETRIC TRIDIAGONAL MATRIX, EISPACK,
C QL METHOD
C***AUTHOR Smith, B. T., et al.
C***DESCRIPTION
C
C This subroutine is a translation of the ALGOL procedure TQLRAT.
C
C This subroutine finds the eigenvalues of a SYMMETRIC
C TRIDIAGONAL matrix by the rational QL method.
C
C On Input
C
C N is the order of the matrix. N is an INTEGER variable.
C
C D contains the diagonal elements of the symmetric tridiagonal
C matrix. D is a one-dimensional REAL array, dimensioned D(N).
C
C E2 contains the squares of the subdiagonal elements of the
C symmetric tridiagonal matrix in its last N-1 positions.
C E2(1) is arbitrary. E2 is a one-dimensional REAL array,
C dimensioned E2(N).
C
C On Output
C
C D contains the eigenvalues in ascending order. If an
C error exit is made, the eigenvalues are correct and
C ordered for indices 1, 2, ..., IERR-1, but may not be
C the smallest eigenvalues.
C
C E2 has been destroyed.
C
C IERR is an INTEGER flag set to
C Zero for normal return,
C J if the J-th eigenvalue has not been
C determined after 30 iterations.
C
C Calls PYTHAG(A,B) for sqrt(A**2 + B**2).
C
C Questions and comments should be directed to B. S. Garbow,
C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C ------------------------------------------------------------------
C
C***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
C Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
C system Routines - EISPACK Guide, Springer-Verlag,
C 1976.
C C. H. Reinsch, Eigenvalues of a real, symmetric, tri-
C diagonal matrix, Algorithm 464, Communications of the
C ACM 16, 11 (November 1973), pp. 689.
C***ROUTINES CALLED PYTHAG, R1MACH
C***REVISION HISTORY (YYMMDD)
C 760101 DATE WRITTEN
C 890831 Modified array declarations. (WRB)
C 890831 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE TQLRAT
C
INTEGER I,J,L,M,N,II,L1,MML,IERR
REAL D(*),E2(*)
REAL B,C,F,G,H,P,R,S,MACHEP
REAL PYTHAG
LOGICAL FIRST
C
SAVE FIRST, MACHEP
DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT TQLRAT
IF (FIRST) THEN
MACHEP = R1MACH(4)
ENDIF
FIRST = .FALSE.
C
IERR = 0
IF (N .EQ. 1) GO TO 1001
C
DO 100 I = 2, N
100 E2(I-1) = E2(I)
C
F = 0.0E0
B = 0.0E0
E2(N) = 0.0E0
C
DO 290 L = 1, N
J = 0
H = MACHEP * (ABS(D(L)) + SQRT(E2(L)))
IF (B .GT. H) GO TO 105
B = H
C = B * B
C .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT ..........
105 DO 110 M = L, N
IF (E2(M) .LE. C) GO TO 120
C .......... E2(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
C THROUGH THE BOTTOM OF THE LOOP ..........
110 CONTINUE
C
120 IF (M .EQ. L) GO TO 210
130 IF (J .EQ. 30) GO TO 1000
J = J + 1
C .......... FORM SHIFT ..........
L1 = L + 1
S = SQRT(E2(L))
G = D(L)
P = (D(L1) - G) / (2.0E0 * S)
R = PYTHAG(P,1.0E0)
D(L) = S / (P + SIGN(R,P))
H = G - D(L)
C
DO 140 I = L1, N
140 D(I) = D(I) - H
C
F = F + H
C .......... RATIONAL QL TRANSFORMATION ..........
G = D(M)
IF (G .EQ. 0.0E0) G = B
H = G
S = 0.0E0
MML = M - L
C .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........
DO 200 II = 1, MML
I = M - II
P = G * H
R = P + E2(I)
E2(I+1) = S * R
S = E2(I) / R
D(I+1) = H + S * (H + D(I))
G = D(I) - E2(I) / G
IF (G .EQ. 0.0E0) G = B
H = G * P / R
200 CONTINUE
C
E2(L) = S * G
D(L) = H
C .......... GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST ..........
IF (H .EQ. 0.0E0) GO TO 210
IF (ABS(E2(L)) .LE. ABS(C/H)) GO TO 210
E2(L) = H * E2(L)
IF (E2(L) .NE. 0.0E0) GO TO 130
210 P = D(L) + F
C .......... ORDER EIGENVALUES ..........
IF (L .EQ. 1) GO TO 250
C .......... FOR I=L STEP -1 UNTIL 2 DO -- ..........
DO 230 II = 2, L
I = L + 2 - II
IF (P .GE. D(I-1)) GO TO 270
D(I) = D(I-1)
230 CONTINUE
C
250 I = 1
270 D(I) = P
290 CONTINUE
C
GO TO 1001
C .......... SET ERROR -- NO CONVERGENCE TO AN
C EIGENVALUE AFTER 30 ITERATIONS ..........
1000 IERR = L
1001 RETURN
END
|
