GPUで並列処理をやってみる

未分類

この間のインタプリタをはじめから・・・で問題がありご指摘を受けました。
ですので次回より修正を行い、今回は GPGPU を使った計算についてやりたいと思います。

今回 nVidia の GeForce 8800 GTX というボードが手に入りましたので、専用の言語 （現在 GF8x 系のみで動作可能） であるCUDAを利用したいと思います。

1.処理の流れ

  • CPUからの命令でGPUメモリのGPUのデータ領域のメモリを確保
  • CPUからGPUへメモリ内容をコピー
  • GPUで演算処理、エラーの有無をチェック
  • GPUからCPUへ出力用メモリ内容をコピー

という流れで処理を行います。
注意しなければならない点としてGPUからCPUにメモリ内容をコピーする際、GPU内部で出力用メモリに書き出しが行われなかった場合、前回の出力結果とまったく同じものが出てきます。 (再起動してもフラッシュされない場合もありました。)

2.実行中の流れ

本家サンプルなどにあるGPU処理の呼び出しでは

関数名 <<< BLOCK_NUM, THREADS_NUM >>> 関数に渡す引

となっています。 BLOCK_NUM は１つずつマルチプロセッサ (今回は１６個) に置かれて並列実行されます。
また、それぞれのブロックの中で THREADS_NUM の回数分スレッドが実行されます。

3.コーディング

とりあえず軽めに 44 行列の掛け算を300回ほど計算してもらいました。 

//main.cu
#include <stdio.h>
#include <time.h>
#include <cutil.h>
#include "kernel.cu"
#define THREADS_NUM 300
#define BLOCKS_NUM 4
void Run(int argc, char** argv)
{
// 初期化
CUT_DEVICE_INIT();
int nSize = 600;
int n;
n = sizeof(Test)nSize;
printf("size = %d\n", n);
// CPU側のデータを準備
// 乗算されるデータ
Test h_sd1;
h_sd1 = (Test)malloc(n / 2);
// 乗算するデータ
Test h_sd2;
h_sd2 = (Test)malloc(n / 2);
// 返り値の代入用
Test h_sd3;
h_sd3 = (Test)malloc(n / 2);
int i;
// 値の準備
for(i=0; i<nSize/2; i++){
for(int j = 0; j < 4; j++)
{
Obj.col[j] = make_float4(rand()%10, rand()%10, rand()%10, rand()%10);
h_sd1[i].col[j] = Obj.col[j];
Obj.row[j] = make_float4(rand()%10, rand()%10, rand()%10, rand()%10);
h_sd2[i].row[j] = Obj.row[j];
Obj.col[j] = make_float4(0, 0, 0, 0);
h_sd3[i] = Obj;
}
}
printf("\n");
// GPU側メモリの準備
// テストデータ
Test d_sd1;
// cudaMalloc(pointer, size)
CUDA_SAFE_CALL(cudaMalloc((void**)&d_sd1, n / 2));
// cudaMemcpy(dist, from, size, CPU->GPU)
CUDA_SAFE_CALL(cudaMemcpy(d_sd1, h_sd1, n / 2, cudaMemcpyHostToDevice));
Test d_sd2;
// cudaMalloc(pointer, size)
CUDA_SAFE_CALL(cudaMalloc((void)&d_sd2, n / 2));
// cudaMemcpy(dist, from, size, CPU->GPU)
CUDA_SAFE_CALL(cudaMemcpy(d_sd2, h_sd2, n / 2, cudaMemcpyHostToDevice));
// 値を得るためのデータ代入用
Test* d_sd3;
CUDA_SAFE_CALL(cudaMalloc((void)&d_sd3, n/2));
// スレッド数
dim3 threads(THREADS_NUM, 1, 1);
// ブロック数
dim3 grid(BLOCKS_NUM, BLOCKS_NUM, 1);
// タイマー処理、msで出ます。
unsigned int hTimer;
CUT_SAFE_CALL( cutCreateTimer(&hTimer) );
CUT_SAFE_CALL( cutResetTimer(hTimer) );
CUT_SAFE_CALL( cutStartTimer(hTimer) );
// GPUにジョブを投げます
test<<< grid, threads >>>(d_sd3, d_sd1, d_sd2, 0, nSize / 2);
// 終了待ち
CUDA_SAFE_CALL( cudaThreadSynchronize() );
CUT_SAFE_CALL( cutStopTimer(hTimer) );
double gpuTime = cutGetTimerValue(hTimer);
// GPUの処理に問題が起きていないかの確認
// メモリ領域をオーバーするなどkernel.cuが実行できなかった場合エラーになる
CUT_CHECK_ERROR("Kernel execution failed");
// 演算結果の取得
CUDA_SAFE_CALL(cudaMemcpy(h_sd3, d_sd3, n / 2, cudaMemcpyDeviceToHost));
FILE fp;
if((fp = fopen("log.txt", "w")) == NULL)
{
exit(-1);
}
for(i=0; i<nSize / 2; i++)
{
for(int j = 0; j < 4; j++){
for(int k = 0; k < 4; k++)
{
fprintf(fp, "%d, %d, %d\n", i, j, k);
fprintf(fp, "x1 = %f, y1 = %f, z1 = %f, w = %f\n",
h_sd1[i].col[j].x, h_sd1[i].col[j].y, h_sd1[i].col[j].z, h_sd1[i].col[j].w);
fprintf(fp, "x2 = %f, y2 = %f, z2 = %f, w = %f\n",
h_sd2[i].row[k].x, h_sd2[i].row[k].y, h_sd2[i].row[k].z, h_sd2[i].row[k].w);
switch(k){
case 0:
fprintf(fp, "x = %f\n",h_sd3[i].col[j].x);
break;
case 1:
fprintf(fp, "y = %f\n", h_sd3[i].col[j].y);
break;
case 2:
fprintf(fp, "z = %f\n",h_sd3[i].col[j].z);
break;
case 3:
fprintf(fp, "w = %f\n",h_sd3[i].col[j].w);
break;
}
} // k
} // j
fprintf(fp, "threadIdx.x = %d, blockIdx.x = %d\n", h_sd3[i].threadx,
h_sd3[i].blockx);
} // i
printf("Time: %f ms\n", gpuTime);
fclose(fp);
// クリーンアップ
free(h_sd1);
free(h_sd2);
free(h_sd3);
CUDA_SAFE_CALL(cudaFree(d_sd1));
CUDA_SAFE_CALL(cudaFree(d_sd2));
CUDA_SAFE_CALL(cudaFree(d_sd3));
}
int main(int argc, char** argv)
{
Run(argc, argv);
// 終了処理
CUT_EXIT(argc, argv);
}

次にGPU内部処理(kernel.cu)です。 

#include <math.h>
struct Test
{
float4 col[4], row[4];
int threadx, thready, threadz;
int blockx, blocky, blockz;
int nSize;  // 今回使用しません。
int grid;   // 今回使用しません。
};
static struct Test Obj;
// GPU側の処理
/
global・・・CPU側からでもGPU側からでも呼び出せる関数（必ずvoid型）
device・・・GPU側からGPUのみ呼び出せる関数（GPUのグローバルメモリに置かれる）
host・・・CPU側からCPUのみ呼び出せる関数（特に指定しない場合全てhostになる）
test関数（出力メモリ, 入力データ1（乗算される行列）, 入力データ2（乗算するデータ）,
スレッド分割これ以上できなくなった場合などに使用する開始点変数, 終了点変数)
/

global void
test(Test fOut, Test fIn1, TestfIn2, int nStart, int nEnd)
{
/ 並列動作の順序
while(スレッドX < THREADNUMS) {
while(ブロックX < BLOCKNUMS) {
while(ブロックY < BLOCKNUMS) {
// 計算処理
ブロックY++;
}
ブロックX++;
　　　　　ブロックY = 0;
}
スレッドX++;
ブロックX = 0;
}
/
// スレッドのインデックス（0-3）
const unsigned int tid = threadIdx.x;
// グリッドXのインデックス（0-299）
const unsigned int bidx = blockIdx.x;
// グリッドYのインデックス（0-299）
const unsigned int bidy = blockIdx.y;
int i = bidx;
int j = bidy;
// 行列積を1要素ずつ計算します、tid 行 i 列の計算です
switch(j) {
case 0:
fOut[tid].col[i].x =
fIn1[tid].col[i].x * fIn2[tid].row[j].x +
fIn1[tid].col[i].y * fIn2[tid].row[j].y +
fIn1[tid].col[i].z * fIn2[tid].row[j].z +
fIn1[tid].col[i].w * fIn2[tid].row[j].w;
break;
case 1:
fOut[tid].col[i].y =
fIn1[tid].col[i].x * fIn2[tid].row[j].x +
fIn1[tid].col[i].y * fIn2[tid].row[j].y +
fIn1[tid].col[i].z * fIn2[tid].row[j].z +
fIn1[tid].col[i].w * fIn2[tid].row[j].w;
break;
case 2:
fOut[tid].col[i].z =
fIn1[tid].col[i].x * fIn2[tid].row[j].x +
fIn1[tid].col[i].y * fIn2[tid].row[j].y +
fIn1[tid].col[i].z * fIn2[tid].row[j].z +
fIn1[tid].col[i].w * fIn2[tid].row[j].w;
break;
case 3:
fOut[tid].col[i].w =
fIn1[tid].col[i].x * fIn2[tid].row[j].x +
fIn1[tid].col[i].y * fIn2[tid].row[j].y +
fIn1[tid].col[i].z * fIn2[tid].row[j].z +
fIn1[tid].col[i].w * fIn2[tid].row[j].w;
break;
}
/ 初期化用です。
fOut[tid].col[i].x = 0;
fOut[tid].col[i].y = 0;
fOut[tid].col[i].z = 0;
fOut[tid].col[i].w = 0;
/
// 最終的に何行何列にいるのか、を出力メモリに書き出します 
fOut[tid].threadx = tid;
fOut[tid].thready = threadIdx.y;
fOut[tid].threadz = threadIdx.z;
fOut[tid].blockx = bidx;
fOut[tid].blocky = bidy;
fOut[tid].blockz = blockIdx.z;
}

　結果の一部を抜き出しチェックしてみました。
x1, y1, z1, wが行列ABのAのカラム、x2, y2, z2, wがBのロウです。
　

x1 = 3.000000, y1 = 3.000000, z1 = 6.000000, w = 1.000000
x2 = 5.000000, y2 = 3.000000, z2 = 2.000000, w = 6.000000
w = 42.000000
299, 3, 0 <- スレッドID, ブロックx, ブロックy
x1 = 3.000000, y1 = 6.000000, z1 = 7.000000, w = 7.000000
x2 = 4.000000, y2 = 5.000000, z2 = 8.000000, w = 4.000000
x = 126.000000
299, 3, 1
x1 = 3.000000, y1 = 6.000000, z1 = 7.000000, w = 7.000000
x2 = 1.000000, y2 = 0.000000, z2 = 9.000000, w = 6.000000
y = 108.000000
299, 3, 2
x1 = 3.000000, y1 = 6.000000, z1 = 7.000000, w = 7.000000
x2 = 8.000000, y2 = 4.000000, z2 = 8.000000, w = 2.000000
z = 118.000000
299, 3, 3
x1 = 3.000000, y1 = 6.000000, z1 = 7.000000, w = 7.000000
x2 = 5.000000, y2 = 3.000000, z2 = 2.000000, w = 6.000000
w = 89.000000
threadIdx.x = 299, blockIdx.x = 0

4.パフォーマンスと実用性

上記程度の計算では実行時間は0.074008ms（メモリ確保の時間は考慮していません。）ほどでCPU(Core2 DUO E6600）実行時間と差がでません。
そこで行列の要素数を1616, 3232, 6464..., 512512としてみたところ256256行列では約5倍程度、512*512行列では約30倍ほど性能に開きがでました。
このことより、並列処理が有効な部分に用いればかなりの性能の上昇が期待できます。ただしメモリ転送時間や、一度に渡すことが出来るデータ量の制限がありますので、計算よりもデータ転送のほうが多い処理はパフォーマンスは落ちてしまいます。 担当：松浦