	program FFT

	parameter (NN=65536)
	integer npts
	common  /INFO/npts
	complex x(0:NN-1)
	common  /CDATA/x

	print *, 'Enter # of points: '
	read  *, npts
	
	print *, 'please wait...'

	call GetData(x,npts)
	call FFTit(x,npts)
	call GlbRev(x,npts)
	call OutResult(x,npts)
	
	stop
	end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C			GetData					C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	subroutine GetData(y,n)
	integer n
	complex y(0:n-1)

	integer npts
	common  /INFO/npts
	real slope
	integer i
	
C	use ramp function now
	slope = 1.0 / (npts-1)
	do i = 0, n-1
	    y(i) = CMPLX(i*slope,0.0)
	end do

	return
	end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C			OutResult				C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	subroutine OutResult(y,n)
	integer n
	complex y(0:n-1)

	integer npts
	common  /INFO/npts

	print *, 'result ...', (y(i), i=0,n-1)

	return
	end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C			FFTit					C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	subroutine FFTit(y,n)
	integer n
	complex y(0:n-1)

	integer npts
	common  /INFO/npts
C	integer timeBegin, timeEnd, mclock
	integer k, nx, incrm

C       DIR$ PAR y(BLK)

	incrm = npts/2
	nx = 2

C        timeBegin = mclock()

	do k = Log2(npts)-1, 0, -1
	    call FFTc2(y,npts,nx,incrm)
	    nx = nx*2
	    incrm = incrm/2
C	    print *,'k=',k,(y(i),i=0,n-1)
	end do

C        timeEnd = mclock()
C        print *,'Node', thisNode,'  Elapsed time:',
C     &          timeEnd-timeBegin,'mS'

	return
	end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C			FFTc2					C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	subroutine FFTc2(y,n,nx,increm)
	integer n, nx, increm
	complex y(0:n-1)

	
	integer evenPtr, oddPtr, minPtr, maxPtr, pairLoc, indx, Permute, i
	complex oddVal, evenVal, term2, cexpon
	real    fac
	parameter (PI = 3.14)

	minPtr 	= 0
	maxPtr	= minPtr + nx / 2
CP	print *,'PE',thisNode,'evenPtr=',evenPtr,'increm=',increm,
CP	    &			'minPtr=',minPtr,'maxPtr=',maxPtr
	do pairLoc = minPtr, maxPtr-1
	  indx    = Permute(pairLoc,nx/2)
	  fac	  = 2.0*PI*indx/nx
	  cexpon  = CMPLX(COS(fac),SIN(fac))
	  evenPtr = pairLoc*2*increm
	  oddPtr  = evenPtr + increm
	  do i = 0, increm-1
	    oddVal  = y(oddPtr)
	    evenVal = y(evenPtr)
  	    term2   = cexpon*oddVal
	    y(oddPtr)  = evenVal - term2
	    y(evenPtr) = evenVal + term2
	    evenPtr = evenPtr + 1
	    oddPtr  = oddPtr + 1
	  end do
	end do

	return
	end

	
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C			GlbRev					C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC#` 

   	subroutine GlbRev(y,n)
	integer n
	complex y(0:n-1)
	parameter	(NN = 65536)
	complex		x(0:NN-1)
	integer Permute

	do i=0,n-1
	  x(i) = y(Permute(i))
	end do
	do i=0,n-1
	  y(i) = x(i)
	end do

	return
	end


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C			Permute					C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	integer function Permute(l,m)
	integer l,m

	integer numBit, k, power2, Log2
	logical IsOdd

	numBit 	= Log2(m)
	k	= l
	power2 	= m / 2
	Permute = 0
	do i = 1, numBit
	    if (IsOdd(k)) Permute = Permute + power2
	    power2 = power2 / 2
	    k	   = k / 2
	end do

	return
	end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C			IsOdd					C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	logical function IsOdd(n)
	integer n

	if (2*(n/2) .EQ. n) then
	    IsOdd = .FALSE.
	else
	    IsOdd = .TRUE.
	endif

	return
	end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C			Log2					C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	integer function Log2(n)
	integer n

	integer k

	k = n
	Log2 = 0
	do i = 1, n
	    if (k .GT. 1) then
	    	Log2 = Log2 + 1
	    	k    = k / 2
	    else
		return
		endif
	end do

	end
