非静力学海洋モデル kinaco の GPU による高速化 平成 28 年度高速化ワークショップ ~ 京 を中核とするHPCI メニーコアを見据えて~ 平成 29 年 3 月 24 日秋葉原 UDXカンファレンス 山岸孝輝 1, 松村義正 2 1 高度情報科学技術研究機構 2 東京大学大気海洋研究所 Ver. 1.1
発表の概要 GPU の基本 ハードの特徴実行モデル プログラミングモデル性能を引き出すための基本 海洋モデルのGPUによる高速化事例まとめ ( 本発表では NVIDIA GPU を前提 ) 2
GPU の基本 3
GPU の構造 NVIDIA Kepler GPU のブロックダイアグラム 単純な構造でリソースの大半は演算 多数のコア 4
ハード性能 ( 演算 ) 低クロックのコアを多数 バルクで高い演算性能 高並列計算コア内パイプライン化あり 複雑な機構は無し インオーダで処理無し : プリフェッチ 分岐予測 etc. 高い演算性能 高並列計算単純な命令スケジューリング 5
ハード性能 ( メモリ キャッシュ ) 高いメモリバンド幅 高レイテンシ GDDR5 の性質 キャッシュ :CPU より弱い ( 階層による ) 演算にリソースを多く割り当てた結果 高いメモリバンド幅 高レイテンシ弱いキャッシュ 6
CPU/GPU 実行モデル CPU 逐次処理 GPU 並列処理 スレッド 0 0 1 2 3 4 3 0 3 1 64 65 66 67 68 94 95 64 65 66 67 68 94 95 1 スレッドが連続データを処理 CPU コア数 = スレッド数 基本は 1 スレッドが 1 データを処理 GPU コア数 << スレッド数 各々のハードウェアの特徴によるもの 7
カーネルの CUDA 化 CPU 逐次処理 GPU 並列処理 組み込み命令でスレッド管理 void func(){ global void func(){ for(i=0;i<n;i++){ JST[i] = UTC[i] + 9; } int i = threadidx.x; JST[i] = UTC[i] + 9; 関数利用時にスレッドの総数 形状を指定 1 スレッドが連続データを処理 1 スレッドが 1 データを処理 CUDA: NVIDIA GPU 用の原語拡張 書くだけなら結構簡単だが スレッドとデータの対応付けは完全に書き手に任されている ハードウェア 実行モデル プログラミングモデルの理解のもとでスレッドとデータの対応をつける必要あり 8
CPU 依存性のある命令 時間 実行 待ち データ読み込み待ちが発生 1 スレッドが逐次で処理 do i = 1, N tmp = data( i ) out(i) = tmp + 2.0 end do プリフェッチ キャッシュ ソフトウェアパイプラインなど様々な手段で対応 9
GPU スレッドの切り替え : 実行 : 待ち ( ストール ) : 待機 (active) 0 1 2 3 ストールしたスレッド 待機していた他のスレッドに入れ替え 他のスレッドにもレジスタを割り当てて 準備をさせておく 外部メモリなどへの待避無しにすぐに切り替え可能 時間 本当は 32 スレッドごとに入れ替わる ( ここでは省略 ) 0 1 2 3 0 コアが常に稼働する状態 10
スレッドの切り替え ( 続き ) 効率的にレイテンシを隠蔽するには スレッドが多い方がよいスレッド間の同期はなるべくとりたくない スレッドを 自由に レイテンシを隠蔽出来れば GPU の利点 高い演算性能 と 広いメモリバンド幅 を活用出来る 11
スレッドの階層とメモリの階層 NVIDIA K20c スレッド レジスタ 65K ブロック ( 例 : 256 スレッド ) L1/ シェアードメモリ 48KB グリッド ( 全スレッド ) L2 1MB メモリ 5GB 各階層でどこまでデータを共有できるかを把握する レジスタと L1/ シェアードメモリは複数のスレッド ( 最大で 2048) でリソースを分け合うことに注意 12
コアレスアクセス 0 1 2 3 4 64 65 66 67 68 94 95 3 0 3 1 インデックスが連続したスレッド群 + アドレスが連続したメモリ領域 0 1 2 64 65 66 67 68 94 95 1 5 ex. ストライド2のアクセス メモリ帯域が半分無駄に!! 両方とも連続するようなアルゴリズムが望ましいまたはシェアードメモリの活用 13
スレッドレベルの並列性 独立なスレッドを複数用意スレッドを切り替えてレイテンシを隠蔽 thread 0 thread 1 thread 2 thread 3 x = x + c x = x + b x = x + a y = y + c y = y + b y = y + a z = z + c z = z + b z = z + a w = w + c w = w + b w = w + a 4 つの独立な命令 Volkov (2010) より 14
命令レベルの並列性 thread 0 w = w + b z = z + b y = y + b x = x + b w = w + a z = z + a y = y + a x = x + a 4 つの独立な命令 独立な命令をスレッドの中で並列処理 Volkov (2010) より 15
スレッドレベルの並列性 vs 命令レベルの並列性 基本はデータ並列性 多数のスレッドを用意する 複数データ /1 スレッドを試してみる どこかに最適解がある どちらが良いかは処理の内容次第余計なデータロードの削減も期待できる 16
ここまでのまとめ GPU のハード性能 Good: 高い演算性能 高バンド幅 Bad: 高レイテンシ 弱いキャッシュ GPU の実行モデルとプログラミングモデル スレッドの切り替えによるレイテンシの隠蔽階層化されたスレッドとメモリ 両者の対応 その他性能を引き出すために重要なこと コアレスアクセス命令レベルの並列性 17
海洋モデルの GPU による高速化事例 18
非静力学海洋モデル kinaco ウェッデル海における南極低層水形成の再現シミュレーション - 3 次元 Navier-Stokes 方程式 移流拡散方程式 - 静水圧近似なし 3 次元の流れを陽に計算 - 等方構造格子 - ポワソン / ヘルムホルツ方程式をマルチグリッド法を用いた前処理付き CG 法で求解 - 詳細は Matsumura and Hasumi (2008) CPU 向けの最適化実績有り ~1km のスケールの運動を詳細に表現 GPU 向けの最適化 : - Yamagishi and Matsumura (2016) - SC15, SC16 Poster session 19
3 種類のカーネル最適化事例を紹介 Case1: コアレスアクセス Case2: 命令レベルの並列性 Case3: スレッドレベルの並列性 20
Case1: 疎行列 - ベクトル積 7 点ステンシル計算に相当 CPU 実装 DO k=1, n3 DO j=1, n2 DO i=1, n1 out(i,j,k) = a(-3,i,j,k) * x(i, j, k-1) & + a(-2,i,j,k) * x(i, j-1,k ) & + a(-1,i,j,k) * x(i-1,j, k ) & + a( 0,i,j,k) * x(i, j, k ) & + a( 1,i,j,k) * x(i+1,j, k ) & + a( 2,i,j,k) * x(i, j+1,k ) & + a( 3,i,j,k) * x(i, j, k+1) END DO END DO 係数行列 END DO 係数の空間的配置 2 3-1 1-2 -3 0 a(-3,i,j,k)~a( 3,i,j,k) -3-2 -1 0 1 2 3 キャッシュライン上に 7 点をまとめる 21
GPU でコアレスアクセス a(-3,i,j,k) a(-3,i+1,j,k) a(-3,i+2,j,k) thread(id) thread(id+1) thread(id+2) GPU の各スレッドがストライドアクセス (7 間隔 ) a(i,j,k,-3) a(i+1,j,k,-3) a(i+2,j,k,-3) a(-3:3,i,j,k) a(i,j,k,-3:3) thread(id) thread(id+2) thread(id+1) 次元の入れ替えでコアレスアクセス 22
Case 2: 命令レベルの並列性 できる限り多くのスレッドを設定 (i, j, k) 3 次元スレッド (i, j, k) 1 スレッド for 1 データ i = threadidx%x + blockdim%x * (blockidx%x-1) j = threadidx%y + blockdim%y * (blockidx%y-1) k = threadidx%z + blockdim%z * (blockidx%z-1) out(i,j,k) = a(i,j,k,-3) * x(i, j, k-1) & + a(i,j,k,-2) * x(i, j-1,k ) & + a(i,j,k,-1) * x(i-1,j, k ) & + a(i,j,k, 0) * x(i, j, k ) & + a(i,j,k, 1) * x(i+1,j, k ) & + a(i,j,k, 2) * x(i, j+1,k ) & + a(i,j,k, 3) * x(i, j, k+1) スレッドを切り替えることでレイテンシを隠蔽 23
命令レベルの並列性を確保 独立な命令をスレッド内で繰り返し (i, j) 2 次元スレッド (i, j) 1 スレッド for 1 カラム i = threadidx%x + blockdim%x * (blockidx%x-1) j = threadidx%y + blockdim%y * (blockidx%y-1) DO k=1, n3 out(i,j,k) = a(i,j,k,-3) * x(i, j, k-1) & + a(i,j,k,-2) * x(i, j-1,k ) & + a(i,j,k,-1) * x(i-1,j, k ) & + a(i,j,k, 0) * x(i, j, k ) & + a(i,j,k, 1) * x(i+1,j, k ) & + a(i,j,k, 2) * x(i, j+1,k ) & + a(i,j,k, 3) * x(i, j, k+1) END DO 命令を重ねてレイテンシ隠蔽 本カーネルでは 2 次元スレッドの設定の方が高速 24
Case 3: スレッドレベルの並列性 物理過程の 1 カーネル 処理の概要 DO i,j,k loop tx(i,j,k) = Fx( A(i,j,k) ) ty(i,j,k) = Fy( B(i,j,k) ) tz(i,j,k) = Fz( C(i,j,k) ) END DO DO i,j,k loop out(i,j,k) = Fout( tx(i,j,k), tx(i-1,j,k), ty(i,j,k), ty(i,j-1,k), tz(i,j,k), tz(i,j,k-1) ) END DO i-1 j-1 k-1 (i,j,k) 実際のカーネルはより複雑 要点だけ抽出 25
CPU tuned DO k loop DO i,j loop tz0(i,j) = tz1(i,j) tx (i,j,k) = Fx( A(i,j,k) ) ty (i,j,k) = Fy( B(i,j,k) ) tz1(i,j) = Fz( C(i,j,k) ) END DO 2 次元データでブロック化キャッシュに DO i,j loop out(i,j,k) = Fout( tx(i,j,k), tx(i-1,j,k), ty(i,j,k), ty(i,j-1,k), tz1(i,j), tz0(i,j) ) END DO END DO k 軸のデータは使い回しする 26
i=threadidx%x; j=threadidx%y DO k loop r_tz0 = r_tz1 r_tx = Fx( A(i,j,k) ) r_ty = Fy( B(i,j,k) ) r_tz1 = Fz( C(i,j,k) ) r_tx_im1 = Fx( A(i-1,j,k) ) r_ty_jm1 = Fy( B(i,j-1,k) ) レジスタで確保 GPU tuned 2 次元でスレッドを生成 out(i,j,k) = Fout( r_tx, r_tx_im1, r_ty, r_ty_jm1, r_tz1, r_tz0 ) END DO 隣接格子上の計算を付加 (i,j,k) i-1 j-1 k-1 1 スレッド K 軸ループ 各スレッドが独立同期が不要 27
Case3 次に 何をやるべきか 複数データ /1 スレッドレジスタ増えるがロードストア 演算が減少シェアードメモリの利用 tx, ty, tz を各スレッドで計算させた後スレッド間で共有スレッド間で同期をとる必要あり以上はプロセッサのリソースとのバランスが重要レジスタ シェアードメモリ増加 割り当てスレッド数減少 レイテンシ隠蔽が非効率化 28
海洋モデル高速化 上記以外の施策 カーネルの分割 融合シェアードメモリによるスレッド間での共有 CPU-GPU 間転送の最小化混合精度演算 ( 連立方程式ソルバ前処理を単精度化 ) 粒子追跡コードにてメモリアクセスの改善テクスチャメモリの活用参考 :Yamagishi and Matsumura (2016), SC(15, 16)Poster Session 29
性能評価事例実験設定 CPU (Fujitsu SPARC64VIIIfx) vs GPU (NVIDIA K20c) 1 CPU vs 1 GPU 傾圧不安定を伴う対流混合実験 Visbeck et al. (1996) 外部強制 : 温度 forcing, コリオリ力 等方構造格子 size: (256, 256, 32) 時間ステップ / 総時間 2min/5hours (150 steps) 32 256 スレッド数 256 3 次元 : 200 万 2 次元 : 6 万
性能比較 経過時間 [s]: CPU vs GPU CPU GPU Speedup all components 174.2 37.3 4.7 Poisson/Helmholtz solver 36.8 10.5 3.5 others 137.4 26.8 5.1 演算, メモリ性能 : CPU vs GPU CPU dp: double precision, sp: single precision GPU Computational performance (GFLOPS) 7.7 42.3(dp)/3.8(sp) GFLOPS/PEAK (%) 6.0 3.6(dp)/0.1(sp) Memory throughput (GB/S) 22.2 114.1 GPU は CPU 比較で 4.7 倍高速 31
まとめ GPU の基本を説明 ハードの特長 実行モデル プログラミングモデル GPU で性能を引き出すための基本 非静力学海洋モデルの高速化事例紹介 コアレスアクセス 命令レベルの並列性 レジスタ活用によるスレッドレベルの並列性 紹介しきれなかった細かい工夫も多数 CPU 比較で 5 倍弱の高速化 32
References NVIDIA 公式資料 CUDA C Best Practices Guide, CUDA C Programming Guide Programming Massively Parallel Processors, Third Edition, David B. Kirk, Wen-mei W. Hwu, Morgan Kaufmann. CUDA C プロフェッショナルプログラミング, John Cheng ら, インプレス. V. Volkov. Better performance at lower occupancy. GPU Technology Conference 2010. Dan Cyca. Essential CUDA Optimization Techniques. 2014. http://on-demand.gputechconf.com/gtc/2014/video/s4702- essential-cuda-optimization-techniques-acceleware-part-4.mp4 Matsumura, Y. and Hasumi, H., 2008. A non-hydrostatic ocean model with a scalable multigrid Poisson solver. Ocean Modelling. Yamagishi, T. and Matsumura, Y., 2016. GPU Acceleration of a Nonhydrostatic Ocean Model with a Multigrid Poisson/Helmholtz Solver. Procedia Computer Science 80, 1658-1669. 出川智啓, GPGPU 実践プログラミング, http://www.toffee.jp/streaming/gpgpu/index.html 33