fftw3の使い方について¶
留意点¶
使用時に留意するのは3つ
include ファイルの言及
コンパイルオプションの指定
ファイル中のsubroutine呼び出し
コンパイルオプション¶
gfortran test.f90 -I/usr/local/include -lfftw3
が必要．-Iオプションはincludeファイルの存在場所を指定する．-lfftw3はライブラリの指定を行っている．
subroutine 呼び出し(fftw3の使い方)¶
integer, parameter :: N=256
integer(8) :: plan
complex(kind(0d0)) :: y(N), fft(N)
call dfftw_plan_dft_1d( plan, N, y, fft, FFTW_FORWARD, FFTW_ESTIMATE )
call dfftw_execute_dft( plan, y, fft )
call dfftw_destroy_plan( plan )
よくわからんのはFFTW_FORWARDとFFTW_ESTIMATEだけど，これは特に何か変数を用意する必要はなさそうで，指定すれば事足りるみたい．
使うべき関数の注意点¶
fftw3には複数のsubroutineが含まれているけど，そのうち，幾つかは危険な香りがするみたい．例えば，fftw3のreferenceには次のような文言がある．
call dfftw_execute(plan)
However, we have had reports that this causes problems with some recent optimizing Fortran compilers. The problem is, because the input/output arrays are not passed as explicit arguments to dfftw_execute, the semantics of Fortran (unlike C) allow the compiler to assume that the input/output arrays are not changed by dfftw_execute. As a consequence, certain compilers end up optimizing out or repositioning the call to dfftw_execute, assuming incorrectly that it does nothing.
dfftw_execute(plan)ではなく，別の，それぞれ用途にフィットした関数を使うのがよろしいようです．
You must use the correct type of execute function, matching the way the plan was created. Complex DFT plans should use dfftw_execute_dft, Real-input (r2c) DFT plans should use use dfftw_execute_dft_r2c, and real-output (c2r) DFT plans should use dfftw_execute_dft_c2r. The various r2r plans should use dfftw_execute_r2r.
複素Fourier変換(最もgeneral)が欲しければ，
call dfftw_execute_dft( plan, in, out )
を使いましょう．
fftw3 in fortran の作法¶
planにはinteger(8) を使いましょう．
integer(8) :: plan
ref. 8.1 Fortran-interface routines( http://www.fftw.org/doc/Fortran_002dinterface-routines.html#Fortran_002dinterface-routines )
fortran fftw3のサンプル¶
- 抜粋 :
In Fortran, you would use the following to accomplish the same thing:
double complex in, out
dimension in(N), out(N)
integer*8 plan
call dfftw_plan_dft_1d(plan,N,in,out,FFTW_FORWARD,FFTW_ESTIMATE)
call dfftw_execute_dft(plan, in, out)
call dfftw_destroy_plan(plan)
To transform a three-dimensional array in-place .... In Fortran, you would use this instead:
double complex arr
dimension arr(L,M,N)
integer*8 plan
call dfftw_plan_dft_3d(plan, L,M,N, arr,arr &
& ,FFTW_FORWARD, FFTW_ESTIMATE)
call dfftw_execute_dft(plan, arr, arr)
call dfftw_destroy_plan(plan)
Note that we pass the array dimensions in the “natural” order in both C and Fortran.
fortranでの"natural"な配列の渡し方 == (L,M,N)で良いみたい．
自作サンプル結果¶
正弦波(モード数=2)の波を入力してFFTする．
1program main
2
3
4
5 implicit none
6
7 integer, parameter :: N=256
8 double precision, parameter :: pi = 4.d0 * atan( 1.d0 )
9 double precision, parameter :: xmax = 2.d0 * pi, xmin = 0.d0
10 double precision, parameter :: amp = 1.d0
11 double precision, parameter :: k = 2.d0
12
13 integer(8) :: plan
14
15 integer :: i, j
16 double precision :: dx
17 double precision :: x(N)
18 complex(kind(0d0)) :: y(N), fft(N)
19 double precision :: power(N), theta(N)
20
21 ! data preparation
22
23 include 'fftw3.f'
24
25 do i=1, N
26 dx = ( xmax - xmin ) / dble(N)
27 x(i) = dx * dble( i - 1 )
28 y(i) = amp * sin( k * x(i) )
29 enddo
30
31 open( 50, file='x.dat', form='formatted' )
32 do i=1, N
33 write(50,*) x(i), y(i)
34 enddo
35 close(50)
36
37 ! FFT
38
39 call dfftw_plan_dft_1d( plan, N, y, fft, FFTW_FORWARD, FFTW_ESTIMATE )
40 call dfftw_execute_dft( plan, y, fft )
41 call dfftw_destroy_plan( plan )
42
43 power = abs( fft ) / dble(N)
44 theta = atan2( aimag(fft), real(fft) )
45
46 ! write
47
48 open( 50, file='fft.dat', form='formatted' )
49 do i=1, N
50 write(50,*) i, power(i), theta(i)
51 enddo
52 close(50)
53
54
55end program main
結果は以下．
multi-threaded での注意点もあるみたいだ．
