# 2D Discrete Fourier Transform using Fortran

I am solving the 2D Wave Equation using Fourier Transform. The Discrete Fourier Transform (and the inverse also) is done inside the kx-loop and ky-loop. But this code runs slow, is there anyway to make it much more efficient? I guess the kx-loop, ky-loop inside the i-loop and j-loop makes it slow. I am also open for external package suggestion.

program fourier_transform_solution
implicit none

integer :: i,j,kx, ky, it
real, parameter :: PI = 3.1415926
integer, parameter :: N=100, N_t=100
integer, parameter :: N_k = N/2
real, dimension(N):: x, y
real, dimension(N,N) :: f, ut_init, u_r
real, dimension(N_t) :: t
complex, dimension(N,N) :: u, u_hat, u_hat_prev, ut_hat, ut_hat_prev
complex, dimension(N) :: wx, wy
complex, parameter :: I_complex = complex(0, 1)
real :: dx, xmin=0, xmax=3*PI, dy, ymin = 0, ymax = 3*PI
real :: d, dt = 0.001

do i=1,N_t,1
t(i) = 0 + (i-1)*dt
end do

dx = (xmax-xmin)/(N-1); dy = (ymax-ymin)/(N-1); d = dx/(2*PI);

do i=1,N_k,1
wx(i) = 0 + -(N_k-i)/(d*(N-1)); wx(N_k+i) = 0 + (i)/(d*(N-1))
wy(i) = 0 + -(N_k-i)/(d*(N-1)); wy(N_k+i) = 0 + (i)/(d*(N-1))
end do

do i=1,N,1
x(i) = xmin+ (i-1)*dx
y(i) = ymin+ (i-1)*dy
end do

do i=1,N,1
do j=1,N,1
f(i,j) = 0.6*exp( -((x(i)- 1.5*PI)**2) + -((y(j)- 1.5*PI)**2)  )
u(i,j) = f(i,j)
u_r(i,j) = f(i,j)
ut_init(i,j) = 0
end do
end do

!Fourier Transform
do i=1,(N),1
do j=1,N,1
u_hat_prev(i,j) = 0
ut_hat_prev(i,j) = 0
do kx=1,N,1
do ky=1,N,1
u_hat_prev(i,j) = u_hat_prev(i,j) + u_r(kx, ky)*exp(-(wx(i)*x(kx) + wy(j)*y(ky))*(I_complex))
ut_hat_prev(i,j) = ut_hat_prev(i, j) + ut_init(kx, ky)*exp(-(wx(i)*x(kx) + wy(j)*y(ky))*(I_complex))
end do
end do
end do
end do

do it=2,N_t,1

do i=1,(N),1
do j=1, N, 1
ut_hat(i, j) = ut_hat_prev(i, j) + dt*(-(wx(i) + wy(j))**2)*(u_hat_prev(i, j))
end do
end do

do i=1,(N),1
do j=1,N,1
u_hat(i,j) = u_hat_prev(i,j) + dt*ut_hat(i,j)
end do
end do

!Inverse Fourier
do i=1,N,1
do j = 1,N,1
u(i, j)=0
do kx=1,N,1
do ky=1,N,1
u(i,j)  = u(i,j) + u_hat(kx, ky)*exp((I_complex)*(wx(kx)*x(i) + wy(ky)*y(j)))
end do
end do
u_r(i,j) = (real(u(i,j)))/(N*N)
ut_hat_prev(i,j) = ut_hat(i,j)
u_hat_prev(i,j) = u_hat(i,j)
end do
end do
print *, it
end do

end program fourier_transform_solution
$$$$
`
• You should never ever compute the discrete Fourier transform from the definition. Instead, use FFT (fast Fourier transform). There are packages like P3DFFT, FFTW3, FFTPACK and some others available. Aug 12, 2021 at 7:11
• @VladimirF which one you recommend? I downloaded the intel one but dont mnow Aug 14, 2021 at 4:18