GPU (OpenCL) 上的 3x3 卷积不应该快得多
Shouldn't be 3x3 convolution much faster on GPU (OpenCL)
我正在学习如何为 GPU 优化代码。我读到了内存位置的重要性。我还看到了一些 tutorials and examples 的 GPU 卷积。基于此，我编写并测试了几个自己的内核。令人惊讶的是，我发现最简单的 naive kerell 是最快的！？ 并且它比 CPU 快 <10 倍。 （是的，我通过 运行 kenrnel 64x 摊销了 upload/download 时间）。
我做错了什么？ 我希望卷积正是 GPU 优化的那种操作。如果我能得到100x speed-up on matrix multiplication，为什么卷积这么慢？
性能 [CPU ticks/pixel]（越低越好）：
- CPU-幼稚 9.5
- GPU-naive 1.64
- GPU-local 2.56
- GPU-local_async 15.10
- GPU-扫描线-私有 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 次）并且逻辑是最小（无 ifs/switches）
测试是在 my laptop with CPU Intel Core i7 6700HQ Skylake and GPU nVidia 960M by running the kernels 64x/frame on floating point array of 256x256 pixels. The code full can be seen here 上完成的。
===========内核代码===========
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] );
}
kernel 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二维全局=(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-幼稚
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 可能）
矩阵multiplication/with 16x16 子矩阵）vs 卷积（17x17 画笔大小）：
矩阵：L2重用率随主矩阵大小增加，或L1重用率随子矩阵大小增加(L1)
- 卷积：所有图像尺寸的总重用率相同，但 L1 使用率随着画笔尺寸的增加而增加（好）
矩阵：每个工作组 16*16*16 次乘法 + 16*16*16 次加法
- 卷积：每个线程 17*17 次加法 + 1 次乘法（不好）
矩阵：统一线程使用，没有if-else，所有本地内存都被重新使用
- 卷积：需要加载比边界（厚度为 16 的幽灵墙）至少 16 像素，这些边界将被相邻工作组重新使用，但这些相邻工作组可能在另一个计算单元中，只需使用 L2 而不是在同一个计算单元上使用 L1（丑陋）
- 这就是为什么我建议异步工作组副本在同一计算单元（和 L1）上使用这些邻居并提高重用率。
矩阵：增加补丁大小也会增加子矩阵乘法中立方幂率的重用（但由于每行补丁较少，因此减少了 L2 重用，这使得总重用像平方功率比）
- 卷积：增加补丁大小会增加平方功率率的重用
矩阵：本地内存必须至少为 2x tile area (sub mat-mat mul)
- 卷积：局部内存必须至少是瓦片面积+鬼墙面积
矩阵：可以在私有内存中做4x4次子乘法（每个元素使用4次）即4x4内存=64加+64乘
- 卷积：将 4x4 加载到私有内存中只做 4 像素计算（对于 3x3 画笔），这意味着 4x4 内存 = 36 add + 4 mul
拥有一个重加法内核为另一个重乘法内核留出了空间，可以同时工作或在同一内核中异步工作。也许如果您将其用于图像处理，也许您可以在其中添加一些 "blend" 或 "resize" 内核，以便它们一起工作？
扫描线版本正在加载 3 个元素，执行 9 add + 1 mul 然后重复，加载的元素停留 3 圈，这意味着它们只被重复使用 3 次并且它的邻居（x 或 y directio）可能不会掉落在邻居线程甚至邻居工作组中。 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 个乘法，所以更加平衡，例如：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
现在将它用于工作组级别的扫描线：ac c u m u l a to or (+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 个周期的多线程 reduce（考虑到计算单元中有足够的线程）。
我正在学习如何为 GPU 优化代码。我读到了内存位置的重要性。我还看到了一些 tutorials and examples 的 GPU 卷积。基于此，我编写并测试了几个自己的内核。令人惊讶的是，我发现最简单的 naive kerell 是最快的！？ 并且它比 CPU 快 <10 倍。 （是的，我通过 运行 kenrnel 64x 摊销了 upload/download 时间）。
我做错了什么？ 我希望卷积正是 GPU 优化的那种操作。如果我能得到100x speed-up on matrix multiplication，为什么卷积这么慢？
性能 [CPU ticks/pixel]（越低越好）：- CPU-幼稚 9.5
- GPU-naive 1.64
- GPU-local 2.56
- GPU-local_async 15.10
- GPU-扫描线-私有 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 次）并且逻辑是最小（无 ifs/switches）
测试是在 my laptop with CPU Intel Core i7 6700HQ Skylake and GPU nVidia 960M by running the kernels 64x/frame on floating point array of 256x256 pixels. The code full can be seen here 上完成的。
===========内核代码===========
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] );
}
kernel 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二维全局=(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-幼稚
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 可能）
矩阵multiplication/with 16x16 子矩阵）vs 卷积（17x17 画笔大小）：
矩阵：L2重用率随主矩阵大小增加，或L1重用率随子矩阵大小增加(L1)
- 卷积：所有图像尺寸的总重用率相同，但 L1 使用率随着画笔尺寸的增加而增加（好）
矩阵：每个工作组 16*16*16 次乘法 + 16*16*16 次加法
- 卷积：每个线程 17*17 次加法 + 1 次乘法（不好）
矩阵：统一线程使用，没有if-else，所有本地内存都被重新使用
- 卷积：需要加载比边界（厚度为 16 的幽灵墙）至少 16 像素，这些边界将被相邻工作组重新使用，但这些相邻工作组可能在另一个计算单元中，只需使用 L2 而不是在同一个计算单元上使用 L1（丑陋）
- 这就是为什么我建议异步工作组副本在同一计算单元（和 L1）上使用这些邻居并提高重用率。
- 卷积：需要加载比边界（厚度为 16 的幽灵墙）至少 16 像素，这些边界将被相邻工作组重新使用，但这些相邻工作组可能在另一个计算单元中，只需使用 L2 而不是在同一个计算单元上使用 L1（丑陋）
矩阵：增加补丁大小也会增加子矩阵乘法中立方幂率的重用（但由于每行补丁较少，因此减少了 L2 重用，这使得总重用像平方功率比）
- 卷积：增加补丁大小会增加平方功率率的重用
矩阵：本地内存必须至少为 2x tile area (sub mat-mat mul)
- 卷积：局部内存必须至少是瓦片面积+鬼墙面积
矩阵：可以在私有内存中做4x4次子乘法（每个元素使用4次）即4x4内存=64加+64乘
- 卷积：将 4x4 加载到私有内存中只做 4 像素计算（对于 3x3 画笔），这意味着 4x4 内存 = 36 add + 4 mul
拥有一个重加法内核为另一个重乘法内核留出了空间，可以同时工作或在同一内核中异步工作。也许如果您将其用于图像处理，也许您可以在其中添加一些 "blend" 或 "resize" 内核，以便它们一起工作？
扫描线版本正在加载 3 个元素，执行 9 add + 1 mul 然后重复，加载的元素停留 3 圈，这意味着它们只被重复使用 3 次并且它的邻居（x 或 y directio）可能不会掉落在邻居线程甚至邻居工作组中。 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 个乘法，所以更加平衡，例如：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 现在将它用于工作组级别的扫描线：ac c u m u l a to or (+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 个周期的多线程 reduce（考虑到计算单元中有足够的线程）。