CUDA(十一) 周斌 CUDA程序深入优化

存储优化

CPU-GPU数据传输最小化，依赖于PCIe总线

• Host<->Device 数据传输带宽远低于 global memory（GPU内部）

        8GB/s（PCIe x16 Gen2）vs 156 GB/s & 515Ginst/s(C2050)

• 减少传输

 

中间数据直接在GPU分配、操作、释放。

有时更适合在GPU进行重复计算。现在CPU上计算一部分，然后传到GPU上去，减少重复计算量

如果没有减少数据传输的话，将CPU代码一直到GPU可能无法提升性能。

• 组团传输，较小整个系统外部的开销

大块传输好于小块：10微秒延迟，8GB/s => 如果数据小于80KB，性能将受延时支配

• 内存传输于计算时间重叠

双缓存double buffering平方缓存，一块做缓存，一块做操作，缓存或者操作结束，相互切换

Coalescing 合并

• Global memory 延时：400-800 cycles：最重要的性能影响因子。带宽虽然很大，但是延时很严重，

在Fermi， global memory 默认缓存于 一级缓存 L1

通过给nvcc指令设置参数“-Xptxas -dlcm=cg” 可以绕过一级缓存L1： 只缓存于二级缓存L2;

如果缓存：warp 的读写请求落到L1 cache line，只需要一次传输

# transaction = # L1 line accessed

如果没有缓存：有一些和并原则

但是传输大小可以减至32字节块

Memory Coalescing合并访存

给定一个矩阵以行优先的方式存储于global memory，对于thread来说合适的访存模式是什么。每个线程访问邻近的内存地址1,2,3,4，相邻的线程获取的时行存储的相邻数据

Coalescing 合并（Too old）

• 如果满足访存合并条件，半个warp 的线程访问Global memory的32，64，或128位宽数据，结果仅需要1或2次传输。half-warp of thread因为存储的数据类型不一样，32bit单精度浮点型=4byte。一个warp32个线程4个byte相当于128字节数据。一个warp需要128个字节数据，将128个字节数据排成一排，每个线程访问各自的位置。
• 依赖compute capability

1.0和1.1有更严格的访问要求

• Float（32-bit）举例

compute capability 1.0 和 1.1

• 第K个线程必须访问同一段里面的第k个字（或者2段连续128字节里面的第k个字），并不是所有线程都要参与。

Coalesces - 1 transaction
Out of sequence - 16 transactions
Misaligned - 16 transactions

混乱的随机访问对GPU的压力巨大

compute capability 1.2 和 higher

• 解决32B，64B和128B的内存段的传输问题
• 更小的传输可以避免浪费带宽

1 transaction - 64B segment

2 transactions - 64B and 32B segments

1 transaction - 128B segment

 

设计规则话的访存模式，避免出现随机访存，避免出现一个线程访问连续的存储空间，访存系统没办法进行合并，会反复的寻址，浪费带宽。

合并举例

• 小型kernel拷贝数据时有效带宽

偏移和步长对性能的影响

• 2款GPUs

GTX 280

              Compute cabability 1.3

              峰值带宽141 GB/s

FX 5600

              Compute capability 1.0

              峰值带宽77 GB/s

Copy with Offset偏移

__global__ void offsetCopy(float *odata, float *idata, int offset)
{
    int xid = blockIdx.x * blockDim.x + threadIdx.x + offset;
    odata[xid] = idata[xid]
}

Copy with Stride跳着步长

__global__ void strideCopy(float *odata, float *idata, int stride)
{
    int xid = (blockIdx.x * blockDim.x + threadIdx.x) * stride;
    odata[xid] = idata[xid];
}

• 按步长访问内存的场合，通常发生在多维问题里面，步长很大（>>18），多维数据经常发生
• 按步长访问global memory 可以通过shared memory避免。通过shared memory 规则化数据，通过shared memory做global memory不规则的访问

shared memory

• 比global memory快上百倍
• 可以通过缓存数据减少global memory访存次数
• 线程可以通过shared memory协作，线程协作

• 用来避免不满足合并条件的访存

读入shared memory重排顺序，从而支持合并寻址。

shared memory架构

• 很多线程访问存储器

因此，存储器被划分为banks线程区块

连续的32-bit访存被分配到连续的banks

• 每个bank每个周期可以响应一个地址

如果有多个bank的话可以同时响应更多地址申请

• 对同一bank进行多个并发访存将导致bank冲突

冲突的访存必须串行执行

No Bank Conflicts	No Bank Conflicts
线性寻址	随机1：1置换
stride == 1

bank编号并不是内存编号，内存地址取余。

2-way Bank Conflicts	8-way Bank Conflicts
线性寻址	线性寻址
stride == 2	stride == 8

 

几路的bank冲突：每个bank同时有几个线程在访问。虽然物理地址不同，但是bank映射后的地址（取余）是一样的。导致shared memory性能下降几倍。

 

shared memory bank 冲突

• shared memory 跟registers一样快，如果没有bank冲突的话，一个数量级的差异
• warp_seriallize profiler分析器的可以反映冲突情况，warp串行化
• 快速情况

half-warp内所有线程访问不同banks，没有冲突。

half-warp内所有线程读取同一地址，没有冲突（广播），不是同一个bank

• 慢速情况

Bank Confilct: half-warp 内多个线程访问同一个bank

访存必须串行化

代价 = 多个线程同时访问一个bank的线程数最大值，倍数性能的下降。

举例：Transpose 矩阵转置

• 每个线程块在矩阵的一个瓦片上操作
• 原始版本存在对global memory按步长访问的情况，

Element transposed by a half-warp of threads

读入一行，写入一列。

对于global memory 是不合并的，读入是合并读入写入不合并。出现读写总有一个是不合并的，按行读就要按列写，按列读就要按行写。

原始Transpose

__global__ void transposeNavie(float *odata, float *idata, int width, int height)
{
    int xIndex = blockIdx.x * TILE_WIDTH + threadIdx.x;//按照block方式合并读入
    int yIndex = blockIdx.y * TILE_WIDTH + threadIdx.y;

    int index_in = xIndex + width * yIndex;//非合并的方式写入
    int index_out = yIndex+ height * xIndex;

    odata[index_out] = idata[index_in];
}

通过shared memory 实现合并

• 先将瓦片的多列元素存入shared memory，再以连续化的数据写入global memory
• 需要同步__syncthreads()因为线程需要用到其他线程存储到shared memory 的数据

Element transposed by a half-warp of threads

也就是先将数据读到shared memory再在shared memory里面转置一下。shared memory读入的时候需要同步的。

__global__ void transposeCoalescee(float *odata, float *idata, int width, int height)
{
    _shared__ float tile[TILE_DIM][TILE_DIM];

    int xIndex = blockIdx.x * TILE_DIM + threadIdx.x;
    int yIndex = blockIdx.y * TILE_DIM + threadIdx.y;//读入的方式与之前的一致
    int index_in = xIndex + (yIndex) * width;

    xIndex = blockIdx.y * TILE_DIM + threadIdx.x;
    yIndex = blockIdx.x * TILE_DIM + threadIdx.y;
    int index_out = xIndex + (yIndex) * width;

    tile[threadIdx.y][threadIdx.x] = idata[index_in];//写入的时候用合并的方式写到shared memory里面去
    
    __syncthreads();

    odata[index_out] = title[threadIdx.x][threadIdx.y];
}

Transpose 存在Bank冲突

• 瓦片内16 * 16个floats 存在于shared memory

列中的数据存于相同的bank

读入瓦片一列数据时存在16-way bank confict

• 解决方案-填充shared memory 数组

__shared__ float title[TILE_DIM][TILE_DIM+1]

反对角线上的数据存在相同的bank

访存的时候又16号地址，例如0和16号地址同时被访问

填充shared memory

m00	m01	m02	m03
m10	m11	m12	m13
m20	m21	m22	m23
m30	m31	m32	m33

读入的时候没有冲突，访问不同的bank

m00	m01	m02	m03	m04
m10	m11	m12	m13	m14
m20	m21	m22	m23	m24
m30	m31	m32	m33	m34

写的时候按列写，产生bank冲突

__shared__ float title[TILE_DIM][TILE_DIM+1]4*5

访存的时候依旧是4*4，第5列是不实际使用的只是用来占位。

每一列的访问时斜对角访问效果，访问落在不同的bank里面，达到不同的访问效果。

横向访问不同bank没有冲突，列项访问同一个bank有冲突，但是作为矩阵转置，总是包含横向和纵向的访问，增加占位之后，横纵向都是落在不同的bank，也就是说填充bank+1可以有效的避免bank冲突

CUDA的texture纹理

• Texture是用于读入数据的一个对象，放在global memory里面，只能读入对象
• 优点，针对游戏或者图像处理

数据被缓存

              特别适用于无法合并访存的场合，小数据凌乱读取

支持过滤

              线性/双线性/三线性插值

             专用硬件

warp模式（针对“越界”寻址）

              裁剪到边缘或者重复，越界坐标被折叠wrap（算术取模），越界坐标被边界最近的值替换clamp。，

一维、二维、三维寻址

             以整数或者归一化小数作为坐标

CUDA的texture纹理类型

• 绑定到线性存储空间

Global memory 地址绑定到texture

仅1D

整数寻址

不支持过滤，不支持寻址模式

• 绑定到CUDA数组

块状线性的CUDA array 绑定到texture

1D, 2D, or 3D 

小数（尺寸归一化）

支持过滤

支持寻址模式（折叠，重复）

• 绑定到pitch linear(CUDA 2.2)

Global memory 地址绑定到texture

2D

小数/整数寻址、过滤和折叠/重复寻址模式，类似于CUDA arrays

CUDA 纹理操作步骤

• Host(CPU) 代码：

分配或者获取内存（global linear/pitch linear, or CUDA array）

创建纹理引用对象

                 目前必须属于文件域

纹理引用绑定到内存或数组

用完以后

                   为纹理引用解绑，释放资源

• Device(kernel)代码：

通过纹理引用获取数据

线性内存问题：text1Deftch()

数组纹理：tex1D() or tex2D() or tex3D()

Pitch线性纹理：tex2D()

__global__ shiftCopy(float *idata, float *odata, int shift)
{
    int xid = blockIdx.x * blockDim.x +threadIdx.x;
    odata[xid] = idata[xid +shift];
}

texture <float> texRef;

__global__ void textureShiftCopy(float *idata, float *odata, int shift)
{
    int xid = blockIdx.x * blockDim.x +threadIdx.x;
    odata[xid] = tex1Dfetch(texRef, xid +shift);
}

 使用texture对小数据访问比较好，性能更好。

总结

• 如果遵循一些简单的原则，GPU硬件在数据可并行计算问题上，可以达到很好的性能：

有效利用并行性

尽可能合并内存

利用shared memory

开发其他存储空间

       Texture

      Constant

减少bank冲突（shared memory使用）

SM资源分割

SM资源动态分配

SM硬件基本的完整单元

有助于调节整个配备，就是启动block，线程应该怎样调度资源。

SM

Thread block slots	(G80 Limits)8
Thread slots	768
Registers	8K registers/ 32 memory
Shared Memory	16K

如果我们有

8 blocks of 96 threads（每个block有96个线程8*96=768个线程）

4 blocks of 192 threads（每个block有192个线程）

But not 8 blocks of 192 threads

如果有（假设每个block含256个线程）

768 threads(3 blocks)，每个线程用10registers

512 threads(2 blocks)，每个线程用11registers  11*768>8k

资源分割是一个多约束的形式。

多用寄存器将减少线程级别的并行，但是由于每一个线程使用的寄存器增多了，有些情况下反而既降低系统性能

performance cliff：增加资源用量后导致并行性急剧下降，导致占用率下降，例如增加寄存器数量。除非为了隐藏global memory访存延迟。

CUDA Occupancy Calculator

CUDA占用率计算器

\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v10.0\tools下的一个excel 计算占用率，表示计算机硬件繁忙程度，如果占用率低，意味着录用率低下。

kernel启动参数配置

Grid Size试探法

• # of blocks > # of SM软件的block数大于硬件是SM数

每个SM至少有一个work-group在执行，不会有人在闲着

• # of blocks / # of SM > 2

多个blocks可以在SM上并发执行

如果一个block在等待同步，启动另外一个block

• # of blocks / # of SM > 100对支持未来设备用于很好的伸缩性

包含grid和block大小。

• 块大小必须是32（warp尺寸）的倍数
• 需要尽量多的warp尽可能的隐藏延时
• 最少：64.通常采用128 or 256视程序而定
• 依赖于问题本身，请多试验！，程序可以改变块大小

Latency Hiding隐藏延时

关键点：

指令按顺序发布

一个线程的任一操作数没有预备好时将阻塞等待

延时在切换线程时被隐藏

总结：

需要足够多的线程来隐藏延时

Occupancy占用率

Occupancy：一个SM里面激活warp与最大可容纳warp数目的比值

最大warp数目：32 in Tesla, 48 in Fermi

SM资源动态分割

某型号GPU

Shared memory按block分割

Registers按线程分割:<=63

Thread block slots:<=8

Thread slots:<=1536

以上任何一个因素都可能是：同一时间在同一SM启动线程的数量的影响因子

多因素约束，资源动态分割

延时隐藏利用率计算

• 假设global memory需要400cycles，我们需要400/2（假定每个算数指令需要两个时钟周期来完成）=200条算数指令来隐藏延时
• 例如，假设代码中针对每个global memory访存含8条独立算数指令（不依赖于global memory的访问），需要200/8 ~ 26 warps（26/48=54% occupancy）
• 注意超过54%，在这个例子里面更高的利用率无法带来性能的提升，因为延时已经隐藏掉了。
• 如果一条指令用到寄存器里面由上一条语句写入的结果，延时大约~24cycles
• 所以我们需要24/2（一个指令需要两个时钟周期来完成）=13warp（大约24个时钟周期，因此要大一个）是来隐藏由于依赖寄存器引起的延时，相应的利用率是：27%

Data Prefetching数据预读

再一次global memory读操作和实际用到这个数据的语句中间，插入独立于以上数据的指令，可以隐藏访存延迟

float m = Md[i];//将数据从存储器读入到寄存器
float f = a * b + c * d;//系统通过切换多个warp按隐藏。编译器会自动的将它插入这个位置。执行指令，不依赖于读内存的操作。
float f2 = m * f;//使用取出的m

从global memory预读数据可以有效的提升独立指令的数量，在global memory读取和使用两者之间。

帮助编译器完成延时掩藏的效果。

例如矩阵相乘中瓦片化的使用shared memory的过程

for (/* ... */)
{
    // Load current tile into shared memory
    __syncthreads();
    // Accumulate dot product
    __syncthreads();
}

使用数据预期，有一个循环，不断的读入瓦片数据到shared memory，然后在shared memory里面做累加。

引入预读操作的瓦片化matrix multiple

// Load first tile into registers把第一个瓦片的数据先读到寄存器

for(/* ... */)
{
    // Deposit registers into shared memory然后从寄存器放入共享存储器
    __syncthreads();
    // Load next tile into registers利用shared memory作为预取空间，将下一个瓦片嘟嘟奥shared memory，由于下一个瓦片的数据和现在做乘累加的数据没有依赖性，所以他们就会填满这个访存的延时隐藏。为下一个循环节点预读。
    // Accumulate dot product//，然后在做现在瓦片的乘累加。这些指令被足够多的线程执行，从而隐藏了下一个瓦片预读内存产生的延时。
    __syncthreads();
}

可以得到100%的性能提升。

指令吞吐量优化

指令优化

如果你发现代码性能受限于指令

      如果不够仔细的话，计算密集型算法很容易受限于带宽

     典型情况，在存储器和执行配置优化完成后，再担心指令优化。指令计算能力很强

Fermi 算术指令吞吐量

Int & fp32: 2cycles

fp64: 2cycles

Fp32 transendental: 8 cycles

Int divide and modulo are expensive

除以2^n，采用“>>n”移位的方式，前提是2的整数次

以2^n求模，采用“&(2^n - 1)”

避免double到float的类型自动转换

添加"f" 到float常量（e.g. 1.0f），因为缺省是double

Fermi 缺省： -tz=false,-pree-div=true,-pree-sqrt=true for IEEE compliance

运行时数学库和固有函数

两种类型的运行时数学库函数

func():

慢但是精度高(5 ulp or less)

Examples: sin(x) , exp(x), pow(x,y)

__func();

快但是精度低(see prog, guide for full details)

Examples: __sin(x), __exp(x), __pow(x,y)

其他固有函数:

__sincos(), __rcp(), ...

详细列表间《CUDA编程指南》附录C.2

-use-fast-math: func() 强制转换为__func()强制将精度高但是慢的函数转换为速度快但是精度低的函数。

循环展开

for (int k = 0; k < BLOCK_SIZE; ++k)//++k循环计数器，计数器更新；<判断，判断有没有到循环的边界，分支
{
    Pvalue += Ms[ty][k] * Ns[k][tx];//访存地址的运算
}

每轮循环包含的指令：

一条浮点数乘法

一条浮点数加法

还有其它运算？编译器还有很多额外的工作

指令混合：

2条浮点预算指令

1条循环分支指令

2地址运算指令

1循环计数器自增指令

• 仅1/3是浮点计算
• 但是期望达到理论1TFLOP(Fermi)
• 考虑loop unrolling

Pvalue +=
    Ms[ty][0] * Ns[0][tx] +
    Ms[ty][1] * Ns[1][tx] +
    ...
    Ms[ty][15] * Ns[15][tx]; //BLock_SIZE = 16

• 不再有循环

不再有循环计数器更新

不再有分支

常量索引 - 不再有地址运算

自动实现：

#pragma unroll BLOCK_SIZE
for (int k = 0; k < BLOCK_SIZE; ++k)
{
    Pvalue += Ms[ty][k] * Ns[k][tx];
}

循环展开有什么缺点？可扩展性不强，如果blocksize变换比较多，变量没发去适应，可能出现错误。