用于形态学膨胀卷积等算法的库达原位记忆竞赛问题

Question

我用CUDA编写了一个扩展内核，当我的输入和输出图像是不同的缓冲区时，它工作得很好，但当我在原地调用内核时，我理解为内存竞争问题，即输入和输出缓冲区指向相同的内存位置。

我试过：

a.利用合作团体，

b.使用互斥和原子加法，但正如本论文和网络上的几个来源所建议的，

c.采用无锁块间同步，提出了同一张纸中的同步方案.

我所有的尝试都失败了，因为：

答:因为我的输入缓冲区是一个const指针，当我不得不将它转换成一个void*参数(这是有意义的)时，我有一个编译错误(这是有意义的)，所以我不能更进一步。

没有工作是因为我面临着一种普遍的行为:我有16x16块，每个块都有32x32个线程。同步块应该将互斥量增加到256个，但是程序块在添加48个原子后会增加。

c不起作用，因为它没有块间同步，尽管我直接从纸上使用的代码对我来说似乎很好。我可以通过添加一些__syncthreads()来稍微提高比赛效果。

这是膨胀函数；

template <typename T>
__global__ void GenericDilate2dImg_knl(const ImageSizeInfo imgSizeInfo,
                                       volatile int* syncArrayIn, volatile int* syncArrayOut, 
                                       const unsigned long localSizeX, const unsigned long localSizeY,
                                       const int borderPolicyType, const T outOfImageValue,
                                       const struct StructuringElementInfo seInfo,
                                       const T* pInBuf, T* pOutBuf)
{
    // Extract sizeX, sizeY, etc. from imgSizeInfo
    SPLIT_SIZES_FROM_STRUCT(imgSizeInfo)

    // Declare the shared buffer pSharedBuf
    extern __shared__ char pSharedMem[];
    T* pSharedBuf = reinterpret_cast<T*>(pSharedMem);

    const unsigned long x = blockDim.x * blockIdx.x + threadIdx.x;
    const unsigned long y = blockDim.y * blockIdx.y + threadIdx.y;
    const unsigned long planIdx = blockDim.z * blockIdx.z + threadIdx.z;

    const unsigned long nbPlans = sizeZ * sizeC * sizeT;

    const unsigned long idx = x + y * sizeX + planIdx * sizeX*sizeY;

    // Copy the input image data into shared memory
    if (x < blockDim.x * gridDim.x && y < blockDim.y * gridDim.y && planIdx < blockDim.z * gridDim.z) {
        copyDataToSharedMemory2d(pInBuf, sizeX, sizeY, planIdx,
                                 localSizeX, localSizeY, 
                                 seInfo._paddingX, seInfo._paddingY,
                                 borderPolicyType, outOfImageValue,
                                 pSharedBuf);
    }

    // Wait to ensure that the copy is terminated
    if (pInBuf == pOutBuf) {
        // Grid synchronization for in-situ case
        //__gpu_sync(gridDim.x * gridDim.y);        // Use a mutex
        __gpu_sync2(1, syncArrayIn, syncArrayOut);  // Use a lock-free barrier
    }
    else
        // The input and ouput buffers point to different data 
        // -> we simply need to synchronize the threads inside the block
        __syncthreads();

    // Compute the convolution for pixels inside the image
    if (x < sizeX && y < sizeY && planIdx < nbPlans) {
        T vMax = 0;
        for (unsigned int curCoefIdx = 0; curCoefIdx < seInfo._nbOffsets; ++curCoefIdx) {
            const unsigned int sx = threadIdx.x + seInfo._paddingX + seInfo._pOffsetsX[curCoefIdx];
            const unsigned int sy = threadIdx.y + seInfo._paddingY + seInfo._pOffsetsY[curCoefIdx];
            const unsigned long sidx = sx + sy * localSizeX;
            const T curVal = pSharedBuf[sidx];
            vMax = (vMax > curVal ? vMax : curVal);
        }

        // Round the result
        pOutBuf[idx] = vMax;
    }
}

我要从全局内存复制到共享内存的功能是：

template <typename T>
__device__ void copyDataToSharedMemory2d(const T* pInBuf,
                                         const unsigned long sizeX, const unsigned long sizeY, const unsigned long planIdx,
                                         const unsigned long localSizeX, const unsigned long localSizeY,
                                         const int paddingX, const int paddingY,
                                         const int borderPolicyType, const T outOfImageValue,
                                         T* pSharedBuf)
{
    const int x = blockDim.x * blockIdx.x + threadIdx.x;
    const int y = blockDim.y * blockIdx.y + threadIdx.y;
    const int localX = threadIdx.x;
    const int localY = threadIdx.y;

    // Fill the shared buffer tile by tile
    // A tile is related to the group size
    const unsigned int groupSizeX = blockDim.x;
    const unsigned int groupSizeY = blockDim.y;

    // For each tile
    for (int offsetY = 0; offsetY < localSizeY; offsetY += groupSizeY) {
        int curLocalY = localY + offsetY;
        int curGlobalY = y + offsetY - paddingY;
        for (int offsetX = 0; offsetX < localSizeX; offsetX += groupSizeX) {
            int curLocalX = localX + offsetX;
            int curGlobalX = x + offsetX - paddingX;

            // If the current coordinate is inside the shared sub-image
            if (curLocalX < localSizeX && curLocalY < localSizeY) {
                const int idx = curLocalX + curLocalY * localSizeX;
                pSharedBuf[idx] = getPixel2d(pInBuf, sizeX, sizeY, curGlobalX, curGlobalY, planIdx, borderPolicyType, outOfImageValue);
            }
        }
    }
}

其中，getPixel2d允许我管理图像之外的数据：

template <typename T>
__device__
T getPixel2d(const T* pInBuf,
             const unsigned long sizeX, const unsigned long sizeY,
             const int x, const int y, const int z,
             const int borderPolicyType, const T outOfImageValue)
{
    int x_inside = x;
    if (x < 0 || x >= sizeX) {
        switch (borderPolicyType) {
        case 0://outside the image, there is a constant value
            return outOfImageValue;
        case 1://outside the image, we propagate the data at the image borders
            if (x < 0)
                x_inside = 0;
            else // x >= sizeX
                x_inside = sizeX - 1;
            break;
        case 2://Miror effect
            if (x < 0)
                x_inside = -(x + 1);
            else // x >= sizeX
                x_inside = sizeX - ((x - sizeX) + 1);
            break;
        }
    }

    // y-coordinate inside the image
    int y_inside = y;
    if (y < 0 || y >= sizeY) {
        switch (borderPolicyType) {
        case 0://outside the image, there is a constant value
            return outOfImageValue;
        case 1://outside the image, we propagate the data at the image borders
            if (y < 0)
                y_inside = 0;
            else // y >= sizeY
                y_inside = sizeY - 1;
            break;
        case 2://Miror effect
            if (y < 0)
                y_inside = -(y + 1);
            else // y >= sizeY
                y_inside = sizeY - ((y - sizeY) + 1);
            break;
        default: break;
        }
    }

    return pInBuf[x_inside + y_inside * sizeX + z * sizeX * sizeY];
}

现在，下面是我的块间同步功能：

// Using a mutex
__device__ volatile int g_mutex;
__device__ void __gpu_sync(int goalVal) {
    //thread ID in a block
    int tid_in_block = threadIdx.x * blockDim.y + threadIdx.y;
    // only thread 0 is used for synchronization
    if (tid_in_block == 0) {
        atomicAdd((int*)&g_mutex, 1);
        printf("[%d] %d Vs %d\n", blockIdx.x * gridDim.y + blockIdx.y, g_mutex, goalVal);
        //only when all blocks add 1 to g_mutex
        //will g_mutex equal to goalVal
        while (g_mutex </*!=*/ goalVal) {
            ;//Do nothing here
        }
    }
    __syncthreads();
}

// Lock-free barrier
__device__ void __gpu_sync2(int goalVal, volatile int* Arrayin, volatile int* Arrayout) {
    // thread ID in a block
    int tid_in_blk = threadIdx.x * blockDim.y + threadIdx.y;
    int nBlockNum = gridDim.x * gridDim.y;
    int bid = blockIdx.x * gridDim.y + blockIdx.y;
    // only thread 0 is used for synchronization
    if (tid_in_blk == 0) {
        Arrayin[bid] = goalVal;
    }
    if (bid == 1) {
        if (tid_in_blk < nBlockNum) {
            while (Arrayin[tid_in_blk] != goalVal) {
                ;//Do nothing here
            }
        }
        __syncthreads();
        if (tid_in_blk < nBlockNum) {
            Arrayout[tid_in_blk] = goalVal;
        }
    }
    if (tid_in_blk == 0) {
        while (Arrayout[bid] != goalVal) {
            ;//Do nothing here
        }
    }
    __syncthreads();
}

我在现场计算得到的图像是：

​

​

我使用了一个11x15结构功能，共享缓冲区的大小是(nbThreadsPerBlock+2*paddindX) * (nbThreadsPerBlock+2*paddindY)。错误的结果(由箭头显示)出现在某些块的顶部，但总是位于相同的位置和相同的值。我希望记忆竞赛的结果会更随机.

是否有更好的方法来管理就地计算，或者有什么理由阻止网格同步工作？

编辑我使用的图像的大小是510x509，我在NVidia QuadroRTX5000上运行我的代码。

Answer 1

我通常会建议这样一个问题的最小可重现性的例子，以及你正在运行的GPU的指示，但我们可能可以不这样做。简而言之，正如你已经发现的那样，你想要做的事情不会可靠地发挥作用。

您已经选择了一个线程策略，在网格中为每个输出点分配一个线程：

    pOutBuf[idx] = vMax;

这是明智和好的。基于此，我想：

我有16x16块，每个块都有32x32个线程。

输入图像为512x512 (每个方向为16x32个线程，每个输出点一个线程)。

正如您已经说过的，您在网格中有256个块(1024个线程中的每个线程)。此外，对于就地情况，我们可以将内核简化为以下伪代码：

__global__ void GenericDilate2dImg_knl(...){
  read_in_image();
  grid_wide_sync();
  write_out_image();
}

因此，要使这种方法工作，read_in_image()步骤必须能够读取整个图像，然后才能进行任何写入。但是，在一般情况下，您的方法将不起作用，而且显然也不适用于您的特定GPU。为了按照上面的方式读取整个图像，我们必须同时将网格中的每个线程块驻留在GPU中的SMs上。所有256个块都需要存放，并在SM上运行。但是GPU并没有为这样的事情提供内在的保证。如果GPU中有24条SMs，每个SMs最多可以容纳2048个线程，那么GPU的“运行”或“瞬时”容量为24*2048个线程，或48个线程块。将没有足够的空间让所有256个线程块运行。您的算法不仅依赖于此，而且所有3种网格同步方法都依赖于这个概念。

在48个“原子加法”之后，您的第二个网格同步方法停止了，这一事实向我提供了上面的示例数字。这是一个看似合理的解释，为什么该方法可能会失败:您的GPU只允许您的48个线程块驻留，而其他208个线程块正在等待，还没有存放在任何SM上，因此不允许它们的任何线程运行。这208个线程块中的线程需要运行，以获取相关的输入数据，并满足网格范围内同步的要求。但是他们没有在运行，因为他们在等待空间来打开一个SM。空间永远不会在SM上打开，因为完整的SMs有等待网格同步点的线程块。所以你陷入僵局了。

这个问题在一般情况下不容易解决。任何网格同步机制，包括协作组，都有一个固有的要求，即所有线程块实际上都可以在您的特定GPU上同时调度。因此，在一般情况下，当我们不知道数据集大小或我们将要运行的GPU时，这个问题是相当困难的。

一种可能的方法是将输入数据集划分为区域，并让内核一次处理一个区域。这可能需要多个网格同步，一个用于处理每个区域中的入/出分区，另一个用于处理内核在跨区域过程中的进度。您还必须小心处理区域边缘。

如果您知道数据集大小和正在运行的GPU的具体情况，另一种可能的方法是确保运行在GPU上的“足够大”来处理数据集大小。例如，一个A100图形处理器可能同时驻留216个块，因此在这种情况下，您可以处理一个稍微小一些的图像大小，可能是14x32=448高度和448个宽度维度。

考虑到针对这个特定示例的就地或现场工作的这些方法需要相当的复杂性，我个人强烈地希望在输出不同于输入的情况下使用这种方法。这种方法可能也会运行得更快。从性能角度来看，网格范围内的同步并不是一个“自由”的结构。