The following routines (based on the original EISPACK library) perform a diagonalization of a real symmetric matrix based on the QL algorithm. The first step is to reduce the matrix to tridiagonal form (TRED2) and then a further routine (TQLI) finds the eigenvalues and eigenvectors of a tridiagonal matrix:
(Fortran 77 version)
PROGRAM D11R4
C Driver for routine TQLI
PARAMETER(NP=10,TINY=1.0E-6)
DIMENSION A(NP,NP),C(NP,NP),D(NP),E(NP),F(NP)
DATA C/5.0,4.0,3.0,2.0,1.0,0.0,-1.0,-2.0,-3.0,-4.0,
* 4.0,5.0,4.0,3.0,2.0,1.0,0.0,-1.0,-2.0,-3.0,
* 3.0,4.0,5.0,4.0,3.0,2.0,1.0,0.0,-1.0,-2.0,
* 2.0,3.0,4.0,5.0,4.0,3.0,2.0,1.0,0.0,-1.0,
* 1.0,2.0,3.0,4.0,5.0,4.0,3.0,2.0,1.0,0.0,
* 0.0,1.0,2.0,3.0,4.0,5.0,4.0,3.0,2.0,1.0,
* -1.0,0.0,1.0,2.0,3.0,4.0,5.0,4.0,3.0,2.0,
* -2.0,-1.0,0.0,1.0,2.0,3.0,4.0,5.0,4.0,3.0,
* -3.0,-2.0,-1.0,0.0,1.0,2.0,3.0,4.0,5.0,4.0,
* -4.0,-3.0,-2.0,-1.0,0.0,1.0,2.0,3.0,4.0,5.0/
DO 12 I=1,NP
DO 11 J=1,NP
A(I,J)=C(I,J)
11 CONTINUE
12 CONTINUE
CALL TRED2(A,NP,NP,D,E)
CALL TQLI(D,E,NP,NP,A)
WRITE(*,'(/1X,A)') 'Eigenvectors for a real symmetric matrix'
DO 16 I=1,NP
DO 14 J=1,NP
F(J)=0.0
DO 13 K=1,NP
F(J)=F(J)+C(J,K)*A(K,I)
13 CONTINUE
14 CONTINUE
WRITE(*,'(/1X,A,I3,A,F10.6)') 'Eigenvalue',I,' =',D(I)
WRITE(*,'(/1X,T7,A,T17,A,T31,A)') 'Vector','Mtrx*Vect.','Ratio'
DO 15 J=1,NP
IF (ABS(A(J,I)).LT.TINY) THEN
WRITE(*,'(1X,2F12.6,A12)') A(J,I),F(J),'div. by 0'
ELSE
WRITE(*,'(1X,2F12.6,E14.6)') A(J,I),F(J),
* F(J)/A(J,I)
ENDIF
15 CONTINUE
WRITE(*,'(/1X,A)') 'press ENTER to continue...'
READ(*,*)
16 CONTINUE
END
SUBROUTINE TRED2(A,N,NP,D,E)
DIMENSION A(NP,NP),D(NP),E(NP)
IF(N.GT.1)THEN
DO 18 I=N,2,-1
L=I-1
H=0.
SCALE=0.
IF(L.GT.1)THEN
DO 11 K=1,L
SCALE=SCALE+ABS(A(I,K))
11 CONTINUE
IF(SCALE.EQ.0.)THEN
E(I)=A(I,L)
ELSE
DO 12 K=1,L
A(I,K)=A(I,K)/SCALE
H=H+A(I,K)**2
12 CONTINUE
F=A(I,L)
G=-SIGN(SQRT(H),F)
E(I)=SCALE*G
H=H-F*G
A(I,L)=F-G
F=0.
DO 15 J=1,L
A(J,I)=A(I,J)/H
G=0.
DO 13 K=1,J
G=G+A(J,K)*A(I,K)
13 CONTINUE
IF(L.GT.J)THEN
DO 14 K=J+1,L
G=G+A(K,J)*A(I,K)
14 CONTINUE
ENDIF
E(J)=G/H
F=F+E(J)*A(I,J)
15 CONTINUE
HH=F/(H+H)
DO 17 J=1,L
F=A(I,J)
G=E(J)-HH*F
E(J)=G
DO 16 K=1,J
A(J,K)=A(J,K)-F*E(K)-G*A(I,K)
16 CONTINUE
17 CONTINUE
ENDIF
ELSE
E(I)=A(I,L)
ENDIF
D(I)=H
18 CONTINUE
ENDIF
D(1)=0.
E(1)=0.
DO 23 I=1,N
L=I-1
IF(D(I).NE.0.)THEN
DO 21 J=1,L
G=0.
DO 19 K=1,L
G=G+A(I,K)*A(K,J)
19 CONTINUE
DO 20 K=1,L
A(K,J)=A(K,J)-G*A(K,I)
20 CONTINUE
21 CONTINUE
ENDIF
D(I)=A(I,I)
A(I,I)=1.
IF(L.GE.1)THEN
DO 22 J=1,L
A(I,J)=0.
A(J,I)=0.
22 CONTINUE
ENDIF
23 CONTINUE
RETURN
END
SUBROUTINE TQLI(D,E,N,NP,Z)
DIMENSION D(NP),E(NP),Z(NP,NP)
IF (N.GT.1) THEN
DO 11 I=2,N
E(I-1)=E(I)
11 CONTINUE
E(N)=0.
DO 15 L=1,N
ITER=0
1 DO 12 M=L,N-1
DD=ABS(D(M))+ABS(D(M+1))
IF (ABS(E(M))+DD.EQ.DD) GO TO 2
12 CONTINUE
M=N
2 IF(M.NE.L)THEN
IF(ITER.EQ.30)PAUSE 'too many iterations'
ITER=ITER+1
G=(D(L+1)-D(L))/(2.*E(L))
R=SQRT(G**2+1.)
G=D(M)-D(L)+E(L)/(G+SIGN(R,G))
S=1.
C=1.
P=0.
DO 14 I=M-1,L,-1
F=S*E(I)
B=C*E(I)
IF(ABS(F).GE.ABS(G))THEN
C=G/F
R=SQRT(C**2+1.)
E(I+1)=F*R
S=1./R
C=C*S
ELSE
S=F/G
R=SQRT(S**2+1.)
E(I+1)=G*R
C=1./R
S=S*C
ENDIF
G=D(I+1)-P
R=(D(I)-G)*S+2.*C*B
P=S*R
D(I+1)=G+P
G=C*R-B
DO 13 K=1,N
F=Z(K,I+1)
Z(K,I+1)=S*Z(K,I)+C*F
Z(K,I)=C*Z(K,I)-S*F
13 CONTINUE
14 CONTINUE
D(L)=D(L)-P
E(L)=G
E(M)=0.
GO TO 1
ENDIF
15 CONTINUE
ENDIF
RETURN
END
(Fortran 90 version)
SUBROUTINE tred2(a,d,e,novectors)
USE nrtype; USE nrutil, ONLY : assert_eq,outerprod
IMPLICIT NONE
REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a
REAL(SP), DIMENSION(:), INTENT(OUT) :: d,e
LOGICAL(LGT), OPTIONAL, INTENT(IN) :: novectors
INTEGER(I4B) :: i,j,l,n
REAL(SP) :: f,g,h,hh,scale
REAL(SP), DIMENSION(size(a,1)) :: gg
LOGICAL(LGT) :: yesvec=.true.
n=assert_eq(size(a,1),size(a,2),size(d),size(e),'tred2')
if (present(novectors)) yesvec=.not. novectors
do i=n,2,-1
l=i-1
h=0.0
if (l > 1) then
scale=sum(abs(a(i,1:l)))
if (scale == 0.0) then
e(i)=a(i,l)
else
a(i,1:l)=a(i,1:l)/scale
h=sum(a(i,1:l)**2)
f=a(i,l)
g=-sign(sqrt(h),f)
e(i)=scale*g
h=h-f*g
a(i,l)=f-g
if (yesvec) a(1:l,i)=a(i,1:l)/h
do j=1,l
e(j)=(dot_product(a(j,1:j),a(i,1:j)) &
+dot_product(a(j+1:l,j),a(i,j+1:l)))/h
end do
f=dot_product(e(1:l),a(i,1:l))
hh=f/(h+h)
e(1:l)=e(1:l)-hh*a(i,1:l)
do j=1,l
a(j,1:j)=a(j,1:j)-a(i,j)*e(1:j)-e(j)*a(i,1:j)
end do
end if
else
e(i)=a(i,l)
end if
d(i)=h
end do
if (yesvec) d(1)=0.0
e(1)=0.0
do i=1,n
if (yesvec) then
l=i-1
if (d(i) /= 0.0) then
gg(1:l)=matmul(a(i,1:l),a(1:l,1:l))
a(1:l,1:l)=a(1:l,1:l)-outerprod(a(1:l,i),gg(1:l))
end if
d(i)=a(i,i)
a(i,i)=1.0
a(i,1:l)=0.0
a(1:l,i)=0.0
else
d(i)=a(i,i)
end if
end do
END SUBROUTINE tred2
SUBROUTINE tqli(d,e,z)
USE nrtype; USE nrutil, ONLY : assert_eq,nrerror
USE nr, ONLY : pythag
IMPLICIT NONE
REAL(SP), DIMENSION(:), INTENT(INOUT) :: d,e
REAL(SP), DIMENSION(:,:), OPTIONAL, INTENT(INOUT) :: z
INTEGER(I4B) :: i,iter,l,m,n,ndum
REAL(SP) :: b,c,dd,f,g,p,r,s
REAL(SP), DIMENSION(size(e)) :: ff
n=assert_eq(size(d),size(e),'tqli: n')
if (present(z)) ndum=assert_eq(n,size(z,1),size(z,2),'tqli: ndum')
e(:)=eoshift(e(:),1)
do l=1,n
iter=0
iterate: do
do m=l,n-1
dd=abs(d(m))+abs(d(m+1))
if (abs(e(m))+dd == dd) exit
end do
if (m == l) exit iterate
if (iter == 30) call nrerror('too many iterations in tqli')
iter=iter+1
g=(d(l+1)-d(l))/(2.0_sp*e(l))
r=pythag(g,1.0_sp)
g=d(m)-d(l)+e(l)/(g+sign(r,g))
s=1.0
c=1.0
p=0.0
do i=m-1,l,-1
f=s*e(i)
b=c*e(i)
r=pythag(f,g)
e(i+1)=r
if (r == 0.0) then
d(i+1)=d(i+1)-p
e(m)=0.0
cycle iterate
end if
s=f/r
c=g/r
g=d(i+1)-p
r=(d(i)-g)*s+2.0_sp*c*b
p=s*r
d(i+1)=g+p
g=c*r-b
if (present(z)) then
ff(1:n)=z(1:n,i+1)
z(1:n,i+1)=s*z(1:n,i)+c*ff(1:n)
z(1:n,i)=c*z(1:n,i)-s*ff(1:n)
end if
end do
d(l)=d(l)-p
e(l)=g
e(m)=0.0
end do iterate
end do
END SUBROUTINE tqli