FFTW (2-D Serial) (Complex to Real to Complex)

!

! Discussion:

!

! In FORTRAN it is the FIRST dimension

! of the complex coefficient array that is "half" the size.

!

implicit none

integer ( kind = 4 ), parameter :: Nx = 5

integer ( kind = 4 ), parameter :: Ny = 3

integer ( kind = 4 ), parameter :: Nh = ( Nx / 2 ) + 1

include "fftw3.f"

integer ( kind = 4 ) i,j

real ( kind = 8 ) ux(Nx,Ny)

complex ( kind = 8 ) ukx(Nh,Ny),ukx_dum(Nh,Ny)

integer ( kind = 8 ) plan_forward,plan_backward

integer,parameter :: seed = 99999999

call srand(seed)

write ( *, '(a)' ) ' '

write ( *, '(a)' ) 'TEST04'

write ( *, '(a)' ) ' Demonstrate FFTW3 on a 2D real array'

write ( *, '(a,i8)' ) ' NX = ', nx

write ( *, '(a,i8)' ) ' NY = ', ny

write ( *, '(a)' ) ' '

write ( *, '(a)' ) ' Transform data to FFT coefficients.'

write ( *, '(a)' ) ' Backtransform FFT coefficients to recover'

write ( *, '(a)' ) ' the data.'

write ( *, '(a)' ) ' Compare recovered data to original data.'

!

! Compute the data.

!

do i = 1,Nh

do j = 1,Ny

ukx(i,j) = (1.0d0,0.0d0)*rand() + (0.0d0,1.0d0)*rand()

enddo

write ( *, '(a)' ) ' '

write ( *, '(a)' ) ' Input Data:'

write ( *, '(a)' ) ' '

do i = 1,Nh

do j = 1,Ny

write(*, *) i,j,ukx(i,j), "ukx"

enddo

!

! Make a plan for the FFT, and backward transform the data.

!

call dfftw_plan_dft_c2r_2d_ (plan_backward, Nx, Ny, ukx, ux, FFTW_ESTIMATE)

call dfftw_execute_ (plan_backward)

write ( *, '(a)' ) ' '

write ( *, '(a)' ) ' Output FFT Coefficients:'

write ( *, '(a)' ) ' '

do i = 1,Nx

do j = 1,Ny

ux(i,j) = ux(i,j)/(float(Nx)*float(Ny))

write(*, *) i,j,ux(i,j), "ux"

enddo

!

! Make a plan for the backward FFT, and recover the original data.

!

call dfftw_plan_dft_r2c_2d_ (plan_forward, Nx, Ny, ux, ukx_dum, FFTW_ESTIMATE)

call dfftw_execute_ (plan_forward)

write ( *, '(a)' ) ' '

write ( *, '(a)' ) ' Recovered input data divided by NX * NY:'

write ( *, '(a)' ) ' '

do i = 1,Nh

do j = 1,Ny

write(*, *) i,j,ukx_dum(i,j), "ukx_dum"

enddo

!

! Discard the information associated with the plans.

!

call dfftw_destroy_plan_ ( plan_forward )

call dfftw_destroy_plan_ ( plan_backward )

! return

end