CUDA SDK例子分析（2）：transpose

Transpose是一个矩阵转置的例子，通过两个功能相同的核函数：transpose()和transpose_naive()展示了shared memory的优势。

Transpose.cu中是host端程序，与上节的template基本相同，请读者自己分析。需要注意的有三处：

1. 在测速之前，首先分别运行了一次transpose()和transpose_naive（），这样可以防止将CUDA的启动时间计入。

2. 上节的grid只有一个block，而本节中，grid和threads都是二维的。

3. 在调用transpose()和transpose_naive()时，<<<>>>中的参数只有两个，与上节的template相比，少了一个shared memory size，原因我们将在下面分析。

首先看看transpose_naive()：

 

__global__ void transpose_naive(float *odata, float* idata, int width, int height)

{

   unsigned int xIndex = blockDim.x * blockIdx.x + threadIdx.x;

   unsigned int yIndex = blockDim.y * blockIdx.y + threadIdx.y;

  

   if (xIndex < width && yIndex < height)

   {

       unsigned int index_in = xIndex + width * yIndex;

       unsigned int index_out = yIndex + height * xIndex;

       odata[index_out] = idata[index_in];

   }

}

xIndex和yIndex是根据线程在block中的位置 threadIdx.x, threadIdx.y 和线程所在 block 在整个 grid 中的标号 blockIdx.x 和 blockIdx.y ，以及 blcok 的长和宽 : blockDim.x ， blockDim.y 来计算线程在整个 grid 的位置 :

 

对矩阵进行转置的过程实际上就是从 index_in 位置读入数据，然后写到 index_out 处

 

以上程序的想法非常自然，似乎一切都是那么理所当然。的确，这个程序可以完成我们需要的矩阵转置功能，然而却忽略了一个重要的问题： GPU 对显存的 coalesced access 。判断是否 coalesced access 的简单原则之一是：对存储器的访问按照 threadIdx.x 连续。例如， index_in 展开以后的值是

width * yIndex + blockDim.x * blockIdx.x + threadIdx.x

     满足按行访问，符合 coalesced reading

而 index_out 的值实际上是：

height * xIndex + blockDim.y * blockIdx.y + threadIdx.y

则是在按行访问，造成了 non-coalesced writes 。

 

在 transpose() 中，这个问题通过使用 shared_memory 得到了解决：

__global__ void transpose(float *odata,

                          float *idata,

                          int width,

                          int height) // 指针外的其他参数，如 width 和 height 传入显卡会被存储到 shared memory 中

{

   __shared__ float block[(BLOCK_DIM+1)*BLOCK_DIM];// 在 template 中， shared 之前有一个 extern ，说明 shared memory size 由外部定义，因此方括号留空，而此处 block 的大小由方括号内的数字决定。

 

   unsigned int xBlock = __mul24(blockDim.x, blockIdx.x);

   unsigned int yBlock = __mul24(blockDim.y, blockIdx.y);

   unsigned int xIndex = xBlock + threadIdx.x;

   unsigned int yIndex = yBlock + threadIdx.y;

   unsigned int index_out, index_transpose;

 

   if (xIndex < width && yIndex < height)      // 保证内存访问不会超过矩阵边界

   {

       // load block into smem

       unsigned int index_in =

           __mul24(width, yIndex) + xIndex;// 线程需要读入的数据在矩阵中的位置

     

       unsigned int index_block =

           __mul24(threadIdx.y, BLOCK_DIM+1) + threadIdx.x;// 线程中要处理的数据在 shared memory 中的位置， __mule24 是快速 int 乘法，第二个参数是 BLOCK_DIM+1 而不是 BLOCK_DIM ，这是为了防止产生 bank conflict 。

     

       // load a block of data into shared memory

       block[index_block] = idata[index_in];// 将数据读入 shared memory

 

       index_transpose = __mul24(threadIdx.x, BLOCK_DIM+1) + threadIdx.y;

    

       index_out = __mul24(height, xBlock + threadIdx.y) +

           yBlock + threadIdx.x; // 向 global 写入的时候，仍然是按照 threadIdx.x 连续访问，保证了 coalesced writing 。

   }

   __syncthreads();

 

   if (xIndex < width && yIndex < height)      // 保证内存访问不会超过矩阵边界

   {

       // write it out (transposed) into the new location

       odata[index_out] = block[index_transpose];

   }

}

根据 nVidia 的注释， transpose 和 transpose_naive 的效率可以相差一个数量级以上！一段不太长的矩阵转置就包含了任务拆分， coalesced access ， bank conflict 等内容，可以说是大有乾坤。要写出高效的 CUDA 程序，还真是需要考虑周全。

 

源自：张舒Blog

来自 “ ITPUB博客 ” ，链接：http://blog.itpub.net/22785983/viewspace-619736/，如需转载，请注明出处，否则将追究法律责任。

转载于:http://blog.itpub.net/22785983/viewspace-619736/