スライド 1

Size: px

Start display at page:

Download "スライド 1"

かずまさたみや
4 years ago
Views:

1 ATI Stream SDK による天文物理計算の高速化会津大学中里直人

2 計算事例 : 重力 N 体計算 No.2

3 プログラム :N 体の重力計算 (1) No.3 既存のアルゴリズムやアプリケーションを CAL で実装するには前提として並列計算可能な問題でなくては利用する意味がない GPU のアーキテクチャにあわせてアルゴリズムを変更する必要あり GPU のメモリに合わせた効率のよいデータ構造を考える必要あり以下例題として重力多体問題の計算を実装した場合を説明

4 プログラム :N 体の重力計算 (2) No.4 計算すべきモノ : 重力加速度とポテンシャル a ( m r i j ( r r j i 2 r j ) 2 ) 3/ 2 p ( r i m r j j 2 2 ) 1/ 2 入力 : 位置ベクトルと質量ベクトル 3 成分 + スカラー 1 出力 : 加速度ベクトルとポテンシャルベクトル 3 成分 + スカラー 1

5 プログラム :N 体の重力計算 (3) No.5 データ構造 (array of structure) 入力 : (x,y,z,w) の 4 成分を 1 粒子にわりあて単精度浮動小数点変数を利用ホスト上 x[n], y[n], z[n], m[n] GPU メモリ上 x[0] y[0] z[0] m[0] x[1] y[1] z[1] m[1] x[2] y[2] z[2] m[2] x[0] y[0] z[0] m[0] x[1] y[1] z[1] m[1] x[2] y[2] z[2] m[2] x[0] y[0] z[0] m[0] x[1] y[1] z[1] m[1] x[2] y[2] z[2] m[2] x[0] y[0] z[0] m[0] x[1] y[1] z[1] m[1] x[2] y[2] z[2] m[2] 出力 : (x,y,z,w) の 4 成分を 1 粒子にわりあてホスト上 ax[n], ay[n], az[n], p[n] GPU メモリ上 ( 入力と同様 ) ax[0] ay[0] az[0] p[0] ax[1] ay[1] z[1] p[1] ax[2] ay[2] az[2] p[2]

6 プログラム :N 体の重力計算 (4) No.6 ループ計算の詳細 1. r i (x,y,z,m) を読み込む (1 回のみ ) 2. r j を読み込む 3. r i とr j の間の相互作用を計算する a 4. 結果を積算する以上の 2 から 4 の繰り返し ( m r i j ( r r j i 2 r j ) 2 ) 3/ 2

7 プログラム :N 体の重力計算 (5) No.7 ループによる計算部分 whileloop ige r88.x, r100.x, r77.x break_logicalnz r88.x sample_resource(0)_sampler(0) r0, r2 sub r5.xyz, r0.xyz, r4.xyz dp4 r6, r5, r5 rsq r7, r6 mul r8, r7, r7.xyz1 mul r8, r8, r7.xyz1 mul r9, r8, r5.xyz1 mad r3, r9, r0.w, r3 add r2.x, r2.x, l1.x iadd r100.x, r100.x, l0.z umod r101.x, r100.x, r77.y if_logicalz r101.x add r2.0y, r2.0y, l1.x endif endloop loop 終了判定相互作用計算 index 更新

8 プログラム :N 体の重力計算 (6) No.8 相互作用計算の注意点 2 2 r i r j dp4 によるの計算 rsq による x -0.5 の計算 r2.xy によるデータ読み込み GPU メモリを 2 次元メモリとして定義しているので r2.xy を正しく update する必要がある x,y の update pattern (domain = {10,10}) 0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 1,1 1,2 1,3 1,4 1,5 1,6 1,7 1,8 1,9 2,0 2,1 2,2 2,3 2,4 2,5 2,6 2,7 2,8 2,9 3,0 3,1 3,2 3,3 3,4 3,5 3,6 3,7 3,8 3,9 GPU 上の粒子メモリの構造

9 プログラム :N 体の重力計算 (7) No.9 IL プログラム全体 il_ps_2_0 dcl_input_interp(linear) v0.xy dcl_output_generic o0 dcl_cb cb0[1] dcl_resource_id(0)_type(2d,unnorm)_fmtx(float)_fmty(float)_fmtz(float)_fmtw(float) dcl_literal l0, 1.0, 0x0, 0x1, 1.0 mov r100.0, r100 mov r2.00, r2 mov r3.0000, r3 ftoi r77.xy, cb0[0].zw mov r5. w, cb0[0].y sample_resource(0)_sampler(0) r4, v0.xy whileloop ige r88.x, r100.x, r77.x break_logicalnz r88.x sample_resource(0)_sampler(0) r0, r2 sub r5.xyz, r0.xyz, r4.xyz dp4 r6, r5, r5 rsq r7, r6 mul r8, r7, r7.xyz1 mul r8, r8, r7.xyz1 mul r9, r8, r5.xyz1 mad r3, r9, r0.w, r3 add r2.x, r2.x, l0.x iadd r100.x, r100.x, l0.z umod r101.x, r100.x, r77.y if_logicalz r101.x add r2.0y, r2.0y, l0.x endif endloop mov o0, r3 ret_dyn end

10 プログラム :N 体の重力計算 (8) No.10 ホストプログラムの流れ 1. 初期化データ読み込み GPUメモリの設定 2. 座標データの送信 3. GPU kernelの実行 4. kernel 終了待ち 5. 結果を回収 6. 積分して新しい座標を得る以上の 2 から 6 を繰り返し実行

11 プログラム :N 体の重力計算 (9) No.11 性能の測定 (Radeon 4870) ボードのピーク性能 : ~ 1.2 TFLOPS 160 (VLIW units) x 5 (FP units) x 2 (FPMAD) x 750 (MHz) ~ 1200 GFLOPS 積和演算 (FPMAD) A = αb + C (2 浮動小数点演算 ) 1 interaction 38 演算で count した実測値 ~ 200 GFLOPS 全命令で積和演算の活用は不可能 4 要素の演算になっていない部分がある

12 プログラム :N 体の重力計算 (10) No.12 VLIW 命令の解析 CAL APIによりVLIW 命令を取得可能 VLIWのスロットが多く埋まるほど性能が高い slot 命令数割合 (%) 計 19 おおよそ 50% は効率が 2 割以下

13 N 体計算の最適化 (1) No.13 高速化の余地今の実装は array of structure を利用個々のkernelが1 粒子を担当 1 粒子読み込み 1 粒子の計算 one iterationで1 相互作用を計算計算数 / メモリの比率が小さい 4 要素のベクトル SIMD 命令が効率的でない Structure of Array にすると個々のkernelが4 粒子を担当 4 粒子読み込みして 4 x 4 = 16 相互作用計算数 / メモリの比率が大幅増加ベクトル SIMD 命令が少し効率的になる

14 N 体計算の最適化 (2) No.14 1 度に 1 相互作用 (AoS 構造の場合 ) 1 度に 4 相互作用 (AoS でアンローリング )

15 N 体計算の最適化 (3) No.15 1 度に 16 相互作用 (SoA でアンローリング ) Structure of Arrayに変更個々のTPが4 粒子を担当イテレーションごとに4 粒子読み込みメモリ読み出しあたりの演算量が増加する

16 N 体計算の最適化 (4) No.16 最適化結果 slot 命令数割合 (%) 計 19 slot 命令数割合 (%) 計 81 VLIW 命令の実行効率が大幅に増加演算性能は約 4 倍 (~ 800 Gflops)

17 N 体計算の最適化 (5) No.17 ループ制御構造の最適化 slot 命令数割合 (%) 計 81 slot 命令数割合 (%) 計 66 命令数の削減と効率の更なる向上さらに 20% ほど性能が向上

18 N 体計算の最適化 (6) No.18 結果のまとめ最新 GPU における性能 ~ 2.6 Tflops

19 計算事例 : 高精度数値積分 No.19

20 高精度演算の必要性 No.20 倍精度では十分ではない問題条件数が非常に大きい (>10 16 ) 行列メッシュを再帰的に分割するAMR 分割数が 50 以上となると倍精度では不足ファインマンループの数値積分二重指数関数型積分公式 ε 算法精度の足りない例 : ~1.1726@ 倍精度 a b f b 6 a 2 (11a 2 b 2 b 6 121b 4 2) 5.5b 8 a 2b

21 高精度演算の実現 No.21 整数演算でエミュレーション 2 進数による浮動小数点演算をそのまま実装仮数部を整数配列で表現単純な実装では分岐命令が多くなる FP 演算でエミュレーション Knuth (1969), Dekkar (1971) 数値をFP 変数の和と考える A ai :( a1 a a b a 1.0, b a n )

22 FP 演算でエミュレーション No.22 四倍精度 (DD) 演算の場合変数 2つの倍精度変数で表現精度仮数部 106 bit, 指数部 11 bit 加算 20 回の倍精度演算演算密度 5.0 演算 /1 語読み出し乗算 23 回の倍精度演算演算密度 5.7 演算 /1 語読み出し演算密度が高いため GPU での計算にむいているキャッシュありの現代の CPU にも向いている?

23 DD 演算の CPU での性能 No.23 CPU での演算性能まとめ加算の場合で Mflops 乗算の場合で Mflops 演算器のレイテンシがボトルネック x86 アーキテクチャでは論理レジスタが少ないためループアンローリングは効かない

24 GPU での DD 計算 (1) No.24 IL により実装 IL は VLIW の機械語に翻訳される以下 VLIW 命令数での結果加算 21 step 乗算 25 step 除算 53 step 性能予測 750 MHz 時秒間 1.2x10 11 個の VLIW 命令加算, 乗算, 除算 : 5.7, 5.2, 2.3 Gflops

25 GPU での DD 計算 (2) No.25 単独演算での VLIW スロットの分布命令 5 slots 4 slots 3 slots 2 slots 1 slots 計加算乗算除算演算器の利用率が低いため演算性能が低めになっている演算が連続するとスロットがより埋まるため演算性能が向上すると予測される

26 高精度演算の応用例 No.26 ファインマンループの積分素粒子衝突実験の検証に必要とされる情報落ちが発生するため倍精度では困難多重積分を100 万組のパラメータについて計算一例では 5.5x10 16 FP operations

27 実性能の評価 (1) No.27 我々が開発中のコンパイラによる結果 ILコードは 555 行 VLIW 命令は481 step 81% は 4または5 slotsが埋まっている命令融合の効果を確認利用レジスタ数は 39 個性能向上の余地がある slot 命令数割合 (%) 計 481

28 実性能の評価 (2) No.28 CPU, GPU, GRAPE-DR において級数の項数を変化させて実機で計算した CPU ~ 80 Mflops GPU ~ Gflops ( 経過時間 sec ) I/Oが高速のためN 依存が小さい倍精度演算器の利用率 ~ Gflops 相当

29 計算事例 : 複雑な N 体計算 No.29

30 実用的な N 体計算 No.30 O(N 2 ) アルゴルリズムは演算精度がフルに必要な場合のみ利用必ずしも高精度は必要ない N ~ 100 k くらいまでに実用的近似的だが高速な計算方法 O(N log N) アルゴリズム Particle-Mesh method (FFT based) Oct-tree method O(N) アルゴリズム Fast-Multipole Method

31 Oct-tree Method No.31 ツリーデータ構造を利用して遠方の粒子集団を多重極モーメントに置き換える n 個を置き換えると演算量を 1/n にできる O(N log N) の演算量 << O(N 2 )

32 ツリーデータ構造 No.32

33 No.33 Program Flow of Oct-tree Method 1. ツリーデータ構造の構築ツリーはノード ( 複数粒子 ) と粒子からなる 2. ツリーを根から走査して個々のノードの多重極を計算 3. それぞれの粒子に対して 1. ツリーを根から走査する 1. もし対象が粒子なら力の計算 2. もし対象がノードなら置換条件をテストして力の計算またはノードの子ノードを走査する

34 Oct-tree on GPU No.34 第 3 ステップを GPU に実装した個々の粒子に対するツリー走査は並列計算可能 GPU での効率よく動作する可能性があるツリー法のベクトル化と同等の手法単純な N 体手法との大きな違いはメモリアクセスが一様からランダムになること重力計算以外にも応用あり Coulomb force short-range MD force Hydrodynamics (SPH) Any algorithm required neighbor particles

35 Flow of Oct-tree on GPU No ツリーデータ構築 2. 多重極計算 3. ツリーデータを GPU メモリに書き込み 4. 個々の粒子について (on GPU) 1. ツリーを走査して展開条件をテスト 2. テスト結果に基づき力の計算または更なるツリー走査 5. 結果を GPU メモリから回収する

36 Oct-tree on GPU: 最適化 No.36 データ配置の最適化粒子の分布を最適化することで性能が向上キャッシュヒット率

37 キャッシュヒット率の高い粒子配置 No.37 Morton-order

38 Performance No.38

Microsoft PowerPoint - CAL1

Microsoft PowerPoint - CAL1 AMD GPU プログラミング N.Nakasato GPU のプログラミング (1)! 2000 年頃以降 " programable shader が搭載されるようになる HLSL, GLSL, Cg 等のシェーディング言語で shader をプログラム可能になった " ただし制限が多かった " General Purpose GPU programming GPU での汎用計算 (GPGPU)!