c++ - CUDA GPU 上的小型二进制搜索

Question

我有一个int64_t 类型的大型设备数组inputValues。该数组的每 32 个元素按升序排序。我有一个未排序的搜索数组 removeValues。

我的目的是在inputValues 中查找removeValues 中的所有元素，并将它们标记为-1。实现这一目标的最有效方法是什么？如果有帮助，我正在使用 3.5 cuda 设备。

我不是在寻找更高级别的解决方案，即我不想使用 thrust 或 cub，但我想使用 cuda 内核来编写它。

我最初的方法是在一个线程 block 中加载共享内存中的每个 32 值。每个线程还从 removeValues 加载单个值，并在共享内存数组上执行独立的二进制搜索。如果找到，则使用 if 条件设置该值。

这种方法不会涉及很多银行冲突和分支机构分歧吗？你认为分支分歧可以通过在实现二分搜索时使用三元运算符来解决吗？即使解决了，银行冲突又如何消除？由于排序数组的大小是 32，是否可以使用 shuffle 指令实现二分查找？这会有帮助吗？

编辑:我添加了一个示例来展示我打算实现的目标。

假设 inputValues 是一个 vector ，其中每 32 个元素被排序:
[2, 4, 6, ..., 64], [95, 97, ..., 157], [1, 3, ..., 63], [...]

此数组的典型大小介于 32*2 到 32*32 之间。值的范围可以从 0 到 INT64_MAX。

removeValues 的一个例子是:
[7, 75, 95, 106]

此数组的典型大小范围从 1 到 1024。

removeValues 操作后将是: [-1, 75, -1, 106]

inputValues 中的值保持不变。

Answer 1

我同意@harrism 的回答(现已删除)和评论。由于我在非推力方法上付出了一些努力，因此我将介绍我的发现。

我尝试使用 __shfl() 天真地在 warp 级别实现二进制搜索，然后在数据集中重复该二进制搜索，将数据集传递给每个 32 元素组。

这很尴尬，但我的代码比 thrust 慢了大约 20 倍(事实上，如果你用 nvprof 仔细计时，它可能比 thrust 更糟糕)。

我使数据大小比问题中提出的稍大一些，因为问题中的数据大小太小以至于时间不准确。

下面是 2 种方法的完整示例:

1. 问题中大致概述了什么，即使用 warp shuffle 创建二分搜索，它可以针对 32 元素有序数组搜索最多 32 个元素。对尽可能多的 32 元素有序数组重复此过程，将整个数据集传递给每个有序数组(希望您现在可以开始看到一些低效率。)

2. 使用推力，与@harrism 概述的基本相同，即对分组数据集进行排序，然后运行矢量化 thrust::binary_search。

例子如下:

$ cat t1030.cu
#include <stdio.h>
#include <assert.h>
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/sort.h>
#include <thrust/binary_search.h>

typedef long mytype;

const int gsize = 32;
const int nGRP = 512;
const int dsize = nGRP*gsize;//gsize*nGRP;

#include <time.h>
#include <sys/time.h>
#define USECPSEC 1000000ULL

unsigned long long dtime_usec(unsigned long long start){

  timeval tv;
  gettimeofday(&tv, 0);
  return ((tv.tv_sec*USECPSEC)+tv.tv_usec)-start;
}

template <typename T>
__device__ T my_shfl32(T val, unsigned lane){
  return __shfl(val, lane);
}

template <typename T>
__device__ T my_shfl64(T val, unsigned lane){
  T retval = val;
  int2 t1 = *(reinterpret_cast<int2 *>(&retval));
  t1.x = __shfl(t1.x, lane);
  t1.y = __shfl(t1.y, lane);
  retval = *(reinterpret_cast<T *>(&t1));
  return retval;
}

template <typename T>
__device__ bool bsearch_shfl(T grp_val, T my_val){
  int src_lane = gsize>>1;
  bool return_val = false;
  T test_val;
  int shift = gsize>>2;
  for (int i = 0; i <= gsize>>3; i++){
    if (sizeof(T)==4){
      test_val = my_shfl32(grp_val, src_lane);}
    else if (sizeof(T)==8){
      test_val = my_shfl64(grp_val, src_lane);}
    else assert(0);
    if (test_val == my_val) return_val = true;
    src_lane += (((test_val<my_val)*2)-1)*shift;
    shift>>=1;
    assert ((src_lane < gsize)&&(src_lane > 0));}
  if (sizeof(T)==4){
    test_val = my_shfl32(grp_val, 0);}
  else if (sizeof(T)==8){
    test_val = my_shfl64(grp_val, 0);}
  else assert(0);
  if (test_val == my_val) return_val = true;
  return return_val;
}

template <typename T>
__global__ void bsearch_grp(const T * __restrict__ search_grps, T *data){

  int idx = threadIdx.x+blockDim.x*blockIdx.x;
  int tid = threadIdx.x;
  if (idx < gsize*nGRP){
    T grp_val = search_grps[idx];
    while (tid < dsize){
      T my_val = data[tid];
      if (bsearch_shfl(grp_val, my_val)) data[tid] = -1;
      tid += blockDim.x;}
  }
}


int main(){

  // data setup
  assert(gsize == 32);  //mandatory (warp size)
  assert((dsize % 32)==0);  //needed to preserve shfl capability
  thrust::host_vector<mytype> grps(gsize*nGRP);
  thrust::host_vector<mytype> data(dsize);
  thrust::host_vector<mytype> result(dsize);
  for (int i = 0; i < gsize*nGRP; i++) grps[i] = i;
  for (int i = 0; i < dsize; i++) data[i] = i;
  // method 1: individual shfl-based binary searches on each group
  mytype *d_grps, *d_data;
  cudaMalloc(&d_grps, gsize*nGRP*sizeof(mytype));
  cudaMalloc(&d_data, dsize*sizeof(mytype));
  cudaMemcpy(d_grps, &(grps[0]), gsize*nGRP*sizeof(mytype), cudaMemcpyHostToDevice);
  cudaMemcpy(d_data, &(data[0]), dsize*sizeof(mytype), cudaMemcpyHostToDevice);
  unsigned long long my_time = dtime_usec(0);
  bsearch_grp<<<nGRP, gsize>>>(d_grps, d_data);
  cudaDeviceSynchronize();
  my_time = dtime_usec(my_time);
  cudaMemcpy(&(result[0]), d_data, dsize*sizeof(mytype), cudaMemcpyDeviceToHost);
  for (int i = 0; i < dsize; i++) if (result[i] != -1) {printf("method 1 mismatch at %d, was %d, should be -1\n", i, (int)(result[i])); return 1;}
  printf("method 1 time: %fs\n", my_time/(float)USECPSEC);
  // method 2: thrust sort, followed by thrust binary search
  thrust::device_vector<mytype> t_grps = grps;
  thrust::device_vector<mytype> t_data = data;
  thrust::device_vector<bool> t_rslt(t_data.size());
  my_time = dtime_usec(0);
  thrust::sort(t_grps.begin(), t_grps.end());
  thrust::binary_search(t_grps.begin(), t_grps.end(), t_data.begin(), t_data.end(), t_rslt.begin());
  cudaDeviceSynchronize();
  my_time = dtime_usec(my_time);
  thrust::host_vector<bool> rslt = t_rslt;
  for (int i = 0; i < dsize; i++) if (rslt[i] != true) {printf("method 2 mismatch at %d, was %d, should be 1\n", i, (int)(rslt[i])); return 1;}
  printf("method 2 time: %fs\n", my_time/(float)USECPSEC);

  // method 3:  multiple thrust merges, followed by thrust binary search



  return 0;
}

$ nvcc -O3 -arch=sm_35 t1030.cu -o t1030
$ ./t1030
method 1 time: 0.009075s
method 2 time: 0.000516s
$

我在 Linux、CUDA 7.5、GT640 GPU 上运行它。显然，不同 GPU 上的性能会有所不同，但如果有任何 GPU 能显着缩小差距，我会感到惊讶。

简而言之，您最好使用像 thrust 或 cub 这样经过良好调整的库。如果你不喜欢推力的整体性，你可以试试 cub。我不知道 cub 是否有二进制搜索，但对整个排序数据集的单个二进制搜索是 not a difficult thing to write ，并且它是所涉及时间的较小部分(对于方法 2 - 可使用 nvprof 或其他计时代码识别)。

由于您的 32 元素分组范围已经排序，我还考虑了使用多个 thrust::merge 操作而不是单一排序的想法。我不确定哪个会更快，但由于 thrust 方法已经比 32 元素随机搜索方法快得多，我认为 thrust(或 cub)是显而易见的选择。