!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  FILE: ser_fft.f
!  DESCRIPTION: see mpi_2dfft.f
!  AUTHOR: William Smith - adapted from George Gusciora's MPI C version
!  LAST REVISED: 
!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  fft -  n-point in-place decimation-in-time FFT of complex vector 
!         "data" using the n/2 complex twiddle factors in "w_common".
!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
      subroutine fft(data, w_common, n, logn)

      complex  data(*), w_common(*)
      integer  n, logn

      complex  temp
      integer  incrvec, i0, i1, i2, nx, ip2
      integer  ivec2
    
! bit-reverse the input vector 

      call bit_reverse (data, n)


! do the first logn-1 stages of the fft */

      i2 = logn
      incrvec = 1

      do 
      ivec2 = incrvec 
      incrvec = incrvec * 2
      if (incrvec .ge. n) exit			! do while (incrvec .lt. n-1)

          i2 = i2 - 1
          do i0 = 0, ivec2-1
              do i1 = 0, n-1, incrvec
                  ip2 = 2 ** i2

                  temp = data(i0+i1 + ivec2 + 1) * w_common(i0*ip2 + 1)
                  data(i0+i1 + ivec2 + 1) = data(i0+i1 + 1) - temp
                  data(i0+i1 + 1) = data(i0+i1 + 1) + temp
              enddo
          enddo
      enddo
      
!
! do the last stage of the fft
!
      do i0 = 0, n/2-1
          temp = data(i0 + n/2 + 1) * w_common(i0 + 1)
          data(i0 + n/2 + 1) = data(i0 + 1) - temp
          data(i0 + 1) = data(i0 + 1) + temp
      enddo

      end

      
!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! bit_reverse - simple (but somewhat inefficient) bit reverse 
!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

      subroutine bit_reverse (a, n)
      complex  a(*)
      integer  n
      
      complex  temp
      integer  i, j, k
      
      j = 0
      do i=0, n-3

          if (i .lt. j) then
                temp = a(j+1)
                a(j+1) = a(i+1)
                a(i+1) = temp
          endif
          k = n / 2

          do while (k .le. j) 
            j = j - k 
            k = k / 2
          enddo
          j = j + k
      enddo

      end