在CUDA中进行1D Min-convolution

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])

(编辑：循环从 1 改为从 0 开始)

如果最小值改成求和，我们可以使用快速傅里叶变换（FFT），但在最小值情况下，不存在这样的类比。 因此，我想通过使用 GPU（CUDA）来尽可能地加快这个简单算法的速度。 我很乐意找到已经存在的代码来实现此功能（或者实现求和情况而不使用 FFT，以便我可以根据自己的需求进行适当的修改），但目前为止我的搜索没有得到什么好结果。 我的用例将涉及大小在 1,000 到 100,000 之间的 a 和 b。

问题：

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

• 如果我要自己实现，从结构上讲，CUDA 核心应该如何设计才能最大程度地提高效率？ 我尝试了一种简单的解决方案，即每个线程计算一个 c[i]，但这似乎不是最佳方法。 在线程块结构和内存访问模式方面有哪些建议？

- dan_x

Answer 1

更快的版本:

__global__ void convAgB(double *a, double *b, double *c, int sa, int sb)
{
    int i = (threadIdx.x + blockIdx.x * blockDim.x);
    int idT = threadIdx.x;
    int out,j;

    __shared__ double c_local [512];

    c_local[idT] = c[i];

    out = (i > sa) ? sa : i + 1;
    j   = (i > sb) ? i - sb + 1 : 1;

    for(; j < out; j++)
    {    
       if(c_local[idT] > a[j] + b[i-j])
          c_local[idT] = a[j] + b[i-j]; 
    }   

    c[i] = c_local[idT];
} 

**Benckmark:**
Size A Size B Size C Time (s)
1000   1000   2000   0.0008
10k    10k    20k    0.0051
100k   100k   200k   0.3436
1M     1M     1M     43,327

旧版本， 对于1000到100000之间的尺寸，我使用这个朴素版本进行了测试：

__global__ void convAgB(double *a, double *b, double *c, int sa, int sb)
{
    int size = sa+sb;

    int idT = (threadIdx.x + blockIdx.x * blockDim.x);
    int out,j;


    for(int i = idT; i < size; i += blockDim.x * gridDim.x)
    {
        if(i > sa) out = sa;
        else out = i + 1;

        if(i > sb) j = i - sb + 1;
        else j = 1;


        for(; j < out; j++)
        {
                if(c[i] > a[j] + b[i-j])
                    c[i] = a[j] + b[i-j];
        }
    }
}

我用一些随机的双精度数字填充了数组a和b，并用999999填充了c（仅用于测试）。使用您的函数（未进行任何修改）验证了c数组（在CPU中）。

我还从内部循环中删除了条件语句，因此它将只测试一次。

我不是100％确定，但我认为以下修改是有道理的。由于你有 i - j >= 0，这与 i >= j相同，这意味着一旦j > i，就永远不会进入这个块'X'（因为j ++）：

if(c[i] > a[j] + b[i-j])
   c[i] = a[j] + b[i-j];

在变量out上计算循环条件，如果i > sa，则意味着当j == sa时循环将结束；如果i < sa，则意味着由于条件i >= j的存在，循环将在i + 1时（更早地）结束。

另一个条件i - j < size(b)意味着当i > size(b) + 1时，将开始执行块“X”，因为j始终等于1。因此，我们可以将j设置为应该开始的值，如下：

if(i > sb) j = i - sb + 1;
else j = 1;

请尝试使用真实的数据数组测试此版本，并给我反馈。欢迎任何改进意见。

编辑：可以实现一项新的优化，但这并没有什么大的区别。

if(c[i] > a[j] + b[i-j])
    c[i] = a[j] + b[i-j];

我们可以通过以下方式消除if语句：

double add;
...

 for(; j < out; j++)
 {
   add = a[j] + b[i-j];
   c[i] = (c[i] < add) * c[i] + (add <= c[i]) * add;
 }

需要：

if(a > b) c = b; 
else c = a; 

这段代码的意思是 c = (a < b) * a + (b <= a) * b。

如果 a > b，那么 c = 0 * a + 1 * b；=> c = b； 如果 a <= b，那么 c = 1*a + 0 *b； => c = a；

**Benckmark:**
Size A Size B Size C Time (s)
1000   1000   2000   0.0013
10k    10k    20k    0.0051
100k   100k   200k   0.4436
1M     1M     1M     47,327

我正在测量从CPU复制到GPU、运行内核和从GPU复制到CPU所需的时间。

GPU Specifications   
Device                       Tesla C2050
CUDA Capability Major/Minor  2.0
Global Memory                2687 MB
Cores                        448 CUDA Cores
Warp size                    32

- dreamcrash

Answer 2

一个可能对大的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 *);

< p > 【免责声明：未经测试或基准测试，使用需谨慎】 < p > 这是单精度浮点数，但同样的方法也适用于双精度浮点数。对于整数，您需要将 C99 的 INFINITY 宏替换为类似 INT_MAX 或 LONG_MAX 的内容，但原则上仍然相同。

- talonmies

---

谢谢！这比我的天真基准在大小为（1000,1000）的问题的1000倍迭代中快了4倍。基准：852.6毫秒；这个：225.3毫秒。 - dan_x

@dan_x 抱歉问了这么多问题，你用的是什么类型的GPU？ - Robert Crovella

@RobertCrovella 很高兴回答。我可以使用几个，但现在主要在 Tesla C1060 上运行。 - dan_x

Answer 3

我使用了你的算法。我认为它会对你有帮助。

const int Length=1000;

__global__ void OneD(float *Ad,float *Bd,float *Cd){
    int i=blockIdx.x;
    int j=threadIdx.x;
    Cd[i]=99999.99;
    for(int k=0;k<Length/500;k++){
        while(((i-j)>=0)&&(i-j<Length)&&Cd[i+k*Length]>Ad[j+k*Length]+Bd[i-j]){
            Cd[i+k*Length]=Ad[j+k*Length]+Bd[i-j];
    }}}

我已经使用每个块500线程，每个网格500块。由于我的设备每个块的线程数限制为512，因此我使用了500线程。我将所有数组的大小都设置为Length(=1000)。

工作原理：

1. i存储块索引，j存储线程索引。 2. 使用for循环，因为线程数少于数组大小。 3. 使用while循环迭代Cd[n]。 4. 我没有使用共享内存，因为我使用了大量块和线程，所以每个块所需的共享内存量较低。

注：如果您的设备支持更多线程和块，请将k<Length/500替换为k<Length/(支持的线程数)。

- Fr34K