CUDA编程练习（十二）双调归并排序

Simple_Bright&Man

于 2024-06-04 20:51:11 发布

阅读量168

点赞数 3

分类专栏： CUDA编程文章标签： c++ python 算法

本文链接：https://blog.csdn.net/SimpleBright_Man/article/details/139442882

版权

CUDA编程专栏收录该内容

14 篇文章 0 订阅

订阅专栏

1 想法

难啃的东西

双调归并的核心思想是两段连续数组经过一个操作，合并为一个，性质不变。说起来简单，实际步骤会比较困难，因为涉及许多边界条件确定，复杂的升序和降序控制。

2 原理解释

2.1 双调归并排序

包含两个部分：如何双调归并，如何排序。

图1中，同字母索引的相互比较交换。

图2中，BM使得该区域变成尖尖对应的排序方式（升序/降序），逐渐变到最后。

Question：排序深度是否是log的平方？

3 代码

多block版本很复杂

3.1 单block版本

唯一需要注意的，是为什么2*threadIdx.x-(threadIdx.x & (stride-1))能将线程映射到对应的位置。你可以自己写一个对应关系看看。

// 交换器
__device__ void comparator(ELEM_TYPE &A, ELEM_TYPE &B, uint dir) {
    ELEM_TYPE t;
    // 交换各自独立
    if ((A>B) == dir) {  // dir=1就是升序
        t=A; A=B; B=t;
    }
}
// 单block双调排序
__global__ void bitonicSort_singleBlock(ELEM_TYPE* d_src, ELEM_TYPE* d_tar, uint n, uint dir) {
    __shared__ ELEM_TYPE temp[SHARED_SIZE];

    int tx = threadIdx.x;
    temp[tx] = d_src[tx];
    temp[tx+BLOCK_DIM] = d_src[tx+BLOCK_DIM];


    // 双调排序至全局双调
    for (uint size=2; size<n; size<<=1) {  // size是，当前要构造的双调规模
        // 确定当前排序次序
        uint curdir = dir ^ ((threadIdx.x & (size/2)) != 0);
        // 双调归并
        for (uint stride=size/2; stride>0; stride>>=1) {
            __syncthreads();
            uint pos = 2*threadIdx.x-(threadIdx.x & (stride-1));
            comparator(temp[pos], temp[pos+stride], curdir);
        }
    }
    // 全局双调后，把它变为单调，最终使用一次归并
    for (uint stride=n/2; stride>0; stride>>=1) {
        __syncthreads();
        uint pos = 2*threadIdx.x-(threadIdx.x & (stride-1));
        comparator(temp[pos], temp[pos+stride], dir);
    }
    __syncthreads();

    // 移动回去
    d_tar[tx] = temp[tx];
    d_tar[tx+BLOCK_DIM] = temp[tx+BLOCK_DIM];
}

3.2 多block版本（包括归并）

这个对应关系过于复杂，尤其是如何将线程对应到它应该执行的大小判断。这里只做代码分析：首先&(n-1)的意思是，向n取余数。comparatorI取值范围就是0~n/2-1。再和size/2取＆，意思就是你以size/2为单位跳跃，你是第奇数个还是偶数个。

为什么分两个函数？因为stride可以小于一个共享内存大小是，不必进行多block之间同步了。

__global__ void bitonicSort_multiBlock_intra(ELEM_TYPE* d_src, ELEM_TYPE* d_tar, uint n) {
    __shared__ ELEM_TYPE temp[SHARED_SIZE];

    uint tidx = blockIdx.x*blockDim.x*2 + threadIdx.x;
    uint tx = threadIdx.x;
    uint dir = blockIdx.x & 1;

    temp[tx] = d_src[tidx];
    temp[tx+BLOCK_DIM] = d_src[tidx+BLOCK_DIM];
    // 双调排序至共享内存双调
    for (uint size=2; size<SHARED_SIZE; size<<=1) {  // size是，当前要构造的双调规模
        // 确定当前排序次序
        uint curdir = dir ^ ((tx & (size/2)) != 0);
        // 双调归并
        for (uint stride=size/2; stride>0; stride>>=1) {
            __syncthreads();
            uint pos = 2*tx-(tx & (stride-1));
            comparator(temp[pos], temp[pos+stride], curdir);
        }
    }

    // 全局双调后，把它变为单调，最终使用一次归并
    for (uint stride=BLOCK_DIM; stride>0; stride>>=1) {
        __syncthreads();
        uint pos = 2*tx-(tx & (stride-1));
        comparator(temp[pos], temp[pos+stride], dir);
    }
    __syncthreads();

    // 移动回去
    d_tar[tidx] = temp[tx];
    d_tar[tidx+BLOCK_DIM] = temp[tx+BLOCK_DIM];
}

// 块间合并
__global__ void bitonicMergeGlobal(ELEM_TYPE* d_src, ELEM_TYPE* d_tar, uint n, uint size, uint stride, uint dir) {
    uint tidx = blockIdx.x*blockDim.x+threadIdx.x;
    uint curdir = dir ^ ( (tidx & (n/2-1) & (size/2)) != 0);
    uint pos = 2*tidx-(tidx&(stride-1));
    ELEM_TYPE A = d_src[pos], B = d_tar[pos+stride];
    comparator(A, B, curdir);
    d_tar[pos] = A;
    d_tar[pos+stride] = B;
}

__global__ void bitonicMergeShared(ELEM_TYPE* d_src, ELEM_TYPE* d_tar, uint n, uint size, uint dir) {
    __shared__ ELEM_TYPE ss[SHARED_SIZE];
    uint tx = threadIdx.x;
    d_src += blockIdx.x*SHARED_SIZE + threadIdx.x;
    d_tar += blockIdx.x*SHARED_SIZE + threadIdx.x;

    ss[tx] = d_src[0];
    ss[tx+SHARED_SIZE/2] = d_src[SHARED_SIZE/2];
    uint comparatorI = (blockIdx.x*blockDim.x+tx) & (n/2-1);
    uint curdir = dir ^ ((comparatorI&(size/2)) != 0);

    for(uint stride=SHARED_SIZE/2; stride>0; stride>>=1) {
        __syncthreads();

        uint pos = 2*tx-(tx & (stride-1));
        comparator(ss[pos], ss[pos+stride], curdir);
    }

    d_tar[0] = ss[tx];
    d_tar[SHARED_SIZE/2] = ss[tx+SHARED_SIZE/2];
}