OpenCL：并行求和n个整数

Question

我需要创建一个OpenCL内核函数，该函数使用并行算法对数组n中的numbers整数求和。

我应该使用类似于以下的算法：

parallel_summation(A):
    # ASSUME n = |A| is a power of 2 for simplicity

    # level 0:
    in parallel do:
      s[i] = A[i]                 for i = 0, 1, 2, ..., n-1

    # level 1:
    in parallel do:
      s[i] = s[i] + s[i+1]        for i = 0, 2, 4, ...

    # level 2:
    in parallel do:
      s[i] = s[i] + s[i+2]        for i = 0, 4, 8, ...

    # level 3:
    in parallel do:
      s[i] = s[i] + s[i+4]        for i = 0, 8, 16, ...

    # ...
    # level log_2( n ):
    s[0] = s[0] + s[n/2]

    return s[0]

所以，我想出了以下内核代码：

kernel void summation(global uint* numbers,
                      global uint* sum,
                      const  uint  n,
                      const  uint  work_group_size,
                      local  uint* work_group_buf,
                      const  uint  num_of_levels) {

    // lets assume for now that the workgroup's size is 16,
    // which is a power of 2.


    int i = get_global_id(0);

    if(i >= n)
        return;

    int local_i = get_local_id(0);

    uint step = 1;
    uint offset = 0;

    for(uint k = 0; k < num_of_levels; ++k) {

        if(k == 0) {

            work_group_buf[local_i] = numbers[i];

        }  else {

            if(local_i % step == 0) {
                work_group_buf[local_i] += work_group_buf[local_i + offset];
            }

        }

        if(offset == 0) {
            offset = 1;
        } else {
            offset *= 2;
        }

        step *= 2;

        barrier(CLK_LOCAL_MEM_FENCE);

    }

     atomic_add(sum, work_group_buf[0]);

}

但是有一个错误，因为我没有收到预期的结果。 numbers是一个包含从1到n的数字的缓冲区。 num_of_levels是log2（每个工作组的工作项数），在我当前的例子中是4（log2（16））。

我究竟做错了什么？

注意：我没有收到任何错误，只是错误的结果。例如，我有一个从0到999999的1000000个元素的数组，这些元素的总和应该是1783293664，但我得到的是1349447424。

Answer 1

我修了几个bug。有一些错误，我错过了这部分s[0] = s[0] + s[n/2]，你可以从这个新版本看到。

kernel void summation(global uint* numbers,
                          global uint* sum,
                          const  uint  n,
                          local  uint* work_group_buf,
                          const  uint  num_of_levels) {

        const int i = get_global_id(0);
        const int local_i = get_local_id(0);

        private uint step = 2;
        private uint offset = 1;


        if(i < n)
            work_group_buf[local_i] = numbers[i];

        barrier(CLK_LOCAL_MEM_FENCE);

        for(uint k = 1; k < num_of_levels; ++k) {

            if((local_i % step) == 0) {
                work_group_buf[local_i] += work_group_buf[local_i + offset];
            }

            offset *= 2;
            step *= 2;

            barrier(CLK_LOCAL_MEM_FENCE);
        }

        work_group_buf[0] += work_group_buf[get_local_size(0) / 2];

        if(local_i == 0)
            atomic_add(sum, work_group_buf[0]);

}

请注意，现在我只在sum中添加到最终work_group_buf只是每个work_group_buf[0]的第一个元素（即local_i == 0），因为该位置将包含工作组中所有元素的总和。

这实际上似乎适用于大小高达32的工作组（功率为2）。换句话说，这个内核似乎只适用于大小为2,4,8,16和32个工作项的工作组。