ShibaFFT とは、Fortran で記述された簡単な混合基数 FFT ライブラリです。 離散フーリエ変換の系列長を 2 のべき乗に限る必要はありません。 系列長の素因数分解に、3,5,7,11,... など、2 以外の素数が含まれていても FFT できます。(ただし、大きな素数が含まれると計算量は O(N^2) になります)
ShibaFFT はシンプルです。1200 行程度の１つのファイルで実装されています。 そのため、提供する機能も１次元の複素離散フーリエ変換のみです。潔いでしょ(^^; Fortran も FORTRAN77 ではなく Fortran90 を採用していますので、読みやすく、 メンテナンスしやすいライブラリになっています。
ShibaFFT は、Fortran で書かれ、OpenMP で並列化されている、 そこそこのパフォーマンスの FFT ライブラリです。 比較的大きな 2 のべき乗の系列長の変換で、 ShibaFFT は、FFTW の FFTW_ESTIMATE におけるパフォーマンスに迫ります。
ShibaFFT は「しば犬 FFT」です。その昔、作者がしば犬を飼っていたからです。え？ 何の意味があるのかって？何の意味もありません。ソフトウェアの名前なんて、 そんなもんです。
ShibaFFT は こちらのページ からダウンロードできます。ライセンスは MIT ライセンスです。
ShibaFFT のパフォーマンスを調べるために、いくつかベンチマークを行いました。 まずは、シングルスレッドでのベンチマークを示しましょう。比較対象は、 FFTPACK5.1D、大浦さんの FFT ライブラリ、FFTE-6.0 です。実行環境は MacBook Pro 15 inch Late 2013(Haswell Core i7 2GHz, 4 core 8 thread) 上の OS X El Capitan / gfortran 5.3 です。OpenMP は Intel OpenMP Runtime を使いました。 浮動小数点数の精度は倍精度です。
------+-----------+-----------------+-----------------+----------------- length|FFTPACK[us]| OOURA [us]| FFTE ST [us]| ShibaFFT [us] ------+-----------+-----------------+-----------------+----------------- 2^( 1)| 0.09| 0.01( 11%)| 0.04( 44%)| 0.02( 22%) 2^( 2)| 0.10| 0.03( 30%)| 0.07( 70%)| 0.04( 40%) 2^( 3)| 0.20| 0.07( 35%)| 0.11( 55%)| 0.10( 50%) 2^( 4)| 0.27| 0.13( 48%)| 0.20( 74%)| 0.19( 70%) 2^( 5)| 0.51| 0.37( 72%)| 0.43( 84%)| 0.45( 88%) 2^( 6)| 0.89| 0.84( 94%)| 0.90(101%)| 0.92(103%) 2^( 7)| 1.89| 1.85( 97%)| 1.97(104%)| 1.96(103%) 2^( 8)| 3.83| 4.08(106%)| 4.46(116%)| 4.44(115%) 2^( 9)| 8.62| 8.83(102%)| 9.84(114%)| 9.87(114%) 2^(10)| 18.24| 19.21(105%)| 21.64(118%)| 21.26(116%) 2^(11)| 47.10| 41.56( 88%)| 52.41(111%)| 49.37(104%) 2^(12)| 100.76| 91.23( 90%)| 110.67(109%)| 109.01(108%) 2^(13)| 222.05| 197.44( 88%)| 243.96(109%)| 243.16(109%) 2^(14)| 488.15| 436.36( 89%)| 596.74(122%)| 531.05(108%) 2^(15)| 1161.46| 919.60( 79%)| 1368.29(117%)| 1209.96(104%) 2^(16)| 2399.87| 2002.60( 83%)| 3039.06(126%)| 2371.88( 98%) 2^(17)| 5094.79| 4328.91( 84%)| 6330.21(124%)| 4982.03( 97%) 2^(18)| 12127.08| 10064.58( 82%)| 13845.31(114%)| 12366.67(101%) 2^(19)| 26975.00| 22770.83( 84%)| 28359.38(105%)| 26206.25( 97%) 2^(20)| 50731.25| 51960.42(102%)| 60627.08(119%)| 54225.00(106%) 2^(21)| 105979.17| 103712.50( 97%)| 131720.83(124%)| 111450.00(105%) 2^(22)| 237116.67| 207466.67( 87%)| 316025.00(133%)| 244700.00(103%) 2^(23)| 708216.67| 465833.33( 65%)| 670633.33( 94%)| 732783.33(103%) 2^(24)| 1572366.67| 1018000.00( 64%)| 1441100.00( 91%)| 1612566.67(102%) ------+-----------+-----------------+-----------------+----------------- FFTE ST = FFTE with Single Thread
length の欄は FFT の系列長を表しています。FFTPACK の欄が FFTPACK5.1D の倍精度版での実行時間です。DFT と IDFT の実行時間の和を表示しています。 OOURA の欄は大浦さんの FFT ライブラリ fftsg.f での実行時間です。FFTE ST の欄は FFTE-6.0 での実行時間です。そして ShibaFFT の欄が ShibaFFT-1.2 での実行時間です。実行時間の単位はマイクロ秒で、初期化の時間は含まれていません。 括弧の中の％付きの数字は FFTPACK の実行時間に対して何パーセントの実行時間になるかを表しています。
大浦さんのライブラリがブッチギリの性能ですね。これはもう、 シングルスレッドなら大浦さんのライブラリを使うしかありません。ShibaFFT は、 まあ、そこそこ善戦してます(^^; 最悪でも FFTPACK の 120％ ほどの処理時間に収まっています。
続いて、マルチスレッドでのベンチマークです。FFTPACK、大浦さんのライブラリは、 シングルスレッドしかサポートしてないので、シングルスレッドのままです。FFTE は OpenMP を有効にすると、系列長が 2^16 で落ちてしまうので(私の環境だけ？)、 代わりに FFT_OPENMP と比較しました。 ShibaFFT も OpenMP によるマルチスレッドバージョンを調べました。 小さな系列長では、FFT_OPENMP の計測が困難だったので 2^10 から計測開始してます。
------+-----------+-----------------+-----------------+----------------- length|FFTPACK[us]| OOURA [us]| FFT_OPENMP[us]| ShibaFFT MT [us] ------+-----------+-----------------+-----------------+----------------- 2^(10)| 18.24| 19.21(105%)| 61.10(334%)| 18.75(102%) 2^(11)| 47.10| 41.56( 88%)| 85.73(182%)| 33.24( 70%) 2^(12)| 100.76| 91.23( 90%)| 150.60(149%)| 51.94( 51%) 2^(13)| 222.05| 197.44( 88%)| 232.23(104%)| 102.52( 46%) 2^(14)| 488.15| 436.36( 89%)| 458.30( 93%)| 194.56( 39%) 2^(15)| 1161.46| 919.60( 79%)| 790.49( 68%)| 439.39( 37%) 2^(16)| 2399.87| 2002.60( 83%)| 1685.81( 70%)| 870.44( 36%) 2^(17)| 5094.79| 4328.91( 84%)| 3282.29( 64%)| 2021.61( 39%) 2^(18)| 12127.08| 10064.58( 82%)| 8893.75( 73%)| 4364.06( 35%) 2^(19)| 26975.00| 22770.83( 84%)| 20394.79( 75%)| 10990.62( 40%) 2^(20)| 50731.25| 51960.42(102%)| 46854.17( 92%)| 19335.42( 38%) 2^(21)| 105979.17| 103712.50( 97%)| 90029.17( 84%)| 41458.33( 39%) 2^(22)| 237116.67| 207466.67( 87%)| 289383.33(122%)| 110700.00( 46%) 2^(23)| 708216.67| 465833.33( 65%)| 1005183.33(141%)| 505700.00( 71%) 2^(24)| 1572366.67| 1018000.00( 64%)| 2334900.00(148%)| 1110933.33( 70%) ------+-----------+-----------------+-----------------+-----------------
うん、ShibaFFT がブッチギリですね。まあ、FFTPACK がシングルスレッドなんで、 当たり前なんですが。それにしても、系列長 2^23, 2^24 の大浦さんのライブラリはすごいですね。シングルスレッドなのに、マルチスレッドの ShibaFFT より速いとは。FFTPACK も 2^10 は神がかりですね。
最後は、2 のべき乗以外の系列長でのベンチマークです。とりあえず 10 のべき乗を調べてみました。大浦さんのライブラリは、 2 のべき乗しかサポートしないので、外しています。
------+-----------+-----------------+-----------------+----------------- length|FFTPACK[us]| FFTE ST [us]| ShibaFFT [us]| ShibaFFT MT [us] ------+-----------+-----------------+-----------------+----------------- 10^(1)| 0.22| 0.17( 77%)| 0.20( 90%)| 0.21( 95%) 10^(2)| 1.58| 2.11(133%)| 2.79(176%)| 2.78(175%) 10^(3)| 21.99| 29.83(135%)| 39.34(178%)| 39.36(178%) 10^(4)| 286.60| 402.63(140%)| 504.97(176%)| 162.50( 56%) 10^(5)| 3634.33| 5924.33(163%)| 5750.00(158%)| 1874.67( 51%) 10^(6)| 43130.00| 77300.00(179%)| 67766.67(157%)| 22933.33( 53%) 10^(7)| 751966.67| 926766.67(123%)| 998800.00(132%)| 554100.00( 73%) ------+-----------+-----------------+-----------------+----------------- FFTE ST = FFTE with Single Thread
実は、ShibaFFT は、系列長 N が N = (2^p)*(3^q)*(5^r) の形になる時だけ速くなるチートが仕込んであります。なので、10 のべき乗の時も、 そこそこの性能を示します。しかし、系列長の素因数分解に、 7 以上の素数が含まれるとガクンと性能が劣化します。 2,3,5 には専用のコードが書かれているのですが、 7 以上は素朴で汎用的なコードが使われるからです。
さて、ShibaFFT の配布物には上記のベンチマークを行うコードが含まれていますが、 当然、FFTPACK、大浦さんのライブラリ、FFTE、FFT_OPENMP などのコードは含まれていません。なのでベンチマークプログラムを動かすには、 それらを入手して適切に配置してやる必要があります。
まず、FFTPACK ですが、以下の URL から入手します。
https://people.sc.fsu.edu/~jburkardt/f_src/fftpack5.1d/fftpack5.1d.html
そして、libfftpack5.1d.a を作り、ShibaFFT-1.2 フォルダに配置します。 FFTPACK のビルドはかなりの変態です。 f90split と言うコマンドと専用のシェルスクリプトでビルドします。 方法は「fftpack インストール f90split」とかで検索してみてください。
次に大浦さんのライブラリですが、以下の URL から fft.tgz を入手します。
http://www.kurims.kyoto-u.ac.jp/~ooura/fft-j.html
そして、含まれる fftsg.f を ShibaFFT-1.2 フォルダに配置します。
次に、FFTE ですが、以下の URL から ffte-6.0.tgz を入手します。
http://www.ffte.jp/
そして、ShibaFFT-1.2 フォルダで展開します。
次に FFT_OPENMPですが、以下の URL から fft_openmp.f90 ファイルを入手します。
https://people.sc.fsu.edu/~jburkardt/f_src/fft_openmp/fft_openmp.html
そして、ShibaFFT-1.2 フォルダに配置します。あとは、
make benchmark
make benchmark-mt
make benchmark10
で、それぞれシングルスレッドベンチマーク、マルチスレッドベンチマーク、系列長が 10 のべき乗のベンチマークが行われます。ruby も使っているので、 デフォルトで無い場合はインストールしといてくださいね。
ShibaFFT を使うには準備が必要です。ShibaFFT-1.2.tar.gz を展開すると、ShibaFFT-1.2 フォルダができます。その中で必要なファイルは ShibaFFT フォルダに入っています。 ShibaFFT フォルダに入って以下のコマンドを実行すると、
make module
OpenMP に対するチューニングが行われたモジュールが生成されます。 Mac OS X で Intel OpenMP Runtime を使う場合は、Makefile.osx の
OMPLIB = -lgomp
#OMPLIB = -L$(HOME)/lib -liomp5
を
#OMPLIB = -lgomp
OMPLIB = -L$(HOME)/lib -liomp5
にしてから実行してください。もちろん、Intel OpenMP Runtime が導入されている必要があります。$HOME/lib に libiomp5.dylib があるとします。 (こちらのページ で導入法が紹介されています)
ShibaFFT の使い方をサンプルコードで見てみましょう。まずは、use で shibafft モジュールを読み込みます。n を処理する系列のサイズとして、allocate で領域を確保します。W の範囲が 0 から n までであることに注意してください。 つまり、W のサイズは n + 1 です。
続いて、shiba_fft_init で W を初期化します。これで変換の準備は整いました。 shiba_fft_fwd で離散フーリエ変換を実行します。x が入力でかつ出力です。 y は作業領域として使われます。
program hello
use shibafft ! ShibaFFT モジュールを使う
implicit none
integer, parameter :: n = 2**10 ! n は離散フーリエ変換の系列長
integer, parameter :: dbl = kind(0d0)
complex(dbl), allocatable :: x(:), y(:), W(:)
! ...
allocate(x(0:n-1), y(0:n-1), W(0:n))
! ...
call shiba_fft_init(n, W) ! DFT の重み W を初期化。
call shiba_fft_fwd(n, x, y, W) ! FFT を実行。x は入出力、y は作業領域
! ...
deallocate(x, y, W)
! ...
end program
カレントフォルダに ShibaFFT フォルダがある場合、 コンパイルするには以下のようにします。
gfortran -O3 -IShibaFFT -c hello.f90
gfortran hello.o ShibaFFT/shibafft.o -lgomp -o hello
あるいは、Mac OS X で Intel OpenMP Runtime を使うなら以下のようにします。
gfortran-mp-5 -O3 -IShibaFFT -c hello.f90
gfortran-mp-5 hello.o ShibaFFT/shibafft.o -L$HOME/lib -liomp5 -o hello
シングルスレッドの場合でも、ShibaFFT 自身は OpenMP ありでコンパイルされているので -lgomp あるいは -liomp5 オプションが必要です。 それが気持ち悪い場合は、ShibaFFT フォルダ内の Makefile.* から -fopenmp を外して make module してください。
以下のように宣言されているとします。
integer, parameter :: n = <離散フーリエ変換の系列長>
complex(kind(0d0)) x(0:n-1), y(0:n-1), W(0:n)
次のサブルーチンが提供されています。
＜初期化ルーチン＞
shiba_fft_init(n, W)
フーリエ変換の重み W を初期化します。変換を行う前に呼び出しておく必要が
あります。
＜シングルスレッドルーチン＞
shiba_fft_fwd(n, x, y, W)
x を離散フーリエ変換します。x が入力でかつ出力です。y は作業領域です。
変換結果は 1/n で正規化されます。
shiba_fft_fwd0(n, x, y, W)
x を離散フーリエ変換します。x が入力でかつ出力です。y は作業領域です。
変換結果は正規化されません。
shiba_fft_fwdu(n, x, y, W)
x を離散フーリエ変換します。x が入力でかつ出力です。y は作業領域です。
変換結果は sqrt(1d0/n) で正規化されます。
shiba_fft_fwdn(n, x, y, W)
x を離散フーリエ変換します。x が入力でかつ出力です。y は作業領域です。
変換結果は 1/n で正規化されます。
shiba_fft_inv(n, x, y, W)
x を逆離散フーリエ変換します。x が入力でかつ出力です。y は作業領域です。
変換結果は正規化されません。
shiba_fft_inv0(n, x, y, W)
x を逆離散フーリエ変換します。x が入力でかつ出力です。y は作業領域です。
変換結果は正規化されません。
shiba_fft_invu(n, x, y, W)
x を逆離散フーリエ変換します。x が入力でかつ出力です。y は作業領域です。
変換結果は sqrt(1d0/n) で正規化されます。
shiba_fft_invn(n, x, y, W)
x を逆離散フーリエ変換します。x が入力でかつ出力です。y は作業領域です。
変換結果は 1/n で正規化されます。
＜マルチスレッドルーチン＞
shiba_pfft_fwd(n, x, y, W)
x を離散フーリエ変換します。x が入力でかつ出力です。y は作業領域です。
変換結果は 1/n で正規化されます。
shiba_pfft_fwd0(n, x, y, W)
x を離散フーリエ変換します。x が入力でかつ出力です。y は作業領域です。
変換結果は正規化されません。
shiba_pfft_fwdu(n, x, y, W)
x を離散フーリエ変換します。x が入力でかつ出力です。y は作業領域です。
変換結果は sqrt(1d0/n) で正規化されます。
shiba_pfft_fwdn(n, x, y, W)
x を離散フーリエ変換します。x が入力でかつ出力です。y は作業領域です。
変換結果は 1/n で正規化されます。
shiba_pfft_inv(n, x, y, W)
x を逆離散フーリエ変換します。x が入力でかつ出力です。y は作業領域です。
変換結果は正規化されません。
shiba_pfft_inv0(n, x, y, W)
x を逆離散フーリエ変換します。x が入力でかつ出力です。y は作業領域です。
変換結果は正規化されません。
shiba_pfft_invu(n, x, y, W)
x を逆離散フーリエ変換します。x が入力でかつ出力です。y は作業領域です。
変換結果は sqrt(1d0/n) で正規化されています。
shiba_pfft_invn(n, x, y, W)
x を逆離散フーリエ変換します。x が入力でかつ出力です。y は作業領域です。
変換結果は 1/n で正規化されています。
ShibaFFT の作者 OK おじさんは、C++ で実装された高速な FFT ライブラリである OTFFT も開発しています。ShibaFFT は OTFFT を Fortran で実装してみたらどうなるだろうかという興味から始まりました。始めて見ると、 これがなかなかに手強い。C++ で使えていた高速化技法の多くが、 Fortran では使えなかったり意味がなかったりしたからです。
例えば、Fortran ではプログラマが陽に AVX を使うことができません。 コンパイルオプションや OpenMP で間接的に使うことを指示することはできますが、 思ったような効果は得られませんでした。
そこで、方針を転換し、そこそこのパフォーマンスを維持しつつ、簡単な実装の FFT ライブラリを作ることにしました。ShibaFFT は、 2 のべき乗の系列長のパフォーマンスは、FFTPACK や FFTE と言った有名どころと比べても遜色ありません。何より実装が簡単なので、混合基数 FFT を勉強したい人の教材として最適です。
そんなわけで、興味がわいたら使ってみてください。