演習II (連続系アルゴリズム) 第2回: GPGPU 須田研究室 M1 本谷 徹 motoya@is.s.u-tokyo.ac.jp 2012/10/19
GPU 画像 動画処理用ハードウェア 低性能なプロセッサがたくさん詰まっている ピーク性能が非常に高い GPUを数値計算に用いるのがGPGPU Graphics Processing Unit General Purpose GPU TSUBAME2.0: GPUスパコン 本演習ではNVIDIA社のGPUを用いる CUDA プログラミングをしてもらいます
性能指標 (1) FLOPS Floating point operations per second 1秒あたりに何回浮動小数点数演算ができるか TOP500 (http://www.top500.org/ ) にスパコンのLINPACK 性能, ピーク性能が載っている ピーク性能 (理論性能): 物理的な性能の限界 性能を測る際の目標 Xeon X5670 - 数十GFLOPS Tesla M2050 - 倍精度: 515 Gflops 単精度: 1.03 Tflops GPU はピーク性能に近づけるのが難しい
性能指標 (2) バンド幅 通信帯域 1秒あたりにどれだけデータが送れるか Tesla M2050 の理論値は 148 GB/sec GPU内のプロセッサ メモリ間のバンド幅
GPU アーキテクチャ SP: SM: Cuda core Multiprocessor
GPU の情報を調べるには? devicequery コマンドを使う csc ではログインノードで実行してもGPUが無い ジョブに書いて投げてください 計算ノードに入ればGPUが見えます GPUの情報がたくさん返ってくる コア数 メモリ量 メインメモリ キャッシュ 共有メモリ Warp cudagetdevcieproperties を呼ぶ
CUDAプログラミング Compute Unified Device Architecture プログラミングモデルはSPMD CPUからGPUを呼び出す GPUでの命令処理はSIMD CPUをホスト GPUをデバイスと呼んで区別 Single Instruction Multiple Data C, C++ の拡張として実装されている 拡張子は *.cu コンパイラはnvcc #include<cuda.h>, #include<cuda_runtime.h> csc では Torque を通して実行してください ノード は 1つ
CUDAのスレッド管理 実際のSP数よりはるかに多いスレッドを走らせることができる スレッドは階層的になっている グリッドにスレッドを詰め込んで その中にあるスレッドをまとめて実行 グリッドの中にはブロックが2次元的に配置される 1辺の長さは最大で65535 各辺の最大は65535 x 65535 (x1) ブロック処理はSMに割り当てられる ブロックの中にはスレッドが3次元的に配置される 1つのブロックの最大スレッド数は512 各辺の最大は512 x 512 x 64 スレッドは平行に処理される スレッド処理はSPに割り当てられる
ベクトル型 dim3 宣言 3次元ベクトル整数型 x,y,z の要素を持つ構造体のようなもの 要素は全て1で初期化されている dim3 hoge(8, 9, 10); dim3 hoge; hoge.x = 8; hoge.y = 9; hoge.z = 10; dim3 hoge = {8, 9, 10}; dim3 に似たものに uint3 (符号無整数ベクトル型)がある uint3 は構造体と同じ扱いをする unit3 と同様のものとしてshort3, float3, float4 などが用意されている float4 のメンバは x,y,z,w
カーネル関数 (1) CPUから直接呼び出され GPUで計算する関数をカーネル関数と呼ぶ 普通の関数呼び出しに <<<>>> をかませるだけ cf. function<<<dg, Db>>>(a, b, c); Dg, Db はdim3 変数 Dg にはグリッドに詰めるブロックの数を渡す Db には1つのブロックに詰めるスレッドの数を渡す 個々のスレッドの処理をカーネル関数に記述する uint3 blockidx, uint3 threadidx, dim3 griddim, dim3 blockdim の４つが それぞれ予め用意されていてカーネル関数内からアクセスできる blockidx の値を読み出して自分のブロック番号を知る threadidx の値を読み出して自分のブロックにおける自分のスレッド番号を知る griddim には上記Dg の値がそのまま入る blockdim には上記Gb の値がそのまま入る
カーネル関数 (2) GPU(カーネル関数) とCPUは非同期に実行 cudathreadsynchronize() で同期を取る カーネル関数どうしは基本的には同期的 時間の測り方 struct timeval t0,t1; gettimeofday(&t0, NULL); KernelFunction<<<Grid, Thread>>>(Args); cudathreadsynchronize(); gettimeofday(&t1, NULL); printf( Elapsed time = %lf\n, (double)(t1.tv_sec t0.tv_sec) + (double)(t1.tv_usec t0.tv_usec) /(1000*1000)); float f; cudaevent_t event[2]; For (i=0; i<2; i++) cudaeventcreate(&event[i], 0); cudaeventrecord(event[0], stream_1); KernelFunction<<<G,T,0,stream_1>>>(Args); cudaeventrecord(event[1],stream_1); cudaeventsynchronize(event[1]); cudaeventelapsedtime(&f, event[0], event[1]); printf( Elapsed time = %lf\n, f*1000)
関数修飾詞 global device カーネル関数に付く修飾詞 ちなみにカーネル関数の返り値は voidである カーネル関数以外でデバイスのみで実行される関数に付く 同じく void型の関数に付く host ホストで実行される関数の修飾詞 省略可
GPU のメモリ階層 名称 アクセス 速度 レジスタ 自分のスレッドのみ 4 cycle 共有メモリ 同じMultiprocessor内にあるスレッド syncthreads() で同期 デバイスメモリ 同じGPU内にある全てのスレッド 500 cycle 同期にカーネル関数再呼び出しが必要 キャッシュされる 定数メモリ 全てのスレッド キャッシュされる ローカルメモリ レジスタが溢れると使われる デバイスメモリ上に置かれる 12 cycle hit: 20 cycle miss: 500 cycle 500 cycle
メモリ操作関数 cudaerror_t cudasuccess cudaerror_t cudamalloc(void **devptr, size_t count) GPUのデバイスメモリに配列を確保する関数 cudaerror_t cudafree(void *devptr) CUDA関数が成功したときに返る値 確保したデバイスメモリの配列を開放する関数 cudaerror_t cudamemcpy(void *dst, const void *src, size_t count, enum cudamemcpukind kind) src からdstにメモリ転送する関数 kind には 以下の3種類を指定 cudamemcpyhosttodevice cudamemcpydevicetohost cudamemcpydevicetodevice (CPUからGPU) (GPUからCPU) (GPUからGPU) CPU から直接デバイスメモリをいじることはできない
メモリ修飾詞 device constant 定数メモリを修飾する 関数外に宣言する 初期化するかcudaMemcpyToSymbol()で転送する shared デバイスメモリを静的に宣言する場合に修飾させる 関数外に宣言する cudamemcpytosymbol()で転送する 共有メモリを修飾する カーネル関数内に宣言する cudaerror_t cudamemcpytosymbol(const char *symbol, const void *src, sizt_t count, size_t offset, enum cudamemcpykind kind) offset = 0, kind = cudamemcpyhosttodevice
デバイス操作関数 cudaerror_t cudagetdevicecount(int *count) cudaerror_t cudagetdevcieproperties( struct cudadeviceprop *prop, int device) デバイスのプロパティがpropに入って返ってくる関数 cudaerror_t cudasetdevice(int device) 使用可能なデバイスの数がcount に入って返ってくる関数 使用するデバイスを指定する関数 使用できるデバイスはCPU １スレッドにつき １個まで cudaerror_t cudagetdevice(int *device) 現在使用しているデバイス番号がdevice に入って返ってくる関数
サンプルコード ベクトルの和計算 global void AddVector_GPU(float *a, float *b, float *c, int size){ int index = blockidx.x * blockdim.x + threadidx.x; if(index<size) c[index] = a[index] + b[index]; } float *dev_a, *dev_b, *dev_c; cudamalloc((void**)&dev_a, sizeof(float)*n); cudamalloc((void**)&dev_b, sizeof(float)*n); cudamalloc((void**)&dev_c, sizeof(float)*n); dim3 Grid, Block; Block.x = 196; Grid.x = N/196 + 1; cudamemcpy(dev_a, a, sizeof(float)*n, cudamemcpyhosttodevice); cudamemcpy(dev_b, b, sizeof(float)*n, cudamemcpyhosttodevice); AddVector_GPU<<<Grid, Block>>>(dev_a, dev_b, dev_c, N); cudathreadsynchronize(); cudamemcpy(c, dev_c, sizeof(float)*n, cudamemcpydevicetohost); cudafree(dev_a); cudafree(dev_b); cudafree(dev_c);
Warp (1) GPUはSIMD型の命令実行を行う 同じ命令を32個のスレッドが同時に実行する Warp はout-of-order 実行 8つのSPが4cycle かけて実行する 同時実行する32個のスレッド単位をWarp(縦糸)と呼ぶ 準備のできたWarpから実行していく int 型の組み込み変数warpSize がカーネル関数 に用意されている もちろん 32
Warp (2) GPUでは条件分岐の分岐先が全て 実行される 有効なスレッドの結果のみが反映 不必要な命令実行が発生する このことをWarp divergence と呼ぶ なるべくWarp divergence は避ける if(a) B; else C;
デバイスメモリの coalescing メモリアクセスは16スレッドが同時に行う グローバルメモリのアクセス及び転送は以下の３つ (Compute capability 1.2 以上) Warp の半分 32バイト境界でアラインメントされた32バイト 64バイト境界でアラインメントされた64バイト 128バイト境界でアラインメントされた128バイト 同時に行われるメモリアクセスはなるべく少なく済ます 同時に発生するメモリアクセスの場所をできるだけ固まらせ かつでき るだけ同じアラインメント内に納めることを coalescing と呼ぶ
共有メモリ (1) GPUでのソフトウェアキャッシュ グローバルメモリより非常に高速 プログラムに明示的に記述する 同じブロック内のスレッドがアクセスできる カーネル関数に shared 修飾詞を付けて宣言 syncthreads() を用いて同期を取る必要がある
共有メモリ (2) 共有メモリも16 スレッドが同時にアクセスする 32個のバンクで構成されている 各バンクに対して1回につき1つのアドレスにアクセス 16スレッドがそれぞれ別のバンクにアクセスしたい 1つでもバンクが被ることをバンクコンフリクトと呼ぶ バンクコンフリクトはなるべく避ける
初級CUDA (6) 共有メモリ (3) Bank conflict が発生しているケース shared float s[n]; s[threadidx.x*2] = f(x);
初級CUDA (7) 共有メモリ (4) Bank conflict が発生していないケース shared float s[n]; s[threadidx.x] = f(x);
ループアンローリング for ループは展開したい 回数が定数なら nvccが展開してくれる #pragma unroll をforループの直前に書くことで 明示的に展開できる #pragma unroll n とすることで n回分展開される #pragma unroll 1 で展開されなくなる #pragma unroll 4 for(i=0; i<n; i++){ y[i] = f(x[i]); } for(i=0;i<n-4;i+=4){ y[i] = f(x[i]); y[i+1] = f(x[i+1]); y[i+2] = f(x[i+2]); y[i+3] = f(x[i+3]); } for(;i<n;i++){ y[i] = f(x[i]); }
第2回課題 差分法で偏微分方程式の数値シミュレーションをせよ 陽解法を使って2次元拡散方程式を解いてください 初期条件: 境界以外全て 1.0 境界条件: 境界において常に0 0 < r < 0.25 計算領域は任意 上限反復回数は100回 まずはCPUで実装してください GPUを用いて高速化をしてください CPU, GPU ともに FLOPS, GB/s をそれぞれ測定してください 上記方程式１回当たり６FLOPSとしてください
補足 (1) 偏微分方程式 偏導関数を含む等式 2種類以上の変数で偏微分されている 空間方向 時間方向など 常微分方程式は特殊な偏微分方程式 境界条件 初期条件が解を決定づける 解析的に解くのは困難もしくは不可能 計算機を用いて数値解を求める 空間 時間それぞれに関して離散化して数値計算する
補足 (2) 偏微分方程式と有限差分法 連続な偏微分方程式を離散化して解く手法 有限差分法は最も簡易かつ明解な方法の1つ 偏微分されている方向に向かって格子を生成する それぞれの格子の点に対して微分を差分に置き換 えて値を求める
補足 (2) 偏微分方程式と有限差分法 連続な偏微分方程式を離散化して解く手法 有限差分法は最も簡易かつ明解な方法の1つ 偏微分されている方向に向かって格子を生成する それぞれの格子の点に対して微分を差分に置き換 えて値を求める
補足 (3) 偏導関数と離散化誤差 例としてx軸方向の2階偏導関数 Taylor 展開を使う 離散化誤差
補足 (4) 2次元拡散方程式とその離散化 (1) 時間発展型の偏微分方程式 熱や流体中のインクの濃度が拡散していく様子を表す これを離散化する とし 未知数を左辺にもってくる
補足 (5) 2次元拡散方程式とその離散化 (2) 今回の課題では陽解法を使ってもらいます 前の時刻の値を用いて次の時刻の値を求めていく方法 今回の2次元拡散方程式の離散化スキームでは自分を 含めて周りの5つの点から次の自分の値を求めている 赤が現在 青が未来
補足 (6) ダブルバッファリング 時間発展型の偏微分方程式を解く過程において 全ての時間の値を保存していくと容量が足りない 時間を更新していくには最低で2つの領域が必要 片方の領域の数値を計算に使い もう片方の領域 に次の時間の値を埋めていく
課題について renzoku-enshu@is.s.u-tokyo.ac.jp に提出 氏名 学籍番号 出題日を明記 レポートを pdf にした上でメールで提出 演習IIであることも明記 Subject は 演習２第 回レポート 提出者氏名 プログラム+レポート 課題の半分以上の提出が単位の取得条件 全８回を予定 締切は各課題出題日から２週間後 http://olab.is.s.u-tokyo.ac.jp/~reiji/enshu2_12.htmlを参照