C       Hand Compiled Express Version.
C     this program computes the discrete Fourier transform of a series.
C     The Fourier transforms are finite and complex.
C     the following is the node program.
C
C     NO IMPLICIT VARIABLE ANYWHERE.
C     system calls:
C        KXWRIT,KXREAD and Fortran math library calls.
C
C     user-defined functions:
C          Getdata: generates complex data for the input sequence using
C                   the ramp function. The input sequence is partitioned
C                   into blocks and the values for each block is computed
C                   by the individual nodes.
C          outResult: each node sends the partial soln to the hostnode.
C          FFTit:   this routine finds the NN point FFT of the input data
C                   array. This calls subroutines FFTcom and FFTc2.
C
C          FFTcom:  this computes the fourier transform
C          FFTc2: this combines the even and odd elements of the data vector
C                 within the same processor.
C

	program FFTNode

	include '../include/express.h'

	include '../include/fft.h'

	include '../include/sysdep.h'

	data 	allNode,npid,gtype,intSize,cmplxSize / -1,0,999,4,8 /

	complex x(0:NN-1)
	common  /CDATA/x
	integer status,src

        CALL KXINIT
        status = KXPARA(ENV)

        status = KXBROD(npts, IHOST, intSize, IALNOD, dummy, gtype)

        hostNode = IHOST
        thisNode = ENV(1)
        numNode = ENV(2)


	numDim	= nodedim()
	numLoc  = npts/numNode
	BB = numLoc

	minLoc	= thisNode*numLoc
	maxLoc 	= minLoc + numLoc
 	numLocLg = Log2(numLoc)

	call GetData(x,numLoc)
	call FFTit(x,numLoc)
	call OutResult(x,numLoc)
	
	stop
	end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C			GetData					C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     this subroutine generates data using the ramp function.
C     input parameters:
C          y - a complex vector partitioned into blocks.
C          B - the block size.
C     output parameters:
C           None. It initialises y vector and returns it.
C     system calls used:
C           CMPLX(real,real)
C     user-defined functions:
C           NONE.
C     computation:
C            y(i:0-n-1) = CMPLX((i+minLoc)*slope,0.0)
C     communication:
C            the computation is done independently by each node.

	subroutine GetData(y,BB)
	integer BB
	complex y(0:BB-1)

	include '../include/express.h'

	include '../include/fft.h'

	include '../include/sysdep.h'

	real slope
	integer i
	
C	use ramp function now
	slope = 1.0 / (npts-1)
C	local i ==> globle i+minLoc
	do i = 0, BB-1
	    y(i) = CMPLX((i+minLoc)*slope,0.0)
	end do

	return
	end

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

	include '../include/express.h'	
        include '../include/fft.h'
	include '../include/sysdep.h'

	status = KXWRIT(y,cmplxSize*BB,hostNode,gtype+thisNode)

	return
	end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C			FFTit					C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     this subroutine finds the NN point FFT of distributed complex
C     data array y.
C     input parameters:
C          x - input distributed data vector.
C          BB - the block size
C     output parameters:
C          NONE.
C     system calls:
C          NONE.
C     user-defined subroutines:
C           FFTcom(y,numLoc,nx), FFTc2(y,numLoc,nx,incrm,i)
C

	subroutine FFTit(x,BB)
	integer BB
	complex x(0:BB-1)

	include '../include/express.h'
	include '../include/fft.h'
	include '../include/sysdep.h'

	integer i, j, k, l, s, nx, incrm, incrm2, depth, t
	integer srcPE1, srcPE2, desPE1, desPE2, offset
	integer Log2, Permute
	integer timeBegin, timeEnd, mclock
	complex term2(0:NN-1), xbuf(0:NN-1)
	real    fac(0:NN-1)
	logical mask(0:NN-1), mask1(0:NN-1), anyTrue
	parameter (PI = 3.14)

	incrm2 = npts
	incrm = npts/2
	nx = 2
	depth = Log2(npts)

        timeBegin = mclock()

	do k = depth-1, 0, -1
C	  calculate mask
 	  do i = 0, BB-1 
	    mask(i) = .FALSE. 
	  end do
	  do i = 0, incrm-1
	    do j = 0, nx/2-1
	      l = i + j*incrm2 +incrm
	      if (thisNode .EQ. l/BB) mask(MOD(l,BB)) = .TRUE. 
	    end do
	  end do

C	  calculate fac & term2
	  do s = 0, BB-1
	    if (mask(s)) then
	      fac(s) = 2.0 * PI * Permute((s+minLoc)/incrm2,nx/2) /nx
	      term2(s) = CMPLX(COS(fac(s)), SIN(fac(s))) * x(s)
	    endif
	  end do

C	  calculate the 1st part of x
C	  send
	  desPE1 = (BB*thisNode + incrm) / BB
	  do i = 0, BB-1 
	    mask1(i) = .FALSE. 
	  end do
	  do i = 0, incrm-1
	    do j = 0, nx/2-1
	      l = i + j*incrm2 + incrm
	      if (desPE1 .EQ. (l/BB)) mask1(MOD(l,BB)) = .TRUE.
	    end do
	  end do

	  anyTrue = .FALSE.
	  t = MIN(BB, incrm)
	  do i = 0, t-1
	    if (mask1(i)) anyTrue = .TRUE.
	  end do
	  if (anyTrue) then
	    do i = 0, t-1 
	      xbuf(i) = x(i+BB-t) 
	    end do
	    status = KXWRIT(xbuf,cmplxSize*t,desPE1,gtype)
	  endif

C	  receive
	  srcPE1 = (BB*thisNode-incrm)/BB
	  anyTrue = .FALSE.
	  do i = 0,t-1
	    if (mask(i)) anyTrue = .TRUE.
	  end do
	  if (anyTrue) then
	    src = NOCARE
	    status = KXREAD(xbuf,cmplxSize*t,src,gtype)
	  endif
	  offset = (thisNode-srcPE1)*BB

C	  compute
	  do s = 0, BB-1
	    if (mask(s)) then
	      if ((s-incrm) .GE. 0) then
		x(s) = x(s-incrm) - term2(s)
	      else
		x(s) = xbuf(s-incrm+offset) - term2(s)
	      endif
	    endif
	  end do
	    
C	  calculate the 2nd part of x
	  do i = 0, BB-1 
	    mask(i) = .FALSE. 
	  end do
	  do i = 0, incrm-1
	    do j = 0, nx/2-1
	      l = i+j*incrm2
	      if (thisNode .EQ. (l/BB)) mask(MOD(l,BB)) = .TRUE.
	    end do
	  end do

C	  send
	  desPE1 = (BB*thisNode - incrm) / BB
	  do i = 0, BB-1 
	    mask1(i) = .FALSE. 
	  end do
	  do i = 0, incrm-1
	    do j = 0, nx/2-1
	      l = i + j*incrm2
	      if (desPE1 .EQ. (l/BB)) mask1(MOD(l,BB)) = .TRUE.
	    end do
	  end do

	  anyTrue = .FALSE.
	  do i = BB-1, BB-t, -1
	    if (mask1(i)) anyTrue = .TRUE.
	  end do
	  if (anyTrue) then
	    status = KXWRIT(term2,cmplxSize*t,desPE1,gtype)
	  endif

C	  receive
	  srcPE1 = (BB*thisNode+incrm)/BB
	  anyTrue = .FALSE.
	  do i = BB-1, BB-t, -1
	    if (mask(i)) anyTrue = .TRUE.
	  end do
	  if (anyTrue) then
	    src = NOCARE
	    status = KXREAD(xbuf,cmplxSize*t,src,gtype)
	  endif
	  offset = (srcPE1-thisNode)*BB

C	  compute
	  do s = 0, BB-1
	    if (mask(s)) then
	      if ((s+incrm) .LE. BB-1) then
		x(s) = x(s) + term2(s+incrm)
	      else
		x(s) = x(s) + xbuf(s+incrm-offset)
	      endif
	    endif
	  end do

	  nx = nx*2
	  incrm2 = incrm
	  incrm = incrm/2
	end do

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

	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
