大規模系での高速フーリエ変換1
高橋大介
daisuke@cs.tsukuba.ac.jp

筑波大学計算科学研究センター

2016/5/26

計算科学技術特論B

1

講義内容• 高速フーリエ変換（FFT）• ブロック三次元FFTアルゴリズム• 並列FFTアルゴリズム

2016/5/26

計算科学技術特論B

2

高速フーリエ変換（FFT）• 高速フーリエ変換（FFT）は，離散フーリエ変換（DFT）を高速に計算するアルゴリズム．
• 理学分野での応用例– 偏微分方程式の解法– 畳み込み，相関の計算– 第一原理計算における密度汎関数法• 工学分野での応用例– スペクトラムアナライザ– CTスキャナ・MRIなどの画像処理– 地上波デジタルテレビ放送や無線LANで用いられている
OFDM（直交周波数多重変調）では，変復調処理にFFTを用いている．
2016/5/26

計算科学技術特論B

3

離散フーリエ変換（DFT）• 離散フーリエ変換（DFT）の定義
n −1

y (k ) = ∑ x( j )ωj =0jkn0 ≤ k ≤ n − 1, ωn = e
2016/5/26

計算科学技術特論B− 2πi / n4行列によるDFTの定式化（1/4）• n = 4 のとき，DFTは以下のように計算できる．

y (0) = x(0)ω + x(1)ω + x(2)ω + x(3)ω
000y (1) = x(0)ω + x(1)ω + x(2)ω + x(3)ω
0

1

20

3y (2) = x(0)ω + x(1)ω + x(2)ω + x(3)ω
02

4y (3) = x(0)ω + x(1)ω + x(2)ω + x(3)ω
0

2016/5/26

3

計算科学技術特論B

6695

行列によるDFTの定式化（2/4）• 行列を用いると，より簡単に表すことができる．

 y (0)  ω ω
 y (1)   01ωω= y (2) ω 0 ω 2  0
3 y (3)  ω ω00

ω2
ω4
ω6
ω0

ω   x(0) 
3 ω   x(1) 6ω  x(2)
9 
ω   x(3) 0• n 回の複素数の乗算と，n( n − 1) 回の複素数の加算が必要．
2

2016/5/26

計算科学技術特論B

6

行列によるDFTの定式化（3/4）
•

ω =ω
jk
n

jk mod nの関係を用いると，以下のようnに書き直すことができる．

11   x ( 0)  y (0)  1 1
 y (1)  1 ω 1 ω 2 ω 3   x(1) 
=202
 y (2) 1 ω ω ω   x(2)
 
3
21 
 y (3)  1 ω ω ω   x(3) 
2016/5/26

計算科学技術特論B

7

行列によるDFTの定式化（4/4）• 行列の分解により，乗算回数を減らすことができる．
 y (0)  1 ω 0
 y ( 2)  2ω1
= y (1)  0 0  y (3)  0 00 0  1
0 0  01 ω 1  13 1 ω  00 ω 0 0   x ( 0) 0 1 0 ω   x(1) 0 ω 2 0   x ( 2) 2 
1 0 ω   x(3) 

これを再帰的に行うと，演算量を O (n log n) にできる（データ数 n は合成数である必要がある）
2016/5/26

計算科学技術特論B

8

DFTとFFTの演算量の比較• DFTの実数演算回数

TDFT = 8n − 2n2• FFTの実数演算回数（ n が2のべきの場合）

TFFT = 5n log 2 n
2016/5/26

計算科学技術特論B

9実数演算回数

DFTとFFTの演算量の比較
1.E+07
9.E+06
8.E+06
7.E+06
6.E+06
5.E+06
4.E+06
3.E+06
2.E+06
1.E+06
0.E+008386560回
DFTFFT

51200回
02016/5/26

123456log_2 n

計算科学技術特論B

7891010

バタフライ演算

x ( 0)+y ( 0)ωx(1)−y (1)

y (0) = x(0) + x(1)y (1) = ω{x(0) − x(1)}2016/5/26

計算科学技術特論B

11

Cooley-Tukey FFTの信号流れ図
x(0)
x(1)
x(2)
x(3)
x(4)
x(5)
x(6)
x (7 )
2016/5/26

ω0
ω2
ω
ω1
ω23
ω
0

計算科学技術特論B

ω02
ω

y (0)
ω 0 y ( 4)y ( 2)0ω y ( 6)y (1)0
ω y (5)y (3)0ωy (7 )
12

FFTのカーネル部分の例
SUBROUTINE FFT2(A,B,W,M,L)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION A(2,M,L,*),B(2,M,2,*),W(2,*)
C
DO J=1,L
WR=W(1,J)
WI=W(2,J)
DO I=1,M
B(1,I,1,J)=A(1,I,J,1)+A(1,I,J,2)
B(2,I,1,J)=A(2,I,J,1)+A(2,I,J,2)
B(1,I,2,J)=WR*(A(1,I,J,1)-A(1,I,J,2))-WI*(A(2,I,J,1)-A(2,I,J,2))
B(2,I,2,J)=WR*(A(2,I,J,1)-A(2,I,J,2))+WI*(A(1,I,J,1)-A(1,I,J,2))
END DOEND DO
RETURN
END
計算科学技術特論B
2016/5/26
13

n = n1n2 に対するFFTアルゴリズム• n = n1n2で与えられるとする．

j = j1 + j2 n1

j1 = 0, 1,  , n1 − 1

j2 = 0, 1,  , n2 − 1

k = k 2 + k1n2

k1 = 0, 1,  , n1 − 1 k 2 = 0, 1,  , n2 − 1• 上記の表現を用いると，DFTの定義式を以下のように書き換えることができる．

 j1k1
j2 k 2j1k 2
y (k 2 , k1 ) = ∑  ∑ x( j1 , j2 )ωn2 ωn1n2 ωn1
j1 = 0  j2 = 0
n1 −1 n2 −1• n 点FFTを n1 点FFTと n2 点FFTに分解している．
2016/5/26

計算科学技術特論B

14

Six-Step FFTアルゴリズム
1. 行列の転置
2. n1 組の n2 点 multicolumn FFT
j1k 2
3. ひねり係数（ω n1n2 ）の乗算
4. 行列の転置
5. n2 組の n1 点 multicolumn FFT
6. 行列の転置
2016/5/26

計算科学技術特論B

15

Six-Step FFTアルゴリズム
n1n2n 点FFTを転置

n2n 回実行n1n1n2
転置

転置n1n2
2016/5/26

計算科学技術特論B

16

Six-Step FFTのプログラム例
SUBROUTINE FFT(A,B,W,N1,N2)
COMPLEX*16 A(*),B(*),W(*)
C
CALL TRANS(A,B,N1,N2)
DO J=1,N1
CALL FFT2(B((J-1)*N2+1),N2)
END DO
DO I=1,N1*N2
B(I)=B(I)*W(I)END DOCALL TRANS(B,A,N2,N1)
DO J=1,N2
CALL FFT2(A((J-1)*N1+1),N1)
END DOCALL TRANS(A,B,N1,N2)
RETURN
END
2016/5/26

N1 x N2行列をN2 x N1行列に転置
N1組のN2点multicolumn FFT

ひねり係数（W）の乗算
N2 x N1行列をN1 x N2行列に転置
N2組のN1点multicolumn FFT
N1 x N2行列をN2 x N1行列に転置

計算科学技術特論B

17

ブロック三次元FFTアルゴリズム

2016/5/26

計算科学技術特論B

18

背景• 多くのFFTアルゴリズムはデータがキャッシュに載っている場合には高い性能を示す．
• しかし，問題サイズがキャッシュより大きくなった場合，著しく性能が低下する．
• FFTアルゴリズムにおける目標としては，– いかにしてメモリアクセスを上位のメモリ階層に封じ込めるか．
つまり，いかにキャッシュミスの回数を減らすか．
– かつ主記憶のアクセス回数をできるだけ減らす．

• 通常の三次元 FFTでは，– 3回のmulticolumn FFT– 3回の行列の転置 → ボトルネックが必要．
2016/5/26

計算科学技術特論B

19

目的• ボトルネックとなる行列の転置は，ブロック化することにより、キャッシュヒット率を上げることが可能．
• しかし，通常の三次元FFTではmulticolumn FFTと行列の転置が分離されている．
– 主記憶からキャッシュにロードされたデータは，行列の転置だけ，もしくはmulticolumn FFTにしか用いられていない．

• 行列の転置を単にブロック化するだけでは，キャッシュ内のデータの再利用性が低い．
2016/5/26

計算科学技術特論B

20

方針• 行列の転置で用いられているキャッシュ内のデータをmulticolumn FFTでも活用し，キャッシュ内のデータの再利用性を高くする．
• 主記憶からキャッシュにいったんロードしたデータは，できるだけキャッシュに載っているようにする．
• キャッシュ内のデータは再利用できるだけ再利用し，本当にいらなくなった時点で主記憶に書き戻す．
• 配列にパディングを用いて，キャッシュ内の貴重なデータが主記憶に逃げないようにする．

2016/5/26

計算科学技術特論B

21

三次元FFT• 三次元離散フーリエ変換（DFT）の定義
n1 −1 n2 −1 n3 −1

y (k1 , k 2 , k3 ) = ∑ ∑ ∑ x( j1 , j2 , j3 )j1 = 0 j2 = 0 j3 = 0ω

j3 k 3n3ω

j2 k 2n2

ω

j1k1
n10 ≤ k r ≤ nr − 1, ω nr = exp( −2πi / nr )2016/5/26

計算科学技術特論B

22

通常の三次元FFTアルゴリズム
n2 n3
n1n2
転置
n3n1

転置

行列の転置が3回必要になる．

n1n3n2転置n2 n3n1
2016/5/26

計算科学技術特論B

23

ブロック三次元FFTアルゴリズム

n2 n3nB転置

n3n1nB転置転置n2n2 n3n2 n3n1転置n1
2016/5/26

計算科学技術特論B

24

ブロック三次元FFTの利点• Cooley-Tukey FFT等の通常のFFTアルゴリズム
では，n = n1n2 n3 とすると，– 演算回数： 5n log 2 n （基数2の場合）– 主記憶アクセス回数： 4n log n2• ブロック三次元FFTでは，– 演算回数： 5n log n2– 主記憶アクセス回数：理想的には 12nn1n2 点のデータがキャッシュに入る場合： 8n

2016/5/26

計算科学技術特論B

25

In-Cache FFTアルゴリズム• Multicolumn FFTにおいて，各column FFTがキャッシュに載る場合のin-cache FFTとして，– Cooley-Tukeyアルゴリズム（ビットリバース処理が必要）– Stockhamアルゴリズム（ビットリバース処理は不要）が考えられる．
• 高い基数のFFTを積極的に用いることにより，メモリアクセス回数を減らす．
– 基数4，8のFFTを組み合わせて使用．
– 基数8よりもメモリアクセス回数の多い，基数4のFFTを最大2ステップに抑えることにより，さらにメモリアクセスを減らす．

2016/5/26

計算科学技術特論B

26

最内側ループにおけるロード＋ストア，乗算および加算回数
基数2

基数4

基数8

ロード＋ストア回数

81632

乗算回数

41232

加算回数

62266

総浮動小数点演算回数（ n log 2 n ）浮動小数点演算命令数

54.25

4.083

103498

1.25

2.125

3.063

浮動小数点演算命令数とロード＋ストア回数の比
2016/5/26

計算科学技術特論B

27

n log n

FFTカーネルを用いた場合のn点FFTの必要演算回数
54.5
43.5
32.5
21.5
10.5
0

基数2
基数4
基数8演算回数

ロード＋ストア回数

命令回数
2016/5/26

計算科学技術特論B

28

並列FFTアルゴリズム

2016/5/26

計算科学技術特論B

29

配列の分割方法（1/4）• MPIで並列化を行う際には，各ノードで配列を分割して持つようにすると，メモリを節約することができる．
• ブロック分割– 連続する領域をノード数で分割

P0 P1 P2 P3P0P1P2P3列方向に分割したブロック分割

行方向に分割したブロック分割

2016/5/26

計算科学技術特論B

30

配列の分割方法（2/4）• サイクリック分割– 1列（または1行）ごとに分割– ブロック分割に比べてロードバランスが取りやすい
0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3

列方向に分割したサイクリック分割
2016/5/26

計算科学技術特論B

31

配列の分割方法（3/4）• ブロックサイクリック分割– 複数列（または複数行）ごとに分割– ブロック分割とサイクリック分割の中間P0 P1 P2 P3 P0 P1 P2 P3

列方向に分割したブロックサイクリック分割
2016/5/26

計算科学技術特論B

32

配列の分割方法（4/4）• 二次元分割– 行方向と列方向の両方を分割– 一次元分割よりも通信量が減ることがある– 二次元のブロック分割，サイクリック分割，ブロックサイクリ
ック分割が考えられる．

P0

P1

P2

P3

二次元ブロック分割
2016/5/26

計算科学技術特論B

33

Cooley-Tukey FFTの並列化
y ( 0)y ( 4)y ( 2)
y (6)
y (1)
y (5)
y (3)
y (7 )x ( 0)
P0 x(1)x ( 2)
P1 x(3)x ( 4)
P2 x(5)x ( 6)
P3 x(7)
2016/5/26

計算科学技術特論B

34

並列Cooley-Tukey FFTの通信量• ノード数を P とすると，並列Cooley-Tukey FFT
では， log 2 P ステージの通信が必要になる．
• 各ステージでは n / P 個の倍精度複素数データの通信（MPI_Send，MPI_Recv）が行われるた
め，合計の通信量は，16nlog 2 P （バイト）TCooley−Tukey =P
2016/5/26

計算科学技術特論B

35

並列Six-Step FFTアルゴリズム
n2

全対全P0
通信を用いたn1 P0 P1 P2 P3 転置 n2n1
n2 P0 P1 P2 P3
2016/5/26

P1 P2 P3

全対全 P0 P1 P2 P3通信を用いた
転置 n

全対全通信を用いた
転置1

計算科学技術特論B

36

全対全通信（MPI_Alltoall）を用いた
行列の転置
P0 P1 P2 P3

P0 P1 P2 P3P0P1P2P3018 16 24
9 17 25 全対全010 12 34 56 72310 18 26 通信
11 19 27

8 10 12 14
9 11 13 154512 20 28
13 21 29

16 18 20 22
17 19 21 236714 22 30
15 23 31

24 26 28 30
25 27 29 31

2016/5/26

234567

計算科学技術特論B

8 9 10 11 12 13 14 15
16 17 18 19 20 21 22 23
24 25 26 27 28 29 30 31

ノード内
転置
37

並列Six-Step FFTのプログラム例
SUBROUTINE PARAFFT(A,B,W,N1,N2,NPU)COMPLEX*16 A(*),B(*),W(*)
C
CALL PTRANS(A,B,N1,N2,NPU)
DO J=1,N1/NPU
CALL FFT2(B((J-1)*N2+1),N2)
END DO
DO I=1,(N1*N2)/NPU
B(I)=B(I)*W(I)END DOCALL PTRANS(B,A,N2,N1,NPU)
DO J=1,N2/NPU
CALL FFT2(A((J-1)*N1+1),N1)
END DOCALL PTRANS(A,B,N1,N2,NPU)
RETURN
END
2016/5/26

N1 x N2行列をN2 x N1行列に転置（MPI_ALLTOALLを使用）(N1/NPU)組のN2点multicolumn FFT

ひねり係数（W）の乗算
N2 x N1行列をN1 x N2行列に転置（MPI_ALLTOALLを使用）(N2/NPU)組のN1点multicolumn FFT
N1 x N2行列をN2 x N1行列に転置（MPI_ALLTOALLを使用）

計算科学技術特論B

38並列Six-Step FFTの通信量• ノード数を P とすると，並列Six-Step FFT
では，3回の全対全通信が必要になる．
2• 全対全通信では，各ノードは n / P 個の倍精度複素数データを，自分以外のP − 1 ノードに送ることになるため，合計の通信量は16n
TSix− Step = 3 ⋅ ( P − 1) ⋅ 2 （バイト）P
2016/5/26

計算科学技術特論B

39

並列Cooley-Tukey FFTと並列
Six-Step FFTの通信量の比較• 並列Cooley-Tukey FFTの通信量16n
log 2 P
TCooley−Tukey =P• 並列Six-Step FFTの通信量16n
TSix− Step = 3 ⋅ ( P − 1) ⋅ 2P
• 両者を比較すると， P > 8 の場合には，並列
Six-Step FFTの方が通信量が少なくなる．

2016/5/26

計算科学技術特論B

40

まとめ（1/2）• ブロックFFTアルゴリズムでは，out-of-core算法の考え方を応用し，オンチップキャッシュを「主記憶」とし，メインメモリを「半導体ディスク」として扱っている．
• キャッシュ内のデータの再利用率を高くすることにより，– キャッシュミスを削減．
– その結果，主記憶のアクセス回数も削減．
• ブロックFFTアルゴリズムは，プロセッサの演算速度と主記憶のアクセス速度との差が大きい場合に，より効果的．

2016/5/26

計算科学技術特論B

41

まとめ（2/2）• 並列FFTアルゴリズムでは，問題の領域をどのように分割するかが鍵となる．
– ブロック分割，サイクリック分割，ブロックサイクリック分割• 並列FFTでは通信部分はAll-to-All通信が中心となるので，並列化は比較的容易．
• 通信量を削減するだけではなく，ブロック化等を用いたメモリアクセスの局所化も重要となる．

2016/5/26

計算科学技術特論B

42