Efficient 3D Inverse Fourier Transform on FORTRAN Code with Chi Array

[Ben Wilson]

Hi, I have a FORTRAN code with an array called Chi that I want to run an inverse FT on. I have defined two spaces X and K which each consist of 3 vectors running across my physical verse and inverse space.

My code (If it works??) is extremely slow and inefficient (see below). What is the best way of doing a 3D IFT computationally?

CODE

!Real & Inverse Spaces
do i1 = 1, 2*No_x +1
x_x(i1) = -L/2.d0 + dx*(i1-1)
x_y(i1) = -L/2.d0 + dx*(i1-1)
x_z(i1) = -L/2.d0 + dx*(i1-1)
xi_x(i1) = (i1-(No_k+1))*(2.d0*pi/L)
xi_y(i1) = (i1-(No_k+1))*(2.d0*pi/L)
xi_z(i1) = (i1-(No_k+1))*(2.d0*pi/L)
end do

chi = blablabla
!chi is a function of xi

!3d IFT
do r1 = 1,2*No_x +1
do r2 = 1,2*No_x +1
do r3 = 1,2*No_x +1
chi_ift(r1,r2,r3) = 0
do i1 = 1,2*No_k +1
do i2 = 1,2*No_k +1
do i3 = 1,2*No_k +1
chi_ift(r1,r2,r3) = chi_ift(r1,r2,r3) + chi(i1,i2,i3)*exp(imagi*xi_z(i3)*x_z(r3))
end do
chi_ift(r1,r2,r3) = chi_ift(r1,r2,r3) + chi(i1,i2,i3)*exp(imagi*xi_y(i2)*x_y(r2))
end do
chi_ift(r1,r2,r3) = chi_ift(r1,r2,r3) + chi(i1,i2,i3)*exp(imagi*xi_x(i1)*x_x(r1))
end do
end do
end do
end do
Reply Reply to All Forward More
Click to Reply, Reply all or Forward

 

[Dale]

Definitely you should use a standard library like FFTPACK. Unless it is a classroom exercise you should not be rewriting FFT algorithms.

 

[gsal]

Well, for popular algorithms, what Dale said.

Other than that, you need to learn to manipulate arrays en masse, a-la-matlab, whenever possible. This tells the compiler that assignments can be done in any order and without having to carefully traverse one index at a time, or switching from evaluating one array and then another and loop again.

You initialization do-loop above, can be made at least 3 times faster by not having to repeat the calculations; an even faster if you forget the loop altogether:

n=2*No_x+1
iarr(1:n) = (/ (i1-1, i1=1,n) /)

x_x(1:n) = -L/2.0 + dx*iarr(1:n)
x_y(1:n) = x_x(1:n)
x_Z(1:n) = x_x(1:n)

xi_x(1:n) = ( iarr(1:n) - No_k )*2.0*pi/L
xi_y(1:n) = xi_x(1:n)
xi_z(1:n) = xi_x(1:n)

 

[Ben Wilson]

Well, for popular algorithms, what Dale said.

Other than that, you need to learn to manipulate arrays en masse, a-la-matlab, whenever possible. This tells the compiler that assignments can be done in any order and without having to carefully traverse one index at a time, or switching from evaluating one array and then another and loop again.

You initialization do-loop above, can be made at least 3 times faster by not having to repeat the calculations; an even faster if you forget the loop altogether:

n=2*No_x+1
iarr(1:n) = (/ (i1-1, i1=1,n) /)

x_x(1:n) = -L/2.0 + dx*iarr(1:n)
x_y(1:n) = x_x(1:n)
x_Z(1:n) = x_x(1:n)

xi_x(1:n) = ( iarr(1:n) - No_k )*2.0*pi/L
xi_y(1:n) = xi_x(1:n)
xi_z(1:n) = xi_x(1:n)

i'll give this a whirl

 

[Ben Wilson]

Definitely you should use a standard library like FFTPACK. Unless it is a classroom exercise you should not be rewriting FFT algorithms.

hey I'm a really crummy programmer. Which subroutine do I use and how do I implement it (given I'm doing a 3d complex to real transform)?