C       Hand Compiled 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        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
	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 - 1
 	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          BB - 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/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
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        BB - the block size
C     output parameters:
C        NONE
C     computation:
C        NONE.
C     communication:
C        send result to the host node.CSEND.

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

	include '../include/fft.h'

	include '../include/sysdep.h'

	call csend(gtype+thisNode,y,cmplxSize*BB,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          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/fft.h'

	include '../include/sysdep.h'

	integer i, j, k, l, s, nx, incrm, incrm2, depth
        integer tmp1, tmp2, inxg, pos
	integer srcPE1, srcPE2, desPE1, desPE2, offset1, offset2
	integer Log2, Permute, lb, ub, lb1, ub1, lb2, ub2, s1, s2
	integer timeBegin, timeEnd, mclock
	complex term2(0:NN-1), fac(0:NN-1), fac1
        complex buf1(0:NN-1), buf2(0:NN-1)
	real    tmp(0:NN-1), tp1
	logical 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 the 1st part of x
C	  send
          lb = incrm
          ub = (nx/2-1)*incrm2+incrm-1+incrm
          desPE1 = (minLoc + incrm) / BB
          if (desPE1 .NE. thisNode) then
            if (INT(lb/BB) .LT. desPE1)  lb1 = 0
            if (INT(lb/BB) .EQ. desPE1)  lb1 = MOD(lb,BB)
            if (INT(lb/BB) .GT. desPE1)  lb1 = BB
            if (INT(ub/BB) .LT. desPE1)  ub1 = -1
            if (INT(ub/BB) .EQ. desPE1)  ub1 = MOD(ub,BB)
            if (INT(ub/BB) .GT. desPE1)  ub1 = BB-1

            anyTrue = .FALSE.
            s1 = MAX(MOD(incrm,BB),lb1)
            s2 = MIN(BB-1,ub1)
            tmp1 = desPE1*BB-incrm
            do s = s1, s2
              if (MOD(s+tmp1, incrm2) .LE. incrm-1) then
	        anyTrue = .TRUE.
              end if
	    end do
	    if (anyTrue) then
	      call csend(gtype,x,cmplxSize*BB,desPE1,npid)
	    endif
          endif

          desPE2 = (maxLoc + incrm) / BB
          if ((srcPE2 .NE. thisNode) .AND. (srcPE2 .NE. srcPE1)) then
            if (INT(lb/BB) .LT. desPE2)  lb2 = 0
            if (INT(lb/BB) .EQ. desPE2)  lb2 = MOD(lb,BB)
            if (INT(lb/BB) .GT. desPE2)  lb2 = BB
            if (INT(ub/BB) .LT. desPE2)  ub2 = -1
            if (INT(ub/BB) .EQ. desPE2)  ub2 = MOD(ub,BB)
            if (INT(ub/BB) .GT. desPE2)  ub2 = BB-1

            anyTrue = .FALSE.
            s1 = MAX(0,lb2)
            s2 = MIN(MOD(-1+incrm,BB),ub2)
            tmp1 = desPE2*BB-incrm
            do s = s1, s2
              if (MOD(s+tmp1,incrm2) .LE. incrm-1) then
	        anyTrue = .TRUE.
              end if
	    end do
	    if (anyTrue) then
	    call csend(gtype,x,cmplxSize*BB,desPE2,npid)
	    end if
          end if

          if (INT(lb/BB) .LT. thisNode)  lb1 = 0
          if (INT(lb/BB) .EQ. thisNode)  lb1 = MOD(lb,BB)
          if (INT(lb/BB) .GT. thisNode)  lb1 = BB
          if (INT(ub/BB) .LT. thisNode)  ub1 = -1
          if (INT(ub/BB) .EQ. thisNode)  ub1 = MOD(ub,BB)
          if (INT(ub/BB) .GT. thisNode)  ub1 = BB-1

C	  calculate fac & term2
          pos = -1
          tmp1 = minLoc-incrm
          do s = lb1, ub1
            if (MOD(s+tmp1,incrm2) .LE. incrm-1) then	
              inxg = INT((s+minLoc)/incrm2)*incrm2+incrm
              if (inxg .NE. pos) then
                pos = inxg
                tp1 = 2.0 * PI * Permute(inxg/incrm2,nx/2) /nx
                fac1 = CMPLX(COS(tp1), SIN(tp1))
              end if
              term2(s) = fac1 * x(s)
	    endif
	  end do

C	  receive
          srcPE1 = (minLoc-incrm)/BB
          if (srcPE1 .NE. thisNode) then
	    anyTrue = .FALSE.
            s1 = MAX(0,lb1)
            s2 = MIN(MOD(-1+incrm,BB),ub1)
            do s = s1, s2
              if (MOD(s+tmp1,incrm2) .LE. incrm-1) then
	        anyTrue = .TRUE.
              end if
	    end do
            if (anyTrue) then
CP            print *,'PE',thisNode,' k=',k, ' recv1'
	      call crecv(gtype,buf1,cmplxSize*BB)
            endif
            offset1 = (thisNode-srcPE1)*BB
          end if
 
          srcPE2 = (maxLoc-incrm)/BB
          if ((srcPE2 .NE. thisNode) .AND. (srcPE2 .NE. srcPE1)) then
	    anyTrue = .FALSE.
            s1 = MAX(MOD(incrm,BB),lb1)
            s2 = MIN(BB-1,ub1)
            do s = s1, s2
              if (MOD(s+tmp1,incrm2) .LE. incrm-1) then
	        anyTrue = .TRUE.
              end if
	    end do
            if (anyTrue) then

              call crecv(gtype,buf2,cmplxSize*BB)
	    endif
            offset2 = (thisNode-srcPE2)*BB
          endif

C	  compute
          tmp1 = minLoc-incrm
          do s = lb1, ub1
            if (MOD(s+tmp1,incrm2) .LE. incrm-1) then
	      if ((s-incrm) .GE. 0) then
		x(s) = x(s-incrm) - term2(s)
	      else 
                if (INT((s+tmp1)/BB) .EQ. srcPE1) then
  		  x(s) = buf1(s-incrm+offset1) - term2(s)
                else
  		  x(s) = buf2(s-incrm+offset2) - term2(s)
                endif
	      endif
	    endif
	  end do
	    
C	  calculate the 2nd part of x
C	  send
          lb = 0
          ub = (nx/2-1)*incrm2+incrm-1
          desPE1 = (minLoc - incrm) / BB
          if (desPE1 .NE. thisNode) then
            if (INT(lb/BB) .LT. desPE1)  lb1 = 0
            if (INT(lb/BB) .EQ. desPE1)  lb1 = MOD(lb,BB)
            if (INT(lb/BB) .GT. desPE1)  lb1 = BB
            if (INT(ub/BB) .LT. desPE1)  ub1 = -1
            if (INT(ub/BB) .EQ. desPE1)  ub1 = MOD(ub,BB)
            if (INT(ub/BB) .GT. desPE1)  ub1 = BB-1
            
            anyTrue = .FALSE.
            s1 = MAX(MOD(NN-incrm,BB),lb1)
            s2 = MIN(BB-1,ub1)
            tmp1 = desPE1*BB
            do s = s1, s2
              if (MOD(s+tmp1,incrm2) .LE. incrm-1) then
	        anyTrue = .TRUE.
              end if
	    end do
            if (anyTrue) then
              call csend(gtype,term2,cmplxSize*BB,desPE1,npid)
  	    endif
          endif

          desPE2 = (maxLoc - incrm) / BB
          if ((srcPE2 .NE. thisNode) .AND. (srcPE2 .NE. srcPE1)) then
            if (INT(lb/BB) .LT. desPE2)  lb2 = 0
            if (INT(lb/BB) .EQ. desPE2)  lb2 = MOD(lb,BB)
            if (INT(lb/BB) .GT. desPE2)  lb2 = BB
            if (INT(ub/BB) .LT. desPE2)  ub2 = -1
            if (INT(ub/BB) .EQ. desPE2)  ub2 = MOD(ub,BB)
            if (INT(ub/BB) .GT. desPE2)  ub2 = BB-1

            anyTrue = .FALSE.
            s1 = MAX(0,lb1)
            s2 = MIN(MOD(-1+NN-incrm,BB),ub1)
            tmp1 = desPE2*BB
            do s = s1, s2
              if (MOD(s+tmp1,incrm2) .LE. incrm-1) then
	        anyTrue = .TRUE.
              end if
	    end do
	    if (anyTrue) then
	      call csend(gtype,term2,cmplxSize*BB,desPE2,npid)
	    endif
          endif

          if (INT(lb/BB) .LT. thisNode)  lb1 = 0
          if (INT(lb/BB) .EQ. thisNode)  lb1 = MOD(lb,BB)
          if (INT(lb/BB) .GT. thisNode)  lb1 = BB
          if (INT(ub/BB) .LT. thisNode)  ub1 = -1
          if (INT(ub/BB) .EQ. thisNode)  ub1 = MOD(ub,BB)
          if (INT(ub/BB) .GT. thisNode)  ub1 = BB-1

C	  receive
	  srcPE1 = (minLoc+incrm)/BB
          if (srcPE1 .NE. thisNode) then
	    anyTrue = .FALSE.
            s1 = MAX(0,lb1)
            s2 = MIN(MOD(-1+NN-incrm,BB),ub1)
            do s = s1, s2
              if (MOD(s+minLoc,incrm2) .LE. incrm-1) then
	        anyTrue = .TRUE.
              end if
	    end do
	    if (anyTrue) then

               call crecv(gtype,buf1,cmplxSize*BB)
	    endif
            offset1 = (srcPE1-thisNode)*BB
          endif

          srcPE2 = (maxLoc+incrm)/BB
          if ((srcPE2 .NE. thisNode) .AND. (srcPE2 .NE. srcPE1)) then
	    anyTrue = .FALSE.
            s1 = MAX(MOD(NN-incrm,BB),lb1)
            s2 = MIN(BB-1,ub1)
            do s = s1, s2
              if (MOD(s+minLoc,incrm2) .LE. incrm-1) then
	        anyTrue = .TRUE.
              end if
	    end do
	    if (anyTrue) then

	      call crecv(gtype,buf2,cmplxSize*BB)
	    endif
            offset2 = (srcPE2-thisNode)*BB
          endif

C	  compute
          tmp1 = minLoc + incrm
           do s = lb1, ub1
            if (MOD(s+minLoc,incrm2) .LE. incrm-1) then
	      if ((s+incrm) .LE. BB-1) then
		x(s) = x(s) + term2(s+incrm)
	      else
                if (INT((s+tmp1)/BB) .EQ. srcPE1) then
  		  x(s) = x(s) + buf1(s+incrm-offset1)
                else
  		  x(s) = x(s) + buf2(s+incrm-offset2)
                endif
	      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
