CUDA SDK例子分析（2）：transpose

洛欣發表於2009-11-16
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/，如需轉載，請註明出處，否則將追究法律責任。

相關文章