c     program rede_cubica_2

      implicit real*8 (a-h,o-z)
      Parameter(L=3, Nat=27, d1viz=1.1d0, V=1.0d0)

      Real*8 d(Nat), e(Nat), a(Nat,Nat), r(3,Nat)
      external tqli, pythag, tred2
      Open(10,File='autovalores')

c      Do i=1,Nat
c         Read(10,*) r(1,i), r(2,i), r(3,i)
c      Enddo

      Size=DBLE(L)
      Iat=1
      Do i=1,L
      Do j=1,L
      Do k=1,L
         r(1,Iat)=DBLE(i-1)
         r(2,Iat)=DBLE(j-1)
         r(3,Iat)=DBLE(k-1)
         Iat=Iat+1
      Enddo
      Enddo
      Enddo

      Do iat=1,Nat
      Print*,r(1,iat),r(2,iat), r(3,iat)
      Enddo


      Do iat=1,Nat
      Do jat=1,Nat
         a(iat,jat)=0.0d0
      Enddo
      Enddo

      Do iat=1,Nat-1
      Do jat=iat+1,Nat
         dx=r(1,iat)-r(1,jat)
         dy=r(2,iat)-r(2,jat)
         dz=r(3,iat)-r(3,jat)

         If(dx.gt.(Size/2)) dx=dx-Size
         If(dy.gt.(Size/2)) dy=dy-Size
         If(dz.gt.(Size/2)) dz=dz-Size
         If(dx.lt.(-Size/2)) dx=dx+Size
         If(dy.lt.(-Size/2)) dy=dy+Size
         If(dz.lt.(-Size/2)) dz=dz+Size

         Dist=dsqrt(dx**2 + dy**2 + dz**2)

         If(Dist.lt.d1viz)Then
            a(iat,jat)=V
            a(jat,iat)=V
         EndIf
      Enddo
      Enddo

      Call tred2(a,Nat,Nat,d,e)
      Call tqli (d,e,Nat,Nat,a)
      Call ordena (d,a,Nat)
      do n=1,Nat
         Print*,d(n)
         Write(10,*)d(n)
      Enddo

      stop

      end
c
c --------------------------------------------------
c
      SUBROUTINE tqli(d,e,n,Nat,a)
c
c        d = na entrada, elementos diagonais;na saida, autovalores
c        e = na entrada, elementos fora da diagonal;
c        a = na entrada é a matriz identidade; saida, as colunas são os autovetores
c
      INTEGER n,Nat
      REAL*8 d(Nat),e(Nat),a(Nat,Nat)
      INTEGER i,iter,k,l,m
      REAL*8 b,c,dd,f,g,p,r,s,pythag

      external pythag

      do 11 i=2,n
        e(i-1)=e(i)
 11   continue

      e(n)=0.d0
      do 15 l=1,n
        iter=0

 1      do 12 m=l,n-1
          dd = dabs(d(m))+ dabs(d(m+1))
          if ((dabs(e(m))+dd).eq.dd) goto 2
 12     continue

        m=n
 2      if(m.ne.l)then
          if(iter.eq.30)pause 'too many iterations in tqli'
          iter=iter+1
          g=(d(l+1)-d(l))/(2.d0*e(l))
          r=pythag(g,1.d0)
          g=d(m)-d(l)+e(l)/(g+dsign(r,g))
          s=1.d0
          c=1.d0
          p=0.d0

          do 14 i=m-1,l,-1
            f=s*e(i)
            b=c*e(i)
            r=pythag(f,g)
            e(i+1)=r
            if(r.eq.0.d0)then
              d(i+1)=d(i+1)-p
              e(m)=0.d0
              goto 1
            endif
            s=f/r
            c=g/r
            g=d(i+1)-p
            r=(d(i)-g)*s+2.d0*c*b
            p=s*r
            d(i+1)=g+p
            g=c*r-b
C     Omit lines from here ...
C
C           do 13 k=1,n
C              f=z(k,i+1)
C              z(k,i+1)=s*z(k,i)+c*f
C              z(k,i)=c*z(k,i)-s*f
C13          continue
C
C     ... to here when finding only eigenvalues.
 14       continue

          d(l)=d(l)-p
          e(l)=g
          e(m)=0.d0
          goto 1

        endif
 15   continue

      return
      END
C
c  ------------------------------------------------------------
c.
      FUNCTION pythag(a,b)
      REAL*8 a,b,pythag
      REAL*8 absa,absb
      absa=dabs(a)
      absb=dabs(b)
      if(absa.gt.absb)then
        pythag=absa*dsqrt(1.d0+(absb/absa)**2)
      else
        if(absb.eq.0.d0)then
          pythag=0.d0
        else
          pythag=absb*dsqrt(1.d0+(absa/absb)**2)
        endif
      endif
      return
      END

c      ---------------------------------------------------------

      SUBROUTINE tred2(a,n,Nat,d,e)
C
      IMPLICIT NONE
C------------Variables--------------------------------------
C -- Input/output variables --
      INTEGER n, Nat
      REAL*8 a(Nat,Nat), d(Nat), e(Nat)
C -- Local variables --
      INTEGER i,j,k,l
      REAL*8 f,g,h,hh,scale
C------------Main-------------------------------------------
C
      DO 18 i = n, 2, -1
        l = i - 1
        h = 0.D0
        scale = 0.D0
        IF(l.GT.1) THEN
          DO 11 k = 1, l
            scale = scale + DABS(a(i,k))
 11       CONTINUE
          IF(scale.EQ.0.D0) THEN
C                            (Skip transformation.)
            e(i)=a(i,l)
          ELSE
            DO 12 k = 1, l
              a(i,k) = a(i,k)/scale
C                  (Use scaled a's for transformation.)
              h = h + a(i,k)**2
C                  (Form sigma in h.)
 12         CONTINUE
            f = a(i,l)
            g = - DSIGN( DSQRT(h), f )
            e(i) = scale*g
            h = h - f*g
C                  (Now h is equation (11.2.4).)
            a(i,l) = f - g
C                  (Store u in the i-th row of a.)
            f = 0.D0
            DO 15 j = 1, l
C  Omit following line if finding only eigenvalues.
              a(j,i) = a(i,j)/h
C                  (Store u/H in i-th column of a.)
              g = 0.D0
C                  (Form an element of A*u in g.)
              DO 13 k = 1, j
                g = g + a(j,k)*a(i,k)
 13           CONTINUE
              DO 14 k = j+1, l
                g = g + a(k,j)*a(i,k)
 14           CONTINUE
              e(j) = g/h
C                  (Form element of p in temporarily
C                                  unused element of e.)
              f = f + e(j) * a(i,j)
 15         CONTINUE
            hh = f/(h+h)
C                  (Form K, equation (11.2.11).)
            DO 17 j = 1, l
C                  (Form q and store in e overwriting p.)
              f = a(i,j)
C                  (Note that e(l)=e(i-1) survives.)
              g = e(j) - hh*f
              e(j) = g
              DO 16 k = 1, j
C                  (Reduce a, equation (11.2.13).)
                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
C
C  Omit following line if finding only eigenvalues.
      d(1) = 0.D0
      e(1) = 0.D0
C
      DO 24 i = 1, n
C            (Begin accumulation of transformation matrices.)
C  Delete lines from here ...
        l = i - 1
        IF(d(i).NE.0.D0) THEN
C            (This block skipped when i=1.)
          DO 22 j = 1, l
            g = 0.D0
            DO 19 k = 1, l
C               (Use u and u/H stored in a to form P*Q.)
              g = g + a(i,k)*a(k,j)
 19         CONTINUE
            DO 21 k = 1, l
              a(k,j) = a(k,j) - g*a(k,i)
 21         CONTINUE
 22       CONTINUE
        ENDIF
C  ... to here when finding only eigenvalues.
C
        d(i) = a(i,i)
C              (This statement remains.)
C
C  Also delete lines from here ...
C       (Reset row and column of a to identity matrix
C                                     for next iteration)
        a(i,i) = 1.D0
        DO 23 j = 1, l
          a(i,j) = 0.D0
          a(j,i) = 0.D0
 23     CONTINUE
C  ... to here when finding only eigenvalues.
C
 24   CONTINUE
C
C-----------------------------------------------------------
      RETURN
      END
C***********************************************************
C***************************************************END*****
C  .
c  ---------------------------------------------------------------
c
      subroutine ordena(d,a,Nat)

	implicit real*8 (a-h,o-z)
	dimension d(Nat), a(Nat,Nat)

	do i=1,Nat-1

	   do j=i+1,Nat

	   if(d(i).gt.d(j)) then

	      aux1  = d(j)
	      d(j)  = d(i)
	      d(i)  = aux1

	      do k=1,Nat

	        aux2   = a(k,j)
	        a(k,j) = a(k,i)
	        a(k,i) = aux2

	      end do


	   end if

	   end do

	end do

	return
	end