在OpenCL中计算部分和

Question

一个一维数据集被分成多个段，每个工作项处理一个段。它从段中读取一定数量的元素。这个数量事先不知道，并且对于每个段来说都不同。

例如：

+----+----+----+----+----+----+----+----+----+     <-- segments
  A    BCD  E    FG  HIJK   L    M        N        <-- elements in this segment

在处理完所有的段落后，应该将元素写入连续的输出内存中，例如：

A B C D E F G H I J K L M N

因此，一个片段中元素的绝对输出位置取决于前面片段中元素的数量。 E 在位置4处，因为片段1包含1个元素（A），而片段2包含3个元素。

---

OpenCL内核将每个片段的元素数量写入本地/共享内存缓冲区，并按以下方式工作（伪代码）

kernel void k(
    constant uchar* input,
    global int* output,
    local int* segment_element_counts
) {
    int segment = get_local_id(0);
    int count = count_elements(&input[segment * segment_size]);

    segment_element_counts[segment] = count;

    barrier(CLK_LOCAL_MEM_FENCE);

    ptrdiff_t position = 0;
    for(int previous_segment = 0; previous_segment < segment; ++previous_segment)
        position += segment_element_counts[previous_segment];

    global int* output_ptr = &output[position];
    read_elements(&input[segment * segment_size], output_ptr);
}

每个工作项都必须使用循环计算部分和，其中具有较大id的工作项执行更多次迭代。

在OpenCL 1.2中，是否有更有效的实现方式（每个工作项计算一个序列的部分和，直到其索引为止）？ OpenCL 2似乎提供了work_group_scan_inclusive_add来解决这个问题。

- tmlen

Answer 1

您可以使用以下类似的方法在 log2(N) 次迭代中完成 N 个部分（前缀）求和：

offsets[get_local_id(0)] = count;
barrier(CLK_LOCAL_MEM_FENCE);

for (ushort combine = 1; combine < total_num_segments; combine *= 2)
{
    if (get_local_id(0) & combine)
    {
        offsets[get_local_id(0)] +=
            offsets[(get_local_id(0) & ~(combine * 2u - 1u)) | (combine - 1u)];
    }
    barrier(CLK_LOCAL_MEM_FENCE);
}

给定分段元素的数量

a     b     c        d

连续迭代将产生：

a     b+a   c        d+c

并且

a     b+a   c+(b+a)  (d+c)+(b+a)

这就是我们想要的结果。

因此，在第一次迭代中，我们将段元素计数分为2组，并在其中求和。然后，我们每次将2个组合并成4个元素，并将第一个组的结果传播到第二个组。我们再次将组扩大到8个，以此类推。

关键观察结果也符合每个段的索引的二进制表示:

0: 0b00  1: 0b01  2: 0b10  3: 0b11

索引0不执行任何求和操作。索引1和3在第一个迭代中（位0 / LSB = 1）执行求和操作，而索引2和3在第二个迭代中（位1 = 1）执行求和操作。这解释了这行代码的含义：

    if (get_local_id(0) & combine)

当然需要解释的另一种声明是：

。

        offsets[get_local_id(0)] +=
            offsets[(get_local_id(0) & ~(combine * 2u - 1u)) | (combine - 1u)];

计算我们想要累加到工作项总和的前缀和所在的索引有一点棘手。子表达式(combine * 2u - 1u)在每次迭代（从1开始）时取值为(2ⁿ-1)：

1 = 0b001
3 = 0b011
7 = 0b111
…

通过按位屏蔽这些位后缀（即 i & ~x）工作项索引，这给出了当前组中第一个项目的索引。然后，(combine - 1u) 子表达式给出了第一半的最后一个项目在当前组内的索引。将两者结合起来，就可以得到要累加到当前段中的项目的总索引。结果中有一个小问题：它向左移动了一位：因此，段 1 需要使用 offsets[0]，以此类推，而段 0 的偏移量当然为 0。您可以通过过度分配偏移数组 1 并对从索引 1 开始的子数组执行前缀求和并将索引 0 初始化为 0，或使用条件语句来解决这个问题。上述代码可能存在基于分析优化的微小优化。

- pmdj