小贝子编程

c语言 - 在 GPU (OpenCL) 上不应该更快地进行 3x3 卷积



我正在学习如何优化GPU的代码。我读到了记忆局部性的重要性。我还看过一些 GPU 卷积的教程和示例。在此基础上，我编写并测试了几个自己的内核。令人惊讶的是，我发现最简单的朴素内核是最快的！？它比 CPU 快 <10 倍。(是的，我通过运行狗窝 64x 来摊销上传/下载时间)。

我做错了什么？我希望卷积只是针对GPU进行优化的操作。如果我可以在矩阵乘法上获得 100 倍的速度，为什么卷积这么慢？

性能 [CPU 时钟周期/像素](越低越好)：
  • CPU 幼稚9.5
  • GPU 天真1.64
  • 本地 GPU2.56
  • GPU local_async15.10
  • GPU 扫描线专用7.35
  • GPU scanline_async15.37

编辑：GPU-scanline_async我在阅读有关async_work_group_copy的建议后稍后做出

我想知道两件事：

  • 内核速度是否受内存带宽或计算能力的限制？从我读到的内容中，我期待记忆。但测试结果显示情况恰恰相反。
    • 内核 GPU本地比GPU 朴素慢，即使它执行的全局内存读取要少得多
    • 通过高斯滤波器系数修改内核(即每个像素添加乘法)使其慢>2倍，尽管它执行相同数量的内存读取
    • 但是，如果它受到处理能力的限制，那么为什么我在 GPU 上的矩阵乘法比在 CPU 上快 100 倍？
  • 为什么内核GPU-scanline-private 这么慢？内存局部性要好得多(每像素仅从全局内存读取 3 次而不是 9 次)，逻辑最小(无 ifs/开关)

测试是在我的笔记本电脑上使用 CPU 英特尔酷睿 i7 6700HQ Skylake 和 GPU nVidia 960M 完成的，方法是在 256x256 像素的浮点数组上运行内核 64x/帧。完整的代码可以在这里看到。

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

=========内核GPU 朴素2D 全局=(256,256) 本地=(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_async2D 全局=(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一维全局=(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一维全局=(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-naïve

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

在矩阵乘法中，每个子矩阵(补丁)用于另一个矩阵中所有行中的所有补丁。如果一个补丁中有 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 重用，因为每行的补丁较少，这使得总重用像平方功率率一样)

    • 卷积：增加补丁大小可提高平方功率速率的重用
      • 率

  • 矩阵：本地内存必须至少为 2 倍的平铺面积(子垫垫)

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

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

    • 卷积：将 4x4 加载到私有内存中不会做任何事情，而只是一个 4 像素计算(对于 3x3 画笔)，这意味着 4x4 内存 = 36 加 + 4 mul

拥有一个加法密集型内核为另一个繁重乘法的内核留出了空间，以并发或在同一内核中异步工作。也许如果您将其用于图像处理，也许您可以在内部添加一些"混合"或"调整大小"内核，以便它们协同工作？

扫描线版本是加载 3 个元素，执行 9 个添加 + 1 mul 然后重复，加载的元素停留 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 个 muls，所以更平衡，比如：1Tflops GPU 将有 500Gflops 用于添加，90 Gflops 用于添加。

朴素版本不使用本地内存，为飞行中的更多波前留出更多空间。本地内存版本实际上打破了L1访问模式，并减少了飞行中的波前。这减少了 VALU 占用。

您可以通过在工作组级别而不是线程级别执行扫描线来减少本地内存使用量。我的意思是这样的：

从内存加载：x x x x x x x xx 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)

它具有 1-D 扫描线，单线程用于 N 个周期或多线程化简用于 LogN 周期(考虑计算单元中有足够的线程)。

相关内容

最新更新


热门标签：

javascript python java c# php android html jquery c++ css ios sql mysql arrays asp.net json python-3.x ruby-on-rails .net sql-server django objective-c excel regex ruby linux ajax iphone xml vba spring asp.net-mvc database wordpress string postgresql wpf windows xcode bash git oracle list vb.net multithreading eclipse algorithm macos powershell visual-studio image forms numpy scala function api selenium