|
|
subroutine tql1(n,d,e,ierr)
|
|
|
c
|
|
|
integer i,j,l,m,n,ii,l1,l2,mml,ierr
|
|
|
real d(n),e(n)
|
|
|
real 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 sqrt(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.0e0
|
|
|
tst1 = 0.0e0
|
|
|
e(n) = 0.0e0
|
|
|
c
|
|
|
do 290 l = 1, n
|
|
|
j = 0
|
|
|
h = abs(d(l)) + abs(e(l))
|
|
|
if (tst1 .lt. h) tst1 = h
|
|
|
c .......... look for small sub-diagonal element ..........
|
|
|
do 110 m = l, n
|
|
|
tst2 = tst1 + abs(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.0e0 * e(l))
|
|
|
r = pythag(p,1.0e0)
|
|
|
d(l) = e(l) / (p + sign(r,p))
|
|
|
d(l1) = e(l) * (p + sign(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.0e0
|
|
|
c2 = c
|
|
|
el1 = e(l1)
|
|
|
s = 0.0e0
|
|
|
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 + abs(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
|
|
|
|
