logo标签列表

在GPU（OpenCL）上，3x3卷积是否应该更快？

cperformanceopenclgpgpuconvolution
7
我正在学习如何为GPU优化代码。我了解到内存局部性的重要性，并看过一些有关GPU卷积的 tutorials 和 examples。基于此，我编写并测试了几个自己的核函数。令人惊讶的是，最简单的 naive kernel 是最快的！？而且它只比CPU快<10倍。（是的，我通过64x运行核心来摊销上传/下载时间）。 我做错了什么？我原本期望卷积是GPU优化的那种操作。如果我可以在矩阵乘法上获得 100倍的加速, 为什么卷积如此缓慢？
  • CPU-naive 9.5
  • GPU-naive 1.64
  • GPU-local 2.56
  • GPU-local_async 15.10
  • GPU-scanline-private 7.35
  • GPU-scanline_async 15.37

编辑：GPU-scanline_async 是我在阅读有关 async_work_group_copy 的建议后制作的。

我想知道两件事：

  • 内核速度是由内存带宽还是计算能力限制的？根据我所读的，我会认为是内存。但测试结果却相反。
    • 内核 GPU-local 比 GPU-naive 更慢，尽管它读取了更少的全局内存
    • 通过高斯滤波器系数修改内核（即每个像素添加乘法）使其变慢了 2 倍以上，尽管它读取的内存数量相同
    • 但如果它受到处理能力的限制，那么为什么我在 GPU 上进行矩阵乘法比在 CPU 上快 100 倍？
  • 为什么内核 GPU-scanline-private 如此缓慢？内存局部性要好得多（每个像素从全局内存中读取的次数仅为 3 而不是 9），逻辑最小（没有 if/switches）

这个测试是在我的笔记本电脑上 进行的，使用了CPU Intel Core i7 6700HQ Skylake和GPU nVidia 960M，运行了每帧64x的浮点数组，大小为256x256像素。代码完整版可以在这里查看。

=========== Kernel codes ===========

kernel GPU-Naive 2D global=(256,256) local=(16,16)

__kernel void blur2D_naive(
    __global float* I, 
    __global float* O
){
    const int ix = get_global_id (0)+1;
    const int iy = get_global_id (1)+1;
    const int nx = get_global_size(0)+2;

    int i = iy * nx + ix;

    // 1.6 ticks/pixel
    O[i] =( I[i-nx-1] + I[i-nx] + I[i-nx+1] +
            I[i   -1] + I[i   ] + I[i   +1] +
            I[i+nx-1] + I[i+nx] + I[i+nx+1] ) * 0.11111111111;
    // modified with gaussian mask 4.9 ticks/pixel
    //O[i] =( 0.0625*I[i-nx-1] + 0.125*I[i-nx] + 0.0625*I[i-nx+1] +
    //        0.125 *I[i   -1] + 0.25 *I[i   ] + 0.125 *I[i   +1] +
    //        0.0625*I[i+nx-1] + 0.125*I[i+nx] + 0.0625*I[i+nx+1] );
}

内核 GPU本地化 2D 全局=(256,256) 本地=(16,16)

#define NBx 18 // tile size including borders [halo] 16+2
#define NBy 18
// seems to be slower than naive method
__kernel void blur2D_local(
    __global float* I, 
    __global float* O
){
    __local float L[NBx*NBy];
    const int2 iG  = (int2)(get_global_id  (0)+1 , get_global_id  (1)+1 );
    const int2 nG  = (int2)(get_global_size(0)+2 , get_global_size(1)+2 );
    const int2 iL  = (int2)(get_local_id   (0)+1 , get_local_id   (1)+1 );
    const int2 nL  = (int2)(get_local_size (0)+2 , get_local_size (1)+2 );
    const int2 iGR = (int2)(get_group_id   (0)   , get_group_id   (1)   );

    // copy boundary pixels to local memory
    switch( get_local_id(1) ){ // some threads copy one more of boundary (halo) pixels
        case 4: 
        switch( get_local_id(0) ){ // copy corner points
            case 0: L[        0      ] = I[ nG.x* get_group_id(1)*get_local_size(1)          + get_group_id(0)*get_local_size(0)         ]; break; // upper-left
            case 1: L[         NBx-1 ] = I[ nG.x* get_group_id(1)*get_local_size(1)          + get_group_id(0)*get_local_size(0)+(NBx-1) ]; break; // upper-right
            case 2: L[ (NBy-1)*NBx   ] = I[ nG.x*(get_group_id(1)*get_local_size(1)+(NBy-1)) + get_group_id(0)*get_local_size(0)         ]; break; // lower-left
            case 3: L[ NBy*    NBx-1 ] = I[ nG.x*(get_group_id(1)*get_local_size(1)+(NBy-1)) + get_group_id(0)*get_local_size(0)+(NBx-1) ]; break; // lower-rigth
        }
        // copy border lines 
        case 0: L[               iL.x    ] = I[ nG.x* get_group_id(1)*get_local_size(1)                   + iG.x                                        ]; break; // top    line
        case 1: L[ NBx*(NBy-1) + iL.x    ] = I[ nG.x*(get_group_id(1)*get_local_size(1)+(NBy-1)         ) + iG.x                                        ]; break; // botton line
        case 2: L[ NBx*iL.x              ] = I[ nG.x*(get_group_id(1)*get_local_size(1)+get_local_id(0) ) +  get_group_id(0)*get_local_size(0)          ]; break; // left   line
        case 3: L[ NBx*iL.x    + (NBx-1) ] = I[ nG.x*(get_group_id(1)*get_local_size(1)+get_local_id(0) ) + (get_group_id(0)*get_local_size(0)+(NBx-1)) ]; break; // right  line
    } // each thread coppied at max. 1 border pixels

    int ig = iG.y*nG.x + iG.x;
    int il = iL.y*nL.x + iL.x;
    L[il] = I[ig];             // each thread copy his pixel to local memory

    barrier(CLK_LOCAL_MEM_FENCE);

    const float renorm = 1.0/9.0;
    O[ig] =( L[il-NBx-1] + L[il-NBx] + L[il-NBx+1] +
             L[il    -1] + L[il    ] + L[il    +1] +
             L[il+NBx-1] + L[il+NBx] + L[il+NBx+1] ) / 9.0;
}

内核 GPU-local_async 二维全局=(256,16) 局部=(16,16)

#define nTiles 16
#define NBx 18
#define NBy 18 
#define copy_tile(event,ig0,I,L) { int ig_=ig0; int il_=0; for(int i=0; i<NBy; i++){   event = async_work_group_copy( L+il_, I+ig_, NBx, event ); ig_+=nx; il_+=NBx; } }
// https://streamcomputing.eu/blog/2014-06-19/using-async_work_group_copy-on-2d-data/
__kernel void blur2D_local_async(
    __global float* I, 
    __global float* O
){
    const int nx = get_global_size(0)+2;        
    __local float LI[NBx*NBy*2];
    int iL0 = 0;
    int iL1 = NBx*NBy;        
    event_t event = 0;
    int ig0 = get_group_id(0)*get_local_size(0);
    copy_tile(event,ig0,I,LI);
    for( int it=0; it<nTiles; it++ ){
        int ig   = ig0 + (get_local_id(1)+1)*nx  + get_local_id(0)+1;
        int il   =       (get_local_id(1)+1)*NBx + get_local_id(0) + iL0;
        ig0     += get_local_size(1)*nx;
        event_t event_ = 0;
        copy_tile(event_,ig0,I,LI+iL1);
        wait_group_events(1, &event);
        //barrier(CLK_LOCAL_MEM_FENCE);
        O[ig] =( LI[il-NBx] + LI[il-NBx+1] + LI[il-NBx+2] +
                 LI[il    ] + LI[il    +1] + LI[il    +2] +
                 LI[il+NBx] + LI[il+NBx+1] + LI[il+NBx+2] ) * 0.11111111111;
        int iLtmp=iL0; iL0=iL1; iL1=iLtmp;
        event = event_;
    }
}

内核 GPU-scanline_private 1D 全局=(256) 局部=(32)

__kernel void blur2D_scanline_priv(
    int nx, int ny,
    __global float* I, 
    __global float* O
){ 
    int ig    = get_global_id(0)+1;
    float3 Lm = (float3)( I[ig-1], I[ig], I[ig+1] );  ig += nx;
    float3 L0 = (float3)( I[ig-1], I[ig], I[ig+1] ); 
    for(int iy=1; iy<(ny-1); iy++ ){
        ig += nx;
        float3 Lp= (float3)( I[ig-1], I[ig], I[ig+1] );  
        O[ig-nx] = 
            ( Lm.x + Lm.y + Lm.z +
              L0.x + L0.y + L0.z +
              Lp.x + Lp.y + Lp.z ) * 0.11111111111;              
        Lm=L0; L0=Lp; 
    }
}

内核 GPU-scanline_async 1D 全局=(256) 局部=(32)

 #define NB 34
__kernel void blur2D_scanline_async(
    int nx, int ny,
    __global float* I, 
    __global float* O
){
    __local float  L[NB*4];
    int i0=0;
    int i1=NB;
    int i2=NB*2;
    int i3=NB*3;
    event_t event = 0;
    int ig0 = get_group_id(0)*get_local_size(0);
    event = async_work_group_copy(  L     , I+ig0, NB, event );    ig0 += nx;
    event = async_work_group_copy(  L+NB  , I+ig0, NB, event );    ig0 += nx;   
    event = async_work_group_copy(  L+NB*2, I+ig0, NB, event );    ig0 += nx;
    const int il = get_local_id(0);
    int ig = get_global_id(0)+1;
    for(int iy=1; iy<(ny-2); iy++ ){
        wait_group_events(1, &event);
        event = async_work_group_copy(  L+i3, I+ig0, NB, event ); ig0 += nx;
        ig += nx;
        O[ig] =  
            ( L[i0+il] + L[i0+il+1] + L[i0+il+2] +
              L[i1+il] + L[i1+il+1] + L[i1+il+2] +
              L[i2+il] + L[i2+il+1] + L[i2+il+2] ) * 0.11111111111;
        __local float *Ltmp;
        int itmp=i0; i0=i1; i1=i2; i2=i3; i3=itmp;
    }
}

内核 CPU原始的

void blur(int nx, int ny, float * I, float * O ){
    float renorm = 1.0/9.0;
    for(int iy=1;iy<ny-1;iy++){ for(int ix=1;ix<nx-1;ix++){
        int i   = iy*nx+ix;
        O[i] =( I[i-nx-1] + I[i-nx] + I[i-nx+1] +
                I[i   -1] + I[i   ] + I[i   +1] +
                I[i+nx-1] + I[i+nx] + I[i+nx+1] ) * renorm;
    } }
}
1个回答
3

在矩阵乘法中，每个子矩阵（补丁）用于另一个矩阵的所有行中的所有补丁。如果补丁中有2x2子矩阵，并且主矩阵为20x20，则每个子矩阵用于乘法10次。GPU通常使用16x16或32x32大小的补丁，这意味着对于2kx2k的乘法，每个16x16补丁至少重复使用128次。

MM reuse = 128

加上子矩阵-子矩阵乘法重用，就足以推动GPU的极限。

在3x3卷积中，不会对整个扫描线或整个图像使用3x3补丁。仅重用其像素。

3x3模板：每个像素被相邻的8个模板重用。

5x5模板：每个像素被相邻的24个模板重用。

为了赶上矩阵乘法，需要

11x11 stencil to have a reuse of 120

这种方法比矩阵乘法更本地化，应该比它获得更多的GFLOPS，但它并没有做相同数量的乘法和加法。

它进行了9次加法+1次乘法。

8个潜在的乘法被丢失了。几乎失去了一半的GFLOPS限制。

您应该尝试异步工作组复制。

  • 加载左上角18x18,
  • 加载顶部18x18并异步计算左上角
  • 加载右上角18x18并异步计算顶部并存储异步左上角
  • 加载右侧18x18并异步计算左上角并存储异步顶部
  • 加载...计算...存储...全部异步，因此可以使用本地内存和主内存(主内存可以利用naive版本，L1也许)

矩阵乘法/具有16x16子矩阵与卷积(17x17笔刷大小)：

  • 矩阵: L2重用率随着主矩阵大小的增加而增加，或者L1重用率随着子矩阵大小的增加而增加（L1）

    • 卷积：所有图像大小的总重用率相同，但是使用刷子尺寸增加L1使用率（好）

  • 矩阵：每个工作组有16 * 16 * 16次乘法+ 16 * 16 * 16次加法

    • 卷积：每个线程有17 * 17次加法+ 1次乘法（不好）

  • 矩阵：均匀的线程使用，无if-else，所有本地内存都被重用

    • 卷积：需要加载至少比边界远16个像素（具有16厚度的幽灵墙），这些像素将由邻居工作组重新使用，但是这些邻居工作组可能在另一个计算单元中，并且只使用L2而不是与其位于同一计算单元以使用L1（丑陋）
      • 这就是为什么我建议异步工作组副本，以在同一计算单元上使用那些邻居（和L1），并增加重用率。

  • 矩阵：增加补丁大小还会以立方体力量速率在子矩阵乘法中增加重用（但由于每行的补丁较少，L2重用率下降，这使得总重用率像平方功率一样减少）

    • 卷积：增加补丁大小按平方功率率增加重用

  • 矩阵：本地内存必须至少是瓦片区域的2倍（子mat-mat mul）

    • 卷积：本地内存必须至少为瓷砖区域+幽灵墙区域

  • 矩阵：可以在私有内存中执行4x4个子子乘法（每个元素使用4次），这意味着4x4内存= 64 add + 64 mul

    • 卷积：将4x4加载到私有内存中只进行4像素计算（对于3x3刷子），这意味着4x4内存= 36 add + 4 mul

拥有一个添加密集型内核可以同时或异步地让另一个乘法密集型内核工作。也许如果您将其用于图像处理，可以添加一些“混合”或“调整大小”的内核，以便它们一起工作？

扫描线版本正在加载3个元素，进行9加法+ 1乘法，然后重复，加载的元素保持3次，这意味着它们仅重用3次，它们的邻居（x或y directio）可能不会落在邻居线程甚至不是邻居工作组中。此外，3个加载与1个存储不平衡。如果内存带宽为100 GB / s，则会使用50GB / s进行加载，15 GB / s进行存储，除非它们来自L1。

您可以通过使用累加器减少加/乘不平衡。

store = (accumulator) * 0.1111111
accumulator+=new vector  // 3 adds
accumulator-=old vecotr  // 3 adds

现在是6个加法和1个乘法，更加平衡了：1Tflops的GPU将有500Gflops用于加法，90Gflops用于乘法。

简单版本不使用本地内存，这样可以在飞行中留出更多波前的空间。本地内存版本实际上会破坏L1访问模式，并减少在飞行中的波前数量。这会降低VALU的占用率。

通过对工作组级别而非线程级别进行扫描线，您可以减少本地内存的使用。我的意思是：

从内存加载： x x x x x x x x x x 对其进行扫描线处理： (从左到右，1-D) a b c d e f g h i j 现在在工作组级别上使用它进行扫描线： a c c u m u l a t o r (+new) (从上到下) z x z x z x z x z x (- old)

calculate frontline 1-d scanline:  30 additions for each new row
calculate wide vector 2-d scanline:30*30 additions
each pixel get 1 value instead of adding 3 values
storing: 16x16 multiplications
much less local memory used, more balanced (~8 add 1 mul)

这个处理方式有1D扫描线，可以选择单线程运行N周期，或者是多线程减少LogN周期（如果在计算单元中拥有足够的线程）。

啊哈，我终于通过阅读这篇文章 https://streamcomputing.eu/blog/2014-06-19/using-async_work_group_copy-on-2d-data/ 理解了你的意思... 我之前不知道这个函数。太棒了！谢谢你的提示！ - Prokop Hapala
@ProkopHapala 我检查了你的内核，似乎每个线程都在发出不同的异步复制。关于异步复制的主要规则是所有线程必须给它相同的参数，并且所有参数都必须命中它。我认为它们做的工作较少。我的意思是使用单个异步复制（可能是其分步版本）将整个本地与内部+幽灵加载。最后，异步命令本身不是障碍，如果内存在一些数据写入后立即使用，则内核可能需要在其之前设置障碍。最后，异步复制也可以从本地到全局，这可以与其他内容异步复制。 - huseyin tugrul buyukisik
地址L+i3-il实际上对于所有线程都是相同的。现在，我修改了内核代码以使其更明显（以防编译器出现混淆）（您可以查看更新版本）。索引i0、i1、i2、i3仅取决于get_group_id和常量。我希望洗牌int itmp=i0; i0=i1; i1=i2; i2=i3; i3=itmp;不会让编译器混淆以识别地址是否相同。至于其他评论-是的，我理解你关于18x18瓷砖的想法，但我想从更小/简单的模式开始。我仍然想知道为什么它比以前版本的scanline慢。 - Prokop Hapala
嗨，我实现了blur2D_local_async，它复制18x18个瓷砖，但仍然比朴素内核慢10倍：-(...我不明白为什么。 async_work_group_copy似乎比手动全局到本地复制（blur2D_local）慢得多（> 5倍）。 - Prokop Hapala
实现是否强制它们串行化？你使用的GPU是什么？ - huseyin tugrul buyukisik
显示剩余11条评论
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接