筑波大学計算科学研究センター
CCS HPCサマーセミナー「並列数値アルゴリズムⅡ」
高橋大介
daisuke@cs.tsukuba.ac.jp
筑波大学計算科学研究センター

1

講義内容
• 高速フーリエ変換（Fast Fourier Transform，FFT）• Cooley-Tukey FFTおよび並列化• Six-Step FFTおよび並列化• Nine-Step FFTおよびブロック化，並列化

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

離散フーリエ変換（DFT）• 離散フーリエ変換（DFT）の定義
n −1y (k ) = ∑ x( j )ωj =0

jkn0 ≤ k ≤ n − 1, ωn = e− 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)ω
0

2

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

3

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) 回の複素数の加算が必要．
26

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

ω =ω
jkn

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

1
1   x ( 0)  y (0)  1 1
 y (1)  1 ω 1 ω 2 ω 3   x(1) 
=202
 y (2) 1 ω ω ω   x(2)
 321 
 y (3)  1 ω ω ω   x(3) 7

行列によるDFTの定式化（4/4）• 行列の分解により，乗算回数を減らすことができる．
0y(0)

ω
1


 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 は合成数である必要がある）
8

DFTとFFTの演算量の比較• DFTの実数演算回数
2TDFT = 8n − 2n
• FFTの実数演算回数（ n が2のべきの場合）

TFFT = 5n log 2 n9実数演算回数

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回
012

34

56
log_2 n7891010

バタフライ演算x(0)+y (0)

ωx(1)

y (1)

y (0) = x(0) + x(1)y (1) = ω{x(0) − x(1)}11Cooley-Tukey FFTの信号流れ図
x(0)
x(1)
x(2)
x(3)
x(4)
x(5)
x(6)
x (7 )

ωω20ω
1ωω23
ω
0ω02
ω

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

FFTのカーネル部分の例
C

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,*)

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
13

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

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

16n
TCooley−Tukey =log 2 P （バイト）P15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 2
j1k 2y (k 2 , k1 ) = ∑  ∑ x( j1 , j2 )ωn2 ωn1n2 ωn1j1 = 0  j2 = 0
n1 −1 n2 −1• n 点FFTを n1 点FFTと n2 点FFTに分解している．

16

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

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

n2n 回実行n1n1n2n2
転置

転置

n118

Six-Step FFTのプログラム例
CSUBROUTINE FFT(A,B,W,N1,N2)
COMPLEX*16 A(*),B(*),W(*)
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)
RETURNEND

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行列に転置

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

P0 P1 P2 P3P0P1P2P3

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

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

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

列方向に分割したサイクリック分割
21

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

列方向に分割したブロックサイクリック分割
22

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

P0P1P2P3二次元ブロック分割
23

並列Six-Step FFTアルゴリズム
n2全対全P0通信を用いた
n1 P0 P1 P2 P3 転置 n2n1n2 P0 P1 P2 P3P1 P2 P3

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

全対全通信（MPI_Alltoall）を用いた
行列の転置
P0 P1 P2 P3P0 P1 P2 P3P0P1P2P30
10
10 1

2 3

4 5

6 78 16 24
9 17 25 全対全2345672310 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 3124 26 28 30
25 27 29 31

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

ノード内転置25

並列Six-Step FFTのプログラム例
C

SUBROUTINE PARAFFT(A,B,W,N1,N2,NPU)COMPLEX*16 A(*),B(*),W(*)
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)
RETURNEND

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を使用）
26

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

並列Cooley-Tukey FFTと並列Six-Step FFTの通信量の比較• 並列Cooley-Tukey FFTの通信量

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

28

Six-Step FFTの問題点• Multicolumn FFTにおいて， n 点の各columnFFTがキャッシュに載らない場合，性能が大きく劣化する．
• 分散メモリ型並列計算機で処理する大規模FFT（例えば，2^24点以上）では，高い性能を発揮できない．

29

三次元表現n が n = n1n2 n3 と分解されるとすると
n1 −1 n2 −1 n3 −1y (k3 , k 2 , k1 ) = ∑ ∑ ∑ x( j1 , j2 , j3 )ω
j1 = 0 j2 = 0 j3 = 0

ω
1/ 3

n

j2 k3
n2 n3

ω

j2 k 2n2

j1k 2nj1k 2
n1n2

j3 k 3n3

ω ω ω

j1k1
n1点column FFTにして，キャッシュミスを低減
30

Nine-Step FFTアルゴリズム
1. 行列の転置
2. n1n2 組の n3 点 multicolumn FFT
j2 k3
3. ひねり係数（ω n2 n3 ）の乗算
4. 行列の転置
5. n1n3 組の n2 点 multicolumn FFT
j1k3
j1k 2
6. ひねり係数（ω n ω n1n2 ）の乗算
7. 行列の転置
8. n2 n3組の n1 点 multicolumn FFT
9. 転置

31

Nine-Step FFTアルゴリズム
n2 n3n1n2転置
n3n1転置

1/ 3n 点FFTを

n22/3n 回実行
n2 n1n3n1n3

転置
転置n2 n3n132

ブロックNine-Step FFTアルゴリズム

n2 n3n3n1nB
n2n2 n1n3nB転置

転置

n2 n3

転置

n2転置

n2 n3

転置 n133

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

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

基数4

基数8

ロード＋ストア回数

81632

乗算回数

41232

加算回数

62266

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

54.25

4.083

103498

1.25

2.125

3.063

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

35

n log n

FFTカーネルを用いた場合のn点FFTの必要演算回数
5
4.5
4
3.5
3
2.5
2
1.5
10.5
0

基数2
基数4
基数8

演算回数

ロード＋ストア回数

命令回数
36

Nine-Step FFTのブロック化• 行列の転置で用いられているキャッシュ内のデータをmulticolumn FFTでも活用し，キャッシュ内のデータの再利用性を高くする．
• 主記憶からキャッシュにいったんロードしたデー
タは，できるだけキャッシュに載っているようにする．
• キャッシュ内のデータは再利用できるだけ再利用し，本当にいらなくなった時点で主記憶に書き戻す．
37

並列Nine-Step FFTアルゴリズム

n2 n3n1 P0 P1 P2 P3

ノード内転置

nBn3
全対全通信

nB

ノード内転置n2n2 n3n2 P0 P1 P2 P3n2 n1n3 P0 P1 P2 P3

ノード内転置n2 n3n138

ブロックNine-Step FFTの利点• Stockham FFT等の通常のFFTアルゴリズムでは，– 演算回数： 5n log 2 n– 主記憶アクセス回数： 4n log n2• ブロックNine-Step FFTでは，– 演算回数： 5n log 2 n– 主記憶アクセス回数：理想的には 12n1/ 3• Nine-Step FFTの部分で n 点FFTにブロック化を行っていることから，提案するブロックNine-StepFFTは，「二重ブロック化」アルゴリズムといえる．
39

4並列Nine-Step FFTにおけるブロックサイズ依存性（dual-core Xeon 2.4GHz PCクラスタ, 32 cores)

NB=2
NB=4
NB=8
NB=16
NB=32
NB=6421log_2 N

2826242220018GFLOPS

340

4

並列一次元FFTの性能（dual-core Xeon 2.4GHz PCクラスタ, N = 2^23xP)

FFTE 4.0

GFLOPS

3FFTE 4.0with AT
FFTW
3.2alpha3210124816
Number of cores3241

並列一次元FFTの実行時間の内訳(dual-core Xeon 2.4GHz PCクラスタ, N=2^23xP)

14Time (sec)1210演算時間
通信時間

86420124816
Number of cores3242並列FFTライブラリの例• 商用の並列数値計算ライブラリ– Intel Cluster MKL（Math Kernel Library）• OpenMP版とMPI版が利用可能– AMD ACML（AMD Core Math Library）• OpenMP版が利用可能• オープンソースの並列FFTライブラリ– FFTW（http://www.fftw.org）• OpenMP版とMPI版が利用可能– FFTE（http://www.ffte.jp）• OpenMP版，MPI版とOpenMP+MPI版が利用可能
43

まとめ• 並列数値計算アルゴリズムとして，FFT（高速フーリエ変換）を取り上げた．
• 問題の領域をどのように分割するかが鍵となる．
– ブロック分割，サイクリック分割，ブロックサイクリッ
ク分割• 並列FFTでは通信部分はAll-to-All通信が中心となるので，並列化は比較的容易．
• 通信量を削減するだけではなく，ブロック化等を用いたメモリアクセスの局所化も重要となる．
44

課題• 離散フーリエ変換（DFT）のO(N^2)のアルゴリ
ズムと，高速フーリエ変換（FFT）のO(NlogN)のアルゴリズムを任意のプログラミング言語で実装し，n=65536の場合の実行時間を比較せよ．

45