> [!abstract] 概要
> GPU の魅力的な FLOPS/ドル比と、現代の GPU の速度が急速に向上していることに着目し、高性能科学計算のために GPU クラスタを使うことを提案する。
> 具体例として、格子ボルツマンモデル(LBM)による並列流体シミュレーションを GPU クラスタ上に実装し、ニューヨーク市タイムズスクエア周辺の空気中汚染物質の拡散をシミュレーションした。
> 30 台の GPU ノードを使うと、480×400×80 の LBM を 0.31 秒/ステップで計算でき、CPU クラスタ実装の 4.6 倍の速度となる。
> LBM のほか、セルオートマトン、偏微分方程式ソルバー、有限要素法など、GPU クラスタの科学計算への応用可能性も論じる。
## 論文情報
- タイトル: GPU Cluster for High Performance Computing
- 著者: [[Zhe Fan]]、[[Feng Qiu]]、[[Arie Kaufman]]、[[Suzanne Yoakum-Stover]]
- 所属: [[Stony Brook University]] Center For Visual Computing and Department of Computer Science
- 媒体: ACM/IEEE Conference on Supercomputing 2004(SC'04)
- 発表日: 2004-11-06〜2004-11-12
- 既存の単一ソースノート: [[papers/2004__SC__GPU Cluster for High Performance Computing|2004__SC__GPU Cluster for High Performance Computing]]
## 概要
本論文は、グラフィックス用のコモディティ GPU を多数接続し、単一 GPU の規模制約をクラスタ化によって越える科学計算基盤を提案する。
GPU の計算性能だけでなく、GPU–CPU 転送、ネットワーク通信、データ構造の表現を含めて実装し、LBM による都市流体シミュレーションで性能を検証する。
## 問題設定
当時の GPU は高い演算性能とメモリ帯域を持つ一方、テクスチャメモリが小さく、ポインタを持つ複雑なデータ構造や複雑な制御フローを効率よく扱えなかった。
単一 GPU に収まらない大規模な科学計算を、低コストの GPU ノード群へ分割し、GPU 計算とノード間通信の両方を含む実行時間を抑えることが課題である。
対象アプリケーションには、規則格子上の局所計算と近傍通信を持つ LBM を選んだ。
## GPU 計算モデル
GPU は頂点とフラグメントをストリームとして処理する SIMD 型のプロセッサである。
グラフィックス処理では、頂点処理、ラスタライズ、フラグメント処理、テクスチャ取得が固定順序のパイプラインとして並列に進む。
本論文の一般目的計算では、データをテクスチャの色成分として配置し、ユーザー定義の Cg フラグメントプログラムで計算し、結果をピクセルバッファから次の計算用テクスチャへコピーする。
![[_attachments/GPU-Cluster-for-High-Performance-Computing/fig01-graphics-pipeline.png]]
(図1。頂点を処理して画面上のフラグメントへ変換し、フラグメント処理でテクスチャを取得して色を計算する GPU パイプライン。)
この方式には、配列や規則格子を 2D テクスチャまたはテクスチャスタックへ写像しやすいという利点がある。
一方、複雑なデータ構造には間接テクスチャが必要であり、動的連結リストのようなポインタ依存構造や複雑な制御フローは扱いにくい。
当時のテクスチャメモリは最大 256 MB 程度で、GeForce FX 5800 Ultra の 128 MB 構成では計算格子に実際に使える容量が 86 MB に制限された。
## GPU クラスタの構成
[[Stony Brook Visual Computing Cluster]] は、グラフィックス・科学計算と大規模ボリュームデータの可視化を兼ねるクラスタである。
計算に使う 32 ノードは 1 Gigabit Ethernet スイッチで接続され、各ノードは Pentium Xeon 2.4 GHz のデュアル CPU、2.5 GB メモリ、[[GeForce FX 5800 Ultra]](128 MB)を備える。
ノードは Windows XP または Linux で起動でき、実験では Windows XP を使った。
VolumePro 1000 と HP Sepia-2A コンポジティングカードも搭載されていたが、これらは本計算の中心ではない。
![[_attachments/GPU-Cluster-for-High-Performance-Computing/fig02-cluster-photo.png]]
(図2。Stony Brook Visual Computing Cluster のラック群。)
ネットワーク転送には MPI を使う。
GPU と PC メモリの間の転送は AGP 8x を経由し、GPU への下り方向はピーク 2.1 GB/s、GPU からの上り方向はピーク 133 MB/s と非対称である。
論文は、この上り方向の遅さが通信全体を制限すると指摘し、当時登場予定だった PCI Express では両方向 4 GB/s と複数 GPU の同一 PC 搭載が可能になると見込んだ。
![[_attachments/GPU-Cluster-for-High-Performance-Computing/fig03-cluster-architecture.png]]
(図3。1 Gigabit Ethernet スイッチ、各ノードのネットワークカード、CPU、PC メモリ、AGP、GPU テクスチャメモリの接続。)
GPU の理論ピークは 1 台あたり 16 Gflops、デュアル CPU は約 10 Gflops であり、32 ノード全体では `(16 + 10) × 32 = 832 Gflops` となる。
GPU だけを追加した効果は `16 × 32 = 512 Gflops` で、GPU 価格は 1 台 399 ドル、32 台で 12,768 ドルだった。
クラスタ全体の価格は約 136,000 ドルであり、GPU 追加分だけで計算ピーク 41.1 Mflops/ドルに相当する。
## 格子ボルツマン法の実装
### LBM 流体モデル
LBM は、規則格子上の流体粒子の分布を更新する計算流体力学モデルである。
D3Q19 格子では各格子点に 19 個の速度方向があり、各方向に対応する分布関数 `f_i` を保持する。
計算は、近傍へ粒子を移す streaming と、局所的に分布を平衡へ近づける collision の 2 段階で進む。
本論文は BGK 衝突モデルを使い、複雑な境界形状を格子リンクと境界面の交差位置で表現する。
![[_attachments/GPU-Cluster-for-High-Performance-Computing/fig04-d3q19-lattice.png]]
(図4。D3Q19 格子の 19 方向。速度分布 `f_i` はリンクベクトル `c_i` に対応する。)
LBM は時間と空間について 2 次精度であり、非圧縮流体の Navier–Stokes 方程式を、時間刻みと格子間隔をゼロに近づけた極限で与える。
規則格子上の局所演算が中心であるため GPU 並列化に適するが、境界情報と近傍データをテクスチャへ配置する必要がある。
### 単一 GPU への写像
19 個の速度分布をそれぞれ 1 つのボリュームとして保持し、4 ボリュームを 1 スタックの 2D テクスチャへ詰める。
19 分布は 5 スタックに収まり、密度と速度も同様にテクスチャへ配置する。
境界リンク情報は全格子ではなく、境界付近の小さな矩形に限定して保存する。
![[_attachments/GPU-Cluster-for-High-Performance-Computing/fig05-texture-layout.png]]
(図5。速度方向ごとのボリュームを 2D テクスチャのスタックへ詰めるデータ配置。)
streaming、collision、境界条件はフラグメントプログラムへ変換する。
各パスでフラグメントプログラムが必要なテクスチャを取得し、新しい格子状態を計算してピクセルバッファへ描画し、次のステップ用にテクスチャへ戻す。
この実装は GeForce FX 5900 Ultra で Pentium IV 2.53 GHz のソフトウェア実装より約 8 倍高速だった単一 GPU 実装を基礎にする。
### GPU クラスタへのスケール
LBM 格子を 3D ブロックのサブドメインへ分割し、各 GPU ノードに 1 サブドメインを割り当てる。
境界サイトの分布関数は毎ステップ隣接ノードへ流れるため、GPU からの読み出し、ネットワーク転送、隣接 GPU への書き込みが必要になる。
サブドメインを立方体に近づけると境界面積と体積の比が小さくなり、通信量を抑えられる。
![[_attachments/GPU-Cluster-for-High-Performance-Computing/fig06-subdomain-decomposition.png]]
(図6。各 GPU が 1 ブロックを担当し、軸方向の近傍と斜め方向の第二近傍へ分布を送る 2D 配置の例。)
通信の中断とスイッチングコストを抑えるため、隣接ノード間の交換を複数ステップへ分ける。
第二近傍との斜め方向のデータは直接交換せず、最近傍ノードを中継する間接パターンとする。
16 ノードの 2D 配置では、列方向と行方向の交換を 4 ステップに分けることで通信相手を整理する。
![[_attachments/GPU-Cluster-for-High-Performance-Computing/fig07-communication-schedule.png]]
(図7。16 ノードの LBM 通信スケジュール。色の異なる矢印が 4 つの交換ステップを表す。)
サブドメインが `N^3` のとき、最近傍へ送るデータ量は `5N^2`、第二近傍へ送るデータ量は `N^2` より小さい次数の項になる。
間接転送は最近傍間のパケットを大きくし、通信パターンを単純化することで、特に 3D 配置でネットワーク性能を改善する。
16 ノード未満では各ステップの `MPI barrier()` が性能を改善するが、16 ノードを超えると同期オーバーヘッドが利益を上回る。
GPU から CPU への読み出しでは、分散したテクスチャ位置の境界データをフラグメントプログラムで 1 テクスチャへ集約し、`glGetTexImage()` の読み出し回数を減らす。
## 性能評価
CPU クラスタ版と GPU クラスタ版は同じクラスタ上に実装し、開発・最適化にそれぞれ約 3 人月を費やした。
公平な比較のため CPU クラスタでは各ノード 1 スレッドだけを計算に使い、ネットワーク通信には別スレッドを使った。
各ノードは `80^3` サブドメインを計算し、GPU 版の時間は計算、GPU–CPU 通信、ネットワーク通信のうち計算と重ならない部分を含む。
### 表1: 1 ステップ実行時間
| ノード数 | CPU 合計(ms) | CPU 計算(ms) | GPU–CPU 通信(ms) | ネットワーク非重複(ms) | ネットワーク合計(ms) | GPU 合計(ms) | 高速化 |
|---:|---:|---:|---:|---:|---:|---:|---:|
| 1 | 1420 | 214 | — | — | — | 214 | 6.64 |
| 2 | 1424 | 216 | 13 | 0 | 38 | 229 | 6.22 |
| 4 | 1430 | 224 | 42 | 0 | 47 | 266 | 5.38 |
| 8 | 1429 | 222 | 50 | 0 | 68 | 272 | 5.25 |
| 12 | 1431 | 230 | 50 | 0 | 80 | 280 | 5.11 |
| 16 | 1433 | 235 | 50 | 0 | 85 | 285 | 5.03 |
| 20 | 1436 | 237 | 50 | 0 | 87 | 287 | 5.00 |
| 24 | 1437 | 238 | 50 | 0 | 90 | 288 | 4.99 |
| 28 | 1439 | 237 | 50 | 11 | 131 | 298 | 4.83 |
| 30 | 1440 | 237 | 50 | 25 | 145 | 312 | 4.62 |
| 32 | 1440 | 237 | 49 | 31 | 151 | 317 | 4.54 |
(表1。各ノードが `80^3` サブドメインを計算する場合の CPU クラスタと GPU クラスタの 1 ステップ実行時間。)
1 ノードでは GPU クラスタが 6.64 倍高速であり、通信ボトルネックが無ければ到達できる理論上限の目安になる。
28 ノードまではネットワーク通信を完全に計算へ重ね合わせられるが、それ以上では非重複部分が現れ、高速化は 30 ノードで 4.62 倍、32 ノードで 4.54 倍へ低下する。
![[_attachments/GPU-Cluster-for-High-Performance-Computing/fig08-communication-time.png]]
(図8。ノード数に対するネットワーク通信時間。青線より下が計算と重なる部分、影の部分が重ならない部分。)
![[_attachments/GPU-Cluster-for-High-Performance-Computing/fig09-speedup.png]]
(図9。GPU クラスタと CPU クラスタの高速化率。)
### 表2: スケール効率
| ノード数 | 計算セル数/秒 | 高速化 | 効率 |
|---:|---:|---:|---:|
| 1 | 2.3M | — | — |
| 2 | 4.3M | 1.87 | 93.5% |
| 4 | 7.3M | 3.17 | 79.3% |
| 8 | 14.4M | 6.26 | 78.3% |
| 12 | 20.9M | 9.09 | 75.8% |
| 16 | 27.4M | 11.91 | 74.4% |
| 20 | 34.0M | 14.78 | 73.9% |
| 24 | 40.7M | 17.70 | 73.8% |
| 28 | 45.9M | 19.96 | 71.3% |
| 30 | 47.0M | 20.43 | 68.1% |
| 32 | 49.2M | 21.39 | 66.8% |
(表2。GPU クラスタの計算セル数、1 ノード比の高速化、並列効率。)
![[_attachments/GPU-Cluster-for-High-Performance-Computing/fig10-efficiency.png]]
(図10。ノード数に対する GPU クラスタの計算効率。32 ノードでは 66.8% となる。)
32 GPU ノードでは `640×320×80 = 15.6M` セルを 0.317 秒/ステップ、49.2M セル/秒で計算した。
これは同時期の IBM SP2 や IBM Power4 の LBM 実装と比較可能な性能であり、GPU クラスタは大幅に低い価格でスーパーコンピュータと競合できると論じる。
固定問題サイズでノード数だけを増やす実験では、サブドメインが小さくなって計算/通信比が下がり、4 ノードから 16 ノードで GPU/CPU 高速化率が 5.3 倍から 2.4 倍へ低下した。
改善策として、より高速なネットワーク、PCI Express、より大きなテクスチャメモリを挙げる。
## ニューヨーク市の拡散シミュレーション
対象地域はマンハッタンの 38th Street から 59th Street、8th Avenue から Park Avenue までで、約 1.66 km×1.13 km、91 ブロック、約 850 棟の建物を含む。
都市モデルを LBM の軸へ合わせ、D3Q19 BGK LBM の `480×400×80` 格子を使った。
格子間隔は約 3.8 m、北東風を右側の速度境界条件として与え、30 GPU ノードで 0.31 秒/ステップを達成した。
1000 ステップ後に LBM の速度分布から得た遷移確率で汚染物質のトレーサ粒子を格子リンクに沿って伝播させた。
![[_attachments/GPU-Cluster-for-High-Performance-Computing/fig11-times-square-area.png]]
(図11。マンハッタン地図上で青い輪郭に囲まれたシミュレーション領域。)
![[_attachments/GPU-Cluster-for-High-Performance-Computing/fig12-velocity-field.png]]
(図12。タイムズスクエア周辺の建物上を流れる速度場。青緑の流線が水平に近い流れ、白い流線が建物越しの鉛直成分を示す。)
![[_attachments/GPU-Cluster-for-High-Performance-Computing/fig13-dispersion-density.png]]
(図13。汚染物質の分散密度をボリュームレンダリングした時刻 1000 のスナップショット。)
シミュレーション結果は当時はオフラインでレンダリングした。
GPU に結果が残っているため、将来は各ノードで直接描画し、DVI 経由のコンポジティングネットワークへ転送するオンライン可視化が可能になると展望する。
## 他の科学計算への適用可能性
LBM と同じく規則格子上の陽解法やセルオートマトンは、サブドメイン分割と境界データ交換によって GPU クラスタへ写像しやすい。
非構造格子では、固定接続を前処理でテクスチャへ配置し、間接テクスチャで隣接点の座標を保持する方法が考えられる。
CPU が複雑なデータ構造を管理し GPU が数値計算を担当する協調実行によって、GPU のデータ構造・制御フロー制約を緩和できる。
![[_attachments/GPU-Cluster-for-High-Performance-Computing/fig14-proxy-points.png]]
(図14。非局所 gather を局所 gather として扱うため、サブドメイン境界に近傍点の proxy points を追加する。)
陰解法の有限差分法や FEM は `Ax = y` の大規模疎行列系を解くため、行列とベクトルを各ノードの local points と neighbor points に分割する。
neighbor points を proxy points として局所ベクトルへ追加し、各反復でネットワークから最新値を取得して局所行列ベクトル積を実行する。
![[_attachments/GPU-Cluster-for-High-Performance-Computing/fig15-matrix-vector-decomposition.png]]
(図15。local points と proxy points から局所行列・ベクトルを構成し、行列ベクトル積の結果を local points に戻す分割。)
## 新規性
- 2004 年当時、単一 GPU 上の小規模な非グラフィックス計算に限られていた GPU 利用を、32 ノードのクラスタへ拡張した。
- GPU ノード内のテクスチャ写像、GPU–CPU 転送、MPI によるノード間通信、通信スケジュールを一つの科学計算実装として扱った。
- LBM を用いた都市規模の空気中汚染物質拡散を、30 GPU ノードで実行し、CPU クラスタとの直接比較を行った。
## 強み / 弱点・課題
### 強み
- GPU の価格性能比を、単一 GPU のデモではなくクラスタ実機で検証した点。
- 通信と計算のオーバーラップ、通信相手を減らすスケジュール、GPU–CPU 読み出しの集約を実装した点。
- LBM だけでなく、陽解法、疎行列解法、FEM、セルオートマトンへ写像原理を一般化した点。
### 弱点・課題
- 1 Gigabit Ethernet と AGP 8x の帯域、テクスチャメモリ容量、フラグメント処理という 2004 年固有の制約に強く依存する。
- CPU 比較は各ノード 1 CPU スレッドであり、SSE 最適化をまだ適用していないため、CPU 実装の比較条件には限界がある。
- 固定問題サイズではノード増加に伴い計算/通信比が低下し、GPU の計算ピークをネットワークが活かしきれない。
- 複雑なデータ構造・制御フローを GPU で直接扱うのは難しく、CPU との協調実行が必要になる。
## 関連
- 概念: [[GPUクラスタ運用]] / [[GPU最適化]] / [[HPCインターコネクトベンチマーク]] / [[格子ボルツマン法]] / [[並列化戦略]]
- エンティティ: [[Zhe Fan]] / [[Feng Qiu]] / [[Arie Kaufman]] / [[Suzanne Yoakum-Stover]] / [[Stony Brook University]] / [[Stony Brook Visual Computing Cluster]] / [[GeForce FX 5800 Ultra]]
- 関連 MOC: [[HPC - MOC]]
## 出典
- [[.raw/papers/GPU-Cluster-for-High-Performance-Computing.pdf]](本文全体、図1〜15、表1〜2)
- [[.raw/papers/GPU-Cluster-for-High-Performance-Computing.txt]](抽出テキスト)