Pro*_*ala 5 c performance gpgpu convolution opencl
我正在学习如何优化GPU的代码.我读到了记忆局部性的重要性.我也看过一些GPU卷积的教程和例子.基于此,我编写并测试了几个自己的内核.令人惊讶的是,我发现最简单的天真kernell是最快的!？它比CPU快10倍.(是的,我通过运行kenrnel 64x分摊上传/下载时间).
我做错了什么？我希望卷积只是GPU优化的那种操作.如果我可以在矩阵乘法上获得100倍的加速,为什么卷积是如此之慢？
性能[CPU滴答/像素](越低越好):编辑:GPU-scanline_async我在阅读了有关的建议之后做了async_work_group_copy
我想知道两件事:
测试是在我的笔记本电脑上使用CPU Intel Core i7 6700HQ Skylake和GPU nVidia 960M在256x256像素的浮点阵列上运行内核64x /帧.完整的代码可以在这里看到.
===========内核代码===========
内核 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天真
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;
} }
}
在矩阵乘法中，每个子矩阵（补丁）用于另一个矩阵中所有行中的所有补丁。如果一个 patch 中有 2x2 子矩阵，并且主矩阵是 20x20，则每个子矩阵使用 10 次进行乘法。GPU 通常使用 16x16 或 32x32 大小的补丁，这意味着，对于 2kx2k 乘法，每个 16x16 补丁至少重复使用 128 次。
MM reuse = 128
并添加子矩阵-子矩阵乘法的重复使用，这足以将GPU推向极限。
在 3x3 卷积中，3x3 patch 不用于整个扫描线或整个图片。仅重复使用其像素。
3x3 模板：每个像素被相邻的 8 个模板重复使用。
5x5 模板：每个像素都被相邻的 24 个模板重复使用。
为了赶上矩阵乘法，它需要
11x11 stencil to have a reuse of 120
它也比矩阵乘法更局部，并且应该比它获得更多的 gflops，但它没有进行等量的乘法和加法。
它执行 9 次加法 + 1 次乘法。
8 个潜在的乘法丢失。近一半的 GFLOPS 限制丢失。
您应该尝试异步工作组副本。
矩阵乘法/带有 16x16 子矩阵）与卷积（17x17 画笔大小）：
矩阵：L2复用率随主矩阵大小而增加，或L1复用率随子矩阵大小（L1）而增加
矩阵：每个工作组 16*16*16 乘法 + 16*16*16 加法
Matrix：统一线程使用，没有if-else，所有本地内存都被重用
矩阵：增加补丁大小也会增加子矩阵乘法中三次幂率的重用（但会减少 L2 重用，因为每行补丁较少，这使得总体重用类似于平方幂率）
矩阵：本地内存必须至少为 2 倍图块区域（sub mat-mat mul）
矩阵：可以在私有内存中进行 4x4 子乘法（每个元素使用 4 次），这意味着 4x4 内存 = 64 add+64 mul
拥有一个繁重的加法内核为另一个繁重的乘法内核同时工作或在同一内核中异步工作留下了空间。也许如果您将其用于图像处理，也许您可以在其中添加一些“混合”或“调整大小”内核，以便它们一起工作？
Scanline 版本正在加载 3 个元素，执行 9 add + 1 mul 然后重复，加载的元素停留 3 圈，这意味着它们仅重复使用 3 次，并且其邻居（x 或 y 方向）可能不会落在邻居线程中，甚至不会落入邻居工作组。此外，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 用于加法，90 Gflops 用于乘法。
天真的版本不使用本地内存，为飞行中的更多波前留出更多空间。本地内存版本实际上打破了 L1 访问模式，并减少了飞行中的波前。这减少了VALU占用。
您可以通过在工作组级别而不是线程级别执行扫描线来减少本地内存使用量。我的意思是这样的：
从内存加载：xxxxxxxxx对其执行扫描线：（从左到右，1-D）abcdefghij现在将其用于工作组级别的扫描线：累加器（+新）（从上到下）zxzxzxzxzx（-旧）
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 周期（考虑计算单元中有足够的线程）。