求兩向量距離的CUDA實現

本程式中，我以a[N]和 b[N]代表兩個向量，其歐氏距離計算的序列Ｃ程式碼如下：

CODE:

dis = 0;
for(int i = 0; i < N; i++){
  dis += (a[i] - b[i]) * (a[i - b[i]);
}

dis = sqrt(dis);其中dis表示兩個向量的距離，sqrt表示求平方根函式。從程式碼可以看出，我們可以使用Ｎ個執行緒，每個執行緒計算a,b一維上資料差的平方，但是由於這個計算量太小，為了加大計算量，我採用了一個執行緒計算多個資料差的平方和，在這步後，我們要將所有執行緒資料加起來以得到最終的結果，考慮到ＣＵＤＡ不支援全域性同步而支援塊內同步，因此我們可以先得到每個塊的和，然後再將每個塊的和儲存到全域性儲存器d_temp中。最後再重新建立一個核心來將ｄ＿ｔｅｍｐ的資料求和，由於此時資料量小，因此只要一個塊就夠了。
        考慮到分支會極大的影響到CUDA程式的效能，因此我們將塊內各執行緒資料加和的方式是折半相加，此時可以利用高速的共享儲存器。折半相加的思想是：假設有n ＝ 2^k個執行緒，那麼第一次相加在0和n/2, 1和n/2+1,……，n/2-1 和n-1之間進行。此後執行緒0至n/2-1擁有的資料和等於之前塊內所有資料和，然後可以進行下一次折半相加。在再次折半相加之間為了保證資料一致性，必須使用__syncthreads()函式在塊內所有執行緒間同步。具體程式碼如下：

CODE:

template static __global__ void reduceMultiBlock(const T* __restrict__ d_a, const T* __restrict__ d_b, const int dim, T* __restrict__ d_temp){
//global thread id
1.unsigned int id = blockDim.x*blockIdx.x + threadIdx.x;
//iterate length
2.unsigned int iter = blockDim.x*gridDim.x;
//for temp distance, length = blockDim.x
3.extern __shared__ T s_dis[];
//use for when dim is shorter than iter
4.s_dis[threadIdx.x] = 0;

5.T temp;
6.T tempDis = (T)0;

7.for(int i = id; i < dim; i += iter){
8.  temp = d_a[i] - d_b[i];       
9.  tempDis += temp*temp;       
10.}

11.s_dis[threadIdx.x] = tempDis;

12.__syncthreads();

13.for(int i = (blockDim.x>>1); i > 0; i >>= 1){
14.  if(threadIdx.x < i)
15.                      s_dis[threadIdx.x] += s_dis[i + threadIdx.x];
16.                __syncthreads();
17.}
if(0 == threadIdx.x)
18.d_temp[blockIdx.x] = s_dis[0];
}為了使程式碼能夠同時應用於整形和浮點型，使用了模板機制。下面來具體分析：
        1.取得當前執行緒的索引，CUDA執行緒有兩種組織，一是網格，一是塊。
        2.取得整個網格內的執行緒數。
        3.動態宣告共享儲存器，其大小由核心呼叫的第三個引數決定，單位是位元組。
        4.初始化共享儲存器為0。
        5.-6.宣告兩個暫存器，並初始化一個為0。
        7.-10.遍歷兩個向量的所有維並將向量對應維差的平方加到對應執行緒的暫存器上，為了保證全域性儲存器的合併訪問，採用了迴圈讀取模式，迴圈讀取是指，假設一個有５個執行緒，１０個數，那個執行緒０就讀第０和第５個數，執行緒１讀第１個和第六個數，以此類推。
        11.-12.將塊內執行緒取得的對應維差的平方和存入共享儲存器並同步，同步的目的是防止有些執行緒可能先執行15.語句而取得了不正確結果。
        13.-17.將塊內共享儲存器的資料全加到共享儲存器陣列的第０個元素上。
        18. 將整個塊的資料和存入全域性儲存器。為什麼要做這一步的原因在於：共享儲存器並不持久，當核心執行完後，其就會被回收，而全域性儲存器是持久的。這類似於ＣＰＵ中記憶體和棧。

        在獲得每個塊的資料和並存入全域性儲存器後，下一步要做的就是將上一核心存入全域性儲存器的資料加起來得到最終的結果。由於內容和上一核心函式差不多，因此就不詳細說了。程式碼如下：

CODE:

template
static __global__ void reduceOneBlock(T* __restrict__ d_temp, const int len){
  unsigned int tid = threadIdx.x;
  T temp = (T)0;
//length == blockDim.x
  extern __shared__ T s_dis[];
//use for when len is shorter than d_temp's length
  s_dis[tid] = (T)0;

  for(int i = tid; i < len; i += blockDim.x){
    temp += d_temp[i];               
  }

  s_dis[tid] = temp;

  __syncthreads();

  for(int i = (blockDim.x>>1); i > 0; i >>= 1){
    if(tid < i)
                        s_dis[tid] += s_dis[i + tid];
   __syncthreads();
     }
     if(0 == tid)
     *d_temp = s_dis[0];
}經過這兩個核心的處理，d_temp[0]就是向量所有維差的平方和了。我們只要將其從視訊記憶體上覆制到記憶體並求其平方根就得到了這兩個向量的距離。

來自 “ ITPUB部落格 ” ，連結：http://blog.itpub.net/23057064/viewspace-676605/，如需轉載，請註明出處，否則將追究法律責任。