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        csend, crecv 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
C          Permute: maps the offset from the beginning of partially
C                   transformed array into the index of the element stored
C                   at that location.
C
C
C         IsOdd(int) : if the input is odd this returns TRUE, else it
C                      returns FALSE.
C         Log2(int) : finds the logarithm to the base 2 of the input.
C         communication:
C            It is through IPSC communication calls.

	program FFTNode

	include '/usr/include/fcube.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

	thisNode = mynode()
	hostNode = myhost()

	call crecv(gtype,npts,intSize)


	numNode	= numnodes()
	numDim	= nodedim()
	numLoc  = npts/numNode
	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          n - dimension of the complex vector.
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,n)
	integer n
	complex y(0:n-1)

	include '../include/fft.h'

	include '../include/sysdep.h'

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

	return
	end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C			OutResult				C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     this subroutine passes the solution computed by the current
C     node to the hostnode.
C     input parameters:
C        y - complex vector partitioned into blocks.
C        n - dimension of the vector.
C     output parameters:
C        NONE
C     computation:
C        NONE.
C     communication:
C        send result to the host node.CSEND.

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

	include '../include/fft.h'

	include '../include/sysdep.h'

	call csend(gtype+thisNode,y,cmplxSize*n,hostNode,npid)

	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          y - input distributed data vector.
C          n - the dimension of the data vector.
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(y,n)
	integer n
	complex y(0:n-1)

	include '../include/fft.h'

	include '../include/sysdep.h'

	integer timeBegin, timeEnd, mclock
	integer i, k, nx, incrm

	incrm = npts/2
	nx = 2

        timeBegin = mclock()

	do k = numDim+numLocLg-1, 0, -1
	    if (k .GE. numLocLg) then
		call FFTcom(y,numLoc,nx)
	    else

		    call FFTc2(y,numLoc,nx,incrm,i)

	    endif
	    nx = nx*2
	    incrm = incrm/2

	end do

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

	return
	end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C			FFTcom					C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     this subroutine computes the fourier transform.
C     input parameters:
C          y: input distributed data vector.
C          n: dimension of the data vector.
C         nx: step
C     output parameters:
C         the output is returned in y itself
C     system calls used:
C         csend, crecv calls. FORTRAN mathlib calls- CMPLX, SIN, COS
C     user-defined function calls:
C         Isodd
C     computation:
C           y(i:0:n-1) = cVal(i)-term2(i) if the current location is odd.
C           y(i:0:n-1) = cVal(i)+term2(i) if the current location is even.
C     communication:
C           send y to neighboring processor if odd(current node - nodeincrm)
C           send y to neighboring processor if even(current node + nodeincrm)
C           receive cVal

	subroutine FFTcom(y,n,nx)
	integer	n,nx
	complex	y(0:n-1)

	include '../include/fft.h'

	include '../include/sysdep.h'

	parameter (PI = 3.14)
	integer nodeIncrm, loc, indx, Permute
	complex cVal(0:NN-1), cexpon, term2(0:NN-1)
	real	fac

	nodeIncrm	= numNode/nx
	loc		= thisNode/nodeIncrm

	if (IsOdd(loc)) then
	    call csend(gtype,y,cmplxSize*n,thisNode-nodeIncrm,npid)
	else
	    call csend(gtype,y,cmplxSize*n,thisNode+nodeIncrm,npid)
	endif
	call crecv(gtype,cVal,cmplxSize*n)


	indx 	= Permute(loc/2,nx/2)
	fac	= 2.0*PI*indx/nx
	cexpon	= CMPLX(COS(fac),SIN(fac))

	if (IsOdd(loc)) then
	  do i = 0, n-1
	    term2(i) = cexpon*y(i)
	    y(i) = cVal(i)-term2(i)
	  end do
	else
	  do i = 0, n-1
	    term2(i) = cexpon*cVal(i)
	    y(i) = y(i) + term2(i)
	  end do
	endif

	return
	end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C			FFTc2					C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     this subroutine combines the even and odd elements of the
C     data vector within the same processor.
C     input parameters:
C         y - the distributed data vector
C         n - number of elements in the data vector.
C        nx - step = 2
C        increm - half of npts
C        start - start location
C     output parameters:
C        the result is returned in vector y itself.

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

	
	include '../include/sysdep.h'

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

	minPtr 	= thisNode * nx / (2*numNode)
	maxPtr	= minPtr + nx / (2*numNode)

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

	    evenPtr = evenPtr + 2*increm
	    oddPtr  = oddPtr  + 2*increm
	  end do
	end do

	return
	end

	
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C			Permute					C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     this function maps the offset from the beginning of the partially
C     transformed array, x, into the index of the element stored at
C     that location.
C     input parameters:
C         l: offset l
C         m: number of elements
C     output parameters:
C

	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
C     this subroutine checks whether the input is even or odd.
C     returns true if n is odd.

	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
C     finds the logarithm to the base 2 of the input and returns the
C     truncated value.

	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
