※ [本文转录自 VideoCard 看板] 作者: a5000ml (咖啡里的海洋蓝) 看板: VideoCard 标题: [分享] CUDA 程式设计(12) -- 速成篇(下) 时间: Wed Nov 26 22:00:46 2008 ============================================================================ 第七招 合并存取 ============================================================================ 【合并存取】(coalesced I/O) 是 CUDA 中最基本且最重要的最佳化手段, 因为 GPU 的计算能力太强, 使得效能瓶颈卡在显示记忆体到 GPU 之间的 I/O 上, 合并存取可让多个显示记忆体的交易合并成一次, 而加速记忆体的存取. 现阶段 GPU 在合并存取上是自动发生的, 以半个 warp 为单位(16 个相邻的执行绪), 如果它们的资料位址是连续的, 就会被合并, 所以使用上很简单, 只要 threadIdx 对齐即可, 它可以合并 4, 8, 16 bytes 的资料, 成为一次 or 两次的交易, 合并成的最大封包长度可为 32, 64, 128 bytes (其中 32 bytes 的封包只有在 版本 1.2 以後支援, 避免位址分散情况下的 overhead)，以下为连续的资料位址 的合并情况 16 x 4 bytes -> 64 bytes 16 x 8 bytes -> 128 bytes 16 x 16 bytes -> 2 x 128 bytes 【范例: 矩阵的 transpose】 ======================================================================== #include <stdio.h> #include <stdlib.h> #include <time.h> #include <math.h> #include <cuda.h> //transpose m*n matrix a to n*m matrix b on host --------------------- void transpose_host(float* b, float* a, int m, int n){ for(int y=0; y<m; y++) for(int x=0; x<n; x++){ b[x*m+y]=a[y*n+x]; } } //transpose naive (读取合并). ---------------------------------------- __global__ void transpose_naive_cr(float* b, float* a, int m, int n){ int x=blockIdx.x*blockDim.x + threadIdx.x; int y=blockIdx.y*blockDim.y + threadIdx.y; if(y<m && x<n){ b[x*m+y]=a[y*n+x]; } } //transpose naive (写入合并). ---------------------------------------- __global__ void transpose_naive_cw(float* b, float* a, int m, int n){ //x,y 的 threadIdx 足标对调. int x=blockIdx.x*blockDim.x + threadIdx.y; int y=blockIdx.y*blockDim.y + threadIdx.x; if(y<m && x<n){ b[x*m+y]=a[y*n+x]; } } //transpose shared (读取 & 写入合并). ------------------------------- __global__ void transpose_shared(float* b, float* a, int m, int n){ //宣告共享记忆体. __shared__ float s[256]; //读取合并. int x=blockIdx.x*blockDim.x + threadIdx.x; int y=blockIdx.y*blockDim.y + threadIdx.y; if(y<m && x<n){ int t=threadIdx.y*blockDim.x + threadIdx.x; int i=y*n+x; s[t]=a[i]; } __syncthreads(); //写入合并 (x,y 的 threadIdx 足标对调). x=blockIdx.x*blockDim.x + threadIdx.y; y=blockIdx.y*blockDim.y + threadIdx.x; if(y<m && x<n){ //共享记忆体中的 threadIdx 足标亦要对调. int t=threadIdx.x*blockDim.y + threadIdx.y; int o=x*m+y; b[o]=s[t]; } } //计算相对误差. ---------------------------------------------------- double rd(float*a, float*b, int size){ double s=0, d=0; for(int k=0; k<size; k++){ double w=a[k]-b[k]; s+=a[k]*a[k]; d+=w*w; } return sqrt(d/s); } //timer functions -------------------------------------------------- time_t timer[10]; void set_timer(int k=0){ timer[k]=clock(); } double get_timer(int k=0){ return (double)(clock()-timer[k])/CLOCKS_PER_SEC; } //测试程式. (输入 m*n 矩阵) ---------------------------------------- void test(int m, int n, int loop=100, int loop_host=10){ int size=m*n; printf("matrix size: %d x %d\n", m,n); //配置主记忆体, 并设定初始值. float *a=new float[size]; float *b=new float[size]; float *c=new float[size]; for(int k=0; k<size; k++){ a[k]=(float)rand()*2/RAND_MAX-1; b[k]=0; } //配置显示记忆体, 载入资料. float *ga, *gb; cudaMalloc((void**)&ga, size*sizeof(float)); cudaMalloc((void**)&gb, size*sizeof(float)); cudaMemcpy(ga, a, size*sizeof(float), cudaMemcpyHostToDevice); cudaMemcpy(gb, b, size*sizeof(float), cudaMemcpyHostToDevice); //网格区块设定. dim3 grid(n/16+1,m/16+1,1); //网格要含盖所有范围, 所以除完要加 1. dim3 block(16,16,1); //区块设定 16x16. printf("grid(%d,%d,%d)\n",grid.x,grid.y,grid.z); printf("block(%d,%d,%d)\n",block.x,block.y,block.z); //测试 transpose_host 函数. set_timer(); for(int k=0; k<loop_host; k++){ transpose_host(b,a,m,n); } double t0=get_timer()/loop_host; printf("host time: %g ms\n",t0*1000); //测试 transpose_naive_cr 函数. cudaMemset(gb,0,size*sizeof(float)); transpose_naive_cr<<<grid,block>>>(gb,ga,m,n); cudaMemcpy(c, gb, size*sizeof(float), cudaMemcpyDeviceToHost); set_timer(); cudaThreadSynchronize(); for(int k=0; k<loop; k++){ transpose_naive_cr<<<grid,block>>>(gb,ga,m,n); } cudaThreadSynchronize(); double t1=get_timer()/loop; printf("naive(r) time: %g ms (%dx)\t error: %g\n",t1*1000, (int)(t0/t1),rd(b,c,size)); //测试 transpose_naive_cw 函数. cudaMemset(gb,0,size*sizeof(float)); transpose_naive_cw<<<grid,block>>>(gb,ga,m,n); cudaMemcpy(c, gb, size*sizeof(float), cudaMemcpyDeviceToHost); set_timer(); cudaThreadSynchronize(); for(int k=0; k<loop; k++){ transpose_naive_cw<<<grid,block>>>(gb,ga,m,n); } cudaThreadSynchronize(); double t2=get_timer()/loop; printf("naive(w) time: %g ms (%dx)\t error: %g\n",t2*1000, (int)(t0/t2),rd(b,c,size)); //测试 transpose_shared 函数. cudaMemset(gb,0,size*sizeof(float)); transpose_shared<<<grid,block>>>(gb,ga,m,n); cudaMemcpy(c, gb, size*sizeof(float), cudaMemcpyDeviceToHost); set_timer(); cudaThreadSynchronize(); for(int k=0; k<loop; k++){ transpose_shared<<<grid,block>>>(gb,ga,m,n); } cudaThreadSynchronize(); double t3=get_timer()/loop; printf("shared(r/w)time: %g ms (%dx)\t error: %g\n",t3*1000, (int)(t0/t3),rd(b,c,size)); //释放记忆体 cudaFree(ga); cudaFree(gb); delete [] a; delete [] b; delete [] c; } //主函数 ------------------------------------------------------------------ int main(){ srand(time(0)); printf("----------------------------\n"); test(2047,4000); printf("----------------------------\n"); test(2048,4000); printf("----------------------------\n"); test(2049,4000); return 0; } 【说明】 ========================================================================== (1) 程式中使用 dim3(16,16,1) 的区块设定, 所以 threadIdx.x=0~15 形成所谓的 「半个 warp」. (2) 在 transpose_naive_cr() 中 int x=blockIdx.x*blockDim.x + threadIdx.x; 包含半个 warp (threadIdx.x), 所以 x 在半个 warp 中连续, 读取 a[y*n+x] 被合并, 形成所谓的【读取合并】(coalesced read) (3) 在 transpose_naive_cw() 中, 把 naive_cr 版本 (x,y) 的 threadIdx 足标对调, 换成 y 在半个 warp 上是连续的, 使得写入 b[x*m+y] 被合并 int x=blockIdx.x*blockDim.x + threadIdx.y; int y=blockIdx.y*blockDim.y + threadIdx.x; 形成所谓的【写入合并】(coalesced write). (4) 在 transpose_shared() 中, 上半周期应用 half warp 到读取位址的连续 int x=blockIdx.x*blockDim.x + threadIdx.x; int y=blockIdx.y*blockDim.y + threadIdx.y; int i=y*n+x; <--- (OOO) + threadIdx.x 形成【读取合并】, 并把资料存於快速的共享记忆体中. 下半周期要对调 threadIdx 足标, 使得 half warp 在写入位址连续 int x=blockIdx.x*blockDim.x + threadIdx.y; int y=blockIdx.y*blockDim.y + threadIdx.x; int o=x*m+y; <--- (XXX) + threadIdx.x 形成【写入合并】 (5) transpose_shared() 中使用共享记忆体的理由是资料的读取和写入使用不同的 执行绪 (threadIdx 足标对调), 所以执行绪之间必需交换资料, 注意足标对调时 共享记忆体中的 threadIdx 足标亦要对调, 以对应到相对的元素 读取 int t=threadIdx.y*blockDim.x + threadIdx.x; 写入 int t=threadIdx.x*blockDim.y + threadIdx.y; 【执行结果】 ========================================================================= 这次测试的机型为 GTX260 vs Intel E8400(3.00GHz), 我们可以看出, 合并存取的另一个好处是减少 bank conflict 对它的影响, 当 host 因为 bank conflict 变慢时, GTX260 还是维持差不多的水准. 平均 shared 版的 transpose 的效能大概是 host 的 20 多倍 ---------------------------- matrix size: 2047 x 4000 grid(251,128,1) block(16,16,1) host time: 47 ms naive(r) time: 8 ms (5x) error: 0 //读取合并 (naive_cr) naive(w) time: 4 ms (11x) error: 0 //写入合并 (naive_cw) shared(r/w)time: 2.1 ms (22x) error: 0 //读写合并 (shared) ---------------------------- matrix size: 2048 x 4000 (这是可以产生 bank conflict 的情况) grid(251,129,1) block(16,16,1) host time: 305 ms naive(r) time: 8.5 ms (35x) error: 0 naive(w) time: 4.1 ms (74x) error: 0 shared(r/w)time: 2.1 ms (145x) error: 0 //因为 host 太逊了 ---------------------------- matrix size: 2049 x 4000 grid(251,129,1) block(16,16,1) host time: 46 ms naive(r) time: 8.1 ms (5x) error: 0 naive(w) time: 4.2 ms (10x) error: 0 shared(r/w)time: 2.2 ms (20x) error: 0 【GT200 的改进】 ========================================================================== 在 G80/G90 系列 (compute 1.0 & 1.1), 合并存取必需要符合许多要求, 例如 半个 warp 中的资料位址要按照 threadIdx 的顺序, 封包的边界位址必需满足 64 bytes 或 128 bytes 的对齐要求, 否则会放弃合并, 产生 16 个记忆体要求. 在 GT200 系列 (compute 1.2 之後), 合并存取有大幅度改进, 容许半个 warp 的 资料位址不按照顺序排列, 不需 16 个执行绪全部连续, 容许部份连续的情况 (此时会进行部份位址合并, 并拆成数个封包发出), 在节省频宽方面也会选择 最有效率的方式分解封包, 并引进 32 bytes 的小型封包, 避免较小的连续区块 仍需使用大的合并封包, 整体记忆体效能比 G80/G90 系列好甚多. 【范例：GTX260 半个 warp 的存取顺序】 ========================================================================== (1) 将范例中的 transpose_naive_cr 的 x 反序 __global__ void transpose_naive_cr(float* b, float* a, int m, int n){ int x=blockIdx.x*blockDim.x + (15-threadIdx.x); int y=blockIdx.y*blockDim.y + threadIdx.y; if(y<m && x<n){ b[x*m+y]=a[y*n+x]; } } (2) 使用人工乱数, 将范例中的 transpose_naive_cw 改为随机顺序 __global__ void transpose_naive_cw(float* b, float* a, int m, int n){ int x=blockIdx.x*blockDim.x + threadIdx.y; int y=blockIdx.y*blockDim.y + (threadIdx.x*7+5)%16; if(y<m && x<n){ b[x*m+y]=a[y*n+x]; } } 半个 warp 的存取顺序对应如下 ----------------------- tid.x (tid.x*7+5)%16 ----------------------- 0 5 1 12 2 3 3 10 4 1 5 8 6 15 7 6 8 13 9 4 10 11 11 2 12 9 13 0 14 7 15 14 ----------------------- 【执行结果】 ============================================================================ 在 GTX260 上执行的时间仍相同, 基本上不太受存取顺序的影响。 ---------------------------- matrix size: 2047 x 4000 grid(251,128,1) block(16,16,1) host time: 47 ms naive(r) time: 7.8 ms (6x) error: 0 //naive_cr 版 naive(w) time: 4.1 ms (11x) error: 0 //naive_cw 版 shared(r/w)time: 2 ms (23x) error: 0 //shared 版 ---------------------------- matrix size: 2048 x 4000 (这是可以产生 bank conflict 的情况) grid(251,129,1) block(16,16,1) host time: 308 ms naive(r) time: 8.5 ms (36x) error: 0 naive(w) time: 4.1 ms (75x) error: 0 shared(r/w)time: 2 ms (153x) error: 0 ---------------------------- matrix size: 2049 x 4000 grid(251,129,1) block(16,16,1) host time: 46 ms naive(r) time: 8.1 ms (5x) error: 0 naive(w) time: 4.1 ms (11x) error: 0 shared(r/w)time: 2 ms (22x) error: 0 ---------------------------- 【补充: warp】 ========================================================================== warp 是 GPU 硬体上一次执行的实际单位, 一个区块可分成数个 warp, 分时在 multiprocessor 中执行, 例如 128 threads 的区块 => 4 warps 200 threads 的区块 => 6+1 warps (多出的未填满仍算 1 个 warp) 详细情形在之後硬体篇会再 review, warp 和 threadIdx 的关系如下 ------------------- warp threadIdx ------------------- 0 0~31 1 32~63 2 64~95 3 96~127 ... ------------------- 半个 warp 是记忆体合并的单位, 也就是一个 warp 的记忆体读写可分成两组 独立进行合并, 和 threadIdx 的关系如下 ----------------------- half warp threadIdx ----------------------- 0 0~15 1 16~31 2 32~47 3 48~63 ... ----------------------- 【补充: threadIdx 顺序】 ========================================================================= 若区块中的执行绪使用 3D 结构安排时, 其顺序是 threadIdx.x 在最里面, 然後是 threadIdx.y, 最外层是 threadIdx.z, 结构类似 C/C++ 的三维阵列 在记忆体中的顺序 threads[z][y][x]; 打平成一维 tid = threadIdx.z * (blockDim.y*blockDim.x) threadIdx.y * (blockDim.x) threadIdx.x; 所以在范例中使用 dim(16,16,1) 区块, 其 threadIdx.x 维度包含 16 个执行绪, 刚好每一个 x 维度形成半个 warp, 整个区块有 8 个 warps. -- 。o O　○。o O　○。o O　○。o O　○。o O　○。o 国网 CUDA 中文教学 DVD 影片 (免费线上版) 请至国网的教育训练网登入 https://edu.nchc.org.tw BT 牌的种子下载点 http://www.badongo.com/file/12156676 --

※ 发信站: 批踢踢实业坊(ptt.cc) ◆ From: 114.45.212.139 --

※ 发信站: 批踢踢实业坊(ptt.cc) ◆ From: 114.45.212.139