GPU 上的广义滑动窗口计算

Posted 2023-02-14

【中文标题】GPU 上的广义滑动窗口计算

【英文标题】：Generalized sliding-window computation on the GPU

【发布时间】：2011-12-01 03:50:21

【问题描述】：

这里有一些 Python 代码在两个 3D 矩阵 X 和 Y 上实现滑动窗口计算。

import numpy

def sliding_dot( X,Y ) :

    assert X.ndim == Y.ndim == 3
    iw,ih,id = X.shape
    fw,fh,fd = Y.shape

    assert id == fd
    assert fw < iw and fh < ih

    ow,oh = iw-fw+1,ih-fh+1
    out = numpy.zeros( [ow,oh] )

    for x in xrange(ow) :
        for y in xrange(oh) :
            window = X[x:x+fw,y:y+fh,:]
            out[x,y] = numpy.dot( window.flatten(),Y.flatten() )

    return out

#################    

A_dims = (640,480,32)
B_dims = (6,6,32)

A = numpy.random.rand(*A_dims)
B = numpy.random.rand(*B_dims)

sliding_dot(A,B)

一般来说，Y 在第一维和第二维上总是比 X 小很多，但在第三维上它们是相等的。

请注意，我们可以将 numpy.dot() 替换为 Y 和窗口的任何函数。这与卷积有点不同，因为 Y 仅沿 X 的第一和第二维度滑动。我正在寻找一种有效的策略来使用 CUDA 有效地实现这种滑动窗口计算。有人想给我一些方向吗？干杯！

更新：您可以在下面的回答中观看我在其他用户的帮助下完成优化过程。

【参考方案1】：

在像 CUDA 这样的架构中，尝试设计一个“通用”实现来适应您可能想要的任何操作将是一个巨大的权衡。对于您的具体点积示例，这是一个典型的归约操作，这是一个非常有用的实现：

__constant__ int ldaX[3];
__constant__ int ldaY[3];
__constant__ int dimX[3];
__constant__ int dimY[3];

template<typename real,int blocksize>
__global__ void sliding_k(const real *X, const real *Y, real *out)

    __shared__ volatile real buffer[blocksize];

    int tid = threadIdx.x;
    int gid = blockIdx.x * gridDim.y + blockIdx.y;

    real value = (real)0;
    int xpos = (blockIdx.y * ldaX[2]) + (blockIdx.x * ldaX[1]);
    int ypos = 0;
    for(int i=0; i<dimY[0]; i++) 
        for(int jk=tid; jk<ldaY[1]; jk+=blocksize) 
            value += X[xpos+jk] * Y[ypos+jk];
        
        xpos += ldaX[1];
        ypos += ldaY[1];
    

    buffer[tid] = value;
    __syncthreads();

# pragma unroll
    for(int i=(tid+32); ((tid<32)&&(i<blocksize)); i+=32)
        buffer[tid] += buffer[i];

    if (tid < 16) buffer[tid] += buffer[tid + 16];
    if (tid < 8)  buffer[tid] += buffer[tid + 8];
    if (tid < 4)  buffer[tid] += buffer[tid + 4];
    if (tid < 2)  buffer[tid] += buffer[tid + 2];
    if (tid == 0) out[gid] = buffer[0] + buffer[1];

您可以将任何类型的归约运算符替换为点积使用的浮点乘加/求和运算，并且代码应该可以正常工作。每个窗口计算由单个块执行。有足够的并行工作来证明在这个窗口大小下每个窗口一个块是合理的。这允许合并的全局内存访问，并且在 Fermi 卡上，大量的 L1 缓存命中。

这里我只在代码中建立了一个假设，即源数组和窗口数组的第三维是相等的。这允许内部两个循环“融合”成一个操作，因为它们共享公共内存布局。使用改进版本的参考代码在 Python 中运行测试工具，主机代码用 PyCUDA 编写，我得到了：

In [15]: %timeit -n3 -r3 out2=sliding_cuda(A,B)
3 loops, best of 3: 49.8 ms per loop

In [16]: %timeit -n3 -r3 out=sliding_dot(A,B)
3 loops, best of 3: 2.18 s per loop

In [17]: (numpy.abs(out2-out)/numpy.abs(out)).max()
Out[17]: 4.2921323635558404e-15

在 3GHz Phenom II 和 GTX470 上运行时，在 635x475 2D 网格上使用 64 个线程块 - 即。使用可分页的主机内存分配，包括模块加载、设置和内存传输，速度提高了大约 50 倍。在不包括内存传输和设置开销的情况下，内核本身比 Python 快大约 100 倍。请注意，这是一个双精度版本 - Python 默认使用双精度浮点运算。

【讨论】：

感谢发帖！抱歉，我还没有机会评估您的解决方案。只是好奇你为什么不使用基于纹理的实现。

只是因为我怀疑这样做会有很大的性能改进。我的基于块的版本完全合并了主矩阵和窗口矩阵的读取，这比通过纹理随机读取要快，而且 Fermi L1 缓存大于纹理缓存，因此命中率可能同样高。我对其他矩阵运算的经验表明绑定到纹理并没有更快。

【参考方案2】：

好吧，这里有一些想法：

您对 numpy.dot 执行 ~640*480 次迭代，它本身处理 6*6*32 个元素。并行化点积几乎不值得：192 个并行线程对于 GPU 来说是不够的，减少 CUDA 是额外的麻烦。因此，IMO，并行化任务的最佳方法是将输出数组的一个元素分配给每个线程。

现在关于内存：输出数组将在全局内存中，没有太多选择。对于输入数据，A 看起来非常适合纹理内存，因为相邻线程访问相邻元素。或者，您可以手动将其“缓存”在共享内存中，但在这种情况下，与简单地使用纹理相比，它看起来并没有多大优势。对于B，共享内存不好，因为它会导致bank冲突，因为当你计算点积时，half-warp中的所有线程都访问同一个B的元素（你可以从不同线程中的不同元素开始求和，但是那（再次）看起来并不乐观）。所以选择是纹理还是常数。我投票支持常量，因为 (a) 常量内存适合设备上所有线程访问的数据，(b) 你不会污染纹理缓存。

以上只是我的猜测，要真正获得良好的性能，您最好尝试不同的变体......

关于您的幼稚实现的更新

for (int Yi = 0; Yi < Ydims[0]; Yi++ )

在这里，您可以在每次迭代时访问全局内存。这是一个巨大的性能杀手。由于您有 3 个维度，因此您最好将 int *Ydims 替换为 int3 Ydims（Xdims 和 outdims 相同）。

out[out_indx] += X[X_indx]*Y[Y_indx];

再一次，一个非常糟糕的主意。创建一个寄存器变量并使用它执行所有操作。在内核结束时只写入一次全局数组。

这些优化是您应该做的第一件事。第二件事是为您制作X 和Y 3D 纹理，因此对它们的访问将被缓存。我想，在这之后 CUDA 会胜过 CPU。

如需进一步优化，您最好阅读CUDA C Best Practices Guide。必读，你会更好地了解如何编写高效的 GPU 代码（现在你的实现太天真了）

【讨论】：

谢谢！尝试了您的建议并将每个输出像素映射到单个线程。没有尝试做任何内存优化。到目前为止，结果喜忧参半。

哇，很棒的帮助！据我所知，内核参数存储在本地内存中，而本地内存在芯片外。有什么办法可以让 outdims、Xdims 和 Ydims 存储到片上内存？

@BrianTheLion 不，内核参数存储在片上共享内存中，通常几乎与寄存器一样快。您可能会混淆 OpenCL 的本地内存，它与 CUDA 的共享相同，而 CUDA 的本地内存实际上只是片外全局内存的一部分。

酷。我现在猜测我的 v0.2 性能是由于我使用 1D 纹理，因此没有获得 2D 优化缓存的好处。

【参考方案3】：

v0.1 - 幼稚的实现

这是我第一次天真地尝试完成这项工作：

__global__ void sliding_dot(float *out, int *outdims, float *X, int *Xdims, float *Y, int *Ydims )

    int i = threadIdx.x + blockDim.x * blockIdx.x;
    int j = threadIdx.y + blockDim.y * blockIdx.y;
    int Y_indx = 0;
    int X_indx = 0;
    if ( i < outdims[0] & j < outdims[1] )
    
        int out_indx = j + i*outdims[1];
        for (int Yi = 0; Yi < Ydims[0]; Yi++ )
        
            for (int Yj = 0; Yj < Ydims[1]; Yj++ )
            
                for (int k = 0; k < Ydims[2]; k++ )
                
                    Y_indx = k + Yj*    Ydims[2] + Yi*    Ydims[2]*Ydims[1];
                    X_indx = k + (j+Yj)*Xdims[2] + (i+Yi)*Xdims[2]*Xdims[1];
                    out[out_indx] += X[X_indx]*Y[Y_indx];
                
            
        
    

到目前为止，结果并不理想。选择块大小 (32,32,1) 和网格尺寸 p,q 使得 p*32 >= outdims[0] 和 q*32 >= outdims[1] ：

method=[ sliding_dot ] gputime=[ 7013.280 ] cputime=[ 18.000 ] occupancy=[ 0.667 ] 
method=[ sliding_dot ] gputime=[ 6945.184 ] cputime=[ 7.000 ] occupancy=[ 0.667 ] 
method=[ sliding_dot ] gputime=[ 6990.816 ] cputime=[ 6.000 ] occupancy=[ 0.667 ] 
method=[ sliding_dot ] gputime=[ 6931.648 ] cputime=[ 6.000 ] occupancy=[ 0.667 ] 

---

v0.2 - texture<float,1>

我希望每个人都能像我一样从中学到很多东西！我遵循@aland 的建议并获得了相当大的加速：

texture<float,1> X;
texture<float,1> Y;

__global__ void dotconv(float *out, int2 outdims, int3 Xdims, int3 Ydims )

    int i = threadIdx.x + blockDim.x * blockIdx.x;
    int j = threadIdx.y + blockDim.y * blockIdx.y;

    if ( i < outdims.x & j < outdims.y )
    
        int out_indx = j + i*outdims.y;
        float total = 0.0f;
        int X_indx = 0;
        int Y_indx = 0;
        for (int Yi=0; Yi<Ydims.x; Yi++ )
        
            for (int Yj=0; Yj<Ydims.y; Yj++ )
            
                for (int k=0; k<Ydims.z; k++ )
                
                    Y_indx = k + Yj*    Ydims.z + Yi*    Ydims.z*Ydims.y;
                    X_indx = k + (j+Yj)*Xdims.z + (i+Yi)*Xdims.z*Xdims.y;
                    total += tex1Dfetch(X,X_indx)*tex1Dfetch(Y,Y_indx);
                
            
        
        out[out_indx] = total;
    

但我们的运行速度仍然没有 CPU 快：

method=[ dotconv ] gputime=[ 2224.928 ] cputime=[ 24.000 ] occupancy=[ 0.667 ] 
method=[ dotconv ] gputime=[ 2222.592 ] cputime=[ 7.000 ] occupancy=[ 0.667 ] 
method=[ dotconv ] gputime=[ 2225.216 ] cputime=[ 10.000 ] occupancy=[ 0.667 ] 
method=[ dotconv ] gputime=[ 2222.752 ] cputime=[ 10.000 ] occupancy=[ 0.667 ] 

---

v0.3 - texture<float,3>

texture<float,3,cudaReadModeElementType> X;
texture<float,3,cudaReadModeElementType> Y;

__global__ void dotconv(float *out, int2 outdims, int3 Xdims, int3 Ydims )

    int i = threadIdx.x + blockDim.x * blockIdx.x;
    int j = threadIdx.y + blockDim.y * blockIdx.y;
    if ( i < outdims.x & j < outdims.y )
    
        int out_indx = j + i*outdims.y;
        float total = 0.0f;
        for (int Yi=0; Yi<Ydims.x; Yi++ )
        
            for (int Yj=0; Yj<Ydims.y; Yj++ )
            
                for (int k=0; k<Ydims.z; k++ )
                
                    total += tex3D(X,k,j+Yj,i+Yi) * tex3D(Y,k,Yj,Yi);   
                
            
        
        out[out_indx] = total;
    

这实际上比 v0.2 慢了一点

method=[ dotconv ] gputime=[ 2403.360 ] cputime=[ 35.000 ] occupancy=[ 0.667 ] 
method=[ dotconv ] gputime=[ 2392.160 ] cputime=[ 15.000 ] occupancy=[ 0.667 ] 
method=[ dotconv ] gputime=[ 2396.448 ] cputime=[ 15.000 ] occupancy=[ 0.667 ] 
method=[ dotconv ] gputime=[ 2398.880 ] cputime=[ 16.000 ] occupancy=[ 0.667 ] 

感谢您的建议！

【讨论】：

在您最快的 v0.2 版本中有很多“低悬的果实”。您当前正在对点积内部循环中的每个 fmad 执行 14 整数运算。这是一个巨大的开销，14 个 iops 中至少有 12 个是冗余的。

【参考方案4】：

您可能想尝试将您的读数与您存储的总和分开。

所以每个内核应该有 3 个部分：

从纹理内存读取，存储到整个块的共享内存

__shared blockX[ Ydims.z ][ Ydims.y ][ Ydims.x ];
__shared blockY[ Ydims.z ][ Ydims.y ][ Ydims.x ];
// NOTE: MAKE EACH THREAD LOAD k ELEMENTs * 2 rather than each thread loading Ydims.X*Y*Z elements
blockX[k][yj][yi] = ...
blockY[k][yj][yi] = ...
__syncthreads(); // <-- critical -- all threads in block must finish
// reading from shared memory before any may use the values.

#pragma 展开您的 for 循环。 这将显着增加您的 ILP，并为您的恒定循环大小减少分支

确保您的共享内存访问适当跨步，否则银行冲突会影响您的性能。

【讨论】：

谢谢！共享内存优化是我今天早上一直在做的事情。我们应该很快就会知道这里的结果。