ベクトルの内積のGPU実装
今回ベクトルの内積の計算をGPUで行う。要素数は65536*4=262144。ベクトルの各要素は1/1,1/2,1/3,1/4,1/5・・・であり、この場合内積の結果が
に収束することが分かっている。ベクトル長を長くするほど正確なπが求まる。(が、実際はfloat精度の桁数の問題ですぐ頭打ちになる）ということで各要素の掛け算のところを並列に行うことを考える。普通に考えると1コアで262144回の計算を2コアなら131072回、4コアなら65536回と並列度をあげればあげるほど1コアあたりの計算量が減って速く計算ができるはずだ。
まぁその通りなのだが、ある程度並列度を上げるとメモリ律速になって思うように速くならないことがある。というかそんなことだらけだ。
なので今回の記事ではそこのところも少し突っ込んでまとめていこうと思う。
GPUの並列数を指定する個所は２つある。ComputeShader側の
[numthreads(x, y, z)]とC#側の
shader.Dispatch(k, X, Y, Z);だ。
この場合x*X*y*Y*z*Z並列で処理が行われる。
・x,y,zのほうで1グループあたりのスレッド数を指定して並列化
・X,Y,Zのほうでグループ数を指定して並列化
ができる。
じゃあ２つの使い分けはというと、これは説明しだすとかなり長くなるので次の記事で。最低限覚えておくべき知識は、「numthreadsには1～1024までの数字しか指定できない」「とりあえずnumthreadsには64の倍数を指定しておくと良い」というところ。
この知識をふまえ早速GPGPUでバリバリ高速化・・・したいところだがそこは抑えて、まずは「1並列」で処理するコードを書いて答えを確認してみる。
1並列
前回の記事のHoge0をusing System.Collections;
using System.Collections.Generic;
using UnityEngine;
using System;
public class Hoge0 : MonoBehaviour
{
public ComputeShader shader;
void Start()
{
int time0 = Gettime();
int N = 262144;//=2^18
float[] host_vec1 = new float[N];
float[] host_vec2 = new float[N];
float[] host_ans = new float[1];
//初期値をセット
for (int i=0;i< N; i++)
{
host_vec1[i] = 1.0f / (i + 1.0f);// 1/1 , 1/2 , 1/3 , 1/4・・・・
host_vec2[i] = 1.0f / (i + 1.0f);// 1/1 , 1/2 , 1/3 , 1/4・・・・
}
ComputeBuffer vec1 = new ComputeBuffer(host_vec1.Length, sizeof(float));
ComputeBuffer vec2 = new ComputeBuffer(host_vec2.Length, sizeof(float));
ComputeBuffer ans = new ComputeBuffer(1, sizeof(float));
int k = shader.FindKernel("CSMain");
// host to device
int time1 = Gettime() - time0;
vec1.SetData(host_vec1);
vec2.SetData(host_vec2);
int time2 = Gettime() - time0;
//引数をセット
shader.SetBuffer(k, "vec1", vec1);
shader.SetBuffer(k, "vec2", vec2);
shader.SetBuffer(k, "ans", ans);
// GPUで計算
int time3 = Gettime() - time0;
shader.Dispatch(k, 1, 1, 1);
int time4 = Gettime() - time0;
// device to host
ans.GetData(host_ans);
int time5 = Gettime() - time0;
//計測時間表示
Debug.Log("CPU→GPU転送前\t" + time1);
Debug.Log("CPU→GPU転送後\t" + time2);
Debug.Log("Dispatch前\t" + time3);
Debug.Log("Dispatch後\t" + time4);
Debug.Log("GPU→CPU転送後\t" + time5);
//結果表示
float calc_pi = Mathf.Sqrt(host_ans[0] * 6.0f);
Debug.Log("π =" + calc_pi.ToString("f10"));
//解放
vec1.Release();
vec2.Release();
ans.Release();
}
// Update is called once per frame
void Update()
{
}
//現在の時刻をms単位で取得
int Gettime()
{
return DateTime.Now.Millisecond + DateTime.Now.Second * 1000
+ DateTime.Now.Minute * 60 * 1000 + DateTime.Now.Hour * 60 * 60 * 1000;
}
}とし
HogeCS0を
#pragma kernel CSMain //Read and Write RWStructuredBuffer<float>vec1; RWStructuredBuffer<float> vec2; RWStructuredBuffer<float> ans; [numthreads(1, 1, 1)] void CSMain(uint3 id : SV_DispatchThreadID) { float s = 0.0; for (int i = 0; i < 262144; i++) {//=2^18 s += vec1[i] * vec2[i]; } ans[0] = s; }
として実行。
πが3.141くらいまで求まった。成功だ！
GPUの計算時間計測の落とし穴
GPUの計算時間についてはどうか。一番最初の時刻を0msとしており、命令が実行される時点での時刻(ms)を右に表示してある。上のはった画像を見るに、CPU→GPUのデータ転送は1ms。GPUの計算はなんと0msで行えているようだ。そしてGPU→CPUの転送に75msかかっている・・・ーーと思うがこれは間違い！
実際はGPUの計算に75msかかっている！Dispatch命令はGPUに向かってタスクを投げるだけであって、GPUの計算が終了するまで待つ命令ではない！！
【イメージ】
だからどれだけGPUのタスクが重かろうと、
// GPUで計算
int time3 = Gettime() - time0;
shader.Dispatch(k, 1, 1, 1);
int time4 = Gettime() - time0;このtime3とtime4の差はほとんどない。個人的な感覚では0.3μsくらいだ。
そして
// device to host
ans.GetData(host_ans);ここでGPUの計算が終わるのを待ち、そのあとGPU→CPUの転送が行われるというわけ。ちなみにGPUの計算が終わるまで待つという命令(OpenCLで言うclFlush)はないのかというと、私が探した限りでは見つけられなかった。あればGPUの計算時間を直接計算できるのだが。今はGetDataで代用するしかないようだ。
データ転送時間の計算
計測時間についてさらに考察を深めてみよう。時間がない人は「ちゃんとした並列計算」まで飛ばしても構わない。CPU→GPU、GPU→CPUのデータ転送はPCI Expressを経由している。
2019年1月現在PCI Express Gen3.0がもっとも普及しておりGPUとの接続では基本PCI Express Gen3.0 x16となっているだろう。これは片方向16GB/sの帯域速度がでる。今回のプログラムではCPU→GPUに262144要素×4byte×2のデータ転送を行っているが、これはたったの2MBである。この転送にかかる時間は理論上0.122msとなり、実際に計測された1msでは遅いくらいだ（オーバーヘッドが含まるからおかしくはない）。
clock単位の計算時間計測
次にGPU処理部分の時間について考える。今回1並列で愚直に計算しており、せっかくRX480には2304個の演算器があるのに2303個は空回りしているという状態。ただ丁度よい機会なので1演算器あたりの性能をここで計算してみよう。さっきのCSMainだが、コードをみると計算の大部分はループ内の[numthreads(1, 1, 1)] void CSMain(uint3 id : SV_DispatchThreadID) { float s = 0.0; for (int i = 0; i < 262144; i++) {//=2^18 s += vec1[i] * vec2[i]; } ans[0] = s; }
である。75msで262144ループなので1ループ当たり286nsかかっている。s += vec1[i] * vec2[i];
大雑把にGPUのコアクロックを1Ghzとすると1clock=1nsなので、この1行に286clockもかかっていることになる！これは遅い！？
実はFP32演算器は1clockでa+b , a-b , a*b , a*b+cの計算を行うことができる。とするとこの行の計算自体は1clockで終わるはずで、vec1[i],vec2[i]のメモリアクセスがどうも原因であることに気づく。つまりメモリアクセスに約280clockかかっていて、その間FP32演算器が「待ちぼうけ」を喰らっているわけだ。なんという無駄か！
ちゃんとした並列計算
1グループあたりのスレッド数を増やして並列度を上げる
ではまずComputeShaderの[numthreads(1, 1, 1)]を
[numthreads(32, 1, 1)]に置き換える。
コアレスアクセス
そして高速化において「連続したアドレスにアクセスすること」が鍵になる。上記の実験にて演算が1clockに対しメモリアクセスが280clockかかっていたことからも、メモリアクセスの高速化がいかに大事かわかる。今まで話に出てきたメモリは「グローバルメモリ」といって、これに効率よくアクセスするには「コアレスアクセス」が必要だ。（グローバルメモリ＝実体はGDDR5とかGDDR6とかHBM2）
[numthreads(32, 1, 1)]
void CSMain(uint3 id : SV_DispatchThreadID)
{
float s = 0.0;
for (int i = 0; i < 8192; i++) {
s += vec1[id.x + i*32] * vec2[id.x + i * 32];
}
ans[id.x] = s;
}これがコアレスアクセスしている32並列のComputeShaderコードになる。1スレッドあたり8192回ループを行っているが、ベクトル長が262144要素で32並列なのでこうなる。ここでvec1,vec2の添え字がid.x+i*32になっているのに注目。
iが0のとき
これがコアレスアクセス。一度に連続したデータをとってくる。
iが1のとき
これもコアレスアクセス。
コアレスアクセスでは、一度のメモリアクセスに必要な時間(≒レイテンシ)は変わらないが、一度にとってこれるデータ量が増えるためTotalで見ると高速になる。この矢印がてんでバラバラな領域に向かっているとパフォーマンスが落ちる。同じ領域内なら矢印がバラバラでもパフォーマンスは落ちない ←ここ重要！
これができているかできていないかで、同じ並列数でも実行速度が全然違ってくる。記事の下のほうで実験してるが4-5倍ほど平気で変わってくる。
Hoge0(C#)も修正
using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using System;
public class Hoge0 : MonoBehaviour
{
public ComputeShader shader;
void Start()
{
int time0 = Gettime();
int p = 32;//並列数
int N = 65536* 4;
float[] host_vec1 = new float[N];
float[] host_vec2 = new float[N];
float[] host_ans = new float[p];
//初期値をセット
for (int i=0;i< N; i++)
{
host_vec1[i] = 1.0f / (i + 1.0f);// 1/1 , 1/2 , 1/3 , 1/4・・・・
host_vec2[i] = 1.0f / (i + 1.0f);// 1/1 , 1/2 , 1/3 , 1/4・・・・
}
ComputeBuffer vec1 = new ComputeBuffer(host_vec1.Length, sizeof(float));
ComputeBuffer vec2 = new ComputeBuffer(host_vec2.Length, sizeof(float));
ComputeBuffer ans = new ComputeBuffer(host_ans.Length , sizeof(float));
int k = shader.FindKernel("CSMain");
// host to device
int time1 = Gettime() - time0;
vec1.SetData(host_vec1);
vec2.SetData(host_vec2);
int time2 = Gettime() - time0;
//引数をセット
shader.SetBuffer(k, "vec1", vec1);
shader.SetBuffer(k, "vec2", vec2);
shader.SetBuffer(k, "ans", ans);
// GPUで計算
int time3 = Gettime() - time0;
shader.Dispatch(k, 1, 1, 1);
int time4 = Gettime() - time0;
// device to host
ans.GetData(host_ans);
int time5 = Gettime() - time0;
//計測時間表示
Debug.Log("CPU→GPU転送前\t" + time1);
Debug.Log("CPU→GPU転送後\t" + time2);
Debug.Log("Dispatch前 \t" + time3);
Debug.Log("Dispatch後 \t" + time4);
Debug.Log("GPU→CPU転送後\t" + time5);
//結果表示
float calc_pi = 0.0f;
for (int i=0;i< p; i++)//最後の最後の結果はCPUで加算
{
calc_pi += host_ans[i];
}
calc_pi = Mathf.Sqrt(calc_pi * 6.0f);
Debug.Log("π =" + calc_pi.ToString("f10"));
//解放
vec1.Release();
vec2.Release();
ans.Release();
}
// Update is called once per frame
void Update()
{
}
//現在の時刻をms単位で取得
int Gettime()
{
return DateTime.Now.Millisecond + DateTime.Now.Second * 1000
+ DateTime.Now.Minute * 60 * 1000 + DateTime.Now.Hour * 60 * 60 * 1000;
}
}そして計算時間は
13-6=7ms。結構縮んだ。
75ms→7msで10倍高速化。
もちろんこの7msの中にGPU→CPU転送時間も含まれているから実際はもう少し短い。
内積の最後の足し算だが、今まではGPU側のans[0]に結果を全部まとめていたが、今回はGPU側でans[0]～ans[31]に各スレッドの結果を代入し、CPU側でans[0]～ans[31]を合計している。
GPU側で最後の32要素の合計をすることもできるがやや難易度が上がるので次の機会とする。
グループ数を増やして並列度を上げる
ではC#側のDispatch関数でshader.Dispatch(k, 128, 1, 1);こうして
128*32=4096並列で実行するとしよう。そうするとans[0]～ans[4095]まで必要になるがそこも修正して計算時間を計ってみると・・・
(コードは省略)
7-6=1ms。は、速い！75ms→7ms→1msと少なくとも75倍高速化できており、時間解像度の問題でもはや正確に測定できてるのか疑うレベル。
今回の高速化ではコアレスアクセスのところはいじってないため、単に立ち上げるスレッド数が32→4096に増えたことによるものだろう。細かくは私もわからないが、多くのスレッドからメモリアクセス要求があったほうが、その帯域を使い切れるといったイメージだ。きっとさっきの32並列のやつだけではその帯域を全然使い切れてなかったということなのではなかろうか。
4096並列のときの計算時間再測定
さて今度は計算時間をちゃんと測定するため、もっと計算の規模を大きくしなければいけなさそうだ。Nを65536*4から512倍の65536*2048にする。Hoge0(C#)
using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using System;
public class Hoge0 : MonoBehaviour
{
public ComputeShader shader;
void Start()
{
int time0 = Gettime();
int bp = 32;//ComputeShaderで指定している並列数
int gp = 128;//Dispatchで指定している並列数
int N = 65536*2048;
float[] host_vec1 = new float[N];
float[] host_vec2 = new float[N];
float[] host_ans = new float[bp*gp];
//初期値をセット
for (int i=0;i< N; i++)
{
host_vec1[i] = 1.0f / (i + 1.0f);// 1/1 , 1/2 , 1/3 , 1/4・・・・
host_vec2[i] = 1.0f / (i + 1.0f);// 1/1 , 1/2 , 1/3 , 1/4・・・・
}
ComputeBuffer vec1 = new ComputeBuffer(host_vec1.Length, sizeof(float));
ComputeBuffer vec2 = new ComputeBuffer(host_vec2.Length, sizeof(float));
ComputeBuffer ans = new ComputeBuffer(host_ans.Length , sizeof(float));
int k = shader.FindKernel("CSMain");
// host to device
int time1 = Gettime() - time0;
vec1.SetData(host_vec1);
vec2.SetData(host_vec2);
int time2 = Gettime() - time0;
//引数をセット
shader.SetBuffer(k, "vec1", vec1);
shader.SetBuffer(k, "vec2", vec2);
shader.SetBuffer(k, "ans", ans);
// GPUで計算
int time3 = Gettime() - time0;
shader.Dispatch(k, gp, 1, 1);
int time4 = Gettime() - time0;
// device to host
ans.GetData(host_ans);
int time5 = Gettime() - time0;
//計測時間表示
Debug.Log("CPU→GPU転送前\t" + time1);
Debug.Log("CPU→GPU転送後\t" + time2);
Debug.Log("Dispatch前 \t" + time3);
Debug.Log("Dispatch後 \t" + time4);
Debug.Log("GPU→CPU転送後\t" + time5);
//結果表示
float calc_pi = 0.0f;
for (int i=0;i< bp * gp; i++)//最後の最後の結果はCPUで加算
{
calc_pi += host_ans[i];
}
calc_pi = Mathf.Sqrt(calc_pi * 6.0f);
Debug.Log("π =" + calc_pi.ToString("f10"));
//解放
vec1.Release();
vec2.Release();
ans.Release();
}
// Update is called once per frame
void Update()
{
}
//現在の時刻をms単位で取得
int Gettime()
{
return DateTime.Now.Millisecond + DateTime.Now.Second * 1000
+ DateTime.Now.Minute * 60 * 1000 + DateTime.Now.Hour * 60 * 60 * 1000;
}
}HogeCS0(ComputeShader)
#pragma kernel CSMain //Read and Write RWStructuredBuffer<float>vec1; RWStructuredBuffer<float> vec2; RWStructuredBuffer<float> ans; [numthreads(32, 1, 1)] void CSMain(uint3 id : SV_DispatchThreadID) { float s = 0.0; for (int i = 0; i < 32768; i++) {//65536*2048/32/128=32768 s += vec1[id.x + i * 4096] * vec2[id.x + i * 4096]; } ans[id.x] = s; }
これで実行すると
2794-2503=261ms
512倍にしているので実際は261/512=0.51msだったことがわかる。
まとめると
1並列：75ms
32並列：7ms
4096並列：0.51ms
となり、最初の1並列プログラムと比較すると149倍高速になったことになる。
ソースはこちら
https://github.com/toropippi/Uinty_GPGPU_SharedMem_AtomicAdd/tree/master/gpgpu_sample2
πが全然収束してないがこれはfloat精度なので仕方ない。double型で計算し直したところ3.141592までは正確に計算できた。
ランダムアクセスは遅い
最後に、コアレスアクセスでない方法を試す。先ほどの4096並列の512倍規模にしていたやつのコードを流用する。今度は添え字がid.x*32768 + iになっている。[numthreads(32, 1, 1)] void CSMain(uint3 id : SV_DispatchThreadID) { float s = 0.0; for (int i = 0; i < 32768; i++) {//65536*2048/32/128=32768 s += vec1[id.x* 32768 + i] * vec2[id.x* 32768 + i]; } ans[id.x] = s; }
iが0のとき、こうなる。iが1のときも同様、離れている個所へのメモリアクセスとなり、これはランダムアクセスといって効率が悪い。
この方式でやると・・・
3713-2544=1169ms
やはり遅くなっている！
261msから1169msなので4.48倍遅い。
まとめ
ベクトルの内積をやるといって、結局メモリの話ばかりしてたような気がする。ただ実際、GPUの計算はグローバルメモリのアクセスがボトルネックになってくるので、パフォーマンスに直結するので大事なことだと思っている。※本ブログで紹介しているソースコードはすべてCC0 パブリックドメインです。
ご自由にお使いください
前
次
UnityでGPGPUその3 ComputeShaderの共有メモリ(shared memory)
