【问题标题】：Which way to order a shared 2D/3D array for parallel reduction over 1 dimension in CUDA/OpenCL?哪种方式可以在 CUDA/OpenCL 中订购共享 2D/3D 阵列以并行减少 1 维？
【发布时间】：2015-02-09 10:11:35
【问题描述】：

总体目标

我要对二分图进行一些简化，由两个用于顶点的密集数组和一个密集数组表示，该数组指定是否存在一条边。比如说，两个数组是 a0[] 和 a1[]，所有的边都像 e[i0][i1] （即从 a0 中的元素到 a1 中的元素）。

有~100+100 个顶点，和~100*100 条边，所以每个线程负责一条边。

任务 1：最大减少

对于 a0 中的每个顶点，我想找到与其连接的所有顶点（在 a1 中）的最大值，然后反过来：将结果分配给数组 b0，对于 a1 中的每个顶点，我想求连接顶点的最大 b0[i0]。

为此，我：

1) 加载到共享内存中

    #define DC_NUM_FROM_SHARED 16
    #define DC_NUM_TO_SHARED 16
    __global__ void max_reduce_down(
            Value* value1
        , Value* max_value_in_connected
        , int r0_size, int r1_size
        , bool** connected
        )
    {
        int id_from;
        id_from = blockIdx.x * blockDim.x + threadIdx.x;
        id_to   = blockIdx.y * blockDim.y + threadIdx.y;
        bool within_bounds = (id_from < r0_size) && (id_to < r1_size);

        //load into shared memory
        __shared__ Value value[DC_NUM_TO_SHARED][DC_NUM_FROM_SHARED]; //FROM is the inner (consecutive) dimension
        if(within_bounds)
            value[threadIdx.y][threadIdx.x] = connected[id_to][id_from]? value1[id_to] : 0;
        else
            value[threadIdx.y][threadIdx.x] = 0;
        __syncthreads();

        if(!within_bounds)
            return;

2) 减少

for(int stride = DC_NUM_TO_SHARED/2; threadIdx.y < stride; stride >>= 1)
{
    value[threadIdx.y][threadIdx.x] = max(value[threadIdx.y][threadIdx.x], dc[threadIdx.y + stride][threadIdx.x]);
    __syncthreads();
}

3) 回写

max_value_connected[id_from] = value[0][threadIdx.x];

任务2：最好的k

类似的问题，但减少仅适用于 a0 中的顶点，我需要从 a1 中的连接中找到 k 个最佳候选者（k 约为 5） .

1) 我用零元素初始化共享数组，除了第一个位置

int id_from, id_to;
id_from = blockIdx.x * blockDim.x + threadIdx.x;
id_to   = blockIdx.y * blockDim.y + threadIdx.y;

__shared Value* values[MAX_CHAMPS * CHAMPS_NUM_FROM_SHARED * CHAMPS_NUM_TO_SHARED]; //champion overlaps
__shared int* champs[MAX_CHAMPS * CHAMPS_NUM_FROM_SHARED * CHAMPS_NUM_TO_SHARED]; // overlap champions


bool within_bounds = (id_from < r0_size) && (id_to < r1_size);
int i = threadIdx.y * CHAMPS_NUM_FROM_SHARED + threadIdx.x;
if(within_bounds)
{
    values[i] = connected[id_to][id_from] * values1[id_to];
    champs[i] = connected[id_to][id_from] ? id_to : -1;
}
else
{
    values[i] = 0;
    champs[i] = -1;
}

for(int place = 1; place < CHAMP_COUNT; place++)
{
    i = (place * CHAMPS_NUM_TO_SHARED + threadIdx.y) * CHAMPS_NUM_FROM_SHARED + threadIdx.x;
    values[i] = 0;
    champs[i] = -1;
}
if(! within_bounds)
    return;
__syncthreads();

2) 减少它

for(int stride = CHAMPS_NUM_TO_SHARED/2; threadIdx.y < stride; stride >>= 1)
{
    merge_2_champs(values, champs, CHAMP_COUNT, id_from, id_to, id_to + stride);
    __syncthreads();
}

3) 写回结果

for(int place = 0; place < LOCAL_DESIRED_ACTIVITY; place++)
    champs0[place][id_from] = champs[place * CHAMPS_NUM_TO_SHARED * CHAMPS_NUM_FROM_SHARED + threadIdx.x];

问题

如何对共享数组中的元素进行排序（转置），以便内存访问更好地使用缓存？在这一点上是否重要，或者我可以从其他优化中获得更多？如果我需要针对任务 2 进行优化，转置边缘矩阵会更好吗？（据我了解，Task 1有对称性，所以没关系）。

附言

我已经延迟展开循环并在加载时进行第一次缩减迭代，因为在我探索更简单的方法之前，我认为这太复杂了。

对于任务 2，最好不要加载零元素，因为数组永远不需要增长，只有在完成 log k 步后才开始缩小。这将使它在共享内存中紧凑 k 倍！但我害怕生成的索引数学。

语法和正确性

不寻常的类型只是 typedef'ed ints/chars/etc - AFAIK，在 GPU 中，尽可能压缩它们是有意义的。我还没有运行代码，不需要检查索引错误。

另外，我正在使用 CUDA，但我也对 OpenCL 的视角感兴趣，因为我认为最好的解决方案应该是相同的，而且我将来也会使用 OpenCL。

【问题讨论】：

标签： algorithm caching cuda opencl reduction

【解决方案1】：

好的，我想我想通了。

我正在考虑的两种选择是在 y 维度上进行归约，并且独立于 x 维度，反之亦然（x 维度是连续的）。在任何情况下，调度程序都能够沿 x 维度将线程组装成 warp，因此可以保证一定的一致性。然而，将相干性扩展到扭曲之外会很棒。此外，由于共享数组的 2D/3D 特性，必须将维度限制为 16 甚至 8。

为确保 warp 内的合并，调度程序必须沿 x 维度组装 warp。

如果减少超过 x 维度，每次迭代后，warp 中的活动线程数将减半。但是，如果减少超过 y 维度，则活动扭曲的数量将减半。

所以，我需要减少超过 y。

除非转置（加载）最慢，否则属于异常情况。

【讨论】：

【解决方案2】：

合并缓冲区读取真的很重要；如果你不这样做，内核可能会慢 32 倍。如果它意味着能够执行它们，则值得进行重新排列通道（当然，重新排列通道也需要合并，但您通常可以利用共享本地内存来执行此操作）。

【讨论】：