【问题标题】:Shouldn't be 3x3 convolution much faster on GPU (OpenCL)GPU 上的 3x3 卷积不应该快得多(OpenCL)
【发布时间】:2017-04-04 20:49:57
【问题描述】:

我正在学习如何为 GPU 优化代码。我读到了内存局部性的重要性。我还看到了一些 GPU 卷积的 tutorialsexamples。基于此,我编写并测试了几个自己的内核。令人惊讶的是,我发现最简单的 naive kernell 是最快的!? 它比 CPU 快

我做错了什么?我希望卷积只是 GPU 优化的那种操作。如果能得到100x speed-up on matrix multiplication,为什么卷积这么慢?

性能 [CPU ticks/pixel](越低越好):
  • CPU-naive 9.5
  • GPU-naive 1.64
  • GPU 本地 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-localGPU-naive 慢,尽管它执行的全局内存读取要少得多
    • 通过高斯滤波器系数修改内核(即为每个像素添加乘法)使其速度变慢了 2 倍,尽管它执行相同数量的内存读取
    • 但如果它受到处理能力的限制,为什么我在 GPU 上的矩阵乘法比在 CPU 上快 100 倍?
  • 为什么内核 GPU-scanline-private 这么慢? 内存局部性要好得多(每个像素从全局内存读取 3 次而不是 9 次)并且逻辑极少(没有 ifs/开关)

测试是在my laptop 上使用 CPU Intel Core i7 6700HQ Skylake 和 GPU nVidia 960M 在 256x256 像素的浮点阵列上以 64x/帧运行内核完成的。 code full can be seen here

===========内核代码===========

内核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-local 2D global=(256,256) local=(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 2D global=(256,16) local=(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 global=(256) local=(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 global=(256) local=(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-naive

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;
    } }
}

【问题讨论】:

    标签: c performance opencl gpgpu convolution


    【解决方案1】:

    在矩阵乘法中,每个子矩阵(补丁)用于另一个矩阵中所有行中的所有补丁。如果补丁中有 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 并计算左上异步并存储顶部异步
    • 加载 .... 计算 ... 存储... 全部异步,因此可以使用本地内存和主内存(主内存将利用幼稚版本,L1 可能)

    矩阵乘法/使用 16x16 子矩阵)与卷积(17x17 画笔大小):

    • 矩阵:L2复用率随着主矩阵大小增加,或者L1复用率随着子矩阵大小(L1)增加

      • 卷积:所有图像大小的总重用率相同,但 L1 使用率随画笔大小增加(良好)
    • 矩阵:每个工作组 16*16*16 次乘法 + 16*16*16 次加法

      • 卷积:17*17 加法 + 每个线程 1 次乘法(错误)
    • 矩阵:统一线程使用,没有if-else,所有本地内存都重用

      • 卷积:需要比边界(厚度为 16 的鬼墙)加载至少 16 个像素,这些边界将由相邻工作组重新使用,但这些相邻工作组可能在另一个计算单元中,并且只使用 L2 而不是打开使用 L1 的相同计算单元(丑陋)
        • 这就是为什么我建议异步工作组副本在同一计算单元(和 L1)上使用这些邻居并提高重用率。
    • 矩阵:增加补丁大小也会增加子矩阵乘法中三次功率率的重用(但由于每行的补丁较少,因此会减少 L2 重用,这使得总重用类似于平方功率率)

      • 卷积:增加补丁大小会增加平方功率率的重用
    • 矩阵:本地内存必须至少是 2x tile 区域(sub mat-mat mul)

      • 卷积:本地内存必须至少是瓦片区域 + 鬼墙区域
    • 矩阵:可以在私有内存中进行 4x4 子乘法(每个元素使用 4 次),这意味着 4x4 内存 = 64 add+64 mul

      • 卷积:将 4x4 加载到私有内存中,除了 4 像素计算(对于 3x3 画笔)没有任何作用,这意味着 4x4 内存 = 36 add + 4 mul

    拥有一个重加法的内核为另一个重乘法内核同时工作或在同一个内核中异步工作留下了空间。也许如果您将其用于图像处理,也许您可​​以在其中添加一些“混合”或“调整大小”内核,以便它们一起工作?


    扫描线版本正在加载 3 个元素,执行 9 add + 1 mul 然后重复,加载的元素停留 3 圈,这意味着它们仅被重复使用 3 次,并且其邻居(x 或 y 方向)可能不会落在邻居中线程甚至邻居工作组。 3 个负载与 1 个存储也是不平衡的。如果内存带宽为 100 GB/s,那么它将使用 50GB/s 进行加载,15 GB/s 用于存储,除非它们来自 L1。

    您可以使用累加器来减少 add/mul 不平衡。

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

    所以现在它是 6 个加法 + 1 个 muls,所以更加平衡:1Tflops GPU 将有 500Gflops 加法,90 Gflops 用于 muls。


    朴素版本不使用本地内存,为更多波前留出了更多空间。本地内存版本实际上打破了 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(-旧)

    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)
    

    这具有一维扫描线,对于 N 个周期是单线程的,对于 LogN 周期是多线程的减少(考虑到计算单元中有足够多的线程)。

    【讨论】:

    • 谢谢,您似乎知识渊博。我必须读几遍才能理解。我不清楚的是:1)您对 异步工作组副本 的建议是否适用于 3x3 卷积,还是适用于更大的模板? 2)您说我的内核很慢,因为它们“不平衡”并且可重用性低(乘法很少,读取很多,加法)。那么我希望像在朴素内核的 gaussian mask mod 中添加乘法不会影响速度?
    • 一般来说,我正在测试这个 3x3 卷积,就像偏微分方程的有限差分求解器等算法的简单模型一样。内存访问模式是相同的——每个像素都需要知道其邻居的值。但是操作会有所不同 - 不同的模板,(带有模板掩码的点积,最少的邻居)。在那种情况下,我猜您对 accumulator 的建议不适用(?)
    • (第一条评论)异步工作组可能适用于更大的块,例如您提到的 18x18。是的,当 mul 单元空闲时,使用它们会填充管道并隐藏。 (第二条评论)如果不仅有加法,那么累加器不可用,但您仍然可以使用多个减法。
    • 啊哈通过阅读streamcomputing.eu/blog/2014-06-19/…,我终于明白了你的意思......我不知道这个功能。很好的提示!谢谢
    • 嗨。我试过 async_work_group_copy (参见 EDIT kernel GPU-scanline_async )。但这并没有帮助,实际上它比以前的 GPU-scanline-private 慢了近 3 倍。你能对此发表评论吗?我希望我不会打扰你太多。
    猜你喜欢
    • 2016-12-31
    • 1970-01-01
    • 1970-01-01
    • 2013-10-21
    • 1970-01-01
    • 2013-10-28
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多