

c      Dear Physics 210 students,
c      Here is a fortran code that calls the library routine `jacobi'
c      to diagonalize a test matrix.  `eigsrt' sorts the eigenvalues.
c      N is the declared dimension of the array holding your matrix.
c      Np is actual dimension of the matrix you want to diagonalize.
c      Note that `jacobi' only works for real symmetric matrices.
c      Richard Scalettar

       program diag2
       implicit none
       integer N,Np
       parameter (Np=64)
       real*8 A(Np,Np),ev(Np)
       real*8 evect(Np,Np)
       integer i,j,nrot

c      SET UP SAMPLE MATRIX

       N=3
       do 20 i=1,N
       do 10 j=1,N
             A(i,j)=0.d0
10     continue
20     continue

       A(2,2)=3.d0
       A(1,1)=7.d0
       A(1,3)=0.5d0
       A(3,1)=0.5d0
       A(3,3)=7.d0

c      GET EIGENVALUES OF A.

       call jacobi(A,N,Np,ev,evect,nrot)
       call eigsrt(ev,evect,N,Np)

c      WRITE EIGENVALUES AND EIGENVECTORS

       do 400 j=1,N
       write (67,*) ' '
       write (67,990) j,ev(j)
       do 300 i=1,N
             write (67,991) i,evect(i,j)
300    continue
400    continue
990    format('eigenvector number  ',i3,'    eigenvalue= ',e16.8)
991    format(i8,2e16.8)

       return
       end



      SUBROUTINE JACOBI(A,N,NP,D,V,NROT)
      implicit double precision(a-h,o-z)
      PARAMETER (NMAX=200)
      double precision A(NP,NP),D(NP),V(NP,NP),B(NMAX),Z(NMAX)
      DO 12 IP=1,N
        DO 11 IQ=1,N
          V(IP,IQ)=0.d0
11      CONTINUE
        V(IP,IP)=1.d0
12    CONTINUE
      DO 13 IP=1,N
        B(IP)=A(IP,IP)
        D(IP)=B(IP)
        Z(IP)=0.d0
13    CONTINUE
      NROT=0
      DO 24 I=1,50
        SM=0.d0
        DO 15 IP=1,N-1
          DO 14 IQ=IP+1,N
            SM=SM+DABS(A(IP,IQ))
14        CONTINUE
15      CONTINUE
        IF(SM.EQ.0.d0)RETURN
        IF(I.LT.4)THEN
          TRESH=0.2d0*SM/dfloat(N)**2
        ELSE
          TRESH=0.d0
        ENDIF
        DO 22 IP=1,N-1
          DO 21 IQ=IP+1,N
            G=100.d0*DABS(A(IP,IQ))
            IF((I.GT.4).AND.(DABS(D(IP))+G.EQ.DABS(D(IP)))
     *         .AND.(DABS(D(IQ))+G.EQ.DABS(D(IQ))))THEN
              A(IP,IQ)=0.d0
            ELSE IF(DABS(A(IP,IQ)).GT.TRESH)THEN
              H=D(IQ)-D(IP)
              IF(DABS(H)+G.EQ.DABS(H))THEN
                T=A(IP,IQ)/H
              ELSE
                THETA=0.5d0*H/A(IP,IQ)
                T=1.d0/(DABS(THETA)+DSQRT(1.d0+THETA**2.d0))
                IF(THETA.LT.0.d0)T=-T
              ENDIF
              C=1.d0/DSQRT(1.d0+T**2.d0)
              S=T*C
              TAU=S/(1.d0+C)
              H=T*A(IP,IQ)
              Z(IP)=Z(IP)-H
              Z(IQ)=Z(IQ)+H
              D(IP)=D(IP)-H
              D(IQ)=D(IQ)+H
              A(IP,IQ)=0.d0
              DO 16 J=1,IP-1
                G=A(J,IP)
                H=A(J,IQ)
                A(J,IP)=G-S*(H+G*TAU)
                A(J,IQ)=H+S*(G-H*TAU)
16            CONTINUE
              DO 17 J=IP+1,IQ-1
                G=A(IP,J)
                H=A(J,IQ)
                A(IP,J)=G-S*(H+G*TAU)
                A(J,IQ)=H+S*(G-H*TAU)
17            CONTINUE
              DO 18 J=IQ+1,N
                G=A(IP,J)
                H=A(IQ,J)
                A(IP,J)=G-S*(H+G*TAU)
                A(IQ,J)=H+S*(G-H*TAU)
18            CONTINUE
              DO 19 J=1,N
                G=V(J,IP)
                H=V(J,IQ)
                V(J,IP)=G-S*(H+G*TAU)
                V(J,IQ)=H+S*(G-H*TAU)
19            CONTINUE
              NROT=NROT+1
            ENDIF
21        CONTINUE
22      CONTINUE
        DO 23 IP=1,N
          B(IP)=B(IP)+Z(IP)
          D(IP)=B(IP)
          Z(IP)=0.d0
23      CONTINUE
24    CONTINUE
      Write(63,*) '50 iterations should never happen'
      RETURN
      END



      SUBROUTINE EIGSRT(D,V,N,NP)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION D(NP),V(NP,NP)
      DO 13 I=1,N-1
        K=I
        P=D(I)
        DO 11 J=I+1,N
          IF(D(J).GE.P)THEN
            K=J
            P=D(J)
          ENDIF
11      CONTINUE
        IF(K.NE.I)THEN
          D(K)=D(I)
          D(I)=P
          DO 12 J=1,N
            P=V(J,I)
            V(J,I)=V(J,K)
            V(J,K)=P
12        CONTINUE
        ENDIF
13    CONTINUE
      RETURN
      END

