7th Symposium	

学際大規模情報基盤共同利用・共同研究拠点 公募型共同研究 平成２７年度採択課題	

jh150047-NA25	

下川辺隆史 (東京工業大学)

大規模GPU/CPU計算に向けた高生産フレームワークの構築と	これを用いた都市気流計算コードの開発	
3.4	  ステンシル計算の実行	

1

研究背景と研究目的		

従来のCPUと比べ、低消費電力で高性能であるGPUによるステンシル計算
は高い性能が得られるものの、それにはアーキテクチャに適した高度な最適化を導入する必要がある｡我々は、通常のC++を記述するだけで、GPUスパコンに必須な最適化をアプリケーションへ導入できる構造格子用GPUフレームワーク開発している。本年度は、これを拡張し、GPU/CPU両対応した高性能・高生産フレームワークを構築する｡アーキテクチャに適した最適化、自動チューニング機構を導入し、高度化する。そして、フレームワークを用い、B/F値の小さい次世代のスパコンに向いたアルゴリズムである格子ボルツマン法によって実装された都市気流計算コードを従来よりも高生産に	  情報処理学会研究報告
IPSJ SIG Technical Report
開発する。開発を通して、GPUスパコンおよびCPUスパコン上で高生産・高性能・大規模な実アプリケーションを開発する技術を確立する。	  
nzü フレームワークを用いユーザ定義関数を全格子点へ適用	  ü 自動チューニングされ実行される（現在、GPUのみ対応）	  viewat(fn,inside,info) = funcf<float>(Diffusion3d(), ptr(f), idx(f),!
ce, cw, cn, cs, ct, cb, cc);!
// f,fnはETArray, insideはRange3Dで関数適用範囲の指定!
// viewat により fn のinside領域内の格子点のみ値を更新!
// infoは自動チューニング用パラメータと結果を保持!
// ptrはETArrayのrawポインタ、idxはインデックス表現 ArrayIndex を渡す!

3.5	  GPUにおける自動チューニング	
CUDAのブロックをxy平面に配置し、各スレッドをある格子点 (i	  ,j)	  に割り当て、z方向にマーチングしながら計算する。x,	  y方向のスレッド数、z方向のマーチングする格子数が自動チューニングさ
表 1 自動チューニングされるパラメータれる。	   Table 1 Auto-tuning parameters.
自動チューニングされるパラメータ	  Number of threads in x direction

4, 8, 16, 32, 64, 128

Number of threads in y direction

1, 2, 4, 8, 16

Number of mesh in z-marching

nzmarching
threads

1, 2, 4, 8, 16

3.6	  
全要素への計算	
表 2 計算格子サイズによる最適な CUDA スレッド数と z 方向のü 全要素への計算は、ステンシル・アプリケーションでしばしば見られる	  マーチング格子数
ü 左辺のETArrayの領域を適用範囲とする（明示的範囲指定は不要）	  Table2 Optimal CUDA thread size and z-marching length forü 式からGPU/CPUコードの生成（式テンプレート）	  mesh sizes.

格子ボルツマン法による都市気流計算の例	  2

高生産フレームワークの概要	  	marchingü 陽解法による直交格子上のステンシル計算を対象	  ü C++/CUDAを用いCPUおよびNVIDIA	  GPU	  に対応	  threadsü フレームワークはC++から利用可能。言語拡張や標準でないプログラミnx

f = 0.125 * (1.0 - cos(k[0] * dx*(pos(axis0) + 0.5)))!
Meshsize - cos(k[1]Threadsz-marching + 0.5)))!
* (1.0* dy*(pos(axis1)* dz*(pos(axis2)+ 0.5)));!
323* (1.0 - cos(k[2](32, 8, 1)
2// 拡散方程式の初期化の例。左辺のf は ETArray。posにより位置情報が渡される。!
364(64, 4, 1)

2

ングモデルを利用しない。	  
à フレームワークおよびユーザコードの可搬性の向上 既存ライ
図 1 z 方向のマーチングを用いた 3 次元計算（xz 断面）
ブラリと既存コードとの連携が可能	
Fig. 1 GPU computation of 3D computational mesh with theü ユーザはステンシル関数を定義し、フレームワークで実行	  256(128, 1, 1)
3.7	  GPU間通信	

3

PBoundaryExchange<Range3D> exchange(inside, rank, connect);!
るため、チューニング・パラメータ決定のための速度低下// 境界でない計算中心領域を示す inside (Range3D型)!
// MPIの自身と隣接ランク情報 rank, connect を渡す!
はアプリケーション全体の実行時間と比較すると無視する
exchange.append(&f); // 任意の型の ETArray を登録!
ことができる。
exchange.transfer(); // 登録された ETArray の境界領域を隣接GPU間で交換!

z-marching (X-Z plane).

フレームワークの実装とプログラミングモデル	  	グラム中で Loop3D::run() 毎にそのオブジェクトを用意
Decomposition

3.1	  フレームワークの構造	

Whole domain

Decomposed subdomains

する必要がある。第三引数以降はこれまで通りファンクタü 複数GPU計算に対応	  Diffusion3d() へ与えられる。
ü フラットMPIによる並列化	  GPUGPU GPU
DeviceAutoTuningLoopEngine3D は、CUDAのブロッü ステンシル計算はGPUあるいはクを二次元に確保し、xy 平面に平行に配置し、ステンシルOpenMP並列によりCPUで実行	  計算関数を実行する。CUDA のブロック内の各スレッドを3.2	  直交格子データ構造	ある (i, j) に割り当て、z 方向にマーチングしながら順に
ü データの大きさと位置を保持した構造	  Boundary exchange
計算を行う。この計算方法の xz 断面を図1 に示す。一般ü 効率的な記述のために、このデータ構造	  ETArray	  を利用する	  に、計算領域に割り当てるスレッド数を多くすると GPU// 直交格子のデータの大きさ、位置!
では性能が向上しやすい。このため、1 スレッドで z 方向
unsigned int length[] = {nx+2*mgnx, ny+2*mgny, nz+2*mgnz};!
int
beginの全格子点を計算するのではなく、[] = {-mgnx, -mgny, -mgnz};!
z 方向でいくつかの領!

325123(128, 1, 1) 2
ü 指定された計算中心領域と	  ETArray	  データ自身の持つ大きさから通信を
8 × 512 × 512(8, 16, 1)

8必要とする境界領域を自動判定	  4

評価実験	
3.3 	
評価実験

4.1	  
計算環境：東京工業大学	  TSUBAME	  DeviceAutoTuningLoopEngine3Dを用いた
Loop3D に 2.5	ü NVIDIA	  Tesla	  K20X	  GPU:	  4224台	
よる最適パラメータ選択機能の有用性を示すため、本フü 各ノード：3	  GPUs,	  2	  Intel	  Xeon	  X5670	  CPU,	  2	  Inﬁniband	  QDR	  レームワークを用い拡散方程式を実装し、性能評価する。

4.2	  結果	
拡散方程式は流体計算等で多く用いられる方程式で，以下フレームワークで実装された拡散方程式の性能評価。左図：単一GPUで自のように表される．
動チューニング機能を用いた場合と用いない場合の実行性能。右図：複数∂fGPU計算の強スケーリング。自動チューニングにより性能向上している。	  
= κ∇2 f(1)∂t域に分割して、CUDA のブロックを割り当てる。z 方向のここで，f は物理変数，κ は拡散係数である．1 方向に 3分割数が dz で、計算領域の格子点数を (nx , ny , nz ) とする

点，3 次元計算では 7 点の格子点を参照する．

Range3D whole(length, begin); // 大きさ、位置を表す!
!

nx// ホスト上、デバイス上に直交格子用データを作成する!
とき、nz /dz の格子点をマーチングすることになる。
ETArray<float, Range3D> f_h(whole, MemoryType::HOST_MEMORY);!
DeviceAutoTuningInfo
は、現状では、自動チューニン
ETArray<float, Range3D>
f_d(whole, MemoryType::DEVICE_MEMORY);!

グするパラメータとして CUDA のブロック内の x と y 方3.3	  ステンシル計算の記述	
向のスレッド数（すなわち、1 つの CUDA ブロックが計ü ある格子点を更新する関数（C++ファンクタ）をユーザが定義	  算する x と y 方向の格子数に等しい）と z 方向にマーチü ArrayIndexによりインデックスを表現	  struct Diffusion3d ングする格子数を指定できる。自動チューニングによっ{ // ユーザ定義のステンシル関数（例は拡散方程式） !
__host__ __device__!て選択されうるスレッド数とマーチングする格子数をfloat operator()(constfloat *f, const ArrayIndex &idx, !
float ce, float表cw,float cn,float cs, float ct, !
1 に示す。
DeviceAutoTuningLoopEngine3Dが Loop3D
float cb, float cc) {!に指定されると、ステンシル計算関数が実行される度にconst float fn = + cc*f[idx.ix()]!
+ ce*f[idx.ix(1,0,0)]+ cw*f[idx.ix(-1,0,0)]!
DeviceAutoTuningInfo
によって新しい自動チューニン+ cn*f[idx.ix(0,1,0)] + cs*f[idx.ix(0,-1,0)]!
グ・パラメータが指定され、その実行時間が計測され、+ ct*f[idx.ix(0,0,1)] + cb*f[idx.ix(0,0,-1)];!
DeviceAutoTuningInfo のオブジェクトへ結果として保持
return fn; // 戻り値は、ある一点の更新する値!
}!
される。DeviceAutoTuningLoopEngine3D は、全てのパ};!

性能測定は，GPU 上での最適化の効果を見るため、単一GPU で行い、NVIDIA Tesla K20X GPU を用いる。Tesla
K20X は単精度計算で 3.95 TFlops のピーク性能を持ち，250 GB/s のメモリバンド幅に達する．
計算格子サイズを 323 から 5123 まで変化させる。また、計算領域のアスペクト比が大きい場合の例として8 × 512 × 512 の計算格子を用いた計算を合わせて行う。本5フレームワークは、カーネル分割による通信を計算で隠蔽まとめと今後の研究計画		するオーバーラップ手法を提供する。この手法では、アスフレームワークに直交格子データ構造を追加し、ステンシル計算、全要素へのペクト比が大きい境界領域とそうでない中心領域に分割さ計算へ対応した。今後は、計算による通信の隠蔽手法、効率的データ入出力、れる。
8 × 512 × 512 の計算格子は、アスペクト比が大きいリダクション計算などを導入する。都市気流計算へフレームワークを適用し、境界領域を想定している。表2 にそれぞれの格子サイズのGPU/CPUスパコン上での高性能な実アプリケーションの開発技術を確立する。	  
場合で、自動チューニングによって選択された最適なパララメータで計測が完了すると、それ以降のステンシル関数

学際大規模情報基盤共同利用・共同研究拠点
第７回シンポジウム	
メータを示す。最適な x 方向のスレッド数は、計算格子サの実行では、実行時間を最短とする最適なパラメータを用

イズの x 方向の長さに大きく依存しており、特にこのパラ2015年 7月 9日，10日 	いる。自動チューニングはプログラムの実行中に行われ、メータを実行時に最適化する必要がある。

THE	  GRAND	  HALL	  （品川）