在GPU上进行广义滑动窗口计算

7
这里有一段Python代码,它在两个三维矩阵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小得多,但在第三维上它们相等。

请注意,我们可以用任何关于Y和窗口的函数替换numpy.dot()。这与卷积略有不同,因为Y仅沿着X的第一和第二维滑动。我正在寻找一种有效的策略,使用CUDA高效地实现这种滑动窗口计算。有人想给我提供一些方向吗?谢谢!

更新:您可以在我的答案中观看我通过其他用户的帮助进行优化过程。

4个回答

4
尝试设计一个“通用”的实现,以适应您可能想要的任何操作,在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默认使用双精度浮点运算。


感謝您的發佈!很抱歉我還沒有機會評估您的解決方案。只是好奇為什麼您沒有選擇基於紋理的實現方法。 - BrianTheLion
仅仅因为我怀疑这样做不会有太大的性能提升。我的基于块的版本已经完全合并了主矩阵和窗口矩阵的读取,这比随机通过纹理读取更快,并且Fermi L1缓存比纹理缓存更大,所以命中率可能同样高。我在其他矩阵操作中的经验表明,绑定到纹理并不更快。 - talonmies

1

好的,这里有一些想法:

您执行了约640*480次numpy.dot迭代,它本身处理了6*6*32个元素。并行化点积几乎不值得:192个并行线程对于GPU来说还不够,并且在CUDA上进行归约会带来额外的麻烦。因此,在我看来,最好的并行化任务的方法是将输出数组的每个元素分配给一个线程。

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

以上只是我的猜测,要实现良好的性能,最好尝试不同的变体...

关于您的天真实现的更新

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

在这里,每次迭代都会访问全局内存。这是一个巨大的性能杀手。由于您有三个维度,最好将您的int *Ydims替换为int3 Ydims(对于Xdimsoutdims也是如此)。
out[out_indx] += X[X_indx]*Y[Y_indx];

再次强调,这是非常糟糕的想法。创建一个寄存器变量并通过它执行所有操作。在内核结束时只需一次性写入全局数组。

这些优化是您应该首先做的事情。其次是将您的 XY 变成 3D 纹理,这样可以缓存对它们的访问。我猜,这样的话 CUDA 将超越 CPU。

为了进一步优化,您最好阅读 CUDA C 最佳实践指南。这是必读的,您会更好地了解如何编写高效的 GPU 代码(现在您的实现过于幼稚)。


谢谢!我尝试了你的建议,将每个输出像素映射到一个单独的线程。目前还没有尝试进行任何内存优化。结果迄今为止是参差不齐的。 - BrianTheLion
哇,太棒了!据我所知,内核参数存储在本地内存中,而本地内存是离芯片的。有没有办法将outdims、Xdims和Ydims存储到芯片内存中? - BrianTheLion
@BrianTheLion 不,内核参数存储在芯片上的共享内存中,这通常几乎与寄存器一样快。您可能会混淆OpenCL'ish本地内存(与CUDA'ish共享相同)和CUDA'ish本地内存,它实际上只是离芯片全局内存的一部分。 - aland
很好。我现在猜测我的v0.2性能是因为我正在使用1D纹理,因此没有得到2D优化缓存的好处。 - BrianTheLion

0

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 - 纹理<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个是冗余的。 - talonmies

0

你可能想尝试将读取、求和和存储分开。

因此,每个内核应该有3个部分:

  1. 从纹理内存中读取,将整个块的数据存储到共享内存中

    __shared blockX[ Ydims.z ][ Ydims.y ][ Ydims.x ];
    __shared blockY[ Ydims.z ][ Ydims.y ][ Ydims.x ];
    // 注意:每个线程加载 k 个元素 * 2,而不是每个线程加载 Ydims.X*Y*Z 个元素
    blockX[k][yj][yi] = ...
    blockY[k][yj][yi] = ...
    __syncthreads(); // <-- 至关重要 -- 所有块中的线程必须完成
    // 从共享内存中读取操作,然后才能使用这些值。
    
  2. #pragma 展开你的 for 循环。
    这将显著提高你的 ILP,并且对于常量循环大小,分支较少。

  3. 确保你的共享内存访问按适当的步长进行,否则银行冲突会影响性能。


谢谢!共享内存优化是我今天早上一直在做的。我们很快就会知道结果。 - BrianTheLion

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接