使用CUDA进行合并排序:小输入数组的高效实现

Question

使用CUDA进行合并排序:小输入数组的高效实现

M. *_*G. 1 arrays sorting merge cuda scientific-computing

我有以下问题:给定两个排序的数组A和B,我必须生成一个排序的数组C,其元素为A和B.

我找到了一些解决这个问题的解决方案,使用CUDA:Merge Path,例如http://www.cc.gatech.edu/~bader/papers/GPUMergePath-ICS2012.pdf

但是,它们的问题是由原始数组的大小给出的,至少有10k个元素.我有一个不同的问题.

我要合并的数组要小得多(最多1000个元素),并且复杂度由要完成的合并次数给出(10次幂的顺序为10,10 ^ 5个大小为1000的数组)相互合并).

他们的部分问题是将原始数组拆分为同等大小的并行处理的部分.我必须合并的数组足够小,完全适合GPU共享内存.

Thrust不是我的第一选择,因为我的过程的输出不是排序数组本身,而是带有元素的计算:所以我认为专门的内核应该比首先排序元素索引更快,然后使用它们进行计算.

该算法的串行版本如下所示:

i=0
j=0
k=0
T=4
while i<N and j<M:
    if A[i]<B[j]:
        start_i = max(0,i-T)
        C[k]=sum(A[start_i:i+1])
        i+=1
    else:
        start_j = max(0,j-T)
        C[k]=sum(B[start_j:j+1])
        j+=1
    k+=1

while i<N:
    start_i = max(0,i-T)
    C[k]=sum(A[start_i:i+1])
    i+=1
    k+=1
while j<M:
    start_j = max(0,j-T)
    C[k]=sum(B[start_j:j+1])
    j+=1
    k+=1

Run Code Online (Sandbox Code Playgroud)

如何利用CUDA功能来解决这个问题？

Answer 1

Rob*_*lla 6

任何CUDA计划的两个最重要的优化目标应该是:

暴露(充分)并行性
有效利用记忆力

在优化过程中肯定会考虑许多其他事项,但这些是首先要解决的两个最重要的事项.

合并操作(与合并排序不完全相同)乍一看是一种固有的串行操作.我们无法正确决定选择哪个项目,从A或B输入数组到下一个输出数组C,直到我们在C中进行了所有先前的选择.在这方面,合并算法(在此实现中) )使得很难公开并行性,问题中链接的论文几乎完全集中在该主题上.

本文描述的算法的目标是将两个输入矢量A和B分解成可以独立工作的多个较小的片段,以暴露并行性.特别是,目标是将所有SM保持在GPU中,并将所有SP保持在SM中.一旦完成了足够的工作分解,每个SP最终都会执行顺序合并(如本文所述):

合并阶段 - 每个核心合并使用与简单顺序合并相同的算法给出的两个子数组.

但是,正如您所指出的,您想要做的事情有些不同.您已经拥有许多阵列,并且希望对这些阵列执行独立的合并操作.由于您的阵列数大约为100,000,因此考虑将每个阵列映射到GPU SP(即线程)是足够的独立工作.这意味着我们就像在论文中一样,在每个核心/ SP /线程上使用简单的顺序合并.所以暴露并行性的问题在你的情况下已经完成(可能已经足够程度了).

在这一点上,我们可以考虑按原样实施.我稍后展示的代码将此作为比较的起点.然而,我们发现性能不是很好,这是因为合并算法从根本上具有依赖于数据的访问序列,因此(更难)安排在GPU上进行合并访问.本文的作者提出通过首先将数据(以合并方式)读入共享内存,然后让算法在共享内存中对其进行处理来缓解此问题,其中对无组织访问的惩罚较少.

我会提出一个不同的方法:

安排顺序合并算法,使A和B的每个元素只需要读一次
以列主要形式安排A,B和C的存储,而不是人们可能考虑的更"自然"的行主存储.这有效地转换了A,B和C向量的存储矩阵.这允许合并访问的改进,因为GPU线程在其各自的A和B向量上通过合并操作导航.它远非完美,但改进很大.

这是一个实现上述想法的工作示例,在每个线程中运行一个简单的顺序合并,每个线程将一个A向量与一个B向量合并:

$ cat t784.cu
#include <stdio.h>
#include <stdlib.h>
#include <thrust/sort.h>
#include <thrust/merge.h>

#define NUM_SETS 100000
#define DSIZE 100
typedef int mytype;

// for ascending sorted data
#define cmp(A,B) ((A)<(B))
#define nTPB 512
#define nBLK 128

#include <time.h>
#include <sys/time.h>
#define USECPSEC 1000000ULL

long long dtime_usec(unsigned long long start){

  timeval tv;
  gettimeofday(&tv, 0);
  return ((tv.tv_sec*USECPSEC)+tv.tv_usec)-start;
}

template <typename T>
__host__ __device__ void smerge(const T * __restrict__  a, const T * __restrict__ b, T * __restrict__ c, const unsigned len_a, const unsigned len_b, const unsigned stride_a = 1, const unsigned stride_b = 1, const unsigned stride_c = 1){

  unsigned len_c = len_a+len_b;
  unsigned nc = 0;
  unsigned na = 0;
  unsigned nb = 0;
  unsigned fa = (len_b == 0);
  unsigned fb = (len_a == 0);
  T nxta = a[0];
  T nxtb = b[0];
  while (nc < len_c){
    if (fa)      {c[stride_c*nc++] = nxta; na++; nxta = a[stride_a*na];}
    else if (fb) {c[stride_c*nc++] = nxtb; nb++; nxtb = b[stride_b*nb];}
    else if (cmp(nxta,nxtb)){
      c[stride_c*nc++] = nxta;
      na++;
      if (na == len_a) fb++;
      else nxta = a[stride_a*na];}
    else {
      c[stride_c*nc++] = nxtb;
      nb++;
      if (nb == len_b) fa++;
      else nxtb = b[stride_b*nb];}}
}



template <typename T>
__global__ void rmtest(const T * __restrict__  a, const T * __restrict__ b, T * __restrict__  c, int num_arr, int len){

  int idx=threadIdx.x+blockDim.x*blockIdx.x;

  while (idx < num_arr){
    int sel=idx*len;
    smerge(a+sel, b+sel, c+(2*sel), len, len);
    idx += blockDim.x*gridDim.x;}
}

template <typename T>
__global__ void cmtest(const T * __restrict__ a, const T * __restrict__ b, T * __restrict__ c, int num_arr, int len, int stride_a, int stride_b, int stride_c){
  int idx=threadIdx.x+blockDim.x*blockIdx.x;
  while (idx < num_arr){
    smerge(a+idx, b+idx, c+idx, len, len, stride_a, stride_b, stride_c);
    idx += blockDim.x*gridDim.x;}
}




template <typename T>
int rmvalidate(T *a, T *b, T *c, int num_arr, int len){

  T *vc = (T *)malloc(2*len*sizeof(T));
  for (int i = 0; i < num_arr; i++){
    thrust::merge(a+(i*len), a+((i+1)*len), b+(i*len), b+((i+1)*len), vc);
#ifndef TIMING
    for (int j = 0; j < len*2; j++)
      if (vc[j] != c[(i*2*len)+j]) {printf("rm mismatch i: %d, j: %d, was: %d, should be: %d\n", i, j, c[(i*2*len)+j], vc[j]); return 0;}
#endif
    }
  return 1;
}

template <typename T>
int cmvalidate(const T *c1, const T *c2, int num_arr, int len){
  for (int i = 0; i < num_arr; i++)
    for (int j = 0; j < 2*len; j++)
      if (c1[i*(2*len)+j] != c2[j*(num_arr)+i]) {printf("cm mismatch i: %d, j: %d, was: %d, should be: %d\n", i, j, c2[j*(num_arr)+i], c1[i*(2*len)+j]); return 0;}
  return 1;
}

int main(){


  mytype *h_a, *h_b, *h_c, *d_a, *d_b, *d_c;
  h_a = (mytype *)malloc(DSIZE*NUM_SETS*sizeof(mytype));
  h_b = (mytype *)malloc(DSIZE*NUM_SETS*sizeof(mytype));
  h_c = (mytype *)malloc(DSIZE*NUM_SETS*sizeof(mytype)*2);
  cudaMalloc(&d_a, (DSIZE*NUM_SETS+1)*sizeof(mytype));
  cudaMalloc(&d_b, (DSIZE*NUM_SETS+1)*sizeof(mytype));
  cudaMalloc(&d_c, DSIZE*NUM_SETS*sizeof(mytype)*2);
// test "row-major" storage
  for (int i =0; i<DSIZE*NUM_SETS; i++){
    h_a[i] = rand();
    h_b[i] = rand();}
  thrust::sort(h_a, h_a+DSIZE*NUM_SETS);
  thrust::sort(h_b, h_b+DSIZE*NUM_SETS);
  cudaMemcpy(d_a, h_a, DSIZE*NUM_SETS*sizeof(mytype), cudaMemcpyHostToDevice);
  cudaMemcpy(d_b, h_b, DSIZE*NUM_SETS*sizeof(mytype), cudaMemcpyHostToDevice);
  unsigned long gtime = dtime_usec(0);
  rmtest<<<nBLK, nTPB>>>(d_a, d_b, d_c, NUM_SETS, DSIZE);
  cudaDeviceSynchronize();
  gtime = dtime_usec(gtime);
  cudaMemcpy(h_c, d_c, DSIZE*NUM_SETS*2*sizeof(mytype), cudaMemcpyDeviceToHost);
  unsigned long ctime = dtime_usec(0);
  if (!rmvalidate(h_a, h_b, h_c, NUM_SETS, DSIZE)) {printf("fail!\n"); return 1;}
  ctime = dtime_usec(ctime);
  printf("CPU time: %f, GPU RM time: %f\n", ctime/(float)USECPSEC, gtime/(float)USECPSEC);
// test "col-major" storage
  mytype *ch_a, *ch_b, *ch_c;
  ch_a = (mytype *)malloc(DSIZE*NUM_SETS*sizeof(mytype));
  ch_b = (mytype *)malloc(DSIZE*NUM_SETS*sizeof(mytype));
  ch_c = (mytype *)malloc(DSIZE*NUM_SETS*sizeof(mytype));
  for (int i = 0; i < NUM_SETS; i++)
    for (int j = 0; j < DSIZE; j++){
      ch_a[j*NUM_SETS+i] = h_a[i*DSIZE+j];
      ch_b[j*NUM_SETS+i] = h_b[i*DSIZE+j];}
  cudaMemcpy(d_a, ch_a, DSIZE*NUM_SETS*sizeof(mytype), cudaMemcpyHostToDevice);
  cudaMemcpy(d_b, ch_b, DSIZE*NUM_SETS*sizeof(mytype), cudaMemcpyHostToDevice);
  gtime = dtime_usec(0);
  cmtest<<<nBLK, nTPB>>>(d_a, d_b, d_c, NUM_SETS, DSIZE, NUM_SETS, NUM_SETS, NUM_SETS );
  cudaDeviceSynchronize();
  gtime = dtime_usec(gtime);
  cudaMemcpy(ch_c, d_c, DSIZE*NUM_SETS*2*sizeof(mytype), cudaMemcpyDeviceToHost);
  if (!cmvalidate(h_c, ch_c, NUM_SETS, DSIZE)) {printf("fail!\n"); return 1;}

  printf("GPU CM time: %f\n", gtime/(float)USECPSEC);
  return 0;
}

$ nvcc -O3 -DTIMING -o t784 t784.cu
$ ./t784
CPU time: 0.030691, GPU RM time: 0.045814
GPU CM time: 0.002784
$

Run Code Online (Sandbox Code Playgroud)

笔记:

当内存组织是行主要时,GPU实际上比天真的单线程CPU代码慢.但对于列主要组织(往往会增加合并访问的机会),GPU代码比我的测试用例的CPU代码快约10倍.对于GPU MergePath 32位整数加速与x86串行合并,这个~10倍加速因子大致在加速因子的范围(~10-20x)内.
使用int与float数据类型相比,CPU时序有显着差异. int似乎更快(在CPU上)所以我在这里显示该版本.(这篇文章中也提到了这种差异.)
-DTIMING添加到编译命令的开关削减了第一个验证功能,以便它只进行CPU合并操作,以进行计时.
基本合并代码模板化,能够处理不同的数据类型,并进行参数化,以便可以在column-major或row major操作中使用.
为了简化演示,我已经免除了CUDA错误检查.但是,每当您遇到CUDA代码时遇到问题,都应该始终使用正确的cuda错误检查.

如何使用推力(正如我在评论中所建议的那样)？应该可以使用thrust :: merge与合适的设备/顺序执行策略,或多或少地模仿我上面所做的. 然而,推力期望向量是连续的,因此,在没有额外复杂性的情况下,它只能用于行主要情况,我们已经看到这种情况严重受到不良内存访问模式的惩罚.应该可以在推力中创建一组置换迭代器,这将允许列主要的跨步访问改善内存场景,但我没有追求.

归档时间：	10 年，3 月前
查看次数：	4649 次
最近记录：	10 年，3 月前