c      program weig(input,output,tape5=input,tape6=output)
       program weig
       implicit double precision(a-h,o-z)
       dimension amtrx(48,3,3),w(40000),rk(3,40000)
       CHARACTER*8 ntransd,intd,inpout
       OPEN(UNIT=8,FILE='ntransd',STATUS='OLD',FORM='FORMATTED')
c      OPEN(UNIT=7,FILE='intd',STATUS='OLD',FORM='FORMATTED')
c      OPEN(UNIT=6,FILE='kpout',STATUS='UNKNOWN',FORM='FORMATTED')
       read(8,*) ntrans
       write(6,*) 'ntran= ' ,ntrans
       read(8,11) (((amtrx(i,j,k),k=1,3),j=1,3),i=1,ntrans)
   10  format(i3)
   11  format(18f4.0)
       call intpnt(wsum,ntrans,amtrx,w,rk)
       stop
       end
       subroutine intpnt(wsum,ntrans,amtrx,w,rk)
       implicit double precision(a-h,o-z)
c
c      subroutine reads and computes k points for integration
c      over the brillouin zone.
c
       dimension amtrx(48,3,3), w(1),rk(3,1),xxk(40000)
c
c      wsum                the sum of the weight factors w(i)
c
       dimension xk(3,40000),modd(40000),rktran(3)
       nrk=0
       nmax=40000
       nrkmax =40000
c
c      read the number of k points to be used in the integration
c      over the brillouin zone.
c      abs(nrk) points will be read.
c      if nrk .le. 0 additional points will be computed.
c
C----- read(5,10) nrk
       nrka = iabs(nrk)
       if (nrk .eq. 0) goto 30
c
c      read parameters used to compute the the integration points.
c      nbandi, ifmini, and ifmaxi are as before the number of 
c      eigenvalues to be computed, and the lowest and highest filled
c      band. nx, ny, and nz are the number of points in the three
c      directions dermined by the lattice wave vectors. sx, sy, and
c      sz shifts the grid of integration points from the origin.
c
 30     continue
     
         write(6,*) ' nx,ny,nz '
c       read(7,*) nx,ny,nz
        read *, nx,ny,nz
        if(nx*ny*nz.gt.nmax)stop ' increase nmax'
         write(6,*) ' sx,sy,sz '
c       read(7,*) sx,sy,sz
        read *, sx,sy,sz
c       nx= 10
c       ny=10
c       nz=10
c       sx=0.d0
c       sy=0.d0
c       sz=0.d0
        if(nx*ny*nz.gt. nmax) stop 1
 40    format(3i5,5x,3i5,3f5.2)
c
c      set up uniform array of k points.
c
       dx = 1.d0/nx
       dy = 1.d0/ny
       dz = 1.d0/nz
       n = 0
       do 50 i=1,nx
       do 50 j=1,ny
       do 50 k=1,nz
       n = n+1
       xk(1,n) = (i-1+sx)*dx
       xk(2,n) = (j-1+sy)*dy
       xk(3,n) = (k-1+sz)*dz
       modd(n) = 1
 50    continue
c
c      reduce to irreducible zone
c
       dw = 1.d0/n
       do 130 i=1,n
       if (modd(i) .eq. 0) goto 130
       nrk = nrk  + 1
       if (nrk .gt. nrkmax) stop 2
       do 60 j=1,3
       rk(j,nrk) = xk(j,i)
 60    continue
       w(nrk) = dw
       if (i .eq. n) goto 130
c
c      operate on rk with the symmetry operations
c
       do 120 j=1,ntrans
       do 80 k=1,3
       rktran(k) = 0.d0
       do 70 l=1,3
       rktran(k) = rktran(k) + amtrx(j,k,l)*rk(l,nrk)
 70    continue
c
c      translate to interval 0-1.
c
       kdel = rktran(k)
       if (rktran(k) .lt. 0) kdel=kdel-1
       rktran(k) = rktran(k) - kdel
 80    continue
c
c      check for points related by symmetry
c
       ip = i+1
       do 120 k=ip,n
       if (modd(k) .eq. 0) goto 120
c
c      both the transformed vector and its inverse
c
       do 90 l=1,3
       if (abs(rktran(l)-xk(l,k)) .gt. 1.0d-5) goto 100
 90    continue
       w(nrk) = w(nrk)+dw
       modd(k) = 0
       goto 120
 100   do 110 l=1,3
       if (abs(1-rktran(l)-xk(l,k)) .gt. 1.0d-5) goto 120
 110   continue
       w(nrk) = w(nrk)+dw
       modd(k) = 0
 120   continue
 130   continue
c140   write(6,150)
 150   format(4h1  i,5x,5hnband,5x,5hifmin,2x,5hifmax,
     1 7x,1hw,19x,2hrk,/)
       wsum = 0.d0
       do 170 i=1,nrk
       wsum = wsum + w(i)
       write(6,160) i,(rk(j,i),j=1,3),w(i)
 160   format(1x,i3,5x,'rk=',3f10.6,4x,'w=',f10.6)
 170   continue
       do 176 i=1,nrk
       w(i)=w(i)*n
 176   continue
       wtsmall = 1.d8
       do 175 i=1,nrk
       if(w(i).lt.wtsmall) wtsmall=w(i)
 175   continue
       nsmall = wtsmall+0.01d0
       
       if (nrk .gt. nrka) write(6,180) nrk-nrka,nx,ny,nz,sx,sy,sz
 180   format(/,5h last,i5,35h k points generated by program from,
     1 11h parameters,/,5x,3hn =,3i5,6x,3hs =,3f6.2)
       if (abs(wsum) .lt. 1.e-8) wsum=0.d0
       write(6,190) wsum
 190   format( ' sum of  weight =', f10.5)
      write(6,*) ' k pts for mixed basis input'
      DO 185 I=1,nrk
 185  write(6,96) rk(1,I),rk(2,I),rk(3,I)
  96  FORMAT(3F10.5)
      write(6,*) ' weights for mixed basis input, before dividing'
      write(6,961) (W(J),J=1,nrk)
  961 FORMAT(16F5.0)
      call div(w,nrk,nsmall)
c
       return
       end
       subroutine div(w,nrk,nsmall)
       implicit double precision(a-h,o-z)
       dimension w(1),rw(40000)
       if(nrk.gt.40000) stop
       do 5 it= 0, nsmall-1
       nfac = nsmall - it
       do 10 i=1,nrk
       rw(i) = w(i)/nfac
       iw = rw(i)
       diff = iw - rw(i)
       if(abs(diff).gt.1.d-6) goto 5 
 10    continue
c --
       write(6,*) ' weights for mixed basis input, after dividing by',nfac
       write(6,961) (rw(J),J=1,nrk)
       stop
 5     continue
  961 FORMAT(16F5.0)
       return
       end