CUDAで行列演算 その2 乗算†
CUDAで行列の乗算を行います．このページで示すコードは CUDA Programming Guideの
3.2.2 Shared Memoryとほとんど同じですのでそちらも参照してください．
行列の定義†
行列の積は和や差と異なり，それぞれの行列のサイズが異なります．
A(N×M) * B(M×L) = C(N×L)
CUDAで行列演算：加減算ではすべての行列のサイズが等しかったので楽でしたが，
行列のサイズが異なると管理が面倒になるので，新しくMatrix構造体を定義します．
1
2
3
4
5
6
| | struct Matrix
{
int width;
int height;
float *elements;
};
|
構造体の定義はC++風で問題ないようです(いちいちtypedefつかわなくてよい)．
widthとheightが行列の列数と行数で，elementsが行列の要素を格納する配列です．
上記の行列Aのサイズだと，width=M，height=Nです．
ホスト側の行列の定義は以下のようにしました．
ただし，行列のサイズは適当です(面倒が少ない用にブロックサイズの倍数にはしてますが)．
1
2
3
4
5
6
7
8
9
10
11
| | Matrix hA, hB, hC;
hA.height = hC.height = 3*BLOCK_SIZE;
hA.width = hB.height = 4*BLOCK_SIZE;
hB.width = hC.width = 5*BLOCK_SIZE;
hA.elements = new float[hA.width*hA.height];
hB.elements = new float[hB.width*hB.height];
hC.elements = new float[hC.width*hC.height];
RandomInit(hA.elements, hA.width*hA.height);
RandomInit(hB.elements, hB.width*hB.height);
|
RandomInitはCUDAで行列演算：加減算のものと同じで，テスト用に乱数を行列の要素に代入する関数です．
デバイス側のメモリ確保とデータ転送も前回と同様に
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
| | Matrix dA, dB, dC;
dA.width = hA.width; dA.height = hA.height;
dB.width = hB.width; dB.height = hB.height;
dC.width = hC.width; dC.height = hC.height;
int size;
size = dA.width*dA.height*sizeof(float);
cutilSafeCall(cudaMalloc((void**)&dA.elements, size));
cutilSafeCall(cudaMemcpy(dA.elements, hA.elements, size, cudaMemcpyHostToDevice));
size = dB.width*dB.height*sizeof(float);
cutilSafeCall(cudaMalloc((void**)&dB.elements, size));
cutilSafeCall(cudaMemcpy(dB.elements, hB.elements, size, cudaMemcpyHostToDevice));
size = dC.width*dC.height*sizeof(float);
cutilSafeCall(cudaMalloc((void**)&dC.elements, size));
|
カーネル†
カーネルの実行は前回と同様です．
1
2
3
4
5
6
7
| | dim3 block(BLOCK_SIZE, BLOCK_SIZE);
dim3 grid((dC.width+block.x-1)/block.x, (dC.height+block.y-1)/block.y);
matrixMul<<< grid, block >>>(dA, dB, dC);
cutilCheckMsg("Kernel execution failed");
|
カーネル名はmatrixMulとしています．
カーネルは行列Cのサイズと同じ数実行され，各カーネルスレッドは行列Cのそれぞれの要素を計算します．
カーネル関数の実装は，
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
| | __global__
void matrixMul(Matrix A, Matrix B, Matrix C)
{
int row = blockIdx.y*blockDim.y+threadIdx.y;
int col = blockIdx.x*blockDim.x+threadIdx.x;
if(row < C.height && col < C.width){
float x = 0.0f;
for(int k = 0; k < A.width; ++k){
x += A.elements[row*A.width+k]*B.elements[k*B.width+col];
}
C.elements[row*C.width+col] = x;
}
}
|
行列Cの各要素をAのrow行とBのcol列からfor文を用いて計算しています．
ホストへの結果の転送と後片付け†
結果をホストに転送し，確保したメモリを解放します．
1
2
3
4
5
6
7
8
9
10
11
12
13
| | size = dC.width*dC.height*sizeof(float);
cutilSafeCall(cudaMemcpy(hC.elements, dC.elements, size, cudaMemcpyDeviceToHost));
cutilSafeCall(cudaFree(dA.elements));
cutilSafeCall(cudaFree(dB.elements));
cutilSafeCall(cudaFree(dC.elements));
delete [] hA.elements;
delete [] hB.elements;
delete [] hC.elements;
|
計算時間†
CUDAで実装したコードの計算時間を計算して，CPU実装と比較してみました．実行環境は以下です．
- GPU : GeForce GTX285
- CPU : Core 2 Duo 3.16GHz
CPUの方が少し古くGPU有利ですが参考データ程度にということで．
時間の計測にはSDKの関数を用いました．
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
| | unsigned int timer = 0;
cutCreateTimer(&timer);
cutStartTimer(timer);
cutilSafeCall(cudaThreadSynchronize());
cutStopTimer(timer);
printf("Processing time: %f (ms) \n", cutGetTimerValue(timer));
cutDeleteTimer(timer);
}
cudaThreadSynchronise()はデバイス側の処理が全て終わるまで待つ関数です．
CPUでの処理用に以下のような関数を用意しました．
#code(C){{
void matrixMulByCPU(const Matrix A, const Matrix B, Matrix C)
{
for(int i = 0; i < C.height; ++i){
for(int j = 0; j < C.width; ++j){
float cval = 0.0f;
for(int k = 0; k < A.width; ++k){
cval += A.elements[i*A.width+k]*B.elements[k*B.width+j];
}
C.elements[i*C.width+j] = cval;
}
}
}
|
行列の行と列のサイズを同じとし，その大きさはブロックサイズを16×16としてその1倍,10倍,100倍として計算してみました．
計算時間は以下の表です．
| 16 | 160 | 1600 |
| GPU | 0.085 | 0.42 | 398 |
| CPU | 0.010 | 12.58 | 21245 |
時間の単位は ms(ミリ秒)です．一回計測しただけなので正確性は低いです．
GPUでの計算部分のみ(カーネル呼出し部のみ)を計測していて，デバイスメモリの確保・解放，デバイスメモリへの転送などは含めていませんが，
それでも行列のサイズが十分大きいときにはGPUでの実装が有用であると考えられます．
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
