応用数理工学特論 第9回 高速フーリエ変換とその並列化

Slides:



Advertisements
Similar presentations
Mathematica による固有値計算の高速化 Eigenvalue calculation speed by Mathematica 情報工学部 06A2055 平塚翔太.
Advertisements

1 高速フーリエ変換 (fast Fourier transform). 2 高速フーリエ変換とは? – 簡単に言うとフーリエ変換を効率よく計算 する方法 – アルゴリズムの設計技法は分割統治法に基 づいている 今回の目的は? – 多項式の積を求める問題を取り上げ、高速 フーリエ変換のアルゴリズムを用いた解法.
HBSP モデル上での 行列積を求めるアルゴリ ム 情報論理工学 吉岡健太.
大規模な三角 Toeplitz 線形方 程式の高速解法とその応用 ○ 安村 修一(法政大学 4 年) 李 磊(法政大学) 日本応用数理学会「行列・固有値の解法とその応用」研究部会 第6回研究会.
ディジタル信号処理 Digital Signal Processing
計算理工学基礎 「ハイパフォーマンスコンピューティングの基礎」
キャッシュ付PRAM上の 並列クィックソートと 並列マージソート
Intel AVX命令を用いた並列FFTの実現と評価
A Q R QR分解とは? → × ◆QR分解 QTQ = I (単位行列) ◆応用例 ◆主な計算方法 n m 今回はこの方法に注目
アルゴリズムイントロダクション第2章 主にソートに関して
プログラミング言語としてのR 情報知能学科 白井 英俊.
近似アルゴリズム 第10章 終了時刻最小化スケジューリング
研究集会 「超大規模行列の数理的諸問題とその高速解法」 2007 年 3 月 7 日 完全パイプライン化シフト QR 法による 実対称三重対角行列の 固有値並列計算 宮田 考史  山本 有作  張 紹良   名古屋大学 大学院工学研究科 計算理工学専攻.
Finger patternのブロック化による 陰的wavelet近似逆行列前処理の 高速化
ファーストイヤー・セミナーⅡ 第8回 データの入力.
スペクトル法による数値計算の原理 -一次元線形・非線形移流問題の場合-
AllReduce アルゴリズムによる QR 分解の精度について
第10回 ソート(1):単純なソートアルゴリズム
P,Q比が変更可能なScaLAPACKの コスト見積もり関数の開発
デジタル信号処理①
デジタル信号処理③
データ構造と アルゴリズム 第二回 知能情報学部 新田直也.
IT入門B2 ー 連立一次方程式 ー.
香川大学工学部 富永浩之 情報数学1 第2-2章 合同式の逆元と応用 香川大学工学部 富永浩之
担当 : 山口 匡 伊藤 祐吾 (TA) 宮内 裕輔 (TA)
岩井 儀雄 コンピュータ基礎演習  ー探索、整列ー 岩井 儀雄
システム開発実験No.7        解 説       “論理式の簡略化方法”.
首都大学東京 都市教養学部数理科学コース 関谷博之
応用数理工学特論 線形計算と ハイパフォーマンスコンピューティング
デジタル信号処理④
電気回路学Ⅱ エネルギーインテリジェンスコース 5セメ 山田 博仁.
(ラプラス変換の復習) 教科書には相当する章はない
電気回路Ⅱ 演習 特別編(数学) 三角関数 オイラーの公式 微分積分 微分方程式 付録 三角関数関連の公式
ML 演習 第 7 回 新井淳也、中村宇佑、前田俊行 2011/05/31.
応用数理工学特論 第5回 計算理工学専攻 張研究室 山本有作.
データ構造と アルゴリズム 第十一回 理工学部 情報システム工学科 新田直也.
2.伝送線路の基礎 2.1 分布定数線路 2.1.1 伝送線路と分布定数線路 集中定数回路:fが低い場合に適用
応用数理工学特論 線形計算と ハイパフォーマンスコンピューティング
7-3.高度な木 (平衡木) AVL木 平衡2分木。回転操作に基づくバランス回復機構により平衡を保つ。 B木
文献名 “Performance Tuning of a CFD Code on the Earth Simulator”
計算理工学基礎 「ハイパフォーマンスコンピューティングの基礎」
応用数理工学特論 第6回 計算理工学専攻 張研究室 山本有作.
応用数理工学特論 第9回 高速フーリエ変換の高性能化手法
スペクトル法の一部の基礎の初歩への はじめの一歩
定兼邦彦 今井浩 東京大学理学系研究科 情報科学専攻
6. ラプラス変換.
デザイン情報学科 メディア情報設計 河原英紀
デザイン情報学科 メディア情報設計 河原英紀
トーリックイデアルの グレブナ基底を求める アルゴリズム – F4およびF5 –
ディジタル信号処理 Digital Signal Processing
変換されても変換されない頑固ベクトル どうしたら頑固になれるか 頑固なベクトルは何に使える?
先進的計算基盤システムシンポジウム SACSIS2007併設企画 マルチコアプログラミングコンテスト 「Cellスピードチャレンジ2007」
プログラミング 4 整列アルゴリズム.
情報処理Ⅱ 第2回:2003年10月14日(火).
生物情報ソフトウェア特論 (2)たたみ込みとハッシュに 基づくマッチング
「データ学習アルゴリズム」 第3章 複雑な学習モデル 報告者 佐々木 稔 2003年6月25日 3.1 関数近似モデル
補講:アルゴリズムと漸近的評価.
電気回路学Ⅱ エネルギーインテリジェンスコース 5セメ 山田 博仁.
アルゴリズムとプログラミング (Algorithms and Programming)
「マイグレーションを支援する分散集合オブジェクト」
電気回路学Ⅱ 通信工学コース 5セメ 山田 博仁.
精密工学科プログラミング基礎 第7回資料 (11/27実施)
4.プッシュダウンオートマトンと 文脈自由文法の等価性
電気回路学Ⅱ 通信工学コース 5セメ 山田 博仁.
精密工学科プログラミング基礎Ⅱ 第2回資料 今回の授業で習得してほしいこと: 配列の使い方 (今回は1次元,次回は2次元をやります.)
応用数理工学特論 線形計算と ハイパフォーマンスコンピューティング
分散メモリ型並列計算機上での行列演算の並列化
Q q 情報セキュリティ 第7回:2005年5月27日(金) q q.
6.2 高速フーリエ変換 (1)FFT(fast Fourier transform)とは
2008年 7月17日 応用数理工学特論 期末発表 鈴木綾華,程飛
Presentation transcript:

応用数理工学特論 第9回 高速フーリエ変換とその並列化 応用数理工学特論 第9回 高速フーリエ変換とその並列化 2006年6月29日 計算理工学専攻 山本有作

今回の講義の目標 (1) FFTの原理と基本的な技法を学ぶ。 信号処理,偏微分方程式の求解など,広い応用範囲 メーカー提供のライブラリ,フリーウェアなどが多数存在 しかし,FFTには用途に応じて様々な変種が存在 実数データのFFT,分散メモリ向けのFFT,など 使いたいタイプのFFTが,ライブラリにあるとは限らない。 FFTの原理と基本的な技法を理解し,必要に応じて既存のソフトウェアを改造して使う力を身に付ける。 FFTに関するテクニカルタームを学ぶ。

今回の講義の目標 (2) FFTを通じて並列処理及び高性能計算の基本的技法を学ぶ。 アルゴリズムの並列化技法 並列アルゴリズムの性能分析法 通信オーバーヘッド削減のための技法 ベクトルプロセッサ,キャッシュマシン向けの高性能化技法 高性能化技法は応用範囲が広い。

もくじ FFTの基礎 単体プロセッサ向けの高速化技法 分散メモリ向けの並列化技法

今回の講義の流れ 離散フーリエ変換 FFTの原理 Cooley-Tukey FFT Stockham FFT 逆FFT FFTの分解 単体プロセッサ向け高性能化技法 1次元FFTの並列化

1. FFTの基礎 離散フーリエ変換 FFTの原理 Cooley-Tukey FFT Stockham FFT

1.1 離散フーリエ変換 (1) DFTの定義 逆DFT 1.1 離散フーリエ変換 (1) DFTの定義 N個の複素数データ a0,a1,… ,aN-1に対し,次の式で定義される c0,c1,… ,cN-1をその離散フーリエ変換(Discrete Fourier Transform)と呼ぶ。 いま,            とおくと,上式は次のように書き直せる。 定義式どおりに計算すると,DFTの計算量はO(N2)である。 逆DFT DFTの逆変換は,次の式により計算できる(確認せよ)。 これを逆DFTと呼ぶ。

1.1 離散フーリエ変換 (2) DFTの意味 FFTの応用 1.1 離散フーリエ変換 (2) DFTの意味 区間 [0, 2π]のN等分点で定義された関数 f (xn) = an を複素指数関数 exp (ikx) (k=0, 1, … , N–1)の重ね合わせに分解 分解の係数 cnを求める演算がDFT 逆に,係数 cn からN等分点での値 an を求める演算が逆DFT FFTの応用 信号処理 微分方程式の解法 畳み込み,相関の計算 多項式の積,多倍長整数の乗算 → πの計算

1.2 FFTの原理 (1) DFTの分解とFFT Nが2のべき乗のとき,DFTの式は,次のように2つの項に分解できる。   ここで,N’ = N/2, ej = a2j, oj = a2j+1。   更に,この式をkに関して前半と後半に分け,   exp(–2πi(k+N/2)/N) = –exp(–2πik/N) を用いると, 従って,N 点のDFTは2つの N/2 点のDFTと,その結果に複素指数関数 exp(–2πik) を掛けて足し合わせる処理に分解できる。 この分解を再帰的に行ってDFTを計算する方法を,高速フーリエ変換(FFT; Fast Fourier Transform)と呼ぶ。 (*)

1.2 FFTの原理 (2) FFTの計算量 FFTを用いてN点のDFTを求めるときの計算量をT(N)とすると,複素指数関数 exp(–2πik) を掛けて足し合わせる処理の計算量は実数の乗算が2N回,加算が3N回の合計5N回だから, T(1) = 0 に注意してこれを解くと(前回の講義参照), 従って,FFTを使うと 5N log2 N の計算量でN点のDFTが計算できる。 T(N) = 2T(N/2) + 5N T(N) = 5N log2 N

1.3 Cooley-Tukey FFT (1) 1次元配列への格納方式 DFTの分解の式(*)において,第1項,第2項に相当するN/2点のDFTをそれぞれ1次元配列の前半,後半に格納する方式を考える。 (*) c0 c1 以上でFFTの原理はわかった。 では,この計算式において,中間結果をどのように配列に格納するかが次の問題。 格納方式により,いろいろなFFTがある。まずはCooley-Tukey FFT そのまま加える。          を掛けて 加える。 cN–1

1.3 Cooley-Tukey FFT (2) Cooley-Tukey FFT 同じ格納方式を,N/2点のDFTのそれぞれに対しても再帰的に適用していくことにより,各ステップでの中間結果の格納場所が定まる。 計算式として(*)を用い,中間結果をこのように1次元配列に格納して計算を進める方法を Cooley-Tukey FFT と呼ぶ(Cooley & Tukey, 1965)。 また,1次元配列への格納形式と計算過程を表現するこのようなグラフをシグナル・フロー・グラフと呼ぶ。 また,各ステップでは2個の要素から2個の要素を計算する処理を繰り返し行う。この処理をバタフライ演算と呼ぶ。 ステップ0 ステップ1 ステップ2 ステップ3 a0 c0 a4 c1 a2 c2 a6 c3 a1 c4 a5 c5 a3 c6 a7 c7 N=8の場合のCooley-Tukey FFT

1.3 Cooley-Tukey FFT (3) Cooley-Tukey FFTの長所 このことを利用すると,ステップ L+1 の中間結果をステップ L の中間結果に上書きでき,配列は1個で済む。 この特徴を持つFFTを,in-place FFT と呼ぶ。 a0 c0 a4 c1 a2 c2 a6 c3 a1 c4 a5 c5 a3 c6 1回のバタフライ演算 a7 c7

1.3 Cooley-Tukey FFT (4) Cooley-Tukey FFTの短所 疑問 入力{aj}が1次元配列中で自然な順序に並ばない。 入力ajのインデックス j を2進数 jp-1…j1j0 (p = log2N),ajの1次元配列中での位置をip-1…i1i0で表すと,格納方式の定義より, したがって, ip-1= j0,ip-2 = j1,…, i0 = jp-1。 すなわち,入力{aj}はビット逆順に並ぶ。 元々の入力が自然な順序に並んでいる場合,並べ替えが必要。 疑問 入力も出力も自然な順に並ぶFFTの計算方式はないか? a0 c0 j0 = 0 なら ip-1 = 0,j0 = 1 なら ip-1 = 1 j1 = 0 なら ip-2 = 0,j1 = 1 なら ip-2 = 1 jp-1 = 0 なら i0 = 0,jp-1 = 1 なら i0 = 1 a4 c1 a2 c2 a6 c3 a1 c4 a5 c5 a3 c6 a7 c7

1.4 Stockham FFT (1) 目標 配列 XL (j, k) の定義 配列 XL (j, k)の性質 入力・出力とも自然な順序で並ぶ(self-sorting)FFTを構成する。 配列 XL (j, k) の定義 αL= 2L ,βL= 2p–L–1 とすると, XL は大きさ 2βL×αLの2次元配列 XL (j, *) は入力データを 2βL個おきに抜き出したαL個の部分列 aj,aj+2βL ,…, aj+2(αL–1 )βL のDFT 配列 XL (j, k)の性質 L = 0のとき X0 (j, 0) = aj,L = p のとき Xp (0, k) = ck すなわち, X0 (j, 0), Xp (0, k) はそれぞれ自然な順序   で並べられた入力データ,出力データと見なせる。 α2 2β2 N = 32,L = 2の ときの X2 (j, k) X0 から始めて X1,X2,…,Xp を順に計算していくことができれば,self-sortingなFFTが構成できる。

1.4 Stockham FFT (2) 配列 XL (j, k)の性質(続き) Stockham FFT Stockham FFTの特徴 DFTの分解の式(*)をXL,XL+1を使って書き直すことにより, XL (j, k) は次の漸化式を満たすことが示せる(確認せよ)。 Stockham FFT この漸化式を用いて, X0 → X1 → … → Xp-1 → Xp を順に計算していく方式を,Stockham FFT と呼ぶ。 Stockham FFTの特徴 Self-sorting である。並べ替えが不要のため,高性能計算に向く。 In-place でない。計算にはサイズ N の配列が2個必要。 XL+1 (j, k)     = XL (j, k) + XL (j+βL, k)・ωNkβL XL+1 (j, k+αL) = XL (j, k) – XL (j+βL, k)・ωNkβL ただし j = 0, 1, … , βL–1, k = 0, 1, … , αL–1

1.4 Stockham FFT (3) Stockham FFT のプログラム 計算の並列性 DO 20,DO 30のループについて完全並列 この2つのループの入れ替えも可能 共有メモリ型計算機での並列化は容易    DO 10 L = 0, p–1      α = 2L      β = 2p–L–1      DO 20 k = 0, α–1        DO 30 j = 0, β–1           XL+1 (j, k)   = XL (j, k) + XL (j+β, k)・ωNkβ           XL+1 (j, k+α) = XL (j, k) – XL (j+β, k)・ωNkβ 30      CONTINUE 20    CONTINUE 10  CONTINUE

1.4 Stockham FFT (4) 第L+1ステップの計算の図解(N =128,L = 3) 全計算ステップの図解(N = 8) XL+1 (j, k) XL+1 (j, k+αL) 2βL XL (j, k) βL 第L段の 演算 XL (j+βL, k) 2αL αL XL (j, k) XL+1 (j, k) 入力 出力

1.5 逆FFTと周波数間引き型FFT (1) 逆FFTの計算方法 (I) 簡便な計算方法 1.1(1)で述べたように,DFT,逆DFTはそれぞれ次の式で計算できる。 従って逆FFTの計算は,FFTの計算式において,ωN = exp(–2πi/N) をωN* (共役複素数)で置き換え,計算結果をNで割ればよい。 簡便な計算方法 次のようにして計算すれば,FFTのプログラムを逆FFTに転用できる。 (1) 入力データ ck を複素共役 ck* にする。 (2) ck* のFFTを計算する。 (3) 結果を複素共役にする。 (4) 1/N をかける。 DFT 逆DFT

1.5 逆FFTと周波数間引き型FFT (2) 逆FFTの計算方法 (II) この方法による逆FFTのプログラム Stockham FFT において,XL から XL+1 を求める式を逆に解き,L に関するループを逆回しにすることによっても,逆FFTは計算できる。 この方法による逆FFTのプログラム    DO 10 L = p–1, 0, –1      α = 2L      β = 2p–L–1      DO 20 k = 0, α–1        DO 30 j = 0, β–1           XL (j, k)   = (XL+1 (j, k) + XL+1 (j, k+α)) / 2           XL (j+β, k) = (XL+1 (j, k) – XL+1 (j, k+α))・ωN–kβ / 2 30      CONTINUE 20    CONTINUE 10  CONTINUE 入力 出力

1.5 逆FFTと周波数間引き型FFT (3) 周波数間引き型FFT 両FFTの使い分け 逆DFTの定義式において,ωN をωN* で置き換え,計算結果に N を掛けると,DFTの定義式となる。 従って,逆FFTの計算方法 (II) において,ωN をωN* で置き換え,計算結果に N を掛けると,再び順方向のFFTを計算するアルゴリズムが得られる。これを周波数間引き型FFTと呼ぶ。 これに対して,(*)式に基づくFFTを時間間引き型FFTと呼ぶ。 両FFTの使い分け 時間間引き型FFTと周波数間引き型FFTは,最内側の式の形が異なるため,プロセッサによっては一方が他方より性能が出やすいことがある。 対象とするプロセッサによって,性能の出やすいほうを選べばよい。

1.6 多次元FFT 2次元DFTの計算方法 多次元FFT Nx×Ny個の複素数データ {ajx, jy}に対し,次の式で定義される {ckx, ky} をその2次元DFTと呼ぶ。 これは,次のように2ステップに分けて計算できる。 ステップ1は Ny 組のデータに対する Nx 点のFFT,ステップ   2は Nx 組のデータに対する Ny 点のFFTとして計算できる。 このとき,計算量は 5 NxNy log2 (NxNy) となる。 多次元FFT 3次元以上のデータに対しても,各方向に対する1次元FFTの繰り返しにより,多次元FFTが計算できる。 x y ステップ1 ステップ2

1.7 FFTの分解 (1) 1次元DFTの2次元への分解 1次元DFTにおいて,N = lm と分解できるとする。このとき,インデックス j,k を と書き直すと,N点DFTの計算式は次のように変形できる。 ここで,内側のΣは l 組のデータに対する m 点DFT,外側のΣは m 組のデータに対する l 点DFTの形をしている。 j = rl + s   (s = 0,…,l–1, r = 0,…,m–1 ) k = pm + q (q = 0,…,m–1, p = 0,…,l–1 )

1.7 FFTの分解 (2) 1次元DFTの2次元への分解(続き) FFTの分解の応用 従って,N = lm 点のDFTの計算は次のように分解できる。 (1) l 組のデータに対する m 点FFT (2) 第 s 組の第 q 要素にひねり係数 exp(–2πiqs / N) を掛ける。 (3) m 組のデータに対する l 点FFT FFTの分解の応用 キャッシュ再利用性の向上,分散メモリ向け並列化を行う際に有効 分解を再帰的に適用することにより,3次元以上への分解も可能 aj arl+s cpm+q r s q s q s q p 4 8 12 1 2 3 1 5 9 13 4 5 6 7 2 6 10 14 8 9 10 11 l 組の m点FFT ひねり 係数 m 組の l点FFT 3 7 11 15 12 13 14 15 ck

2. 単体プロセッサ向けの高性能化技法 ベクトルプロセッサとその高性能化技法 キャッシュマシンとその高性能化技法 最内側ループ長の増加 2. 単体プロセッサ向けの高性能化技法 ベクトルプロセッサとその高性能化技法 キャッシュマシンとその高性能化技法 最内側ループ長の増加 バンクコンフリクトの削減 レジスタ再利用性の向上 キャッシュ再利用性の向上

2.1 ベクトルプロセッサとその高性能化技法 ベクトルプロセッサの特徴 (擬似)ベクトルプロセッサの例 高性能化技法 2.1 ベクトルプロセッサとその高性能化技法 ベクトルプロセッサの特徴 ベクトルレジスタ/演算器による長いベクトルの   高速演算 バンク構成による高い主メモリのスループット (擬似)ベクトルプロセッサの例 NEC SX-7(地球シミュレータ),富士通 VPP5000 日立 SR2201,SR8000,Intel Pentium 4 高性能化技法 最内側ループ長の増加 ループ交換,ループ融合 バンクコンフリクトの削減 2の大きなべき乗でのストライドアクセスの回避 演算器 ベクトルレジスタ 主メモリ 16 32 1 17 33 2 18 34 15 31 47 バンク0 1 2 15 日立 SR8000

2.2 キャッシュマシンとその高性能化技法 キャッシュマシンの特徴 キャッシュマシンの例 高性能化技法 2.2 キャッシュマシンとその高性能化技法 キャッシュマシンの特徴 レジスタ – キャッシュ – 主メモリという   階層的な記憶装置 レジスタ内のデータは高速に演算可能 主メモリ – キャッシュのスループット小 キャッシュマシンの例 Intel Pentium III,IBM Power 4 AMD Athlon,DEC Alpha,など 高性能化技法 レジスタ再利用性の向上 ループ展開 キャッシュ再利用性の向上 ブロッキング 演算器 レジスタ 8~128本程度 キャッシュ 数K~数MB スループット小 主メモリ

2.3 最内側ループ長の増加 Stockham FFT の特徴 最内側ループ長の増加 最内側のループ長が N/2,N/4,…,1 と変化 2.3 最内側ループ長の増加 Stockham FFT の特徴 最内側のループ長が N/2,N/4,…,1 と変化 最内側ループ長の増加 α>βとなった時点で DO 20 と DO 30 のループを交換することにより,最内側ループ長は常に O(N 1/2)以上に保てる。    DO 10 L = 0, p–1      α = 2L      β = 2p–L–1      DO 20 k = 0, α–1        DO 30 j = 0, β–1           XL+1 (j, k)   = XL (j, k) + XL (j+β, k)・ωNkβ           XL+1 (j, k+α) = XL (j, k) – XL (j+β, k)・ωNkβ 30      CONTINUE 20    CONTINUE 10  CONTINUE

2.4 バンクコンフリクトの削減 整合寸法の調整によるバンクコンフリクト削減 2.4 バンクコンフリクトの削減 整合寸法の調整によるバンクコンフリクト削減 k方向が連続アドレスの場合,Stockham FFTでの XL (j, k) と XL (j +βL, k) のアドレス差は αL*βL = N/2 2次元配列 XL の整合寸法 NA を調節することにより,これを NA*βL としてバンクコンフリクトを削減できる。 αL XL (j, k) βL βL XL+1 (j, k) XL+1 (j, k+αL) XL (j+βL, k) αL 余分な列 NA

2.5 レジスタ再利用性の向上 (1) 4基底FFT Stockham FFTにおいて, 2ステップ分の変換を一度に行うことで,ロード/ストアが削減でき,レジスタ再利用性が向上する。 これは,4個の要素に対して演算を行うので,4基底FFTと呼ぶ。 今まで説明してきたFFTは,2基底FFTと呼ぶ。 同様の方式により,8基底FFTも構成できる。 第Lステップ の変換 第L+1ステップ の変換 XL (j, k) XL+1 (j, k) XL+2 (j, k)

2.5 レジスタ再利用性の向上 (2) 4基底FFTの計算式 2.5 レジスタ再利用性の向上 (2) 4基底FFTの計算式 αL = 2L, βL = 2 p–L–2 とおくと,カーネルは次のようになる。 この段階ではロード/ストア削減のみ。演算量は2基底と同じ。 第Lステップの変換 XL+1 (j, k) = XL (j, k) + XL (j+2βL, k)・ω2kβL XL+1 (j, k+αL) = XL (j, k) – XL (j+2βL, k)・ω2kβL XL+1 (j +βL, k) = XL (j +βL, k) + XL (j+3βL, k)・ω2kβL XL+1 (j +βL, k+αL) = XL (j +βL, k) – XL (j+3βL, k)・ω2kβL 第L+1ステップの変換 XL+2 (j, k) = XL+1 (j, k) + XL+1 (j+βL, k)・ωkβL XL+2 (j, k+2αL) = XL+1 (j, k) – XL+1 (j+βL, k)・ωkβL XL+2 (j, k +αL) = XL+1 (j, k +αL) + XL+1 (j+βL, k +αL)・ω ( k+αL ) βL XL+2 (j, k+3αL) = XL+1 (j, k +αL) – XL+1 (j+βL, k +αL)・ω ( k +αL ) βL

2.5 レジスタ再利用性の向上 (3) 4基底FFTにおける演算量削減 2.5 レジスタ再利用性の向上 (3) 4基底FFTにおける演算量削減 第 L+1ステップでの三角関数乗算を,すべて第Lステップに持ってくる。 これにより,ωの乗算回数を削減可能 第Lステップの変換 XL+1 (j, k) = XL (j, k) + XL (j+2βL, k)・ω2kβL XL+1 (j, k+αL) = XL (j, k) – XL (j+2βL, k)・ω2kβL YL+1 (j +βL, k) = XL (j +βL, k)・ωkβL + XL (j+3βL, k)・ω3kβL YL+1 (j +βL, k+αL) = XL (j +βL, k)・ω ( k +αL ) βL – XL (j+3βL, k)・ω ( 3k +αL ) βL = – iXL (j +βL, k)・ωkβL + iXL (j+3βL, k)・ω3kβL ωαLβL = exp (–πi/2) = –i を用いて変形 第L+1ステップの変換 XL+2 (j, k) = XL+1 (j, k) + YL+1 (j+βL, k) XL+2 (j, k+2αL) = XL+1 (j, k) – YL+1 (j+βL, k) XL+2 (j, k +αL) = XL+1 (j, k +αL) + YL+1 (j+βL, k +αL) XL+2 (j, k+3αL) = XL+1 (j, k +αL) – YL+1 (j+βL, k +αL)

2.5 レジスタ再利用性の向上 (4) 4基底FFT,8基底FFTの効果 実数加算/乗算 ロード/ストア Byte/Flop値 2.5 レジスタ再利用性の向上 (4) 4基底FFT,8基底FFTの効果 2基底に比べ,レジスタへのロード/ストア回数,演算量ともに減少   実数加算/乗算 ロード/ストア  Byte/Flop値 2基底      6 / 4    4 / 4    6.4 4基底    22 / 12    8 / 8    3.76  (2基底では 24 / 16) 8基底     66 / 32   16 / 16    2.61  (2基底では 72 / 48) Byte/Flop値: 演算1回当たりに何回のロード/ストアが必要かを示す指標

2.6 キャッシュ再利用性の向上 (1) Stockham FFTのメモリアクセスパターン FFTの分解の利用 効果 2.6 キャッシュ再利用性の向上 (1) Stockham FFTのメモリアクセスパターン 各ステップで,配列 XL (j, k) の全要素(N 個)をアクセス キャッシュサイズを M とするとき,M < N ならば,毎ステップで主メモリに対し,O(N) 回のアクセスが発生 FFT全体では O(N log2 N)回 FFTの分解の利用 いま,M = N 1/2 と仮定する。 このとき,1.7 (2)で示したように, N点のFFTは次のように分解できる。 (1)  M = N 1/2 組のデータに対する M 点FFT (2) 第 s 組の第 q 要素にひねり係数 exp(–2πiqs / N) を掛ける。 (3)  M 組のデータに対する M 点FFT 効果 処理 (1),(3) において,各組のFFTはキャッシュ上で行える。 主メモリのアクセス回数は,(1),(2),(3) でそれぞれO(N) 回 FFT全体では O(N)回

2.6 キャッシュ再利用性の向上 (2) N =16,M = 4の場合の図解 一般の場合 aj arM+s cpM+q 2.6 キャッシュ再利用性の向上 (2) N =16,M = 4の場合の図解 一般の場合 N ~ M r ならば,分解を r–1回再帰的に行い,1つのFFTのサイズを M程度にすればよい。 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 r s q p M 組の M点FFT ひねり 係数 aj ck arM+s cpM+q 1組のFFT キャッシュに乗る部分

3. 分散メモリ向けの並列化技法 2次元FFTの並列化 1次元FFTの並列化

3.1 2次元FFTの並列化 (1) ブロック分割による並列化 計算量とデータ通信量 y方向にデータをブロック分割し,x方向のFFTを実行 その後,all-to-all broadcast によりデータを再分散(転置) x方向にデータをブロック分割し,y方向のFFTを実行 計算量とデータ通信量 プロセッサ数が p のとき, プロセッサあたりの計算量は 5 NxNy log2 (NxNy) / p プロセッサあたりのデータ通信量は NxNy / p,通信回数は p–1回 x y PU0 PU1 PU2 PU3 1 2 3 転置 X方向の変換を行うとき,y方向については独立であることを利用。

3.1 2次元FFTの並列化 (2) サイクリック分割による並列化 計算量とデータ通信量 x方向にFFTを行う処理は,y方向については完全独立 y方向のデータ分割は,どんな形式でもよい。 計算量とデータ通信量 ブロック分割の場合と同じ。 x y PU0 PU1 PU2 PU1 PU0 PU1 転置 PU2 PU3 PU0 1 2 3 1 2 3

3.2 1次元FFTの並列化 (1) FFTの分解の利用 N = NxNy と分解し,N点のFFTを次のように分解 aj ajxNy+jy (1)  Ny 組のデータに対する Nx 点FFT (y方向にデータ分割) (2) 第 jy 組の第 kx 要素にひねり係数 exp(–2πi kx jy / N) を掛ける。 (3)  all-to-all broadcast によりデータを再分散(転置) (4)  Nx 組のデータに対する Ny 点FFT (x方向にデータ分割) この並列FFTアルゴリズムを,転置アルゴリズムと呼ぶ。 aj ajxNy+jy ckxNx+ky jx jy kx jy kx kx 4 8 12 1 2 3 1 5 9 13 ky ky 4 5 6 7 2 6 10 14 8 9 10 11 Ny 組の Nx点FFT ひねり 係数 + 転置 Nx 組の Ny点FFT 3 7 11 15 12 13 14 15 ck 1組のFFT プロセッサ分割の境界

3.2 1次元FFTの並列化 (2) 計算量とデータ通信量 通信のオーバーヘッド プロセッサ数が p のとき, プロセッサあたりの計算量は 5 N log2 N / p プロセッサあたりのデータ通信量は N / p,通信回数は p–1回 通信のオーバーヘッド SR8000の通信/演算の性能比 1ノードの演算性能は8GFLOPS 通信速度は,複素数で0.0625Gword/s N = 230のとき,通信と演算に掛かる時間はほぼ同程度 通信オーバーヘッドを削減する方法は? → 通信の隠蔽