CUDA编程指南阅读笔记(六)

4.CUDA C语言编程接口

接上文。

4.3 CUDA C Runtime

4.3.3 共享内存(Shared Memory)

共享内存是CUDA设备中非常重要的一个存储区域,有效地使用共享内存可以充分利用CUDA设备的潜能,极大提升程序性能。那么,共享内存有哪些特点呢?
1、共享内存(shared Memory)是集成在GPU处理器芯片上的(on-chip),因此相比于存在于显存颗粒中的全局内存(global Memory)和本地内存(local Memory),它具有更高的传输带宽,一般情况下,共享内存的带宽大约是全局内存带宽的7-10倍。
2、共享内存的容量很小。根据NVIDIA官方文档的说法,在计算能力1.x的设备中,每一个流多处理器(Streaming Multiprocessor)上的共享内存容量为16KB。对于计算能力2.x、3.0及3.5的设备该参数为48KB。因此共享内存是稀有资源。
3、共享内存在物理上被划分为很多块,每一块被称为一个存储体(bank)。在同一时刻,CUDA设备可以同时访问多个存储体。因此,如果一次针对共享内存的访存操作需要读取n个地址,而这n个地址恰好分布在n个不同的存储体(bank)中,那么只需要一个存取周期就可以完成n个地址的访存任务了。对于计算能力1.x的设备,共享内存被平均划分为16个存储体。而对于计算能力2.x、3.0及3.5的设备此参数为32。在共享内存中,相邻两块32bit的数据分别属于相邻的两个存储体。存储体每两个时钟周期可以传输32位数据。
4、共享内存既可以静态分配,也可以动态分配。
从共享内存的这些特点中我们可以看出,它实际上相当于一个程序员可以操控的缓存(cache),下面,我们使用矩阵乘法的例子来说明如何有效使用共享内存。
首先,我们使用最直观的方法来完成矩阵乘法C = A x B:读取A的每一行和B的每一列,顺次完成计算任务。矩阵乘法的示意图如下所示:
CUDA编程指南阅读笔记(六)

下面是矩阵乘法的CUDA C主要实现代码:
  1. //Matricesarestoredinrow-majororder:
  2. //M(row,col)=*(M.elements+row*M.width+col)
  3. typedefstruct{
  4. intwidth;
  5. intheight;
  6. float*elements;
  7. }Matrix;
  8. //Threadblocksize
  9. #defineBLOCK_SIZE16
  10. //Forwarddeclarationofthematrixmultiplicationkernel
  11. __global__voidMatMulKernel(constMatrix,constMatrix,Matrix);
  12. //Matrixmultiplication-Hostcode
  13. //MatrixdimensionsareassumedtobemultiplesofBLOCK_SIZE
  14. voidMatMul(constMatrixA,constMatrixB,MatrixC){
  15. //LoadAandBtodevicememory
  16. Matrixd_A;
  17. d_A.width=A.width;d_A.height=A.height;
  18. size_tsize=A.width*A.height*sizeof(float);
  19. cudaMalloc(&d_A.elements,size);
  20. cudaMemcpy(d_A.elements,A.elements,size,cudaMemcpyHostToDevice);
  21. Matrixd_B;
  22. d_B.width=B.width;d_B.height=B.height;
  23. size=B.width*B.height*sizeof(float);
  24. cudaMalloc(&d_B.elements,size);
  25. cudaMemcpy(d_B.elements,B.elements,size,cudaMemcpyHostToDevice);
  26. //AllocateCindevicememory
  27. Matrixd_C;
  28. d_C.width=C.width;d_C.height=C.height;
  29. size=C.width*C.height*sizeof(float);
  30. cudaMalloc(&d_C.elements,size);
  31. //Invokekernel
  32. dim3dimBlock(BLOCK_SIZE,BLOCK_SIZE);
  33. dim3dimGrid(B.width/dimBlock.x,A.height/dimBlock.y);
  34. MatMulKernel<<<dimGrid,dimBlock>>>(d_A,d_B,d_C);
  35. //ReadCfromdevicememory
  36. cudaMemcpy(C.elements,d_c.elements,size,cudaMemcpyDeviceToHost);
  37. //Freedevicememory
  38. cudaFree(d_A.elements);
  39. cudaFree(d_B.elements);
  40. cudaFree(d_C.elements);
  41. }
  42. //MatrixmultiplicationkernelcalledbyMatMul()
  43. __global__voidMatMulKernel(MatrixA,MatrixB,MatrixC){
  44. //EachthreadcomputesoneelementofC
  45. //byaccumulatingresultsintoCvalue
  46. floatCvalue=0;
  47. introw=blockIdx.y*blockDim.y+threadIdx.y;
  48. intcol=blockIdx.x*blockDim.x+threadIdx.xl
  49. for(inte=0;e<A.width;++e)
  50. Cvalue+=A.elements[row*A.width+e]*B.elements[e*B.width+col];
  51. C.elements[row*C.width+col]=Cvalue;
  52. }
可以看出,为了计算矩阵C的任何一个元素,程序都需要从全局内存(global memory)中获得矩阵A的一行和矩阵B的一列。因此,完成这一计算矩阵A被读取了B.width次,矩阵B被读取了A.height次。
现在我们来使用共享内存(shared memory)实现矩阵乘法。假设矩阵C可以被划分为若干个较小的子方阵Csub,我们使用一个线程块(thread block)来负责某一子方阵的计算,线程块中的每一个线程(thread)正好负责子方阵Csub中一个元素的计算。这样划分后,任何一个结果子方阵Csub'(尺寸为block_size * block_size)都是与该方阵具有相同行索引的尺寸为A.width * block_size的A的子矩阵Asub和与该方阵具有相同列索引的尺寸为block_size * B.height的B的子矩阵Bsub相乘所得到。
为了匹配设备的计算资源,两个子矩阵Asub和Bsub被划分为尽可能多的分离的维度为block_size的子方阵,Csub的值便是这些子矩阵相乘后相加所得到的结果。子矩阵乘法的执行顺序都是首先将它们从全局内存(global memory)拷贝到共享内存(shared memory)(线程块中的每一个线程正好负责方阵一个元素的拷贝),然后由线程自己完成相应元素的计算任务,利用寄存器存储局部结果,最后将寄存器的内容与新得到的计算结果依此累加起来得到最终运算结果并将其传输到全局内存(global memory)中。
通过使用这种分治的计算策略,共享内存得到了很好的利用,采用这种方案计算完成时全局内存中矩阵A被访问的次数为B.width / block_size,矩阵B被访问的次数为A.height / block_size,很明显,这为我们节省了非常多的全局内存带宽。优化后的矩阵计算示意图如下所示:
CUDA编程指南阅读笔记(六)
为了提升计算效率,我们为类型Matrix增加了一个成员变量stride。__device__函数用来获得和设置子矩阵的元素。下面是优化后的代码:
  1. //Matricesarestoredinrow-majororder;
  2. //M(row,col)=*(M.elements+row*M.stride+col)
  3. typedefstruct{
  4. intwidth;
  5. intheight;
  6. intstride;
  7. float*elements;
  8. }Matrix;
  9. //Getamatrixelement
  10. __device__floatGetElement(constMatrixA,introw,intcol){
  11. returnA.elements[row*A.stride+col];
  12. }
  13. //Setamatrixelement
  14. __device__voidSetElement(MatrixA,introw,intcol,floatvalue){
  15. A.elements[row*A.stride+col]=value;
  16. }
  17. //GettheBLOCK_SIZExBLOCK_SIZEsub-matrixAsubofAthatis
  18. //locatedcolsub-matricestotherightandrowsub-matricesdown
  19. //fromtheupper-leftcornerofA
  20. __device__MatrixGetSubMatrix(MatrixA,introw,intcol){
  21. MatrixAsub;
  22. Asub.width=BLOCK_SIZE;
  23. Asub.height=BLOCK_SIZE;
  24. Asub.stride=A.stride;
  25. Asub.elements=&A.elements[A.stride*BLOCK_SIZE*row+BLOCK_SIZE*col];
  26. returnAsub;
  27. }
  28. //Threadblocksize
  29. #defineBLOCK_SIZE16
  30. //Forwarddeclarationofthematrixmultiplicationkernel
  31. __global__voidMatMulKernel(constMatrix,constMatrix,Matrix);
  32. //Matrixmultiplication-Hostcode
  33. //MatrixdimensionsareassumedtobemultiplesofBLOCK_SIZE
  34. voidMatMul(constMatrixA,constMatrixB,MatrixC){
  35. //LoadAandBtodevicememory
  36. Matrixd_A;
  37. d_A.width=d_A.stride=A.width;
  38. d_A.height=A.height;
  39. size_tsize=A.width*A.height*sizeof(float);
  40. cudaMalloc(&d_A.elements,size);
  41. cudaMemcpy(d_A.elements,A.elements,size,cudaMemcpyHostToDevice);
  42. Matrixd_B;
  43. d_B.width=d_B.stride=B.width;
  44. d_B.height=B.height;
  45. size=B.width*B.height*sizeof(float);
  46. cudaMalloc(&d_B.elements,size);
  47. cudaMemcpy(d_B.elements,B.elements,size,cudaMemcpyHostToDevice);
  48. //AllocateCindevicememory
  49. Matrixd_C;
  50. d_C.width=d_C.stride=C.width;
  51. d_C.height=C.height;
  52. size=C.width*C.height*sizeof(float);
  53. cudaMalloc(&d_C.elements,size);
  54. //Invokekernel
  55. dim3dimBlock(BLOCK_SIZE,BLOCK_SIZE);
  56. dim3dimGrid(B.width/dimBlock.x,A.height/dimBlock.y);
  57. MatMulKernel<<<dimGrid,dimBlock>>>(d_A,d_B,d_C);
  58. //ReadCfromdevicememory
  59. cudaMemcpy(C.elements,d_C.elements,size,cudaMemcpyDeviceToHost);
  60. //Freedevicememory
  61. cudaFree(d_A.elements);
  62. cudaFree(d_B.elements);
  63. cudaFree(d_C.elements);
  64. }
  65. //MatrixmultiplicationkernelcalledbyMatMul()
  66. __global__voidMatMulKernel(MatrixA,MatrixB,MatrixC){
  67. //Blockrowandcolumn
  68. intblockRow=blockIdx.y;
  69. intblockCol=blockIdx.x;
  70. //Eachthreadblockcomputesonesub-matrixCsubofC
  71. MatrixCsub=GetSubMatrix(C,blockRow,blockCol);
  72. //EachthreadcomputesoneelementofCsub
  73. //byaccumulatingresultsintoCvalue
  74. floatCvalue=0;
  75. //ThreadrowandcolumnwithinCsub
  76. introw=threadIdx.y;
  77. intcol=threadIdx.x;
  78. //Lookoverallthesub-matricesofAandBthatarerequiredtocomputeCsub
  79. //Multiplyeachpairofsub-matricestogetherandaccumulatetheresults
  80. for(intm=0;m<(A.width/BLOCK_SIZE);++m){
  81. //Getsub-matrixAsubofA
  82. MatrixAsub=GetSubMatrix(A,blockRow,m);
  83. //Getsub-matrixBsubofB
  84. MatrixBsub=GetSubMatrix(B,m,blockCol);
  85. //SharedmemoryusedtostoreAsubandBsubrespectively
  86. __shared__floatAs[BLOCK_SIZE][BLOCK_SIZE];
  87. __shared__floatBs[BLOCK_SIZE][BLOCK_SIZE];
  88. //LoadAsubandBsubfromdevicememorytosharedmemory
  89. //Eachthreadloadsoneelementofeachsub-matrix
  90. As[row][col]=GetElement(Asub,row,col);
  91. Bs[row][col]=GetElement(Bsub,row,col);
  92. //Synchronizetomakesurethesub-matricesareloaded
  93. //beforestartingthecomputation
  94. __syncthreads();
  95. //MultiplyAsubandBsubtogether
  96. for(inte=0;e<BLOCK_SIZE;++e)
  97. Cvalue+=As[row][e]*Bs[e][col];
  98. //Synchronizetomakesurethattheprecedingcomputationisdonebefore
  99. //loadingtwonewsub-matricesofAandBinthenextiteration
  100. __syncthreads();
  101. }
  102. //WriteCsubtodevicememory
  103. //Eachthreadwritesoneelement
  104. SetElement(Csub,row,col,Cvalue);
  105. }

同一个block中的线程在95行的for循环中获取到的Asub,Bsub,Csub是一样的,每个线程就负责Csub内元素的计算

http://blog.csdn.net/csgxy123/article/details/10018531

4.CUDA C语言编程接口

接上文。

4.3 CUDA C Runtime

4.3.3 共享内存(Shared Memory)

共享内存是CUDA设备中非常重要的一个存储区域,有效地使用共享内存可以充分利用CUDA设备的潜能,极大提升程序性能。那么,共享内存有哪些特点呢?
1、共享内存(shared Memory)是集成在GPU处理器芯片上的(on-chip),因此相比于存在于显存颗粒中的全局内存(global Memory)和本地内存(local Memory),它具有更高的传输带宽,一般情况下,共享内存的带宽大约是全局内存带宽的7-10倍。
2、共享内存的容量很小。根据NVIDIA官方文档的说法,在计算能力1.x的设备中,每一个流多处理器(Streaming Multiprocessor)上的共享内存容量为16KB。对于计算能力2.x、3.0及3.5的设备该参数为48KB。因此共享内存是稀有资源。
3、共享内存在物理上被划分为很多块,每一块被称为一个存储体(bank)。在同一时刻,CUDA设备可以同时访问多个存储体。因此,如果一次针对共享内存的访存操作需要读取n个地址,而这n个地址恰好分布在n个不同的存储体(bank)中,那么只需要一个存取周期就可以完成n个地址的访存任务了。对于计算能力1.x的设备,共享内存被平均划分为16个存储体。而对于计算能力2.x、3.0及3.5的设备此参数为32。在共享内存中,相邻两块32bit的数据分别属于相邻的两个存储体。存储体每两个时钟周期可以传输32位数据。
4、共享内存既可以静态分配,也可以动态分配。
从共享内存的这些特点中我们可以看出,它实际上相当于一个程序员可以操控的缓存(cache),下面,我们使用矩阵乘法的例子来说明如何有效使用共享内存。
首先,我们使用最直观的方法来完成矩阵乘法C = A x B:读取A的每一行和B的每一列,顺次完成计算任务。矩阵乘法的示意图如下所示:
CUDA编程指南阅读笔记(六)

下面是矩阵乘法的CUDA C主要实现代码:
  1. //Matricesarestoredinrow-majororder:
  2. //M(row,col)=*(M.elements+row*M.width+col)
  3. typedefstruct{
  4. intwidth;
  5. intheight;
  6. float*elements;
  7. }Matrix;
  8. //Threadblocksize
  9. #defineBLOCK_SIZE16
  10. //Forwarddeclarationofthematrixmultiplicationkernel
  11. __global__voidMatMulKernel(constMatrix,constMatrix,Matrix);
  12. //Matrixmultiplication-Hostcode
  13. //MatrixdimensionsareassumedtobemultiplesofBLOCK_SIZE
  14. voidMatMul(constMatrixA,constMatrixB,MatrixC){
  15. //LoadAandBtodevicememory
  16. Matrixd_A;
  17. d_A.width=A.width;d_A.height=A.height;
  18. size_tsize=A.width*A.height*sizeof(float);
  19. cudaMalloc(&d_A.elements,size);
  20. cudaMemcpy(d_A.elements,A.elements,size,cudaMemcpyHostToDevice);
  21. Matrixd_B;
  22. d_B.width=B.width;d_B.height=B.height;
  23. size=B.width*B.height*sizeof(float);
  24. cudaMalloc(&d_B.elements,size);
  25. cudaMemcpy(d_B.elements,B.elements,size,cudaMemcpyHostToDevice);
  26. //AllocateCindevicememory
  27. Matrixd_C;
  28. d_C.width=C.width;d_C.height=C.height;
  29. size=C.width*C.height*sizeof(float);
  30. cudaMalloc(&d_C.elements,size);
  31. //Invokekernel
  32. dim3dimBlock(BLOCK_SIZE,BLOCK_SIZE);
  33. dim3dimGrid(B.width/dimBlock.x,A.height/dimBlock.y);
  34. MatMulKernel<<<dimGrid,dimBlock>>>(d_A,d_B,d_C);
  35. //ReadCfromdevicememory
  36. cudaMemcpy(C.elements,d_c.elements,size,cudaMemcpyDeviceToHost);
  37. //Freedevicememory
  38. cudaFree(d_A.elements);
  39. cudaFree(d_B.elements);
  40. cudaFree(d_C.elements);
  41. }
  42. //MatrixmultiplicationkernelcalledbyMatMul()
  43. __global__voidMatMulKernel(MatrixA,MatrixB,MatrixC){
  44. //EachthreadcomputesoneelementofC
  45. //byaccumulatingresultsintoCvalue
  46. floatCvalue=0;
  47. introw=blockIdx.y*blockDim.y+threadIdx.y;
  48. intcol=blockIdx.x*blockDim.x+threadIdx.xl
  49. for(inte=0;e<A.width;++e)
  50. Cvalue+=A.elements[row*A.width+e]*B.elements[e*B.width+col];
  51. C.elements[row*C.width+col]=Cvalue;
  52. }
可以看出,为了计算矩阵C的任何一个元素,程序都需要从全局内存(global memory)中获得矩阵A的一行和矩阵B的一列。因此,完成这一计算矩阵A被读取了B.width次,矩阵B被读取了A.height次。
现在我们来使用共享内存(shared memory)实现矩阵乘法。假设矩阵C可以被划分为若干个较小的子方阵Csub,我们使用一个线程块(thread block)来负责某一子方阵的计算,线程块中的每一个线程(thread)正好负责子方阵Csub中一个元素的计算。这样划分后,任何一个结果子方阵Csub'(尺寸为block_size * block_size)都是与该方阵具有相同行索引的尺寸为A.width * block_size的A的子矩阵Asub和与该方阵具有相同列索引的尺寸为block_size * B.height的B的子矩阵Bsub相乘所得到。
为了匹配设备的计算资源,两个子矩阵Asub和Bsub被划分为尽可能多的分离的维度为block_size的子方阵,Csub的值便是这些子矩阵相乘后相加所得到的结果。子矩阵乘法的执行顺序都是首先将它们从全局内存(global memory)拷贝到共享内存(shared memory)(线程块中的每一个线程正好负责方阵一个元素的拷贝),然后由线程自己完成相应元素的计算任务,利用寄存器存储局部结果,最后将寄存器的内容与新得到的计算结果依此累加起来得到最终运算结果并将其传输到全局内存(global memory)中。
通过使用这种分治的计算策略,共享内存得到了很好的利用,采用这种方案计算完成时全局内存中矩阵A被访问的次数为B.width / block_size,矩阵B被访问的次数为A.height / block_size,很明显,这为我们节省了非常多的全局内存带宽。优化后的矩阵计算示意图如下所示:
CUDA编程指南阅读笔记(六)
为了提升计算效率,我们为类型Matrix增加了一个成员变量stride。__device__函数用来获得和设置子矩阵的元素。下面是优化后的代码:
  1. //Matricesarestoredinrow-majororder;
  2. //M(row,col)=*(M.elements+row*M.stride+col)
  3. typedefstruct{
  4. intwidth;
  5. intheight;
  6. intstride;
  7. float*elements;
  8. }Matrix;
  9. //Getamatrixelement
  10. __device__floatGetElement(constMatrixA,introw,intcol){
  11. returnA.elements[row*A.stride+col];
  12. }
  13. //Setamatrixelement
  14. __device__voidSetElement(MatrixA,introw,intcol,floatvalue){
  15. A.elements[row*A.stride+col]=value;
  16. }
  17. //GettheBLOCK_SIZExBLOCK_SIZEsub-matrixAsubofAthatis
  18. //locatedcolsub-matricestotherightandrowsub-matricesdown
  19. //fromtheupper-leftcornerofA
  20. __device__MatrixGetSubMatrix(MatrixA,introw,intcol){
  21. MatrixAsub;
  22. Asub.width=BLOCK_SIZE;
  23. Asub.height=BLOCK_SIZE;
  24. Asub.stride=A.stride;
  25. Asub.elements=&A.elements[A.stride*BLOCK_SIZE*row+BLOCK_SIZE*col];
  26. returnAsub;
  27. }
  28. //Threadblocksize
  29. #defineBLOCK_SIZE16
  30. //Forwarddeclarationofthematrixmultiplicationkernel
  31. __global__voidMatMulKernel(constMatrix,constMatrix,Matrix);
  32. //Matrixmultiplication-Hostcode
  33. //MatrixdimensionsareassumedtobemultiplesofBLOCK_SIZE
  34. voidMatMul(constMatrixA,constMatrixB,MatrixC){
  35. //LoadAandBtodevicememory
  36. Matrixd_A;
  37. d_A.width=d_A.stride=A.width;
  38. d_A.height=A.height;
  39. size_tsize=A.width*A.height*sizeof(float);
  40. cudaMalloc(&d_A.elements,size);
  41. cudaMemcpy(d_A.elements,A.elements,size,cudaMemcpyHostToDevice);
  42. Matrixd_B;
  43. d_B.width=d_B.stride=B.width;
  44. d_B.height=B.height;
  45. size=B.width*B.height*sizeof(float);
  46. cudaMalloc(&d_B.elements,size);
  47. cudaMemcpy(d_B.elements,B.elements,size,cudaMemcpyHostToDevice);
  48. //AllocateCindevicememory
  49. Matrixd_C;
  50. d_C.width=d_C.stride=C.width;
  51. d_C.height=C.height;
  52. size=C.width*C.height*sizeof(float);
  53. cudaMalloc(&d_C.elements,size);
  54. //Invokekernel
  55. dim3dimBlock(BLOCK_SIZE,BLOCK_SIZE);
  56. dim3dimGrid(B.width/dimBlock.x,A.height/dimBlock.y);
  57. MatMulKernel<<<dimGrid,dimBlock>>>(d_A,d_B,d_C);
  58. //ReadCfromdevicememory
  59. cudaMemcpy(C.elements,d_C.elements,size,cudaMemcpyDeviceToHost);
  60. //Freedevicememory
  61. cudaFree(d_A.elements);
  62. cudaFree(d_B.elements);
  63. cudaFree(d_C.elements);
  64. }
  65. //MatrixmultiplicationkernelcalledbyMatMul()
  66. __global__voidMatMulKernel(MatrixA,MatrixB,MatrixC){
  67. //Blockrowandcolumn
  68. intblockRow=blockIdx.y;
  69. intblockCol=blockIdx.x;
  70. //Eachthreadblockcomputesonesub-matrixCsubofC
  71. MatrixCsub=GetSubMatrix(C,blockRow,blockCol);
  72. //EachthreadcomputesoneelementofCsub
  73. //byaccumulatingresultsintoCvalue
  74. floatCvalue=0;
  75. //ThreadrowandcolumnwithinCsub
  76. introw=threadIdx.y;
  77. intcol=threadIdx.x;
  78. //Lookoverallthesub-matricesofAandBthatarerequiredtocomputeCsub
  79. //Multiplyeachpairofsub-matricestogetherandaccumulatetheresults
  80. for(intm=0;m<(A.width/BLOCK_SIZE);++m){
  81. //Getsub-matrixAsubofA
  82. MatrixAsub=GetSubMatrix(A,blockRow,m);
  83. //Getsub-matrixBsubofB
  84. MatrixBsub=GetSubMatrix(B,m,blockCol);
  85. //SharedmemoryusedtostoreAsubandBsubrespectively
  86. __shared__floatAs[BLOCK_SIZE][BLOCK_SIZE];
  87. __shared__floatBs[BLOCK_SIZE][BLOCK_SIZE];
  88. //LoadAsubandBsubfromdevicememorytosharedmemory
  89. //Eachthreadloadsoneelementofeachsub-matrix
  90. As[row][col]=GetElement(Asub,row,col);
  91. Bs[row][col]=GetElement(Bsub,row,col);
  92. //Synchronizetomakesurethesub-matricesareloaded
  93. //beforestartingthecomputation
  94. __syncthreads();
  95. //MultiplyAsubandBsubtogether
  96. for(inte=0;e<BLOCK_SIZE;++e)
  97. Cvalue+=As[row][e]*Bs[e][col];
  98. //Synchronizetomakesurethattheprecedingcomputationisdonebefore
  99. //loadingtwonewsub-matricesofAandBinthenextiteration
  100. __syncthreads();
  101. }
  102. //WriteCsubtodevicememory
  103. //Eachthreadwritesoneelement
  104. SetElement(Csub,row,col,Cvalue);
  105. }