ベクトルの加算2011年1月8日 | |||||||||||||||||||
ファイル内容ベクトルどうしの加算は、各要素ごとに独立して行うことができるため、GPU の利用に高速化できることが期待できます。 ベースの C プログラムaddVector.cu
void addVector(const float *a, const float *b, float *c)
{
int i;
for(i = 0; i < N; i++){
c[i] = a[i] + b[i];
}
}
int main(void)
{
float *a;
float *b;
float *c;
...
a = (float *)malloc(N*sizeof(float));
b = (float *)malloc(N*sizeof(float));
c = (float *)malloc(N*sizeof(float));
...
addVector(a, b, c);
...
free(a);
free(b);
free(c);
return 0;
}
これを CUDA 用に書き換えて行きましょう。 CUDA 用のプログラムaddVector2.cu
...
__global__ void addVector(const float *a, const float *b, float *c)
{
int i = blockDim.x*blockIdx.x + threadIdx.x;
c[i] = a[i] + b[i];
}
int main(void)
{
float *a;
float *b;
float *c;
float *a_d;
float *b_d;
float *c_d;
...
cudaMallocHost(&a, N*sizeof(float));
cudaMallocHost(&b, N*sizeof(float));
cudaMallocHost(&c, N*sizeof(float));
cudaMalloc(&a_d, N*sizeof(float));
cudaMalloc(&b_d, N*sizeof(float));
cudaMalloc(&c_d, N*sizeof(float));
...
cudaMemcpy(a_d, a, N*sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(b_d, b, N*sizeof(float), cudaMemcpyHostToDevice);
...
addVector<<<N/BLOCK_SIZE, BLOCK_SIZE>>>(a_d, b_d, c_d);
...
cudaMemcpy(c, c_d, N*sizeof(float), cudaMemcpyDeviceToHost);
...
cudaFreeHost(a);
cudaFreeHost(b);
cudaFreeHost(c);
cudaFree(a_d);
cudaFree(b_d);
cudaFree(c_d);
return 0;
}
まずベクトル 3 つ分のメモリをホスト側にページロックで確保しています。そして同じサイズのベクトルを今度はデバイス側で確保しています。ベクトル a, b に値をセットして、そのデータをデバイスに転送し、addVector() でベクトルの加算を行った後、デバイス側にあるベクトル c のデータをホストに転送しています。 ポイントは addVector() の定義の仕方と呼び出し方です。
__global__ void addVector(const float *a, const float *b, float *c)
{
int i = blockDim.x*blockIdx.x + threadIdx.x;
c[i] = a[i] + b[i];
}
"__global__" というのは、デバイス側にありホストから呼び出すことができる関数であるという意味の指定です。つまり addVector() はデバイスの処理を記述したカーネルのコードなのです。 ベクトルの要素の和についてのループがなくなっています。CPU で計算する場合は 1 つの CPU がベクトルの要素を 1 つずつ繰り返し加算することになりますが、GPU では GPU 内のたくさんのプロセッサがそれぞれベクトルの 1 要素について和を計算するだけで済みます。 空っぽの部屋の床の雑巾がけを想像してみるとわかりやすいかもしれません。一人で行う場合は部屋を繰り返し往復する必要があります。部屋の幅いっぱいの人間がいれば、端から端まで一度行くだけで済んでしまいます。前者が CPU のイメージ、後者が GPU のイメージです。 i の計算における blockDim.x はブロックにおけるスレッド数、blockIdx.x はブロックの ID、threadIdx.x はスレッドの ID です。実際に計算を行うのはスレッドです。addVector() は個々のスレッドに対して呼び出されます。ここではスレッドの所属するブロックの ID とその中のスレッド自身の ID から、スレッドが担当するベクトルの要素を割り当てているのです。 addVector<<<N/BLOCK_SIZE, BLOCK_SIZE>>>(a_d, b_d, c_d); addVector() の呼び出しには <<< ... >>> という謎の記号が付いています。この記号の中ではブロック数 (グリッドサイズ) と各ブロックのスレッド数 (ブロックサイズ) が指定されています。BLOCK_SIZE はブロックがもつスレッドの数で、任意の値です。任意とはいってもデバイスによって最大値の制限はあります (ブロック数についても同様)。N はベクトルのサイズで、上はそれぞれ BLOCK_SIZE 個のスレッドをもつ N/BLOCK_SIZE 個のブロックについて addVecotr() を呼び出します、という意味です。 上の雑巾がけの例でいえば、数名ずつで班になっており、それぞれの班に対して雑巾がけの担当区域が決められ、班の中でそれぞれの人に雑巾がけの担当位置が割り当てられる、というイメージです。 マップドメモリを使うCUDA プログラミングではいちいちホスト - デバイス間の転送処理を記述しなければなりませんが、マップドメモリというものを使うとその必要がなくなります。 addVector3.cu ... cudaSetDeviceFlags(cudaDeviceMapHost); cudaHostAlloc(&a, N*sizeof(float), cudaHostAllocMapped); cudaHostAlloc(&b, N*sizeof(float), cudaHostAllocMapped); cudaHostAlloc(&c, N*sizeof(float), cudaHostAllocMapped); cudaHostGetDevicePointer(&a_d, a, 0); cudaHostGetDevicePointer(&b_d, b, 0); cudaHostGetDevicePointer(&c_d, c, 0); ... addVector<<<N/BLOCK_SIZE, BLOCK_SIZE>>>(a_d, b_d, c_d); ... cudaThreadSynchronize(); ... cudaFreeHost(a); cudaFreeHost(b); cudaFreeHost(c); ... マップドメモリを使うと、ホストとデバイスで (見た目上) 共通のメモリ領域を使うことができます。とはいえホスト - デバイス間のデータ転送はハード的には必要であり、それについては内部で勝手に処理してくれます。cudaSetDeviceFlags() でマップドメモリを使えるようにフラグをセットしています。cudaHostAlloc() の cudaHostAllocMapped でホスト側にマップドメモリを確保し、cudaHostGetDevicePointer() でデバイス側のポインタを取得しています。cudaFreeHost() でマップドメモリを解放していますが、デバイス側の領域は確保していないので解放の必要がありません。cudaThreadSynchronize() はスレッドの同期をとるための関数で、マップドメモリにデータの同期のタイミングを教えるために入れています。 CUBLAS を使うCUDA には "CUBLAS" という CUDA 用の BLAS ライブラリが用意されています。その中にベクトルの加算を行う関数があるので、それを使ってみましょう。 addVector4.cu
...
#include <cublas.h>
...
int main(void)
{
float *a;
float *b;
float *c;
float *a_d;
float *b_d;
...
cublasInit();
cudaMallocHost(&a, N*sizeof(float));
cudaMallocHost(&b, N*sizeof(float));
cudaMallocHost(&c, N*sizeof(float));
cublasAlloc(N, sizeof(float), (void **)&a_d);
cublasAlloc(N, sizeof(float), (void **)&b_d);
...
cublasSetVector(N, sizeof(float), a, 1, a_d, 1);
cublasSetVector(N, sizeof(float), b, 1, b_d, 1);
...
cublasSaxpy(N, 1., a_d, 1, b_d, 1);
...
cublasGetVector(N, sizeof(float), b_d, 1, c, 1);
...
cudaFreeHost(a);
cudaFreeHost(b);
cudaFreeHost(c);
cublasFree(a_d);
cublasFree(b_d);
cublasShutdown();
return 0;
}
CUBLAS を使ったコードのコンパイルにはライブラリのリンクの指定が必要です。 $ nvcc -o addVector4 addVector4.cu -lcublas CUBLAS を使うために cublasInit() で初期化、cublasShutdown() で完了処理を行っています。cublasAlloc() で CUBLAS 用の領域を確保、cublasSetVector()/cublasGetVector() でデータの設定と取得 (データ転送の代わり)、cublasFree() で領域を解放しています。cublasSaxpy() ベクトルの和を求めています。a + b を計算し、結果は b に格納されます (2 番目の引数は a にかかる係数)。 速度の比較上で見た 4 つのプログラムについて実行速度を比較します。処理に数秒の時間がかかるように処理を繰り返すようにしています。 マシンスペック
計算結果
注意: CUDA を使用するプログラムは、たぶん初期化のためだと思うが、たまに数秒余計に時間がかかることがある。 GPU 計算により CPU 計算のおよそ 12 倍の高速化が実現できています。マップドメモリはコードは少しわかりやすくなるものの、結構な時間を要しています。CUBLAS がお手製 GPU コードより時間がかかっているのは、おそらく cublasSaxpy() が純粋にベクトルの和を行う関数ではなく係数の乗算を含んでいるからだろうと思います。それでも CPU 計算の 7 倍は高速になっており、GPU コードについて悩む必要がない分、CUBLAS を使うのも悪くはないと思います。 | |||||||||||||||||||
| PENGUINITIS | |||||||||||||||||||