CUDA形态学膨胀等算法原地运算内存竞争问题咨询

问题描述 #

我在CUDA中编写了形态学膨胀核函数，当输入输出图像使用不同缓冲区时运行正常，但当输入输出指向同一块内存（即原地in-situ运算）时，出现了内存竞争问题。

已尝试的解决方案及失败原因 #

• 方案a：使用协作组（cooperative groups）
  失败原因：因输入缓冲区为const指针，转型为void*参数时报编译错误，无法继续推进
• 方案b：参考相关论文及多个网络资料的方案，使用互斥锁+原子加法
  失败原因：出现异常表现：使用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)
{
    // 从imgSizeInfo提取sizeX、sizeY等参数
    SPLIT_SIZES_FROM_STRUCT(imgSizeInfo)

    // 声明共享缓冲区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;

    // 将输入图像数据拷贝到共享内存
    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);
    }

    // 等待拷贝完成
    if (pInBuf == pOutBuf) {
        // 原地运算场景下的网格同步
        //__gpu_sync(gridDim.x * gridDim.y);        // 使用互斥锁
        __gpu_sync2(1, syncArrayIn, syncArrayOut);  // 使用无锁屏障
    }
    else
        // 输入输出缓冲区指向不同数据，仅需同步块内线程
        __syncthreads();

    // 计算图像内像素的卷积
    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);
        }

        // 写入结果
        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;

    // 逐tile填充共享缓冲区，tile与线程组大小关联
    const unsigned int groupSizeX = blockDim.x;
    const unsigned int groupSizeY = blockDim.y;

    // 遍历每个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 (curLocalX < localSizeX && curLocalY < localSizeY) {
                const int idx = curLocalX + curLocalY * localSizeX;
                pSharedBuf[idx] = getPixel2d(pInBuf, sizeX, sizeY, curGlobalX, curGlobalY, planIdx, borderPolicyType, outOfImageValue);
            }
        }
    }
}

像素读取函数（处理边界） #

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:// 图像外使用常量值
            return outOfImageValue;
        case 1:// 图像外使用边界像素值
            if (x < 0)
                x_inside = 0;
            else // x >= sizeX
                x_inside = sizeX - 1;
            break;
        case 2:// 镜像效果
            if (x < 0)
                x_inside = -(x + 1);
            else // x >= sizeX
                x_inside = sizeX - ((x - sizeX) + 1);
            break;
        }
    }

    // 处理y坐标边界
    int y_inside = y;
    if (y < 0 || y >= sizeY) {
        switch (borderPolicyType) {
        case 0:// 图像外使用常量值
            return outOfImageValue;
        case 1:// 图像外使用边界像素值
            if (y < 0)
                y_inside = 0;
            else // y >= sizeY
                y_inside = sizeY - 1;
            break;
        case 2:// 镜像效果
            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];
}

块间同步函数 #

// 使用互斥锁实现
__device__ volatile int g_mutex;
__device__ void __gpu_sync(int goalVal) {
    // 块内线程ID
    int tid_in_block = threadIdx.x * blockDim.y + threadIdx.y;
    // 仅使用0号线程完成同步
    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);
        // 等待所有块完成累加，g_mutex等于goalVal时结束等待
        while (g_mutex < goalVal) {
            ;// 空转
        }
    }
    __syncthreads();
}

// 无锁屏障实现
__device__ void __gpu_sync2(int goalVal, volatile int* Arrayin, volatile int* Arrayout) {
    // 块内线程ID
    int tid_in_blk = threadIdx.x * blockDim.y + threadIdx.y;
    int nBlockNum = gridDim.x * gridDim.y;
    int bid = blockIdx.x * gridDim.y + blockIdx.y;
    // 仅使用0号线程完成同步
    if (tid_in_blk == 0) {
        Arrayin[bid] = goalVal;
    }
    if (bid == 1) {
        if (tid_in_blk < nBlockNum) {
            while (Arrayin[tid_in_blk] != goalVal) {
                ;// 空转
            }
        }
        __syncthreads();
        if (tid_in_blk < nBlockNum) {
            Arrayout[tid_in_blk] = goalVal;
        }
    }
    if (tid_in_blk == 0) {
        while (Arrayout[bid] != goalVal) {
            ;// 空转
        }
    }
    __syncthreads();
}

问题现象 #

原地运算得到的结果图像如下：

使用了11x15的结构元素，共享缓冲区大小为(nbThreadsPerBlock+2*paddindX) * (nbThreadsPerBlock+2*paddindY)。错误结果（箭头标注处）出现在部分块的顶部，位置和数值均固定，不符合内存竞争通常出现的随机错误表现。
运行环境：图像大小为510x509，运行显卡为NVidia Quadro RTX 5000。

问题疑问 #

是否有更优的原地运算实现方案，或是当前网格同步失效的原因是什么？ #

解答 #

网格同步失效的核心原因 #

你自定义的块间同步机制不符合CUDA的执行模型：NVIDIA GPU的块调度是无序的，且RTX 5000（图灵架构）的SM最多同时驻留几十个块，你启动的256个块无法全部同时在SM上运行，已经运行的块会在空转等待未调度的块时直接造成死锁，这就是互斥锁方案累加48次就阻塞的直接原因（48接近你显卡单SM驻留块上限×SM总数）。
无锁同步方案同样存在该问题，且你使用的无锁代码依赖块ID为1的块先被调度，一旦该块后调度，所有其他块都会直接卡死，同步完全失效。

正确的原地膨胀实现方案 #

无需复杂的网格同步，直接按以下步骤修改即可：

1. 原地运算场景下，在核函数启动前申请一块和输入图像等大的临时全局内存
2. 将核函数的输出先写入临时缓冲区，所有块计算完成后，再用cudaMemcpy将临时缓冲区的数据拷贝回原输入缓冲区，释放临时内存即可。
   如果需要完全避免临时内存开销，可以使用流并分阶段处理图像块，保证相邻块计算时不会互相覆盖需要的原始数据，但开发复杂度更高，常规场景下临时缓冲区方案性价比最高。

针对协作组方案的编译问题 #

const指针转型报错是语法问题，协作组的网格同步功能不需要修改输入缓冲区指针，你只需将协作组同步逻辑加在共享内存拷贝完成后、结果写入全局内存前即可，注意使用协作组网格同步需要用cudaLaunchCooperativeKernel启动核函数，且需要确认显卡支持该特性（你的RTX 5000支持）。 #

内容的提问来源于stack exchange，提问作者ractiv