logo标签列表

在CUDA内核中搜索有序数组

arrayscuda
4
我正在编写一个CUDA核函数，每个线程都必须完成以下任务：假设我有一个已排序的无符号整数数组a（第一个为0），其长度为n，存储在共享内存中，每个线程都必须找到数组索引i，使得a[i] ≤ threadIdx.x且a[i + 1] > threadIdx.x。
一种简单的解决方案可能是：
for (i = 0; i < n - 1; i++)
    if (a[i + 1] > threadIdx.x) break;

但我想这不是最优的解决方式...有没有人能提供更好的建议呢？

第一个条件不是多余的吗？循环终止做了同样的事情。 - talonmies
我这样做是为了避免在 i == n - 1 时通过 a[i + 1] 访问 a 的边界外。 - Filippo Bistaffa
我理解了，我的意思是如果你迭代(0,n-1)，第一个条件就不是必需的。 - talonmies
是的，我明白你的意思，我会用你的改进更新代码。 - Filippo Bistaffa
我在下面发布了一个非常简单且有效的解决方案。 - Drout
3个回答
5

像Robert一样，我认为二分查找肯定比朴素循环更快--二分查找的操作次数上限是O(log(n))，而循环则是O(N)。

非常简单的实现方式：

#include <iostream>
#include <climits>
#include <assert.h>

__device__  __host__
int midpoint(int a, int b)
{
    return a + (b-a)/2;
}

__device__ __host__
int eval(int A[], int i, int val, int imin, int imax)
{

    int low = (A[i] <= val);
    int high = (A[i+1] > val);

    if (low && high) {
        return 0;
    } else if (low) {
        return -1;
    } else {
        return 1;
    }
}

__device__ __host__
int binary_search(int A[], int val, int imin, int imax)
{
    while (imax >= imin) {
        int imid = midpoint(imin, imax);
        int e = eval(A, imid, val, imin, imax);
        if(e == 0) {
            return imid;
        } else if (e < 0) {
            imin = imid;
        } else {         
            imax = imid;
        }
    }

    return -1;
}


__device__ __host__
int linear_search(int A[], int val, int imin, int imax)
{
    int res = -1;
    for(int i=imin; i<(imax-1); i++) {
        if (A[i+1] > val) {
            res = i;
            break;
        }
    }

    return res;
}

template<int version>
__global__
void search(int * source, int * result, int Nin, int Nout)
{
    extern __shared__ int buff[];
    int tid = threadIdx.x + blockIdx.x*blockDim.x;

    int val = INT_MAX;
    if (tid < Nin) val = source[threadIdx.x];
    buff[threadIdx.x] = val;
    __syncthreads();

    int res;
    switch(version) {

        case 0:
        res = binary_search(buff, threadIdx.x, 0, blockDim.x);
        break;

        case 1:
        res = linear_search(buff, threadIdx.x, 0, blockDim.x);
        break;
    }

    if (tid < Nout) result[tid] = res; 
}

int main(void)
{
    const int inputLength = 128000;
    const int isize = inputLength * sizeof(int);
    const int outputLength = 256;
    const int osize = outputLength * sizeof(int);

    int * hostInput = new int[inputLength];
    int * hostOutput = new int[outputLength];
    int * deviceInput;
    int * deviceOutput;

    for(int i=0; i<inputLength; i++) {
        hostInput[i] = -200 + 5*i;
    }

    cudaMalloc((void**)&deviceInput, isize);
    cudaMalloc((void**)&deviceOutput, osize);

    cudaMemcpy(deviceInput, hostInput, isize, cudaMemcpyHostToDevice);

    dim3 DimBlock(256, 1, 1);
    dim3 DimGrid(1, 1, 1);
    DimGrid.x = (outputLength / DimBlock.x) + 
                ((outputLength % DimBlock.x > 0) ? 1 : 0); 
    size_t shmsz = DimBlock.x * sizeof(int);

    for(int i=0; i<5; i++) {
        search<1><<<DimGrid, DimBlock, shmsz>>>(deviceInput, deviceOutput, 
                inputLength, outputLength);
    }

    for(int i=0; i<5; i++) {
        search<0><<<DimGrid, DimBlock, shmsz>>>(deviceInput, deviceOutput,
                inputLength, outputLength);
    }

    cudaMemcpy(hostOutput, deviceOutput, osize, cudaMemcpyDeviceToHost);

    for(int i=0; i<outputLength; i++) {
        int idx = hostOutput[i];
        int tidx = i % DimBlock.x;
        assert( (hostInput[idx] <= tidx) && (tidx < hostInput[idx+1]) );
    } 
    cudaDeviceReset();

    return 0;
}

相较于循环，这个方法可以提升大约五倍的速度：
>nvprof a.exe
======== NVPROF is profiling a.exe...
======== Command: a.exe
======== Profiling result:
 Time(%)      Time   Calls       Avg       Min       Max  Name
   60.11  157.85us       1  157.85us  157.85us  157.85us  [CUDA memcpy HtoD]
   32.58   85.55us       5   17.11us   16.63us   19.04us  void search<int=1>(int*, int*, int, int)
    6.52   17.13us       5    3.42us    3.35us    3.73us  void search<int=0>(int*, int*, int, int)
    0.79    2.08us       1    2.08us    2.08us    2.08us  [CUDA memcpy DtoH]

我相信有些聪明的人可以做得比这更好。但也许这至少给你一些想法。

点赞，这个答案比我的更好，原因至少有两个。它是一个完整的实例，并进行了时间比较，另外它更准确地跟踪了问题陈述，该问题陈述仅限于 threadIdx.x 范围内，并且因此可以利用共享内存。 - Robert Crovella
在《编程珠玑》一书中，Jon Bentley描述了一种非常棒的二分查找优化方法，可能与此相关。它涉及到对中点的仔细初始化，然后使用逐渐减小的2的幂进行条件更新。请参见http://www.cs.bell-labs.com/cm/cs/pearls/search.c中的“binarysearch4”。 - ArchaeaSoftware
1

有人能提供更好的建议吗？

一种蛮力的方法是让每个线程执行二分查找（在 threadIdx.x + 1 上）。

// sets idx to the index of the first element in a that is 
// equal to or larger than key

__device__ void bsearch_range(const int *a, const int key, const unsigned len_a, unsigned *idx){
  unsigned lower = 0;
  unsigned upper = len_a;
  unsigned midpt;
  while (lower < upper){
    midpt = (lower + upper)>>1;
    if (a[midpt] < key) lower = midpt +1;
    else upper = midpt;
    }
  *idx = lower;
  return;
  } 

__global__ void find_my_idx(const int *a, const unsigned len_a,  int *my_idx){
  unsigned idx = (blockDim.x * blockIdx.x) + threadIdx.x;
  unsigned sp_a;
  int val = idx+1;
  bsearch_range(a, val, len_a, &sp_a);
  my_idx[idx] = ((val-1) < a[sp_a]) ? sp_a:-1;
}

这段代码是在浏览器中编写的，未经测试。它是从一段工作代码中修改而来的。如果你有问题，我可以重新检查它。我不建议在没有缓存的设备上使用这种方法（如cc 1.x设备）。
实际上，这是在完整的唯一1D线程索引上进行搜索（blockDim.x * blockIdx.x + threadIdx.x + 1）。你可以将val更改为任何你喜欢的值。
如果你打算启动的线程数大于你的my_idx结果向量的长度，你还可以添加适当的线程检查。
我想象中可能有一种更聪明的方法，可能类似于前缀和。
-2

这是目前最好的算法。它被称为：LPW索引搜索

__global__ void find_position_lpw(int *a, int n)
{
    int idx = threadIdx.x;

    __shared__ int aux[ MAX_THREADS_PER_BLOCK /*1024*/ ];

    aux[idx] = 0;

    if (idx < n)
        atomicAdd( &aux[a[idx]], 1); // atomics in case there are duplicates

    __syncthreads();

    int tmp;

    for (int j = 1; j <= MAX_THREADS_PER_BLOCK / 2; j <<= 1)
    {
        if( idx >= j ) tmp = aux[idx - j];
        __syncthreads();
        if( idx >= j ) aux[idx] += tmp;
        __syncthreads();        
    }

    // result in "i"
    int i = aux[idx] - 1;

    // use "i" here...
    // ...
}
2
这段代码实际上并没有做原问题要求的事情——输入被假定为已经在共享内存中。而且它引入了一个原问题中不存在的限制：你假设所有输入值都在区间[0, MAX_THREADS_PER_BLOCK]上，并且如果不是这样，你的代码将会失败。这不是原问题中所描述的。如果我把这个代码转化成一个设备函数，并且与其他答案中的朴素二分查找进行基准测试，它的速度会慢两倍，同时每个块所需的共享内存量也会增加一倍。我看不出什么"最佳"的地方。 - talonmies
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接