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
|
subroutine tql1(n,d,e,ierr)
c
integer i,j,l,m,n,ii,l1,l2,mml,ierr
double precision d(n),e(n)
double precision c,c2,c3,dl1,el1,f,g,h,p,r,s,s2,tst1,tst2,pythag
c
c this subroutine is a translation of the algol procedure tql1,
c num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
c wilkinson.
c handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
c
c this subroutine finds the eigenvalues of a symmetric
c tridiagonal matrix by the 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 e contains the subdiagonal elements of the input matrix
c in its last n-1 positions. e(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 e 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 e(i-1) = e(i)
c
f = 0.0d0
tst1 = 0.0d0
e(n) = 0.0d0
c
do 290 l = 1, n
j = 0
h = dabs(d(l)) + dabs(e(l))
if (tst1 .lt. h) tst1 = h
c .......... look for small sub-diagonal element ..........
do 110 m = l, n
tst2 = tst1 + dabs(e(m))
if (tst2 .eq. tst1) go to 120
c .......... e(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
l2 = l1 + 1
g = d(l)
p = (d(l1) - g) / (2.0d0 * e(l))
r = pythag(p,1.0d0)
d(l) = e(l) / (p + dsign(r,p))
d(l1) = e(l) * (p + dsign(r,p))
dl1 = d(l1)
h = g - d(l)
if (l2 .gt. n) go to 145
c
do 140 i = l2, n
140 d(i) = d(i) - h
c
145 f = f + h
c .......... ql transformation ..........
p = d(m)
c = 1.0d0
c2 = c
el1 = e(l1)
s = 0.0d0
mml = m - l
c .......... for i=m-1 step -1 until l do -- ..........
do 200 ii = 1, mml
c3 = c2
c2 = c
s2 = s
i = m - ii
g = c * e(i)
h = c * p
r = pythag(p,e(i))
e(i+1) = s * r
s = e(i) / r
c = p / r
p = c * d(i) - s * g
d(i+1) = h + s * (c * g + s * d(i))
200 continue
c
p = -s * s2 * c3 * el1 * e(l) / dl1
e(l) = s * p
d(l) = c * p
tst2 = tst1 + dabs(e(l))
if (tst2 .gt. tst1) 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
|
