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

1.5K Views

May 10, 24

スライド概要

第6回 5月23日 大規模系での高速フーリエ変換1
高速フーリエ変換(FFT)の基礎的なアルゴリズムを、大規模系に適用する方法や問題点について講義する。またブロック三次元FFTアルゴリズム、並列FFTアルゴリズムについて紹介する。

シェア

またはPlayer版

埋め込む »CMSなどでJSが使えない場合

関連スライド

各ページのテキスト
1.

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

2.

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

3.

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

4.

離散フーリエ変換(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

5.

行列による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

6.

行列による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

7.

行列による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

8.

行列による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

9.

DFTとFFTの演算量の比較 • DFTの実数演算回数 2 TDFT = 8n − 2n • FFTの実数演算回数 (n が2のべきの場合) TFFT = 5n log 2 n 2024/5/23 計算科学技術特論B 9

10.

実数演算回数 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

11.

バタフライ演算 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

12.

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

13.

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

14.

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

15.

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

16.

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

17.

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

18.

ブロック三次元FFTアルゴリズム 2024/5/23 計算科学技術特論B 18

19.

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

20.

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

21.

方針 • 行列の転置で用いられているキャッシュ内のデータを multicolumn FFTでも活用し,キャッシュ内の データの再利用性を高くする. • 主記憶からキャッシュにいったんロードしたデータは, できるだけキャッシュに載っているようにする. • キャッシュ内のデータは再利用できるだけ再利用し, 本当にいらなくなった時点で主記憶に書き戻す. • 配列にパディングを用いて,キャッシュ内の貴重な データが主記憶に逃げないようにする. 2024/5/23 計算科学技術特論B 21

22.

三次元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

23.

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

24.

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

25.

ブロック三次元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

26.

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

27.

最内側ループにおける ロード+ストア,乗算および加算回数 基数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

28.

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

29.

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

30.

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

31.

配列の分割方法(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

32.

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

33.

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

34.

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

35.

並列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

36.

並列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

37.

全対全通信(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

38.

並列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

39.

並列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

40.

並列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

41.

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

42.

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