tqli.f90

1       SUBROUTINE tqli(d,e,n,np,z)
2       USE RAD_BIAS
3       IMPLICIT NONE
4       INTEGER, INTENT(IN) ::  n,np
5 !      REAL(KIND=LONG), INTENT(OUT) :: d(np),e(np),z(np,np)
6       REAL(KIND=LONG) :: d(np),e(np),z(np,np)
7 !    USES pythag
8       INTEGER i,iter,k,l,m
9       REAL(KIND=LONG) :: b,c,dd,f,g,p,r,s,pythag
10       do i=2,n
11         e(i-1)=e(i)
12       end do
13       e(n)=0.
14       do l=1,n
15         iter=0
16 1       do m=l,n-1
17           dd=abs(d(m))+abs(d(m+1))
18           if (abs(e(m))+dd.eq.dd) goto 2
19         end do
20         m=n
21 2       if(m.ne.l)then
22           if(iter.eq.30)pause 'too many iterations in tqli'
23           iter=iter+1
24           g=(d(l+1)-d(l))/(2.*e(l))
25           r=pythag(g,1.)
26           g=d(m)-d(l)+e(l)/(g+sign(r,g))
27           s=1.
28           c=1.
29           p=0.
30           do i=m-1,l,-1
31             f=s*e(i)
32             b=c*e(i)
33             r=pythag(f,g)
34             e(i+1)=r
35             if(r.eq.0.)then
36               d(i+1)=d(i+1)-p
37               e(m)=0.
38               goto 1
39             endif
40             s=f/r
41             c=g/r
42             g=d(i+1)-p
43             r=(d(i)-g)*s+2.*c*b
44             p=s*r
45             d(i+1)=g+p
46             g=c*r-b
47 !     Omit lines from here ...
48             do k=1,n
49               f=z(k,i+1)
50               z(k,i+1)=s*z(k,i)+c*f
51               z(k,i)=c*z(k,i)-s*f
52             end do
53 !     ... to here when finding only eigenvalues.
54           end do
55           d(l)=d(l)-p
56           e(l)=g
57           e(m)=0.
58           goto 1
59         endif
60       end do
61       return
62       END SUBROUTINE tqli