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
|
subroutine tqlrat(n,d,e2,ierr)
c
integer i,j,l,m,n,ii,l1,mml,ierr
double precision d(n),e2(n)
double precision b,c,f,g,h,p,r,s,t,epslon,pythag
c
c this subroutine is a translation of the algol procedure tqlrat,
c algorithm 464, comm. acm 16, 689(1973) by reinsch.
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.
c
c d contains the diagonal elements of the input matrix.
c
c e2 contains the squares of the subdiagonal elements of the
c input matrix in its last n-1 positions. e2(1) is arbitrary.
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 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 for dsqrt(a*a + b*b) .
c
c questions and comments should be directed to burton s. garbow,
c mathematics and computer science div, argonne national laboratory
c
c this version dated august 1983.
c
c ------------------------------------------------------------------
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.0d0
t = 0.0d0
e2(n) = 0.0d0
c
do 290 l = 1, n
j = 0
h = dabs(d(l)) + dsqrt(e2(l))
if (t .gt. h) go to 105
t = h
b = epslon(t)
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 = dsqrt(e2(l))
g = d(l)
p = (d(l1) - g) / (2.0d0 * s)
r = pythag(p,1.0d0)
d(l) = s / (p + dsign(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.0d0) g = b
h = g
s = 0.0d0
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.0d0) 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.0d0) go to 210
if (dabs(e2(l)) .le. dabs(c/h)) go to 210
e2(l) = h * e2(l)
if (e2(l) .ne. 0.0d0) 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
|
