FFTWで３次元高速フーリエ変換してみた。

前回の例を元に、3次元に拡張してみた。
koideforest.hatenadiary.com

program fftw3_test_3d
implicit none

integer,       parameter :: WP = selected_real_kind( 8 )
real( WP ),    parameter :: PI = 3.141592653589793
complex( WP ), parameter :: I  = ( 0, 1 )

integer :: i1, i2, i3


! setting variables

integer,     parameter :: N = 2**7
real( WP ) , parameter :: L = 100.

real( WP )    :: dx
complex( WP ) :: temp_c

real( WP ), dimension( 3 ) :: k0, temp_rv
real( WP ), dimension( N ) :: x, k


! data

complex( WP ), dimension( N, N, N ) :: input0, input, output


! plan

integer( WP ) :: plan


include 'fftw3.f'

dx = L / N

do i1 = 1, N
  x( i1 ) = - L/2 + ( i1 - 1 ) * dx
enddo

! Gaussian wave packet
!k0 = (/ 0, 0, 5 /)
k0 = (/ 0, 0, 0 /)
do i3 = 1, N
  do i2 = 1, N
    do i1 = 1, N
      temp_rv = (/ x( i1 ), x( i2 ), x( i3 ) /)
      input( i1, i2, i3 ) = &
        &  exp( - dot_product( temp_rv, temp_rv ) + I * dot_product( k0, temp_rv ) )
    enddo
  enddo
enddo

input0( :, :, : ) = input( :, :, : )

call dfftw_plan_dft_3d( plan, N, N, N, input, output, &
  &                     FFTW_FORWARD, FFTW_ESTIMATE )
call dfftw_execute_dft( plan, input, output )
call dfftw_destroy_plan( plan )

do i1 = 1, N/2
  k( i1             ) =   i1 - 1
  k( N - ( i1 - 1 ) ) = - i1
enddo
k( : ) = ( 2 * PI / L ) *  k( : )

open( 10, file = 'fftw3_forward.dat', status = 'unknown' )
do i3 = 1, N
  i2 = 1  ! ky = 0
  do i1 = 1, N
    temp_c = output( i1, i2, i3 )
    write( 10, '(3e17.5e3)' ) k( i1 ), k( i3 ), abs( temp_c )
  enddo
enddo
close( 10 )

call dfftw_plan_dft_3d( plan, N, N, N, output, input, &
  &                     FFTW_BACKWARD, FFTW_ESTIMATE )
call dfftw_execute_dft( plan, output, input )
call dfftw_destroy_plan( plan )

input( :, :, : ) = input( :, :, : ) / (N**3)

open( 10, file = 'fftw3_backward.dat', status = 'unknown' )
do i3 = 1, N
  i2 = N/2  ! y =0
  do i1 = 1, N
    temp_c = input( i1, i2, i3 ) - input0( i1, i2, i3 )
    write( 10, '(3e17.5e3)' ) x( i1 ), x( i3 ), abs( temp_c )
  enddo
enddo
close( 10 )

end program fftw3_test_3d

プロットはどこか二次元平面を選ばないと出来ないので、逆空間では、実空間ではとなる平面を選択して出力。
で、出力結果が 3 MB 程度だったが、自分のMacBook Airではでのコンパイルでエラー（多分メモリが足りない）が返って来たので、三次元はやはり重たい。