求两向量距离的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/，如需转载，请注明出处，否则将追究法律责任。

转载于:http://blog.itpub.net/23057064/viewspace-676605/