CUDA 中的一维最小卷积
Posted
技术标签:
【中文标题】CUDA 中的一维最小卷积【英文标题】：1D Min-convolution in CUDA 【发布时间】：2012-10-21 01:39:15 【问题描述】：我有两个数组，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）
如果 min 是一个和，我们可以使用快速傅里叶变换 (FFT)，但在 min 的情况下，没有这样的模拟。相反，我想通过使用 GPU (CUDA) 使这个简单的算法尽可能快。我很乐意找到执行此操作的现有代码（或在没有 FFT 的情况下实现 sum case 的代码，以便我可以根据我的目的对其进行调整），但到目前为止我的搜索还没有找到任何好的结果。我的用例将涉及大小在 1,000 到 100,000 之间的 a 和 b。
问题：
是否已经存在有效执行此操作的代码？
如果我要自己实现这个，在结构上，CUDA 内核应该如何看起来才能最大限度地提高效率？我尝试了一个简单的解决方案，其中每个 c[i] 都由一个单独的线程计算，但这似乎不是最好的方法。关于如何设置线程块结构和内存访问模式的任何提示？
【问题讨论】：
【参考方案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 *);
[免责声明：未经测试或基准测试，使用风险自负]
这是单精度浮点，但同样的想法应该适用于双精度浮点。对于整数，您需要将 C99 的 INFINITY 宏替换为 INT_MAX 或 LONG_MAX 之类的东西，但原理保持不变。
【讨论】：
谢谢！通过大小问题（1000,1000），这比我在 1000 倍上的幼稚基线快约 4 倍。基线：852.6 毫秒；这：225.3 毫秒 @dan_x 抱歉所有问题，您在哪种 GPU 上运行？ @RobertCrovella 很乐意回答。我可以访问一些，但现在我主要在 Tesla C1060 上运行。【参考方案2】：更快的版本：
__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 上计算了循环条件 if i > sa，这意味着循环将在j == sa 时完成，如果i < sa 这意味着循环将在i + 1 上完成（更早），因为条件i >= j。
另一个条件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
如果 a > b 那么 c = 0 * a + 1 * b; => c = b; 如果 a 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
【讨论】：
【参考方案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 Threads per 块。而且，500 阻止 per 网格。因为，我的设备中 per 块的线程数限制为 512，我使用了 500 线程。我已将所有数组的大小设为Length (=1000)。
工作：
i 存储块索引，j 存储线程索引。
for 循环用于线程数小于数组大小。
while 循环用于迭代Cd[n]。
我没有使用共享内存，因为我占用了很多块和线程。因此，每个块所需的共享内存量很低。
PS：如果您的设备支持更多线程和块，请将k<Length/500 替换为k<Length/(supported number of threads)
【讨论】：
以上是关于CUDA 中的一维最小卷积的主要内容，如果未能解决你的问题，请参考以下文章
使用 CUDA Thrust 确定每个矩阵列中的最小元素及其位置