CUDA中的1D Min-convolution

Question

我有两个数组,a和b,我想计算"min convolution"来产生结果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)

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

问题:

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

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

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使用了一些随机的双数字和c999999(仅用于测试).我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了循环条件if i > sa,这意味着循环将在完成时j == sa,如果i < sa这意味着循环将i + 1由于条件而完成(之前)i >= j.

另一个条件i - j < size(b)意味着你将开始执行块'X',i > size(b) + 1因为start 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

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 *);

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

这是单精度浮点,但同样的想法应该适用于双精度浮点.对于整数,你就需要更换C99 INFINITY的东西,如宏INT_MAX或LONG_MAX,但原理是一样的,否则.