      program eigen
c
      integer*4 np,n,mp
      parameter (np=534)
      parameter (mp=60000)
      real*8 a(np,np),d(np),v(np,np),vd(mp)
      integer*4 jv(mp),iv(np+1)

      open(1,file='v.dat')
      do i=1,mp
       read(1,*,end=111)vd(i),jv(i)
      enddo
  111 close(1)
      open(1,file='iv.dat')
      do i=1,np+1
       read(1,*,end=112)iv(i)
      enddo
  112 close(1)
      n=534
      do i=1,n
       do j=iv(i),iv(i+1)-1
        a(i,jv(j))=vd(j)
        a(jv(j),i)=a(i,jv(j))
       enddo
      enddo

      write(6,*)'Call JACOBI now...'
      call jacobi(a,n,np,d,v,nrot)
      write(6,*)'JACOBI finished.'

      do i=1,n
       write(2,*)d(i)
      enddo

      end


      SUBROUTINE JACOBI(A,N,NP,D,V,NROT)
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      PARAMETER (NMAX=534)
      DIMENSION 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+ABS(A(IP,IQ))
14        CONTINUE
15      CONTINUE
        IF(SM.EQ.0.D0)RETURN
        IF(I.LT.4)THEN
          TRESH=0.2D0*SM/N**2
        ELSE
          TRESH=0.D0
        ENDIF
        DO 22 IP=1,N-1
          DO 21 IQ=IP+1,N
            G=100.D0*ABS(A(IP,IQ))
            IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP)))
     *         .AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN
              A(IP,IQ)=0.D0
            ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN
              H=D(IQ)-D(IP)
              IF(ABS(H)+G.EQ.ABS(H))THEN
                T=A(IP,IQ)/H
              ELSE
                THETA=0.5D0*H/A(IP,IQ)
                T=1.D0/(ABS(THETA)+SQRT(1.D0+THETA**2))
                IF(THETA.LT.0.D0)T=-T
              ENDIF
              C=1.D0/SQRT(1.D0+T**2)
              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
      PAUSE '50 iterations should never happen'
      RETURN
      END

