第6回 配信講義 計算科学技術特論B（2024）

大規模系での高速フーリエ変換1 高橋大介 [email protected] 筑波大学計算科学研究センター 2024/5/23 計算科学技術特論B 1

講義内容 • 高速フーリエ変換（FFT） • ブロック三次元FFTアルゴリズム • 並列FFTアルゴリズム 2024/5/23 計算科学技術特論B 2

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

離散フーリエ変換（DFT） • 離散フーリエ変換（DFT）の定義 n −1 y (k ) = ∑ x( j )ω j =0 jk n 0 ≤ k ≤ n − 1, ωn = e 2024/5/23 計算科学技術特論B − 2πi / n 4

行列によるDFTの定式化（1/4） • n = 4 のとき，DFTは以下のように計算できる． y (0) = x(0)ω + x(1)ω + x(2)ω + x(3)ω 0 0 0 y (1) = x(0)ω + x(1)ω + x(2)ω + x(3)ω 0 1 2 0 3 y (2) = x(0)ω + x(1)ω + x(2)ω + x(3)ω 0 2 4 y (3) = x(0)ω + x(1)ω + x(2)ω + x(3)ω 0 2024/5/23 3 計算科学技術特論B 6 6 9 5

行列によるDFTの定式化（2/4） • 行列を用いると，より簡単に表すことができる．  y (0)  ω ω  y (1)   0 1 ω ω  =   y (2) ω 0 ω 2   0  3  y (3)  ω ω 0 0 ω 2 ω 4 ω 6 ω 0 2 ω   x ( 0)   3  ω   x(1)  6 ω  x(2)  9  ω   x(3)  0 • n 回の複素数の乗算と，n( n − 1) 回の複素数 の加算が必要． 2024/5/23 計算科学技術特論B 6

行列によるDFTの定式化（3/4） • ω =ω jk n jk mod n n の関係を用いると，以下のよう に書き直すことができる． 1 1   x(0)   y (0)  1 1  y (1)  1 ω 1 ω 2 ω 3   x(1)    =  2 0 2  y (2) 1 ω ω ω   x(2)     3 2 1   y (3)  1 ω ω ω   x(3)  2024/5/23 計算科学技術特論B 7

行列によるDFTの定式化（4/4） • 行列の分解により，乗算回数を減らすことができる．  y (0)  1 ω 0  y (2)   2 ω 1  =  y (1)  0 0     y (3)  0 0 0 0  1  0 0  0 1 ω 1  1 3  1 ω  0 0 ω 0 0   x(0)   0  1 0 ω   x(1)  0 ω 2 0   x(2)  2  1 0 ω   x(3)  これを再帰的に行うと，演算量を O (n log n) に できる（データ数 n は合成数である必要がある） 2024/5/23 計算科学技術特論B 8

DFTとFFTの演算量の比較 • DFTの実数演算回数 2 TDFT = 8n − 2n • FFTの実数演算回数 （n が2のべきの場合） TFFT = 5n log 2 n 2024/5/23 計算科学技術特論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+00 8386560回 DFT FFT 51200回 0 2024/5/23 1 2 3 4 5 6 log_2 n 計算科学技術特論B 7 8 9 10 10

バタフライ演算 x(0) y (0) + ω x(1) y (1) y (0) = x(0) + x(1) y (1) = ω{x(0) − x(1)} 2024/5/23 計算科学技術特論B 11

Cooley-Tukey FFTの信号流れ図 x(0) x(1) x(2) x(3) x(4) x(5) x(6) x (7 ) 2024/5/23 ω ω2 0 ω0 1 ω ω2 ω3 計算科学技術特論B ω0 ω2 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 DO END DO RETURN END 計算科学技術特論B 2024/5/23 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 2 j1k 2 y (k 2 , k1 ) = ∑  ∑ x( j1 , j2 )ωn2 ωn1n2 ωn1 j1 = 0  j2 = 0  n1 −1 n2 −1 • 点FFTを n1 点FFTと n2 点FFTに分解している． 2024/5/23 計算科学技術特論B 14

Six-Step FFTアルゴリズム 1. 行列の転置 2. n1 組の n2 点 multicolumn FFT j1k 2 3. ひねり係数（ω n1n2 ）の乗算 4. 行列の転置 5. n2 組の n1 点 multicolumn FFT 6. 行列の転置 2024/5/23 計算科学技術特論B 15

Six-Step FFTアルゴリズム n1 n2 n 点FFTを 転置 n2 n 回実行 n1 n1 n2 転置 転置 n1 n2 2024/5/23 計算科学技術特論B 16

Six-Step FFTのプログラム例 C SUBROUTINE 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 DO CALL TRANS(B,A,N2,N1) DO J=1,N2 CALL FFT2(A((J-1)*N1+1),N1) END DO CALL TRANS(A,B,N1,N2) RETURN END 2024/5/23 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アルゴリズム 2024/5/23 計算科学技術特論B 18

背景 • 多くのFFTアルゴリズムはデータがキャッシュに載って いる場合には高い性能を示す． • しかし，問題サイズがキャッシュより大きくなった場合， 著しく性能が低下する． • FFTアルゴリズムにおける目標としては， – いかにしてメモリアクセスを上位のメモリ階層に封じ込めるか． つまり，いかにキャッシュミスの回数を減らすか． – かつ主記憶のアクセス回数をできるだけ減らす． • 通常の三次元 FFTでは， – 3回のmulticolumn FFT – 3回の行列の転置 → ボトルネック が必要． 2024/5/23 計算科学技術特論B 19

目的 • ボトルネックとなる行列の転置は，ブロック化する ことにより、キャッシュヒット率を上げることが可能． • しかし，通常の三次元FFTではmulticolumn FFT と行列の転置が分離されている． – 主記憶からキャッシュにロードされたデータは， 行列の転置だけ，もしくはmulticolumn FFTにしか 用いられていない． • 行列の転置を単にブロック化するだけでは， キャッシュ内のデータの再利用性が低い． 2024/5/23 計算科学技術特論B 20

方針 • 行列の転置で用いられているキャッシュ内のデータを multicolumn FFTでも活用し，キャッシュ内の データの再利用性を高くする． • 主記憶からキャッシュにいったんロードしたデータは， できるだけキャッシュに載っているようにする． • キャッシュ内のデータは再利用できるだけ再利用し， 本当にいらなくなった時点で主記憶に書き戻す． • 配列にパディングを用いて，キャッシュ内の貴重な データが主記憶に逃げないようにする． 2024/5/23 計算科学技術特論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 3 n3 ω j2 k 2 n2 ω j1k1 n1 0 ≤ k r ≤ nr − 1, ω nr = exp( −2πi / nr ) 2024/5/23 計算科学技術特論B 22

通常の三次元FFTアルゴリズム n2 n3 n1n2 転置 n3 n1 転置 行列の転置が3回 必要になる． n1n3 n2 転置 n2 n3 n1 2024/5/23 計算科学技術特論B 23

ブロック三次元FFTアルゴリズム n2 n3 nB 転置 n3 n1 nB 転置 転置 n2 n2 n3 n2 n3 n1 転置 n1 2024/5/23 計算科学技術特論B 24

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

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

最内側ループにおける ロード＋ストア，乗算および加算回数 基数2 基数4 基数8 ロード＋ストア回数 8 16 32 乗算回数 4 12 32 加算回数 6 22 66 総浮動小数点演算回数 （ n log 2 n ） 浮動小数点演算命令数 5 4.25 4.083 10 34 98 1.25 2.125 3.063 浮動小数点演算命令数と ロード＋ストア回数の比 2024/5/23 計算科学技術特論B 27

n log n FFTカーネルを用いた場合の n点FFTの必要演算回数 5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 基数2 基数4 基数8 演算回数 ロード＋ストア回数 命令回数 2024/5/23 計算科学技術特論B 28

並列FFTアルゴリズム 2024/5/23 計算科学技術特論B 29

配列の分割方法（1/4） • MPIで並列化を行う際には，各ノードで配列を分割し て持つようにすると，メモリを節約することができる． • ブロック分割 – 連続する領域をノード数で分割 P0 P1 P2 P3 P0 P1 P2 P3 列方向に分割したブロック分割 行方向に分割したブロック分割 2024/5/23 計算科学技術特論B 30

配列の分割方法（2/4） • サイクリック分割 – 1列（または1行）ごとに分割 – ブロック分割に比べてロードバランスが取りやすい 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 列方向に分割したサイクリック分割 2024/5/23 計算科学技術特論B 31

配列の分割方法（3/4） • ブロックサイクリック分割 – 複数列（または複数行）ごとに分割 – ブロック分割とサイクリック分割の中間 P0 P1 P2 P3 P0 P1 P2 P3 列方向に分割したブロックサイクリック分割 2024/5/23 計算科学技術特論B 32

配列の分割方法（4/4） • 二次元分割 – 行方向と列方向の両方を分割 – 一次元分割よりも通信量が減ることがある – 二次元のブロック分割，サイクリック分割，ブロックサイクリ ック分割が考えられる． P0 P1 P2 P3 二次元ブロック分割 2024/5/23 計算科学技術特論B 33

Cooley-Tukey FFTの並列化 P0 P1 P2 P3 y ( 0) y (4) y (2) y (6) y (1) y (5) y (3) y (7 ) x(0) x(1) x(2) x(3) x(4) x(5) x(6) x (7 ) 2024/5/23 計算科学技術特論B 34

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

並列Six-Step FFTアルゴリズム n2 n1 全対全 P0 通信を 用いた P0 P1 P2 P3 転置 n2 n1 n2 P0 P1 P2 P3 2024/5/23 P1 P2 P3 全対全 P0 P1 P2 P3 通信を 用いた 転置 n 全対全 通信を 用いた 転置 1 計算科学技術特論B 36

全対全通信（MPI_Alltoall）を用いた 行列の転置 P0 P1 P2 P3 P0 P1 P2 P3 P0 P1 P2 P3 0 1 0 1 0 1 2 3 4 5 6 7 8 16 24 9 17 25 全対全 2 3 4 5 6 7 2 3 10 18 26 通信 11 19 27 8 10 12 14 9 11 13 15 4 5 12 20 28 13 21 29 16 18 20 22 17 19 21 23 6 7 14 22 30 15 23 31 24 26 28 30 25 27 29 31 2024/5/23 計算科学技術特論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のプログラム例 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 DO CALL PTRANS(B,A,N2,N1,NPU) DO J=1,N2/NPU CALL FFT2(A((J-1)*N1+1),N1) END DO CALL PTRANS(A,B,N1,N2,NPU) RETURN END 2024/5/23 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 2024/5/23 計算科学技術特論B 39

並列Cooley-Tukey FFTと並列 Six-Step FFTの通信量の比較 • 並列Cooley-Tukey FFTの通信量 16n TCooley−Tukey = log 2 P P • 並列Six-Step FFTの通信量 16n TSix− Step = 3 ⋅ ( P − 1) ⋅ 2 P • 両者を比較すると， P > 8 の場合には，並列 Six-Step FFTの方が通信量が少なくなる． 2024/5/23 計算科学技術特論B 40

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

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