このような質問では、何も追加したり変更したりすることなく、誰かがコンパイルして実行できる完全なコードを提供することをお勧めします。一般的に言えば、SOはこれを期待しています。質問はパフォーマンスに関するものでもあるため、実際のタイミング測定方法も完全なコードに含める必要があります。
エラーの修正:
コードには少なくとも 2 つのエラーがあり、そのうちの 1 つは @Jez がすでに指摘しています。この「部分削減」ステップの後:
if ( blockSize >= 128 ) {
if ( threadIdx.x < 64 ) {
sum[threadIdx.x] += sum[threadIdx.x + 64];
}
}
__syncthreads();残りの作業に進む前に が必要です。上記の変更により、カーネルが私の単純なホスト実装に一致する再現可能な結果を生成することができました。また、スレッドブロック全体で同じように評価されない次のような条件付きコードがあるため:
if ( threadIdx.x < 32 ) {
__syncthreads()条件付きコード ブロック内にステートメントを含めることは違法です。
if ( blockSize >= 64 ) { sum[threadIdx.x] += sum[threadIdx.x + 32]; } __syncthreads();
（そして、同じことを行う後続の行についても同様です）。そのため、修正することをお勧めします。これを解決する方法はいくつかありますが、そのうちの 1 つは、volatile型付きポインターを使用して共有データを参照するように切り替えることです。現在ワープ内で操作しているため、volatile修飾子はコンパイラに必要なことを強制します。
volatile UINT *vsum = sum;
if ( threadIdx.x < 32 ) {
if ( blockSize >= 64 ) vsum[threadIdx.x] += vsum[threadIdx.x + 32];
if ( blockSize >= 32 ) vsum[threadIdx.x] += vsum[threadIdx.x + 16];
if ( blockSize >= 16 ) vsum[threadIdx.x] += vsum[threadIdx.x + 8 ];
if ( blockSize >= 8 ) vsum[threadIdx.x] += vsum[threadIdx.x + 4 ];
if ( blockSize >= 4 ) vsum[threadIdx.x] += vsum[threadIdx.x + 2 ];
if ( blockSize >= 2 ) vsum[threadIdx.x] += vsum[threadIdx.x + 1 ];
}
CUDA並列リダクションのサンプル コードと関連する pdfは、良いレビューになるかもしれません。
タイミング/パフォーマンス分析:
たまたま GT 640、cc3.5 デバイスを持っています。それを実行bandwidthTestすると、デバイス間の転送で約 32GB/秒が観察されます。この数値は、デバイス カーネルがデバイス メモリにアクセスしているときに達成可能な帯域幅の合理的なおおよその上限を表しています。また、ベースのタイミングを追加cudaEventして、シミュレートされたデータを使用して、あなたが示したサンプル コードを作成すると、5GB/s ではなく、約 16GB/s のスループットが観察されます。したがって、実際の測定手法はここで役立つ情報になります (実際、カーネルのタイミングとタイミングの違いを分析するには、おそらく完全なコードが必要です)。
問題は残りますが、それを改善できますか？(約 32GB/秒が上限であると仮定します)。
あなたの質問：
バイトにアクセスする方法でバンク競合が発生する可能性はありますか? もしそうなら、私は衝突を避けることができますか?
カーネルは実際にはバイトを 32 ビット量 ( uchar4) として効果的にロードし、各スレッドは隣接する連続した 32 ビット量をロードしているため、カーネルにバンク競合アクセスの問題はないと思います。
私の reinterpret_cast の使い方は正しいですか?
はい、正しいようです (上記の修正を加えた以下のサンプル コードは、カーネルによって生成された結果が単純なホスト関数の実装と一致することを検証します)。
8ビットの符号なし計算を行うためのより良い方法はありますか?
この場合、@njuffa が指摘したように、SIMD 組み込み関数は単一の命令でこれを処理できます (__vsadu4()以下のサンプル コードを参照)。
他にどのような最適化を行うことができますか (私は完全な初心者なので、多くのことを想定しています)。
@MichalHosala によって提案された cc3.0 ワープシャッフル リダクション メソッドを使用します。
__vsadu4()@njuffaによって提案されたバイト数の処理を簡素化および改善するために、SIMD 組み込み関数を利用します。
データベースのベクトル データを再編成して、列優先のストレージに格納します。これにより、通常の並列リダクション方法 (項目 1 で述べた方法であっても) を省き、1 つのスレッドがベクトル比較全体を計算する単純な for ループ読み取りカーネルに切り替えることができます。これにより、カーネルはこの場合 (cc3.5 GT640) のデバイスのメモリ帯域幅にほぼ達することができます。
以下は、3 つの実装を示すコードと実行例です。元の実装 (さらに、正しい結果を生成するための上記の「修正」)、上記のリストの項目 1 と 2 を含めるように変更した opt1 カーネル、上記のリストの 2 と 3 を利用するアプローチであなたのものを置き換える opt2 カーネル。私の測定によると、あなたのカーネルは GT640 の帯域幅の約半分である 16GB/s を達成し、opt1 カーネルは約 24GB/s で実行され (増加は上記の項目 1 と 2 からほぼ同じ割合で発生します)、opt2 カーネルは、データ再編成を使用して、ほぼ全帯域幅 (36GB/秒) で実行されます。
$ cat t574.cu
#include <stdio.h>
#include <stdlib.h>
#define THREADS_PER_BLOCK 128
#define VECTOR_SIZE 1024
#define NUM_DB_VEC 100000
typedef unsigned char BYTE;
typedef unsigned int UINT;
typedef unsigned int uint32_t;
template<UINT blockSize> __global__ void reduction_sum_abs( BYTE* query_vector, BYTE* db_vector, uint32_t* result )
{
extern __shared__ UINT sum[];
UINT db_linear_index = (blockIdx.y*gridDim.x) + blockIdx.x ;
UINT i = threadIdx.x;
sum[threadIdx.x] = 0;
int* p_q_int = reinterpret_cast<int*>(query_vector);
int* p_db_int = reinterpret_cast<int*>(db_vector);
while( i < VECTOR_SIZE/4 ) {
/* memory transaction */
int q_int = p_q_int[i];
int db_int = p_db_int[db_linear_index*VECTOR_SIZE/4 + i];
uchar4 a0 = *reinterpret_cast<uchar4*>(&q_int);
uchar4 b0 = *reinterpret_cast<uchar4*>(&db_int);
/* sum of absolute difference */
sum[threadIdx.x] += abs( (int)a0.x - b0.x );
sum[threadIdx.x] += abs( (int)a0.y - b0.y );
sum[threadIdx.x] += abs( (int)a0.z - b0.z );
sum[threadIdx.x] += abs( (int)a0.w - b0.w );
i += THREADS_PER_BLOCK;
}
__syncthreads();
if ( blockSize >= 128 ) {
if ( threadIdx.x < 64 ) {
sum[threadIdx.x] += sum[threadIdx.x + 64];
}
}
__syncthreads(); // **
/* reduce the final warp */
if ( threadIdx.x < 32 ) {
if ( blockSize >= 64 ) { sum[threadIdx.x] += sum[threadIdx.x + 32]; } __syncthreads();
if ( blockSize >= 32 ) { sum[threadIdx.x] += sum[threadIdx.x + 16]; } __syncthreads();
if ( blockSize >= 16 ) { sum[threadIdx.x] += sum[threadIdx.x + 8 ]; } __syncthreads();
if ( blockSize >= 8 ) { sum[threadIdx.x] += sum[threadIdx.x + 4 ]; } __syncthreads();
if ( blockSize >= 4 ) { sum[threadIdx.x] += sum[threadIdx.x + 2 ]; } __syncthreads();
if ( blockSize >= 2 ) { sum[threadIdx.x] += sum[threadIdx.x + 1 ]; } __syncthreads();
}
/* copy the sum back to global */
if ( threadIdx.x == 0 ) {
result[db_linear_index] = sum[0];
}
}
__global__ void reduction_sum_abs_opt1( BYTE* query_vector, BYTE* db_vector, uint32_t* result )
{
__shared__ UINT sum[THREADS_PER_BLOCK];
UINT db_linear_index = (blockIdx.y*gridDim.x) + blockIdx.x ;
UINT i = threadIdx.x;
sum[threadIdx.x] = 0;
UINT* p_q_int = reinterpret_cast<UINT*>(query_vector);
UINT* p_db_int = reinterpret_cast<UINT*>(db_vector);
while( i < VECTOR_SIZE/4 ) {
/* memory transaction */
UINT q_int = p_q_int[i];
UINT db_int = p_db_int[db_linear_index*VECTOR_SIZE/4 + i];
sum[threadIdx.x] += __vsadu4(q_int, db_int);
i += THREADS_PER_BLOCK;
}
__syncthreads();
// this reduction assumes THREADS_PER_BLOCK = 128
if (threadIdx.x < 64) sum[threadIdx.x] += sum[threadIdx.x+64];
__syncthreads();
if ( threadIdx.x < 32 ) {
unsigned localSum = sum[threadIdx.x] + sum[threadIdx.x + 32];
for (int i = 16; i >= 1; i /= 2)
localSum = localSum + __shfl_xor(localSum, i);
if (threadIdx.x == 0) result[db_linear_index] = localSum;
}
}
__global__ void reduction_sum_abs_opt2( BYTE* query_vector, UINT* db_vector_cm, uint32_t* result)
{
__shared__ UINT qv[VECTOR_SIZE/4];
if (threadIdx.x < VECTOR_SIZE/4) qv[threadIdx.x] = *(reinterpret_cast<UINT *>(query_vector) + threadIdx.x);
__syncthreads();
int idx = threadIdx.x + blockDim.x*blockIdx.x;
while (idx < NUM_DB_VEC){
UINT sum = 0;
for (int i = 0; i < VECTOR_SIZE/4; i++)
sum += __vsadu4(qv[i], db_vector_cm[(i*NUM_DB_VEC)+idx]);
result[idx] = sum;
idx += gridDim.x*blockDim.x;}
}
unsigned long compute_host_result(BYTE *qvec, BYTE *db_vec){
unsigned long temp = 0;
for (int i =0; i < NUM_DB_VEC; i++)
for (int j = 0; j < VECTOR_SIZE; j++)
temp += (unsigned long) abs((int)qvec[j] - (int)db_vec[(i*VECTOR_SIZE)+j]);
return temp;
}
int main(){
float et;
cudaEvent_t start, stop;
BYTE *h_qvec, *d_qvec, *h_db_vec, *d_db_vec;
uint32_t *h_res, *d_res;
h_qvec = (BYTE *)malloc(VECTOR_SIZE*sizeof(BYTE));
h_db_vec = (BYTE *)malloc(VECTOR_SIZE*NUM_DB_VEC*sizeof(BYTE));
h_res = (uint32_t *)malloc(NUM_DB_VEC*sizeof(uint32_t));
for (int i = 0; i < VECTOR_SIZE; i++){
h_qvec[i] = rand()%256;
for (int j = 0; j < NUM_DB_VEC; j++) h_db_vec[(j*VECTOR_SIZE)+i] = rand()%256;}
cudaMalloc(&d_qvec, VECTOR_SIZE*sizeof(BYTE));
cudaMalloc(&d_db_vec, VECTOR_SIZE*NUM_DB_VEC*sizeof(BYTE));
cudaMalloc(&d_res, NUM_DB_VEC*sizeof(uint32_t));
cudaMemcpy(d_qvec, h_qvec, VECTOR_SIZE*sizeof(BYTE), cudaMemcpyHostToDevice);
cudaMemcpy(d_db_vec, h_db_vec, VECTOR_SIZE*NUM_DB_VEC*sizeof(BYTE), cudaMemcpyHostToDevice);
cudaEventCreate(&start); cudaEventCreate(&stop);
// initial run
cudaMemset(d_res, 0, NUM_DB_VEC*sizeof(uint32_t));
cudaEventRecord(start);
reduction_sum_abs<THREADS_PER_BLOCK><<<NUM_DB_VEC, THREADS_PER_BLOCK, THREADS_PER_BLOCK*sizeof(int)>>>(d_qvec, d_db_vec, d_res);
cudaEventRecord(stop);
cudaDeviceSynchronize();
cudaEventSynchronize(stop);
cudaEventElapsedTime(&et, start, stop);
cudaMemcpy(h_res, d_res, NUM_DB_VEC*sizeof(uint32_t), cudaMemcpyDeviceToHost);
unsigned long h_result = 0;
for (int i = 0; i < NUM_DB_VEC; i++) h_result += h_res[i];
printf("1: et: %.2fms, bw: %.2fGB/s\n", et, (NUM_DB_VEC*VECTOR_SIZE)/(et*1000000));
if (h_result == compute_host_result(h_qvec, h_db_vec)) printf("Success!\n");
else printf("1: mismatch!\n");
// optimized kernel 1
cudaMemset(d_res, 0, NUM_DB_VEC*sizeof(uint32_t));
cudaEventRecord(start);
reduction_sum_abs_opt1<<<NUM_DB_VEC, THREADS_PER_BLOCK>>>(d_qvec, d_db_vec, d_res);
cudaEventRecord(stop);
cudaDeviceSynchronize();
cudaEventSynchronize(stop);
cudaEventElapsedTime(&et, start, stop);
cudaMemcpy(h_res, d_res, NUM_DB_VEC*sizeof(uint32_t), cudaMemcpyDeviceToHost);
h_result = 0;
for (int i = 0; i < NUM_DB_VEC; i++) h_result += h_res[i];
printf("2: et: %.2fms, bw: %.2fGB/s\n", et, (NUM_DB_VEC*VECTOR_SIZE)/(et*1000000));
if(h_result == compute_host_result(h_qvec, h_db_vec)) printf("Success!\n");
else printf("2: mismatch!\n");
// convert db_vec to column-major storage for optimized kernel 2
UINT *h_db_vec_cm, *d_db_vec_cm;
h_db_vec_cm = (UINT *)malloc(NUM_DB_VEC*(VECTOR_SIZE/4)*sizeof(UINT));
cudaMalloc(&d_db_vec_cm, NUM_DB_VEC*(VECTOR_SIZE/4)*sizeof(UINT));
for (int i = 0; i < NUM_DB_VEC; i++)
for (int j = 0; j < VECTOR_SIZE/4; j++)
h_db_vec_cm[(j*NUM_DB_VEC)+i] = *(reinterpret_cast<UINT *>(h_db_vec + (i*VECTOR_SIZE))+j);
cudaMemcpy(d_db_vec_cm, h_db_vec_cm, NUM_DB_VEC*(VECTOR_SIZE/4)*sizeof(UINT), cudaMemcpyHostToDevice);
cudaMemset(d_res, 0, NUM_DB_VEC*sizeof(uint32_t));
cudaEventRecord(start);
reduction_sum_abs_opt2<<<64, 512>>>(d_qvec, d_db_vec_cm, d_res);
cudaEventRecord(stop);
cudaDeviceSynchronize();
cudaEventSynchronize(stop);
cudaEventElapsedTime(&et, start, stop);
cudaMemcpy(h_res, d_res, NUM_DB_VEC*sizeof(uint32_t), cudaMemcpyDeviceToHost);
h_result = 0;
for (int i = 0; i < NUM_DB_VEC; i++) h_result += h_res[i];
printf("3: et: %.2fms, bw: %.2fGB/s\n", et, (NUM_DB_VEC*VECTOR_SIZE)/(et*1000000));
if(h_result == compute_host_result(h_qvec, h_db_vec)) printf("Success!\n");
else printf("3: mismatch!\n");
return 0;
}
$ nvcc -O3 -arch=sm_35 -o t574 t574.cu
$ ./run35 t574
1: et: 6.34ms, bw: 16.14GB/s
Success!
2: et: 4.16ms, bw: 24.61GB/s
Success!
3: et: 2.83ms, bw: 36.19GB/s
Success!
$
いくつかのメモ:
- 上記のコード、特にカーネルは、テスト ケースをセットアップした方法で、cc3.0 以降用にコンパイルする必要があります。これは、1 つの 1D グリッドで 100,000 個のブロックを作成しているためです。たとえば、cc2.0 デバイスでそのまま実行することはできません。
- 特に異なるデバイスで実行する場合、opt2 カーネルに対して、グリッドとブロックのパラメーターを変更することで、若干の追加調整が可能な場合があります。私はこれらを 64 と 512 に設定していますが、これらの値は重要ではありません (ただし、ブロックは VECTOR_SIZE/4 スレッド以上である必要があります)。これは、アルゴリズムがグリッド ストライディング ループを使用してベクトル セット全体をカバーするためです。GT640 には SM が 2 つしかないため、この場合、デバイスをビジー状態に保つには 64 個のスレッドブロックで十分です (おそらく 32 個でも問題ありません)。これらを変更して、より大きなデバイスで最大のパフォーマンスを得ることができます。