如何避免CUDA哈希表并行值比较中的竞争条件

在二维(x,y)空间中有大量不同尺寸的点，由结构体(x,y,size)数组表示。需要将该空间下采样至指定网格，将x、y取整至最近网格单元，保留每个单元中尺寸最大的点。

思路是使用哈希表存储由二维索引映射而来的一维索引，哈希表结构为(key:一维索引, value:最大尺寸)，再基于该哈希表生成输出。以下是内核中寻找最大尺寸的代码：

//nx: is the grid size in x
//dx/dy: is the cell size in x/y
//points: is the structure keeping the 2d points with the size as (x,y,size)
//minBounds: has the minimum values of x and y

int i = blockDim.x * blockIdx.x + threadIdx.x;
std::unordered_map<uint64_t, float> indexSizeMap;

// Find point index in the grid and flatten it to 1d index
auto xi = static_cast<uint64_t>(round((points[i].x - minBounds.x) / dx));
auto yi = static_cast<uint64_t>(round((points[i].y - minBounds.y) / dy));
uint64_t index1d = to1dIndex(xi, yi, nx + 1); //mapping 2d index to 1d index

// If the index in the voxel grid is vacant, fill it. Otherwise keep the larger size
if (indexSizeMap.count(index1d))
{
   if (indexSizeMap[index1d] < points[i].size)
   {
      indexSizeMap[index1d] = points[i].size;
   }
}
else
{
   indexSizeMap[index1d] = points[i].size;
}

这段代码存在竞争条件，当两个线程同时填充同一索引时会导致结果不确定：

if (indexSizeMap.count(index1d))
{
   if (indexSizeMap[index1d] < points[i].size)
   {
      indexSizeMap[index1d] = points[i].size;
   }
}
else
{
   indexSizeMap[index1d] = points[i].size;
}

查阅了竞争条件相关文档，大多是关于加减操作的处理方法，但未找到此类场景的解决方案。请问如何避免该竞争条件，或有其他实现方案建议？

方法1：原子操作实现线程安全的最大值更新 #

针对网格单元的最大值更新需求，可以用原子操作保证多线程下的操作原子性，避免竞争：

CUDA内核场景 #

如果是CUDA环境，直接使用atomicMax函数（CUDA 9+支持float类型的原子max）：

// 预先分配全局数组存储每个网格单元的最大size，数组大小为总网格数
__device__ float* gridMaxSizes;

// 线程内计算索引后执行原子更新
atomicMax(&gridMaxSizes[index1d], points[i].size);

若环境不支持float类型的atomicMax，可将float转换为uint32_t，利用整数原子max实现：

__device__ uint32_t* gridMaxSizesAsUint;

uint32_t current = __float_as_uint(points[i].size);
atomicMax(&gridMaxSizesAsUint[index1d], current);
// 后续使用时再转换回float：__uint_as_float(gridMaxSizesAsUint[index1d])

CPU多线程场景 #

用std::atomic结合CAS（比较并交换）实现自定义max原子操作：

#include <atomic>

void atomicMax(std::atomic<float>& target, float value) {
    float current = target.load();
    // 循环直到更新成功，或当前值已不小于目标值
    while (current < value && !target.compare_exchange_weak(current, value)) {}
}

// 线程内调用
atomicMax(gridMaxSizes[index1d], points[i].size);

方法2：预先分配网格数组替代哈希表 #

哈希表本身的多线程竞争处理复杂度高，不如预先计算所有可能的网格单元，分配数组存储每个单元的最大size（初始设为极小值），这种方式效率更高且更易保证线程安全：

// 预先计算总网格数：nx_total = nx + 1，ny_total根据y方向边界计算
size_t totalGridCount = nx_total * ny_total;
float* gridMaxSizes = new float[totalGridCount];
// 初始化所有值为负无穷，确保第一个点能覆盖初始值
std::fill(gridMaxSizes, gridMaxSizes + totalGridCount, -FLT_MAX);

// 线程内代码
int i = blockDim.x * blockIdx.x + threadIdx.x;
auto xi = static_cast<uint64_t>(round((points[i].x - minBounds.x) / dx));
auto yi = static_cast<uint64_t>(round((points[i].y - minBounds.y) / dy));
uint64_t index1d = to1dIndex(xi, yi, nx_total);

// 调用原子操作更新最大值（参考方法1的实现）
atomicMax(&gridMaxSizes[index1d], points[i].size);

方法3：先分组再单线程处理 #

若网格单元数量不多，可先将点按网格索引分组，再单线程处理每组的最大值，完全避免多线程竞争：

// 生成包含网格索引和对应size的配对数组
std::vector<std::pair<uint64_t, float>> indexSizePairs(points.size());
for (int i = 0; i < points.size(); ++i) {
    auto xi = static_cast<uint64_t>(round((points[i].x - minBounds.x) / dx));
    auto yi = static_cast<uint64_t>(round((points[i].y - minBounds.y) / dy));
    uint64_t index1d = to1dIndex(xi, yi, nx + 1);
    indexSizePairs[i] = {index1d, points[i].size};
}

// 按网格索引排序，将同索引的点集中在一起
std::sort(indexSizePairs.begin(), indexSizePairs.end());

// 遍历排序后的数组，提取每组的最大值
std::unordered_map<uint64_t, float> indexSizeMap;
if (indexSizePairs.empty()) return;

uint64_t currentIndex = indexSizePairs[0].first;
float maxSize = indexSizePairs[0].size;
for (size_t j = 1; j < indexSizePairs.size(); ++j) {
    if (indexSizePairs[j].first == currentIndex) {
        maxSize = std::max(maxSize, indexSizePairs[j].second);
    } else {
        indexSizeMap[currentIndex] = maxSize;
        currentIndex = indexSizePairs[j].first;
        maxSize = indexSizePairs[j].second;
    }
}
indexSizeMap[currentIndex] = maxSize;

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