行列の乗算

2011年1月8日

ファイル

03-multMatrix.tar.gz

内容

行列の乗算を GPU で行ってみます。

ベースの C プログラム

matMult.cu

...
void multMatrix(const float *a, const float *b, float *c)
{
	int i, j, k;

	for(i = 0 ; i < N; i++){
		for(j = 0; j < N; j++){
			float v;
			v = 0.;
			for(k = 0; k < N; k++){
				v += a[i*N + k]*b[k*N + j];
			}
			c[i*N + j] = v;
		}
	}
}

int main(void)
{
	float *a;
	float *b;
	float *c;
	
	...

	a = (float *)malloc(N*N*sizeof(float));
	b = (float *)malloc(N*N*sizeof(float));
	c = (float *)malloc(N*N*sizeof(float));

	...

	multMatrix(a, b, c);
	
	...
	
	free(a);
	free(b);
	free(c);

	return 0;
}

CUDA 用のプログラム

multMatrix2.cu

...
__global__ void multMatrix(const float *a, const float *b, float *c)
{
	int i = blockDim.x*blockIdx.x + threadIdx.x;
	int j = blockDim.y*blockIdx.y + threadIdx.y;
	int k;
	float v;

	v = 0.;
	for(k = 0; k < N; k++){
		v += a[i*N + k]*b[k*N + j];
	}

	c[i*N + j] = v;
}

int main(void)
{
	float *a;
	float *b;
	float *c;
	float *a_d;
	float *b_d;
	float *c_d;
	dim3 grid;
	dim3 block;
	
	...

	cudaMallocHost(&a, N*N*sizeof(float));
	cudaMallocHost(&b, N*N*sizeof(float));
	cudaMallocHost(&c, N*N*sizeof(float));

	cudaMalloc(&a_d, N*N*sizeof(float));
	cudaMalloc(&b_d, N*N*sizeof(float));
	cudaMalloc(&c_d, N*N*sizeof(float));
	
	...

	cudaMemcpy(a_d, a, N*N*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(b_d, b, N*N*sizeof(float), cudaMemcpyHostToDevice);

	grid = dim3(N/BLOCK_SIZE, N/BLOCK_SIZE);
	block = dim3(BLOCK_SIZE, BLOCK_SIZE);
	multMatrix<<<grid, block>>>(a_d, b_d, c_d);

	cudaMemcpy(c, c_d, N*N*sizeof(float), cudaMemcpyDeviceToHost);

	...

	cudaFreeHost(a);
	cudaFreeHost(b);
	cudaFreeHost(c);

	cudaFree(a_d);
	cudaFree(b_d);
	cudaFree(c_d);

	return 0;
}

処理の手順はベクトルの加算の場合と同じです。

__global__ void multMatrix(const float *a, const float *b, float *c)
{
	int i = blockDim.x*blockIdx.x + threadIdx.x;
	int j = blockDim.y*blockIdx.y + threadIdx.y;
	int k;
	float v;

	v = 0.;
	for(k = 0; k < N; k++){
		v += a[i*N + k]*b[k*N + j];
	}

	c[i*N + j] = v;
}

ベクトルの加算の場合は各ベクトルの要素ごとにスレッドを割り当てればよかったのですが、今回は行列です。そして乗算なのでアクセスする行列の要素はそれぞれバラバラです。今回は乗算結果を格納する行列の要素ごとにスレッドを割り当てます。スレッドの割り当て方ですが、ブロックやスレッドは 2 次元的に管理することも可能なので、今回はそれを使います。たとえばスレッド ID は (x, y) という形になり、それぞれの値は threadIdx.x や threadIdx.y で取得できます。

	dim3 grid;
	dim3 block;
	
	...
	
	grid = dim3(N/BLOCK_SIZE, N/BLOCK_SIZE);
	block = dim3(BLOCK_SIZE, BLOCK_SIZE);
	multMatrix<<<grid, block>>>(a_d, b_d, c_d);

多次元的にブロックやスレッドを管理する場合は dim3 という型を使います。名前の通り 本来 3 つの値を持ちます。

	grid = dim3(N/BLOCK_SIZE, N/BLOCK_SIZE);
	block = dim3(BLOCK_SIZE, BLOCK_SIZE);

dim3() で dim3 型の値を作ることが可能で、2 次元の場合は 3 番目の引数を省略できます。上ではブロックを BLOCK_SIZE×BLOCK_SIZE 個のスレッドに分け、グリッドを N/BLOCK_SIZE×N/BLOCK_SIZE 個のブロックに分けるように指示しています。

大人数で巨大な絵を描く場面を想像してみましょう。絵を格子状に分け、それぞれの区域を各班に割り当てます。それぞれの班の中でまた区域を格子状に分け、それぞれの区域を担当するメンバーを決定します。そういうイメージです。

あなたが班 (1, 2) のメンバー (3, 5) だと想像しましょう。各班の大きさを 16×16 とすると、あなたは棚 C(19, 37) の製品作りを担当するこになります。さて、C(19, 37) の値を作るためにあなたは何をする必要があるでしょうか。向こうに A の棚と B の棚が見えます。C(19, 37) の製品を作るためには A の 19 行の工具と B の 37 列の工具が必要だとあなたは考えます。あなたは A 棚へ近づき、 その 19 行目の横一列の棚を眺めます。おや、19 行はみな空っぽです。どうやら C(19, 28) の担当の人が使用中のようです。しかたがない、あなたは C(19, 28) 担当の人が仕事を終えるのを待たなければなりません。

ベクトルの和ではそれぞれの要素の計算はみな独立していたのですが、行列の積では計算結果を入れる行列の要素間で使用するデータが重なっています。したがって上記のようにあるスレッドが待ちぼうけを食わされる可能性があります。これがネックとなって、行列の積のような計算では GPU を使ってもベクトルの和のような劇的な効果は得られないと予想できます。

CUBLAS を使う

CUBLAS の行列の乗算の関数を使ってみましょう。

multMatrix3.cu

...
#include <cublas.h>

...

int main(void)
{
	float *a;
	float *b;
	float *c;
	float *a_d;
	float *b_d;
	float *c_d;
	
	...
	
	cublasInit();

	cudaMallocHost(&a, N*N*sizeof(float));
	cudaMallocHost(&b, N*N*sizeof(float));
	cudaMallocHost(&c, N*N*sizeof(float));

	cublasAlloc(N*N, sizeof(float), (void **)&a_d);
	cublasAlloc(N*N, sizeof(float), (void **)&b_d);
	cublasAlloc(N*N, sizeof(float), (void **)&c_d);

	...

	cublasSetVector(N*N, sizeof(float), a, 1, a_d, 1);
	cublasSetVector(N*N, sizeof(float), b, 1, b_d, 1);

	cublasSgemm('n', 'n', N, N, N, 1., a_d, N, b_d, N, 0., c_d, N);

	cublasGetVector(N*N, sizeof(float), c_d, 1, c, 1);

	...

	cudaFreeHost(a);
	cudaFreeHost(b);
	cudaFreeHost(c);

	cublasFree(a_d);
	cublasFree(b_d);
	cublasFree(c_d);

	cublasShutdown();

	return 0;
}

cublasSgemm() が行列の積を計算する関数です。

行列の格納方法を変える

お手製コードと CUBLAS では行列の格納方法が異なっています (a(i,j) に対して前者は a[i*N + j]、後者は a[j*N + i])。配列を CUBLAS 型に変えたときに性能に変化が出るか調べてみましょう。

multMatrix4

...
__global__ void multMatrix(const float *a, const float *b, float *c)
{
	int i = blockDim.x*blockIdx.x + threadIdx.x;
	int j = blockDim.y*blockIdx.y + threadIdx.y;
	int k;
	float v;

	v = 0.;
	for(k = 0; k < N; k++){
		v += a[k*N + i]*b[j*N + k];
	}

	c[j*N + i] = v;
}
...

速度の比較

上で見た 4 つのプログラムについて実行速度を比較します。

GPU 計算により計算速度が CPU 計算のおよそ 2 倍になりました。予想通り、ベクトルの加算の場合の "12 倍" に比べると "劇的" とは言えないものの、そこそこの高速化を実現できています。行列の格納方法を CUBLAS 型にすると少し速度が上がりました。その理由は、multMatrix() での行列 C への値の代入におけるメモリアクセスが連続的になったからではないかと思います。CUBLAS がそれより少し遅いのは、cublasSgemm() に係数の掛け算の処理が含まれるからでしょう。

CUDA プログラミング 目次へ