N上のN個のポイント間の距離の逆数の合計を計算するカーネルを作成しようとしています。Cのシリアルコーダは次のようになります。
average = 0;
for(int i = 0; i < Np; i++){
for(int j = i + 1; j < Np; j++){
average += 1.0e0f/sqrtf((rx[i]-rx[j])*(rx[i]-rx[j]) + (ry[i]-ry[j])*(ry[i]-ry[j]));
}
}
average = average/(float)N;
ここで、rxとryはそれぞれx座標とy座標です。
乱数ジェネレーターを使用するカーネルを介してポイントを生成します。カーネルの場合、4k（8k）ポイントでブロックあたり128（256）スレッドを使用しました。その上で、すべてのスレッドが内側の上の内側のループを実行し、次のように結果がreducesum関数に渡されます。
ポイントの生成：
__global__ void InitRNG ( curandState * state, const int seed ){
int tIdx = blockIdx.x*blockDim.x + threadIdx.x;
curand_init (seed, tIdx, 0, &state[tIdx]);
}
__global__
void SortPoints(float* X, float* Y,const int N, curandState *state){
float rdmn1, rdmn2;
unsigned int tIdx = blockIdx.x*blockDim.x + threadIdx.x;
float range;
if(tIdx < N){
rdmn1 = curand_uniform(&state[tIdx]);
rdmn2 = curand_uniform(&state[tIdx]);
range = sqrtf(0.25e0f*N*rdmn1);
X[tIdx] = range*cosf(2.0e0f*pi*rdmn2);
Y[tIdx] = range*sinf(2.0e0f*pi*rdmn2);
}
}
割引：
__device__
float ReduceSum2(float In){
__shared__ float data[BlockSize];
unsigned int tIdx = threadIdx.x;
data[tIdx] = In;
__syncthreads();
for(unsigned int i = blockDim.x/2; i > 0; i >>= 1){
if(tIdx < i){
data[tIdx] += data[tIdx + i];
}
__syncthreads();
}
return data[0];
}
カーネル：
__global__
void AvgDistance(float *X, float *Y, float *Avg, const int N){
int tIdx = blockIdx.x*blockDim.x + threadIdx.x;
int bIdx = blockIdx.x;
float x , y;
float d = 0.0f;
if(tIdx < N){
for(int i = tIdx + 1; i < N ; i++){
x = X[tIdx] - X[i];
y = Y[tIdx] - Y[i];
d += 1.0e0f/(sqrtf(x*x + y*y));
}
__syncthreads();
Avg[bIdx] = ReduceSum2(d);
}
}
カーネルは次のように構成および起動されます。
dim3 threads(BlockSize,BlockSize);
dim3 blocks(ceil(Np/threads.x),ceil(Np/threads.y));
InitRNG<<<blocks.x,threads.x>>>(d_state,seed);
SortPoints<<<blocks.x,threads.x>>>(d_rx,d_ry,Np,d_state);
AvgDistance<<<blocks.x,threads.x,threads.x*sizeof(float)>>>(d_rx,d_ry,d_Avg,Np);
最後に、データをホストにコピーして戻し、残りの合計を実行します。
Avg = new float[blocks.x];
CHECK(cudaMemcpy(Avg,d_Avg,blocks.x*sizeof(float),cudaMemcpyDeviceToHost),ERROR_CPY_DEVTOH);
float average = 0;
for(int i = 0; i < blocks.x; i++){
average += Avg[i];
}
average = average/(float)Np;
4kポイントの場合、OK！結果は次のとおりです。
Average distance between points (via Kernel) = 108.615
Average distance between points (via CPU) = 110.191
この場合、合計が異なる順序で実行され、両方の結果が互いに異なる可能性があります、私にはわかりません...
しかし、8kになると、結果はまったく異なります。
Average distance between points (via Kernel) = 153.63
Average distance between points (via CPU) = 131.471
私には、カーネルとシリアルコードの両方が同じように書かれているように見えます。浮動小数点数のCUDA計算の精度に不信感を抱く原因は何ですか。これは意味がありますか？または、一部のスレッドがXとYから同時に同じデータをロードすると、グローバルメモリへのアクセスによって競合が発生しますか？または、私がカーネルを書いた方法が何らかの形で「間違っている」（つまり、両方の結果が互いに異なる原因になっていることをしているのでしょうか？）。