CUDA C：我需要修正我的串行程序的偏差

Question

所以我有一个非常基本的二维元胞自动机程序。棘手的部分，或者至少我转向学习 GPU 编程的原因是，我的矩阵往往至少有 1 亿个元素。我已经将我的串行程序的我可能称之为天真的实现（第一稿）变成了 CUDA C。它工作得很好，但我认为我看到了一种进一步改进它的方法。

出于测试目的，我有必要在矩阵的每次迭代之间跟踪活细胞计数（在我的实现中，活细胞用 1 表示，死细胞用 0 表示）。在我的第一稿中，主机负责在设备完成更新矩阵后对活细胞计数的 N^2 个嵌套循环。

这里有一些伪代码可以让您更好地了解（以下代码本身位于循环中，在任意次数的迭代后终止）：

cuda_kern<<<dimGrid, dimBlock>>>( d_grid );
cudaMemcpy( h_grid, d_grid, matrixSize*sizeof(int), cudaMemcpyDeviceToHost );

*liveCell = 0;
for(a=0; a<gridLength; a++)
    for(b=0; b<gridLength; b++)
        if(h_grid[a][b]==1) (*liveCell)++;

printf("live cell count is %d\n", *liveCell);

如您所见，我手头有两个相当昂贵的操作。首先，我必须在每次迭代时在设备和主机之间传输我的矩阵（上面称为网格）。其次，我有 N^2 操作来组合/计算活细胞计数矩阵。我相信我可以通过在每次迭代中仅在主机和设备之间传输活细胞计数来消除这两个步骤。事实上，我已经非常接近编码了，但正如问题标题所示，与我的串行程序和我最新的 CUDA 程序确认的实时计数存在轻微的“随机”偏差。

这是我的第二个内核草案中的一些伪代码：

__global__ void cuda_kern( liveCell, grid )
{
    i = threadIdx.x + blockIdx.x * blockDim.x
    j = threadIdx.y + blockIdx.y * blockDim.y

    if( i < gridLength AND j < gridLength )
    {
        __shared__ int temp[number of threads in block] 
        x, y, count, changeCounter, sum = 0

        // nested for loops here that update int variable 'count'
        // for those familiar with CA it is a basic neighborhood analysis
        for x to arbitrary neighborhood range
            for y to arbitrary neighborhood range
            {
                //count neighbors, that is update 'count' variable
            }

        // __syncthreads()

        if( grid[i*gridLength+j] == 0 AND count == 5 ) 
        { // NOTE: 'count' above could be any arbitrary integer 1 - 8  
            grid[i*gridLength+j] = 1
            changeCounter += 1
        }
        else if( grid[i*gridLength+j] ==  1 AND count >= 5 )
        {
            grid[i*gridLength+j] = 0
        }
        else 
        {
            grid[i*gridLength+j] = 1
            changeCounter += 1
        }

        //__syncthreads()

        temp[threadIdx.x] = changeCounter 

        //__syncthreads()

        if (threadIdx.x == 0)
            for( i = 0; i < N; i++ ) 
            { 
                sum += temp[i]
                //__syncthreads()
            }

        atomicAdd( liveCell, sum )

    } // end of if ( i < gridLength and j < gridLength )   
    return;
} // end of kernel

上述伪代码说明： 我试图忽略不重要的部分，例如元胞自动机函数背后的大部分逻辑。在我的程序中，所有变量都是整数，而矩阵本身是整数类型。我已将尝试同步线程的地方留在 cmets 中（无济于事）。

我的逻辑解释： 如前所述，在讨论我的第一个实现草案时，我想避免在设备和主机之间传输矩阵的昂贵操作以及每次迭代计算主机上活细胞的 N^2 成本。使用变量“changeCounter”可以看到这样做的关键。本质上，每当一个单元格被指定为“1”或活动时，“changeCounter”应该加一。我正在尝试使用由当前线程索引的 shared 变量“temp”来存储“changeCounter”的值。一旦块中的所有线程都完成后，我会尝试将“temp”数组压缩为单个变量“sum”，然后通过原子操作“atomicAdd”将其添加到 liveCell。

我遇到的问题是我使用这种方法的结果（虽然更快）与我原来的草稿结果不匹配，而且每次我在同一个矩阵上运行程序时，每次迭代都会有轻微的偏差活细胞计数。

更新 1k x 1k 矩阵上的示例输出：

GPU 实施的 1st Draft 每次都会产生以下活细胞计数。这些是我希望在我的第二稿中产生的结果。

initial live cell count: 393592 
itr. 0  live cell count: 364118
itr. 1  live cell count: 315417
itr. 2  live cell count: 300413 
itr. 3  live cell count: 284503 

GPU 实施方案 2 次草案的活细胞计数每次都略有不同，但接近上述值。以下是三个单独的运行，可让您了解一下。

运行 A：

initial live cell count: 393592
itr. 0  live cell count: 372402
itr. 1  live cell count: 324114
itr. 2  live cell count: 309580
itr. 3  live cell count: 291393 

运行 B：

initial live cell count: 393592
itr. 0  live cell count: 374139
itr. 1  live cell count: 323948
itr. 2  live cell count: 307214
itr. 3  live cell count: 292582 

运行 C：

initial live cell count: 393592
itr. 0  live cell count: 372391
itr. 1  live cell count: 323105
itr. 2  live cell count: 308295
itr. 3  live cell count: 292512 

每次都会在主机上计算初始活细胞计数，这解释了为什么它在所有示例中都是一致的。

更新结束

感谢您的宝贵时间。如果您需要更多信息，请告诉我。

Answer 1

你不应该这样做吗：

    if (threadIdx.x == 0) {
        for( i = 0; i < N; i++ ) 
        { 
            sum += temp[i]
            //__syncthreads()
        }
        atomicAdd( liveCell, sum )// could probably just update the memory without atomic
    }

正如你所说，块中的所有线程都计算它们的总和，然后你用线程 0 总结所有线程的总和，然后你应该将它添加到指定的指针中。如果每个线程都执行 atomicAdd 会有所不同，您可能应该只对线程 0 执行此操作，并且您可能可以删除 atomic 操作。

另外，在那个循环上做归约：

for (i=blockDim.x/2; i>0; i>>=1) {
     if(threadIdx.x < i) {
         temp[threadIdx.x] += temp[threadIdx.x + i];
     }
    __syncthreads();
}
if(threadIdx.x == 0) {
    sum = temp[0];
    liveCell += sum;
}