使用 managedCuda 汇总数组中的元素

Summing up elements in array using managedCuda

问题描述

我试图让一个汇总数组所有元素的内核工作。内核旨在以每块 256 个线程和任意数量的块启动。传入的数组长度为a,永远是512的倍数,实际上是#blocks * 512。内核的一个块应该总结'its' 512个元素(256个线程可以总结512 个元素使用此算法),将结果存储在 out[blockIdx.x] 中。 out 中值的最终总和,以及块的结果,将在主机上完成。
该内核最多适用于 6 个块,即最多 3072 个元素。但是用超过 6 个块启动它会导致第一个块计算出比其他块(即 out = {572, 512, 512, 512, 512, 512, 512})严格更大的错误结果,这个错误结果是可重现的,多次执行的错误值是相同的。
我想这意味着我的代码中某处存在结构错误,这与 blockIdx.x 有关,但唯一的用途是计算 blockStart,并且这个接缝是正确的计算,也对于第一个块。
我验证了我的主机代码是否为内核计算了正确的块数并传入了一个正确大小的数组。那不是问题。
当然,我在 Whosebug 上阅读了很多类似的问题,但是 none 似乎描述了我的问题(参见 here or here
内核是通过 managedCuda (C#) 调用的,我不知道这是否是一个问题。

硬件

我用的MX150规格如下:

代码

内核

__global__ void Vector_Reduce_As_Sum_Kernel(float* out, float* a)
{   
int tid = threadIdx.x;
int blockStart = blockDim.x * blockIdx.x * 2;
int i = tid + blockStart;

int leftSumElementIdx =  blockStart + tid * 2;

a[i] = a[leftSumElementIdx] + a[leftSumElementIdx + 1];

__syncthreads();

if (tid < 128) 
{
    a[i] = a[leftSumElementIdx] + a[leftSumElementIdx + 1];
}

__syncthreads();

if(tid < 64)
{
    a[i] = a[leftSumElementIdx] + a[leftSumElementIdx + 1];
}

__syncthreads();

if (tid < 32)
{
    a[i] = a[leftSumElementIdx] + a[leftSumElementIdx + 1];
}

__syncthreads();

if (tid < 16)
{
    a[i] = a[leftSumElementIdx] + a[leftSumElementIdx + 1];
}

__syncthreads();

if (tid < 8)
{
    a[i] = a[leftSumElementIdx] + a[leftSumElementIdx + 1];
}

__syncthreads();

if (tid < 4)
{
    a[i] = a[leftSumElementIdx] + a[leftSumElementIdx + 1];
}

__syncthreads();

if (tid < 2)
{
    a[i] = a[leftSumElementIdx] + a[leftSumElementIdx + 1];
}

__syncthreads();

if (tid == 0)
{
    out[blockIdx.x] = a[blockStart] + a[blockStart + 1];
}
}

内核调用

//Get the cuda kernel
//PathToPtx and MangledKernelName must be replaced
CudaContext cntxt = new CudaContext();
CUmodule module = cntxt.LoadModule("pathToPtx");    
CudaKernel vectorReduceAsSumKernel = new CudaKernel("MangledKernelName", module, cntxt);

//Get an array to reduce
float[] array = new float[4096];
for(int i = 0; i < array.Length; i++)
{
    array[i] = 1;
}

//Calculate execution info for the kernel
int threadsPerBlock = 256;
int numOfBlocks = array.Length / (threadsPerBlock * 2);

//Memory on the device
CudaDeviceVariable<float> m_d = array;
CudaDeviceVariable<float> out_d = new CudaDeviceVariable<float>(numOfBlocks);

//Give the kernel necessary execution info
vectorReduceAsSumKernel.BlockDimensions = threadsPerBlock;
vectorReduceAsSumKernel.GridDimensions = numOfBlocks;

//Run the kernel on the device
vectorReduceAsSumKernel.Run(out_d.DevicePointer, m_d.DevicePointer);

//Fetch the result
float[] out_h = out_d;

//Sum up the partial sums on the cpu
float sum = 0;
for(int i = 0; i < out_h.Length; i++)
{
    sum += out_h[i];
}

//Verify the correctness
if(sum != 4096)
{
    throw new Exception("Thats the wrong result!");
}

更新:

非常有用且唯一的答案确实解决了我所有的问题。谢谢!问题是不可预见的竞争条件。

重要提示:

managedCuda 的作者在评论中指出所有 NPPs 方法确实已经在 managedCuda 中实现 (using ManagedCuda.NPP.NPPsExtensions;)。我不知道这一点,我猜很多人都在读这个问题。

您没有正确地将每个块将处理总数组中的 512 个元素的想法合并到您的代码中。根据我的测试,您至少需要进行 2 处更改才能解决此问题:

  1. 在内核中,你错误地计算了每个块的起点:

    int blockStart = blockDim.x * blockIdx.x;
    

    因为 blockDim.x 是 256,但是每个块处理 512 个元素,所以你必须将它乘以 2。(你计算 leftSumElementIdx 时乘以 2 并没有考虑这个 - - 因为它只是乘以 tid).

  2. 在您的主机代码中,您的块数计算不正确:

    vectorReduceAsSumKernel.GridDimensions = array.Length / threadsPerBlock;
    

    array.Length 的值为 2048,threadsPerBlock 的值为 256,这将创建 8 个块。但是正如您已经指出的那样,您的意图是启动块 (2048/512)。所以你需要将分母乘以2:

    vectorReduceAsSumKernel.GridDimensions = array.Length / (2*threadsPerBlock);
    

另外,你的缩小扫描模式坏了。它依赖于 warp 执行顺序,以提供正确的结果,并且 CUDA 不指定 warp 执行顺序。

为了说明原因,让我们举一个简单的例子。我们只考虑一个线程块,数组的起点全为1,就像你初始化它一样。

现在,warp 0 由线程 0-31 组成。你的归约扫描操作是这样的:

a[i] = a[leftSumElementIdx] + a[leftSumElementIdx + 1];

因此 warp 0 中的每个线程将收集另外两个值并将它们相加,并存储它们。线程 31 将取值 a[62]a[63] 并将它们相加。如果 a[62]a[63] 的值仍然是 1,如初始化的那样,那么这将按预期工作。但是 a[62]a[63] 的值由 warp 1 写入 ,由线程 32-63 组成。因此,如果 warp 1 在 warp 0 之前执行(完全合法),那么您将得到不同的结果。这是全局内存竞争条件。这是因为您的输入数组既是中间结果的来源又是目标,__syncthreads() 不会为您解决这个问题。它不会强制扭曲以任何特定顺序执行。

一种可能的解决方案是修复扫描模式。在任何给定的缩减周期中,让我们有一个扫描模式,其中每个线程写入和读取在该周期期间未被任何其他线程触及的值。内核代码的以下改编实现了这一点:

__global__ void Vector_Reduce_As_Sum_Kernel(float* out, float* a)
{
  int tid = threadIdx.x;
  int blockStart = blockDim.x * blockIdx.x * 2;
  int i = tid + blockStart;

  for (int j = blockDim.x; j > 0; j>>=1){
    if (tid < j)
      a[i] += a[i+j];

    __syncthreads();}

  if (tid == 0)
  {
    out[blockIdx.x] = a[i];
  }
}

对于一般目的的归约,这仍然是一个非常慢的方法。这个 tutorial 涵盖了如何编写更快的归约。而且,正如已经指出的那样,managedCuda 可能有一些方法可以完全避免编写内核。