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.コーディング
とりあえず軽めに 4*4 行列の掛け算を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, Test*fIn2, 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が行列A*Bの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）実行時間と差がでません。
そこで行列の要素数を16*16, 32*32, 64*64..., 512*512としてみたところ256*256行列では約5倍程度、512*512行列では約30倍ほど性能に開きがでました。
このことより、並列処理が有効な部分に用いればかなりの性能の上昇が期待できます。ただしメモリ転送時間や、一度に渡すことが出来るデータ量の制限がありますので、計算よりもデータ転送のほうが多い処理はパフォーマンスは落ちてしまいます。
担当：松浦
>2 の「実行中の流れ」の文章が全体的に意味不明です
かなり混乱してしまう書き方でしたので以下のように修正しました。
（修正前）CPUでCPUのデータ領域とGPUのデータ領域のメモリを確保
（修正後）CPUからの命令でGPUメモリのGPUのデータ領域のメモリを確保
>エラーがあるかないかはどうやってチェックするんですか？ そもそも何のエラー？
>CPU 側の処理で、GPU の処理が終わるのをどこで待ってるんですか？
ソースコードのコメントに追加しました。
>関数に適切に分割されてないのも気になります。
小規模であったので、、、かなりひとつの関数が長くなってしまいました。
コメント (2)
2 の「実行中の流れ」の文章が全体的に意味不明です。
CPU から GPU へメモリ内容をコピーするということなので、GPU にメモリがあると思われるのですが、本当に CPU で GPU のデータ領域のメモリを確保することができるのですか？ 確保するのは GPU なのでは？
さらに、エラーがあるかないかはどうやってチェックするんですか？ そもそも何のエラー？
ソースコードについてもっと詳しく書いてください。 関数に適切に分割されてないのも気になります。 高速化のためなのでしょうか？ CPU 側の処理で、GPU の処理が終わるのをどこで待ってるんですか？
投稿者: 齋藤 | 2007年12月24日 11:51