sorting - OpenCL中元素的排序和计数

Question

我想创建一个OpenCL内核，它可以对数百万个ulong进行排序和计数。
有一种适合我需要的特定算法，还是应该去散列表？

要清楚，给出以下输入：

[42, 13, 9, 42]

我想得到这样的输出：

[(9,1), (13,1), (42,2)]

我的第一个想法是修改计数排序-已经进行计数以进行排序-但由于范围广泛，因此需要太多内存。 Bitonic或Radix排序加上某种元素计数方法可能是一种方法，但我想念一种快速的元素计数方法。有什么建议吗？

额外说明：


我正在使用NVIDIA Tesla K40C GPU和Terasic DE5-Net FPGA进行开发。到目前为止，主要目标是使其能够在GPU上运行，但我也对可能非常适合FPGA的解决方案感兴趣。
我知道ulong范围内的某些值未使用，因此我们可以使用它们来标记无效的元素或重复项。
我想使用CPU中的多个线程来消耗GPU的输出，因此希望避免使用任何需要一些后处理（在主机端是这样）的解决方案，该解决方案在输出周围散布着数据。

Answer 1

此解决方案需要两次通过两次音调检查，以对重复项进行计数并删除重复项（将其移至数组末尾）。 Bitonic排序为O(log(n)^2)，因此它将以时间复杂度2(log(n)^2)运行，除非您在循环中运行它，否则应该不会有问题。

为每个元素创建一个简单的结构，以包括重复项的数量，如果该元素已作为重复项添加，则类似于：

// Note: If you are worried about space, or know that there 
// will only be a few duplicates for each element, then 
// make the count element smaller
typedef struct {
  cl_ulong value;
  cl_ulong count : 63;
  cl_ulong seen  : 1; 
} Element; 

算法：

您可以从创建比较功能开始，该功能将重复项移到末尾，如果要将重复项添加到元素的总数中，则对重复项进行计数。这是比较功能背后的逻辑：


如果一个元素是重复元素，而另一个元素不是，则返回非重复元素较小（与值无关），这会将所有重复元素移到末尾。
如果元素是重复项，并且右元素未标记为重复（seen=0），则将右元素的计数添加到左元素的计数，并将右元素设置为重复（seen=1）。这具有将具有特定值的元素的总计数移动到具有该值的数组中最左边的元素的作用，并将所有具有该值的重复项移到末尾。
否则，返回值较小的元素较小。


比较功能如下所示：

bool compare(const Element* E1, const Element* E2) {
  if (!E1->seen && E2->seen) return true;  // E1 smaller
  if (!E2->seen && E1->seen) return false; // E2 smaller

  // If the elements are duplicates and the right element has 
  // not yet been "seen" by an element with the same value 
  if (E1->value == E2->value && !E2->seen) {
    E1->count += E2->count;   
    E2->seen = 1;
    return true;
  }

  // They aren't duplicates, and either
  // neither has been seen, or both have
  return E1->value < E2->value;
}

双音排序具有特定的结构，可以通过图表很好地说明。在该图中，每个元素由一个三元组(a,b,c)引用，其中a = value，b = count和c = seen。

每个图都显示了阵列上的一次双音阶排序（垂直线表示元素之间的比较，而水平线则右移至双音阶排序的下一个阶段）。使用该图以及上面的比较功能和逻辑，您应该能够使自己确信，这满足了要求。

运行1：

运行2：

运行2结束时，所有元素均按值排列。带有seen = 1的副本位于末尾，带有seen = 0的副本位于正确的位置，并且count是具有相同值的其他元素的数量。

实现方式：

这些图用颜色编码以说明双音排序的子过程。我将蓝色块称为一个阶段（图中的每次运行中都有三个阶段）。通常，每次运行都会有ceil(log(N))个阶段。每个阶段都包含多个绿色块（我将这些out-in块称为比较的形状），和红色块（由于它们之间的距离而称为这些constant块）要比较的元素保持不变）。

从图中可以看出，out-in块大小（每个块中的元素）从2开始，并且在每次通过中翻倍。每次通过的constant块大小从out-in块大小的一半开始（在第二（蓝色块）阶段，四个红色块中的每个都有2个元素，因为绿色块的大小为4）相中红色块的每个连续垂直线的一半。同样，一个阶段中constant（红色）块的连续垂直线数始终与索引为0的阶段数相同（相位0的红色块的垂直线为0，相位红色的块的垂直线为1）。 1和第2阶段红色块的2条垂直线-每条垂直线都是调用该内核的迭代）。

然后，您可以为out-in传递创建内核，然后constant传递，然后从主机端调用内核（因为您需要不断同步，这是一个缺点，但是您仍然应该看到顺序性能有了很大的提高）实施）。

从宿主的角度来看，总体上的音调排序可能看起来像：

cl_uint num_elements = 4; // Set number of elements
cl_uint phases       = (cl_uint)ceil((float)log2(num_elements));
cl_uint out_in_block_size = 2;
cl_uint constant_block_size;

// Set the elements_buffer, which should have been created with
// with clCreateBuffer, as the first kernel argument, and the 
// number of elements as the second kernel argument
clSetKernelArg(out_in_kernel, 0, sizeof(cl_mem), (void*)(&elements_buffer));
clSetKernelArg(out_in_kernel, 1, sizeof(cl_uint), (void*)(&num_elements));
clSetKernelArg(constant_kernel, 0, sizeof(cl_mem), (void*)(&elements_buffer));
clSetKernelArg(constant_kernel, 1, sizeof(cl_uint), (void*)(&num_elements));

// For each pass
for (unsigned int phase = 0; phase < phases; ++phase) { 
  //  -------------------- Green Part ------------------------ //

  // Set the out_in_block size for the kernel
  clSetKernelArg(out_in_kernel, 2, sizeof(cl_int), (void*)(&out_in_block_size));

  // Call the kernel - command_queue is the clCommandQueue 
  // which should have been created during cl setup 
  clEnqueNDRangeKernel(command_queue    , // clCommandQueue  
                       out_in_kernel    , // The kernel
                       1                , // Work dim = 1 since 1D array
                       NULL             , // No global offset
                       &global_work_size,
                       &local_work_size ,
                       0                ,
                       NULL             ,
                       NULL);
  barrier(CLK_GLOBAL_MEM_FENCE); // Synchronise

  // ---------------------- End Green Part -------------------- //

  // Set the block size for constant blocks based on the out_in_block_size
  constant_block_size = out_in_block_size / 2;

  //  --------------------  Red Part  ------------------------ // 

  for (unsigned int i  0; i < phase; ++i) {
    // Set the constant_block_size as a kernel argument
    clSetKernelArg(constant_kernel, 2, sizeof(cl_int), (void*)(&constant_block_size));

    // Call the constant kernel
    clEnqueNDRangeKernel(command_queue    , // clCommandQueue  
                         constant_kernel  , // The kernel
                         1                , // Work dim = 1 since 1D array
                         NULL             , // No global offset
                         &global_work_size,
                         &local_work_size ,
                         0                ,
                         NULL             ,
                         NULL);
  barrier(CLK_GLOBAL_MEM_FENCE); // Synchronise

  // Update constant_block_size for next iteration
  constant_block_size /= 2; 
  }
  // ------------------- End Red Part ---------------------- //
}

然后内核就像（您还需要将struct typedef放入内核文件中，以便OpenCL编译器知道'Element'是什么）：

__global void out_in_kernel(__global Element* elements, unsigned int num_elements, unsigned int block_size) {
  const unsigned int idx_upper = // index of upper element in diagram.
  const unsigned int idx_lower = // index of lower element in diagram

  // Check that both indices are in range (this depends on thread mapping)
  if (idx_upper is in range && index_lower is in range) {
    // Do the comparison
    if (!compare(elements + idx_upper, elements + idx_lower) {
      // Swap the elements 
    }
  }
}

constant_kernel外观相同，但是线程映射（如何确定idx_upper和idx_lower）将有所不同。通常，有很多方法可以将线程映射到元素上以模仿图（请注意，所需的线程数是元素总数的一半，因为每个线程可以进行一次比较）。

另一个考虑因素是如何使线程映射具有通用性（这样，如果您有许多元素（不是二的幂），算法就不会中断）。