arrays - CUDA中的一维最小卷积

Question

我有两个数组，a和b，我想计算“最小卷积”以产生结果c。简单的伪代码如下所示:

for i = 0 to size(a)+size(b)
    c[i] = inf
    for j = 0 to size(a)
        if (i - j >= 0) and (i - j < size(b))
            c[i] = min(c[i], a[j] + b[i-j])

(编辑:更改循环从0而不是1开始)

如果最小值是总和，则可以使用快速傅立叶变换(FFT)，但在最小值情况下，没有这样的模拟。相反，我想通过使用GPU(CUDA)来使这种简单的算法尽可能快。我很乐意找到执行此操作的现有代码(或实现不使用FFT的求和情况的代码，以便我可以根据自己的目的进行调整)，但是到目前为止，我的搜索没有得到任何好的结果。我的用例将涉及大小在1,000到100,000之间的a和b。

问题:

是否已经存在有效执行此操作的代码？

如果我要在结构上自己实现这一点，那么CUDA内核应如何看待以最大程度地提高效率？我尝试了一个简单的解决方案，其中每个c [i]由一个单独的线程计算，但这似乎不是最好的方法。关于如何设置线程块结构和内存访问模式的任何技巧？

Answer 1

对于较大的a和b可能有用的替代方法是对c中的每个输出条目使用一个块。使用块可以实现内存合并，这对于限制内存带宽的操作非常重要，并且可以使用相当有效的共享内存减少来将每个线程的部分结果组合为最终的每个块结果。最好的策略可能是每个MP启动与同时运行的块一样多的块，并使每个块发出多个输出点。这消除了与以相对较低的总指令数启动和退出许多块相关的一些调度开销。

有关如何完成此操作的示例:

#include <math.h>

template<int bsz>
__global__ __launch_bounds__(512)
void minconv(const float *a, int sizea, const float *b, int sizeb, float *c)
{
    __shared__ volatile float buff[bsz];
    for(int i = blockIdx.x; i<(sizea + sizeb); i+=(gridDim.x*blockDim.x)) {
        float cval = INFINITY;
        for(int j=threadIdx.x; j<sizea; j+= blockDim.x) {
            int t = i - j;
            if ((t>=0) && (t<sizeb))
                cval = min(cval, a[j] + b[t]);
        }
        buff[threadIdx.x] = cval; __syncthreads();
        if (bsz > 256) {
            if (threadIdx.x < 256) 
                buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+256]);
            __syncthreads();
        }
        if (bsz > 128) {
            if (threadIdx.x < 128) 
                buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+128]); 
            __syncthreads();
        }
        if (bsz > 64) {
            if (threadIdx.x < 64) 
                buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+64]);
            __syncthreads();
        }
        if (threadIdx.x < 32) {
            buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+32]);
            buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+16]);
            buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+8]);
            buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+4]);
            buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+2]);
            buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+1]);
            if (threadIdx.x == 0) c[i] = buff[0];
        }
    }
}

// Instances for all valid block sizes.
template __global__ void minconv<64>(const float *, int, const float *, int, float *);
template __global__ void minconv<128>(const float *, int, const float *, int, float *);
template __global__ void minconv<256>(const float *, int, const float *, int, float *);
template __global__ void minconv<512>(const float *, int, const float *, int, float *);

[免责声明:未经测试或基准测试，使用风险自负]

这是单精度浮点，但相同的想法也适用于 double 浮点。对于整数，您需要将C99 INFINITY宏替换为 INT_MAX或 LONG_MAX，但是其原理保持不变。