!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! Serial - Two-Dimensional Fast Fourier Transform - Fortran Version
! FILE:  ser_2dfft.f
! OTHER FILES:   ser_fft.f  timing_fgettod.c
! DESCRIPTION: Algorithm: Perform a 1D FFT on each row, transpose, then perform a 1D
!   FFT on each row again.  
!
!   Input is an 512x512 complex matrix. Output is an 512x512
!   complex matrix that overwrites the input matrix.  The input matrix is
!   initialized with a point source.  
! 
!   A straightforward unsophisticated 1D FFT kernel is used.  It is
!   sufficient to convey the general idea, but be aware that there are
!   better 1D FFTs available on many systems.
! 
! AUTHOR: William Smith - adapted from C version by George Gusciora
! LAST REVISED: 6/9/96 Blaise Barney
!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

      program    ser_2dfft

      external   timing_fgettod

      integer    MAXN
      real       PI
      parameter  (MAXN = 512)
      parameter  (PI = 3.141592)

      complex    a(MAXN,MAXN)		! input matrix 
      complex    w_common(MAXN/2)	! twiddle factors 
      complex    temp			! temp complex for swapping

      integer    n           		! FFT size 
      integer    logn        		! log base 2 of n 
      integer    nx          		! used to compute logn 
      real       rn    			! real value of n
      integer    i,j         		! index variables 

      integer    errors,sign 		! used for error checking 

      integer    times1(2), times2(2)	! times from timing_fgettod
      integer    fsecs    		! time spent doing 1D FFT's 
      integer    tsecs   		! time spent doing the transpose 
      integer    secs     		! total time of the 2D FFT 
      integer    flops       		! total number of floating point ops 
      real       mflops    		! Mflops/s 
      real       rfsecs    		! float fft seconds
      real       rtsecs    		! float transpose seconds
      real       rsecs    		! float total seconds

!
! initialize the input matrix with a centered point source
!
      n = 512
      logn = 9

      do i=1, n
          do j=1, n
              a(i,j) = (0.0, 0.0)
          enddo
      enddo

      rn = REAL(n)
      a(n/2+1,n/2+1) = (rn, rn)



!
! precompute the complex constants (twiddle factors) for the 1D FFTs */
!
      do i = 1, n/2-1
          j = i - 1
          w_common(i) = (cos((2.0*PI*j)/rn), -sin((2.0*PI*j)/rn))
      enddo

  
!* Get first time - seconds and microseconds
      call timing_fgettod(%REF(times1))

!* 
! first do a set of row FFTs
!*

      do i = 1, n
          call fft(a(1,i), w_common, n, logn)
      enddo

!* Get ending time - in seconds and microseconds
!* Calculate time for 1D ffts
      call timing_fgettod(%REF(times2))

      fsecs = ((times2(1) - times1(1))*1000000) + times2(2) -
     &            times1(2)




!* Get first time - seconds and microseconds
      call timing_fgettod(%REF(times1))

!*
! then transpose the matrix
!*
      
      do i = 1, n
          do j = i+1, n
              temp   = a(i,j)
              a(i,j) = a(j,i)
              a(j,i) = temp
          enddo
      enddo
!* Get ending time - in seconds and microseconds
!* Calculate time for 1D ffts
      call timing_fgettod(%REF(times2))

      tsecs = ((times2(1) - times1(1))*1000000) + times2(2) -
     &            times1(2)




!* Get first time - seconds and microseconds
      call timing_fgettod(%REF(times1))

!*
! then do another set of row FFTs
!*
      
      do i = 1, n
          call fft(a(1,i), w_common, n, logn)
      enddo

!* Get ending time - in seconds and microseconds
!* Calculate time for 1D ffts
      call timing_fgettod(%REF(times2))

      fsecs = fsecs + ((times2(1) - times1(1))*1000000) +
     &        times2(2) - times1(2)



!* Get first time - seconds and microseconds
      call timing_fgettod(%REF(times1))

!*
! transpose the matrix back
!*
      
      do i = 1, n
          do j = i+1, n
              temp   = a(i,j)
              a(i,j) = a(j,i)
              a(j,i) = temp
          enddo
      enddo

!* Get ending time - in seconds and microseconds
!* Calculate time for 1D ffts
      call timing_fgettod(%REF(times2))

      tsecs = tsecs + ((times2(1) - times1(1))*1000000) +
     &        times2(2) - times1(2)



!*
!* summarize the 2d FFT performance 
!*
      print *, n, 'x', n, '2D FFT'
      secs = fsecs + tsecs
      rfsecs = real(fsecs)
      rtsecs = real(tsecs)
      rsecs  = real(secs)
      flops  = (n*n*logn)*10
      mflops = real(flops/1000000.0)
      mflops = mflops/(real(secs)/1000000.)
      print *, ' '
      print *, '1D FFTs   :', rfsecs/1000000., ' sec'
      print *, 'transpose :', rtsecs/1000000., ' sec'
      print *, 'total     :', rsecs/1000000., ' sec'
      print *, ' '
      print *, mflops, 'Mflops '



      end