Expression Templateを使った ベクトル演算のCUDAによる 実装と評価 みずほ情報総研 二田晴彦 GPUシンポジウム 2010 2010年10月19日

発表内容 背景と目的 CUDA概要 Expression Templateによるベクトル演 算ライブラリの実装 評価 まとめ Copyright© 2010 Mizuho Information & Research Institute, Inc.

CUDAでプログラミングする際、ホス トコードとデバイスコードが分離して いる 1-1. 背景と問題点 CUDAでプログラミングする際、ホス トコードとデバイスコードが分離して いる 「プログラミングのハードルが高い」 ベクトル計算に関して、ホストコード のみの記述でGPU上でのベクトル演算 が可能なCUBLASがあるが・・・ 「記述が直観的でない」 「計算内容によっては、メモリ転送の無 駄が生じる」 Copyright© 2010 Mizuho Information & Research Institute, Inc.

一回目の加算の結果x+yが一時的に GPUのグローバルメモリにおかれる 1-2. CUBLASの記述例 ベクトル v = x + y + z 一回目の加算の結果x+yが一時的に GPUのグローバルメモリにおかれる int main() { // CPUメモリ確保(値のセット) float* x = (float*)malloc(N * sizef(float)); … // ベクトルをGPUにコピー cublasSetVector(N, sizeof(float), x, 1, xd, 1); // GPUで計算 cublasSaxpy(N, 1.0f, yd, 1, xd, 1); // x += y cublasSaxpy(N, 1.0f, zd, 1, xd, 1); // x += z // 結果をCPUに cublasGetVector(N, sizeof(float), xd, 1, v, 1); } Copyright© 2010 Mizuho Information & Research Institute, Inc.

ホストコードだけのプログラミングで メモリ転送に無駄のない演算が可能な ベクトル演算フレームワークを実装し 速度面での評価を行う 手段 1-3. 目的 ホストコードだけのプログラミングで メモリ転送に無駄のない演算が可能な ベクトル演算フレームワークを実装し 速度面での評価を行う 手段 Expression Template(式テンプレート)を使用 する 「式の構造をテンプレートで保持し、必要になった際に、計算する手法」 Copyright© 2010 Mizuho Information & Research Institute, Inc.

ベクトル計算(ホストコード) 1-4. 実装したフレームワーク デバイスコードは書く必要なし 計算部分(1.2f * xGPU + yGPU)は、手で書いたデ バイスコードと同等のものがコンパイル時に生 成される(特殊なツール不要 nvccを使うだけ) int main() { 　// CPUメモリ確保(値のセット) float* x = (float*)malloc(N * sizef(float)); … // GPUにコピー GPU::Vector<N> xGPU = x; yGPU = y; vGPU = v; // GPUで計算(v = 1.2f * x + y) vGPU = 1.2f * xGPU + yGPU; // CPUに持ってくる vGPU.CopyToHost(v); Copyright© 2010 Mizuho Information & Research Institute, Inc.

ベクトル計算(ホストコード) この実現のためにC++の機能を多用 1-4. 実装したフレームワーク デバイスコードは書く必要なし 計算部分(1.2f * xGPU + yGPU)は、手で書いたデ バイスコードと同等のものがコンパイル時に生 成される(特殊なツール不要 nvccを使うだけ) int main() { 　// CPUメモリ確保(値のセット) float* x = (float*)malloc(N * sizef(float)); … // GPUにコピー GPU::Vector<N> xGPU = x; yGPU = y; vGPU = v; // GPUで計算(v = 1.2f * x + y) vGPU = 1.2f * xGPU + yGPU; // CPUに持ってくる vGPU.CopyToHost(v); __global__ Kernel(float* v, float* x, float* y) { unsigned int i = threadIdx.x; v[i] = 1.2f * x[i] + y[i]; } この実現のためにC++の機能を多用 Copyright© 2010 Mizuho Information & Research Institute, Inc.

クラスの実体を渡す 2-1. CUDAがサポートするC++ 二つは同じ計算をす るコード 出力されるPTXも同じ になる class Mul { public: __device__ float operator[](size_t i) { return 2.0f * b[i];} float* b; }; template <class O> __global__ void kernel( float* a, O op) { int idx = threadIdx.x; a[idx] = op[idx]; } // ホストコード Mul mul; mul.b = b; kernel<<<1, 64>>>(a, mul); __global__ void kernel( float* a ,float* b) { int idx = threadIdx.x; a[idx] = 2.0f * b[idx]; } // ホストコード kernel<<<1, 64>>>(a, b); 二つは同じ計算をす るコード 出力されるPTXも同じ になる ver. 3.2RCで確認 Copyright© 2010 Mizuho Information & Research Institute, Inc.

2-2. 二つのコードの違い 「クラスの実体」をkernelに渡すやり方 は、kernelは1つだけあらかじめ書いてお き、ホストコードの書き方次第で様々な 計算が可能になる →ホストの記述で様々なデバイスコード class Sub {…}; // a = b – 2; template <class O> __global__ void kernel( float* a, O op) { int idx = threadIdx.x; a[idx] = op[idx]; } // ホストコード Mul mul; mul.b = b; kernel<<<1, 64>>>(a, mul); Sub sub; sub.b = b; kernel<<<1, 64>>>(a, sub); Copyright© 2010 Mizuho Information & Research Institute, Inc.

ソースコードに書いたkernelは一つでも実 際には、2つのPTXのエントリーができる (それぞれ乗算、減算に最適なコード) 以下のようなコードを書いた場合・・・ ソースコードに書いたkernelは一つでも実 際には、2つのPTXのエントリーができる (それぞれ乗算、減算に最適なコード) // ホストコード Mul mul; mul.b = b; kernel<<<1, 64>>>(a, mul); Sub sub; sub.b = b; kernel<<<1, 64>>>(a, sub); .entry _Z6kernelI5MulEvPfT_ ( .param .u32 __cudaparm__Z6kernelI5MulEvPfT__a, .param .align 4 .b8 __cudaparm__Z6kernelI5MulEvPfT__op[4]) .entry _Z6kernelI5SubEvPfT_ ( .param .u32 __cudaparm__Z6kernelI5SubEvPfT__a, .param .align 4 .b8 __cudaparm__Z6kernelI5SubEvPfT__op[4]) Copyright© 2010 Mizuho Information & Research Institute, Inc.

3-1. Expression Template前段階 書いたkernelは1つでもホスト側のコード によりGPU上で様々な計算が可能 例1：ベクトル v = x + y template <class L, class R> class Plus { public: Plus(const L& a, const R& b) : a_(a), b_(b) {} __device__ float operator[](size_t i) const { return a_[i] + b_[i]; } private: const L a_; const R b_; }; Kernel<<<1, 16>>>(v, Plus<float*, float*>(x, y)); Copyright© 2010 Mizuho Information & Research Institute, Inc.

3-2. Expression Template前段階 例2：ベクトル v = x + y + z ツリー構造で式を表現 ホストコードから任意の加算が可 能になるが、記述がすごく面倒く さい Kernel<<<1, 16>>>(v, Plus<Plus<float*, float*>, float*>( Plus<float*, float*>( x, y ), z ) ); v z x y Copyright© 2010 Mizuho Information & Research Institute, Inc.

GPUグローバルメモリにあるベクトル を表すGPU::Vectorクラスを生成 3-3. Expression Template 簡単な記述を可能にするために・・・ GPUグローバルメモリにあるベクトル を表すGPU::Vectorクラスを生成 +演算子をオーバーロードして、その 中では計算を行わずに、Plusオブジェ クトを返却 式の構造を保持することに GPU::Vectorに代入が行われる際に、 kernelに式構造を渡し、式を展開 →直観的な記述が可能に Copyright© 2010 Mizuho Information & Research Institute, Inc.

→ Expression Templateにより、ホスト コードからデバイスコードが生成され る 3-4. 今回実装した演算 ベクトル 基本的な演算 加算、スカラー倍 リダクション ノルム ステンシル計算 3ポイントのステンシル → Expression Templateにより、ホスト コードからデバイスコードが生成され る Copyright© 2010 Mizuho Information & Research Institute, Inc.

評価方法 実験諸元 4-1. 性能評価 種々の演算に対して、CUBLASとの速度比 較や出力されるPTXを確認する OS: Windows XP Pro (32bit) GPU: GeForce GTX 260 CUDA Toolkit: 3.2 RC コンパイラ: Visual Studio 2005 Pro 計算精度：単精度浮動小数点数 Copyright© 2010 Mizuho Information & Research Institute, Inc.

計算式 4-2. 2項のベクトル計算 v = a * x + y CUBLASホストコード 作成したフレームワークでの(ET)ホストコード float a = 2.0f; cublasAlloc(N, sizeof(float), (void**)&xd); cublasAlloc(N, sizeof(float), (void**)&yd); cublasSetVector(N, sizeof(float), x,1, xd, 1); cublasSetVector(N, sizeof(float), y,1, yd, 1); cublasSaxpy(N, a, xd, 1, yd, 1); float a = 2.0f; GPU::Vector<N> xd = x, yd = y, vd = v; vd = a * xd + yd; Copyright© 2010 Mizuho Information & Research Institute, Inc.

速度比較(CPU <-> GPUメモリ転送除) 4-3. 2項のベクトル計算の速度 速度比較(CPU <-> GPUメモリ転送除) ETとCUBLASが同等の速度 高速 Copyright© 2010 Mizuho Information & Research Institute, Inc.

計算式 4-4. 3項のベクトル計算 v = a * x + b * y + z CUBLASホストコード 作成したフレームワークでの(ET)ホストコード float a = 2.0f, b = 3.0f; // メモリ確保と転送 … cublasSaxpy(N, a, xd, 1, zd, 1); // zd += a * xd cublasSaxpy(N, b, yd, 1, zd, 1); // zd += b * yd float a = 2.0f, b = 3.0f; GPU::Vector<N> xd = x, yd = y, zd = z, vd = v; vd = a * xd + b * yd + zd; Copyright© 2010 Mizuho Information & Research Institute, Inc.

速度比較(CPU <-> GPUメモリ転送除) 4-5. 3項のベクトル計算の速度 速度比較(CPU <-> GPUメモリ転送除) ETがCUBLASよりも高速 メモリ転送の無駄が少ないため 高速 Copyright© 2010 Mizuho Information & Research Institute, Inc.

ホストコードに書かれたvd = a * xd + b * yd + zd;から生成されたPTX ld.param.u32 %r5, [paramop+24]; add.u32 %r6, %r5, %r4; ld.global.f32 %f1, [%r6+0]; // z[idx]をグローバルから ld.param.u32 %r7, [paramop+16]; add.u32 %r8, %r7, %r4; ld.global.f32 %f2, [%r8+0]; // y[idx]をグローバルから ld.param.f32 %f3, [paramop+12]; // bのロード mul.f32 %f4, %f2, %f3; // t1 = b * y[idx]の計算 ld.param.u32 %r9, [paramop+4]; add.u32 %r10, %r9, %r4; ld.global.f32 %f5, [%r10+0]; // x[idx]をグローバルから ld.param.f32 %f6, [paramop+0]; // aのロード mad.f32 %f7, %f5, %f6, %f4; // t2 = a * x[idx] + t1 add.f32 %f8, %f1, %f7; // v = t2 + z[idx] ld.param.u32 %r11, [a]; add.u32 %r12, %r11, %r4; st.global.f32 [%r12+0], %f8; // 結果をグローバルへストア 手で書いたkernelと同等のコードが生成 Copyright© 2010 Mizuho Information & Research Institute, Inc.

計算式 4-7. ベクトルのノルム計算 v = Norm(a * x + y) CUBLASホストコード 作成したフレームワークでの(ET)ホストコード // メモリ確保と転送 … cublasSaxpy(N, a, xd, 1, yd, 1); // yd += a * xd float r = cublasSnrm2(N, yd, 1); // r = Norm(yd) GPU::Vector<N> xd = x, yd = y; CPU::Scalar r; r = CalcNorm(a * xd + yd); Copyright© 2010 Mizuho Information & Research Institute, Inc.

速度比較(CPU <-> GPUメモリ転送除) 4-8.ベクトルのノルム計算の速度 速度比較(CPU <-> GPUメモリ転送除) ETがCUBLASよりも高速 メモリ転送の無駄が少ないため 高速 Copyright© 2010 Mizuho Information & Research Institute, Inc.

ラプラス方程式 差分化した計算式(反復計算) 4-9. ラプラス方程式(1次元) 作成したフレームワークでの(ET)ホストコード // GPUのグローバルメモリの準備(uはCPU上のメモリ) GPU::Vector<DIM> uNew = u, uCur = u; GPU::Vector<DIM> uTmp; // 計算ループ(反復計算) for (int i = 0; i < 100; ++i) { uNew = GPU::StencilOf3<F(0.5), 0, F(0.5)>(uCur); // ポインタスワップ uTmp = uNew; uNew = uCur; uCur = uTmp; } Copyright© 2010 Mizuho Information & Research Institute, Inc.

uNew = GPU::StencilOf3<F(0.5), 0, F(0.5)>(uCur)から生成されたPTX // globalメモリからシェアードメモリに読み込み … ld.shared.f32 %f5, [%r6+8]; // uCur[i+1]をsharedからロード mov.f32 %f6, 0f3f000000; // 定数0.5 mul.f32 %f7, %f5, %f6; // t1 = 0.5 * uCur[i+1] ld.shared.f32 %f8, [%r6+0]; // uCur[i-1]をsharedからロード mov.f32 %f9, 0f3f000000; // 定数0.5 mad.f32 %f2, %f8, %f9, %f7; // t2 = 0.5 * uCur[i-1] + t1 ld.param.u32 %r20, [a]; add.u32 %r21, %r20, %r7; st.global.f32 [%r21+0], %f2; // 結果をglobalメモリへ 手で書いたkernelと同等のコードが生成 Copyright© 2010 Mizuho Information & Research Institute, Inc.

波動方程式 差分化した計算式(反復計算) 4-11. 波動方程式(1次元) 作成したフレームワークでの(ET)ホストコード // GPUのベクトルの準備 GPU::Vector<N> uOld = u[0], uCur = u[1], uNew = u[2]; GPU::Vector<N> uTmpGpu; // 計算ループ(反復計算) for (int i = 2; i < M; ++i) { uNew = (2.0f * uCur) + (-1.0f * uOld) + alpha * GPU::StencilOf3<F(1), F(-2), F(1)>(uCur); uTmp = uOld; uOld = uCur; // ポインタ入れ替え uCur = uNew; uNew = uTmp; } Copyright© 2010 Mizuho Information & Research Institute, Inc.

まとめ 今後の課題 5. まとめ ベクトル演算のフレームワークを実装 ホストコードの直観的な記述のみで様々 なベクトルの計算がGPU上で可能に CUDAを使ったほかのアプリケーションでもExpression Templateを使うことで、簡単な記述も可能？ 今後の課題 行列計算への対応 複数GPUの計算の対応 Copyright© 2010 Mizuho Information & Research Institute, Inc.

本発表に関する問い合わせ先 みずほ情報総研 情報・コミュニケーション部 二田 晴彦（ふただ はるひこ） 本発表以外のGPU最適化等に関する相談も受け付けております Copyright© 2010 Mizuho Information & Research Institute, Inc.