      program prepv
      implicit none
      integer*4 nh,ndim,lx,ir,il,iup,idow,itest,ip,iq,k,jn,i,j,nh2,info,what
      parameter(lx=20,nh=400,ndim=nh)
      integer*4 ivic(nh,4)
      real*8 dnrm2,rnupc
      real*8 matrix(4*nh,4*nh)
      integer n,lda,lwmax,lwork
      real*8 w(4*nh)
      parameter (lwmax=100000)
      real*8 work(lwmax)
      real*8, allocatable :: aux(:)
      ALLOCATE( aux( 4*nh * ( 4*nh + 1 ) / 2 ) ) 
      

      what=1
!!    1 ESSL upper triang
!!    2 ESSL lower triang

      do i=1,ndim
       itest=i/lx
       if(itest*lx.eq.i)then
        ir=lx
       else
        ir=0
       endif
       if(itest*lx.eq.i-1)then
        il=lx
       else
        il=0
       endif
       if(i.gt.ndim-lx)then
        iup=ndim
       else
        iup=0
       endif
       if(i.le.lx)then
        idow=ndim
       else
        idow=0
       endif
       ivic(i,1)=i+1-ir ! right
       ivic(i,3)=i-1+il ! left
       ivic(i,2)=i+lx-iup ! up
       ivic(i,4)=i-lx+idow ! down
      enddo


      nh2=2*nh

      do ip=1,4*nh
       do iq=1,4*nh
        matrix(ip,iq)=0.d0
       enddo
      enddo

      do k=1,2
       do ip=1,nh
        jn=ivic(ip,k)
        if(jn.ne.0) then
              matrix(ip,jn)=matrix(ip,jn)-1.d0
              matrix(jn,ip)=matrix(jn,ip)-1.d0
              matrix(ip+nh2,jn+nh2)=matrix(ip+nh2,jn+nh2)-1.d0
              matrix(jn+nh2,ip+nh2)=matrix(jn+nh2,ip+nh2)-1.d0
        endif
       enddo
      enddo
     
      N=4*nh
      LDA=4*nh
      LWORK=-1

      do ip=1,4*nh
       do iq=ip,4*nh
         if(matrix(ip,iq).ne.matrix(iq,ip))then
                 write(*,*)"errore!!",ip,iq
                 stop
         endif
       enddo
      enddo

      if(what.eq.1)then
       k = 0 !! lower triangle -- packing matrix into a linear array 
       DO j = 1, 4*nh   
        DO i = j,4*nh
         k = k + 1
         aux(k) = matrix(i,j)
        END DO
       END DO 
       CALL DSPEV(1,aux,w,matrix,4*nh,4*nh,work,3*4*nh) 
      endif
      if(what.eq.2)then
       k=0 !! upper triangle -- packing matrix into a linear array
       do j=1,4*nh
        do i=1,j
         k=k+1
         aux(k)=matrix(i,j)
        enddo
       enddo
       CALL DSPEV(1,aux,w,matrix,4*nh,4*nh,work,3*4*nh)
      endif
!!      CALL DSYEV('V','L',N,matrix,LDA,W,WORK,LWORK,INFO)
!!      LWORK = MIN( LWMAX, INT( WORK( 1 ) ) )
!!      CALL DSYEV('V','L',N,matrix,LDA,W,WORK,LWORK,INFO)

      write(6,*) 'Eigenval'
      do i=1,4*nh
       write(77,*) i,w(i)
      enddo

      do ip=1,4*nh
       rnupc=0.d0
       do iq=1,4*nh
        rnupc=rnupc+matrix(iq,ip)*matrix(iq,ip)
       enddo
       rnupc=dsqrt(rnupc)
       write(6,*) 'Autovector mod ->',ip,rnupc
      enddo

      write(6,*) 'Check orthonormality'
      do ip=1,4*nh
       do iq=1,4*nh
        rnupc=0.d0
        do k=1,4*nh
         rnupc=rnupc+matrix(k,ip)*matrix(k,iq)
        enddo
        if(ip.ne.iq.and.dabs(rnupc).gt.1.d-8)write(6,*)'Not orthogonal!',ip,iq,rnupc
       enddo
      enddo
      STOP
      end 
