gpu排序__gpu

單機版的雙調排序可以參考 http://blog.csdn.net/sunmenggmail/article/details/42869235

還是這張圖片

基於cuda的雙調排序的思路是：

為每一個元素提供一個線程，如果大於1024個元素，還是提供1024個線程，這是因為__syncthreads只能作為block內的線程同步，而一個block最多有1024個線程，如果元素個數大於1024則每個線程可能就要負責一個以上的元素的比較

就上圖而言，一個矩形代表一次多線程的比較，那麼此圖僅需要6次比較，就可以有右邊的輸出。

#include <vector>#include <algorithm>#include <iostream>#include <time.h>#include <sys/time.h>#include <string.h>#include <math.h>#include <stdlib.h>#include <stdio.h>using namespace std;#define CHECK_EQ1(a,b) do { \if ((a) != (b)) { \ cout <<__FILE__<<" : "<< __LINE__<<" : check failed because "<<a<<"!="<<b<<endl;\cout << cudaGetErrorString(a) <<endl;\exit(1);\}\} while(0)#define CUDA_CHECK(condition)\do {\cudaError_t error = condition;\CHECK_EQ1(error, cudaSuccess);\} while(0)static __device__ __forceinline__ unsigned int __btflo(unsigned int word){    unsigned int ret;    asm volatile("bfind.u32 %0, %1;" : "=r"(ret) : "r"(word));//return the index of highest non-zero bit in a word; for example, 00000110, return 2    return ret;}//for > 1024__global__ void bigBinoticSort(unsigned int *arr, int len, unsigned int *buf) {unsigned len2 = 1 << (__btflo(len-1u) + 1);//unsigned int MAX = 0xffffffffu;unsigned id = threadIdx.x;if (id >= len2) return;unsigned iter = blockDim.x;for (unsigned i = id; i < len2; i += iter) {if (i >= len) {buf[i-len] = MAX;}}__syncthreads();int count = 0;for (unsigned k = 2; k <= len2; k*=2) {for (unsigned j = k >> 1; j > 0; j >>= 1) {for (unsigned i = id; i < len2; i += iter) {unsigned swapIdx = i ^ j;if (swapIdx > i) {unsigned myelem, other;if (i < len) myelem = arr[i];else myelem = buf[i-len];if (swapIdx < len) other = arr[swapIdx];else other = buf[swapIdx-len];bool swap = false;if ((i & k)==0 && myelem > other) swap = true;if ((i & k) == k && myelem < other) swap = true;if (swap) {if (swapIdx < len) arr[swapIdx] = myelem; else buf[swapIdx-len] = myelem;if (i < len) arr[i] = other;else buf[i-len] = other;}}}__syncthreads();}}}//for <= 1024__global__ void binoticSort(unsigned int *arr, int len) {__shared__  unsigned int buf[1024];buf[threadIdx.x] = (threadIdx.x < len ? arr[threadIdx.x] : 0xffffffffu);__syncthreads();for (unsigned k = 2; k <= blockDim.x; k*=2) {//buid k elements ascend or descendfor (unsigned j = k >> 1; j > 0; j >>= 1) {//merge longer binotic into shorter binoticunsigned swapIdx = threadIdx.x ^ j;unsigned myelem = buf[threadIdx.x];unsigned other = buf[swapIdx];__syncthreads();unsigned ascend = k * (swapIdx < threadIdx.x);unsigned descend = k * (swapIdx > threadIdx.x);//if I is front, swap is back; ascend = 0, descend = k//if I is back, swap is front; ascend = k, descend = 0;bool swap = false;if ((threadIdx.x & k) == ascend) {if (myelem > other) swap = true;}if ((threadIdx.x & k) == descend) {if (myelem < other) swap = true;}if (swap) buf[swapIdx] = myelem;__syncthreads();}}if (threadIdx.x < len) arr[threadIdx.x] = buf[threadIdx.x];}template<class T>inline void printVec(T *vec, int len) {for (int i = 0; i < len; ++i) cout <<vec[i] << "\t";cout << endl;}template<class T>inline void printVecg(T *gvec, int len) {T *vec = (T*)malloc(sizeof(T)*len);CUDA_CHECK(cudaMemcpy(vec,gvec,sizeof(T)*len,cudaMemcpyDeviceToHost));printVec(vec,len);free(vec);}void lineSize(int N, int &nblocks, int &nthreads) {if (N <= 1024) {nthreads = (N + 32 - 1)/32*32;//}else {nblocks = (N + 1024 -1)/1024;}}bool validate(unsigned *gvec, int len) {unsigned *vec = (unsigned*)malloc(sizeof(unsigned)*len);CUDA_CHECK(cudaMemcpy(vec,gvec,sizeof(unsigned)*len,cudaMemcpyDeviceToHost));for(int i = 1; i < len; ++i) {if (vec[i] <= vec[i-1]) return false;}return true;}inline int roundUpPower2(int v) {v--;v |= v >> 1;v |= v >> 2;v |= v >> 4;v |= v >> 8;v |= v >> 16;v++;return v;}int main(int argc, char *argv[]) {if (argc != 2) {cout << "len \n";return;}int len = atoi(argv[1]);unsigned int *arr = (unsigned int*)malloc(sizeof(unsigned int)*len);for (int i = 0; i < len; ++i) arr[i] = i;srand((unsigned int)time(NULL));for (int i = len; i >= 2; --i) {int j = rand() % i;swap(arr[i-1], arr[j]);}unsigned* debug;CUDA_CHECK(cudaMalloc((void**)&debug, sizeof(unsigned)*1000));unsigned int* darr, *buf;CUDA_CHECK(cudaMalloc((void**)&darr, sizeof(unsigned int)*len));CUDA_CHECK(cudaMalloc((void**)&buf, sizeof(unsigned int)*len));CUDA_CHECK(cudaMemcpy(darr, arr, sizeof(unsigned int)*len, cudaMemcpyHostToDevice));bigBinoticSort<<<1,1024>>>(darr,len, buf);CUDA_CHECK(cudaPeekAtLastError());CUDA_CHECK(cudaDeviceSynchronize());if (validate(darr, len))cout << "yes\n";elsecout << "no\n";return 1;}

演算法有兩個雙調排序實現，一個用於小於1024個元素，用到了共用記憶體加快訪問速度，但是如果真要排序1024以下的元素，建議還是用cpu版本的快排吧，gpu的在速度上並沒有明顯的優勢，甚至還比cpu慢

如果大於1024元素，就採用另一種方法。這種方法的缺點也是很明顯的，就是不管再多的元素，只能用一個block進行計算，而一個block最多隻能用1024個線程，估計在一萬個元素以內的話，這個方法是gpu上最快的。

經過本人測試，包括thrust的sort(基數排序)， 只有元素數量超過5000個，gpu上的排序演算法才有明顯的優勢。10萬左右的元素，gpu上的排序演算法比cpu有一百倍的提速。

下面會介紹在gpu上進行快速排序。gpu快速排序可以處理非常大的資料，但是會有遞迴深度的限制，當超過遞迴深度時，就可以調用上面所講的雙調排序進行處理。測試表明，速度比thrust還是快一點

gpu上的快排主要參考範例 NVIDIA_CUDA-6.5_Samples/6_Advanced/cdpAdvancedQuicksort

快排只處理大於1024個元素的數組，然後將其分隔為左右兩個子數組，如果子數組長度大於1024則繼續動態遞迴調用快排，如果小於1024則動態調用雙調排序。如果快排的遞迴深度已經超過最大遞迴深度（cuda最大嵌套深度64，但是還受限於每一級所使用的記憶體大小），則直接調用雙調排序。

這段程式的最精彩的地方在於分隔函數

將數組按照warp大小進行分隔，每個warp處理32個元素，通過全域的atomicAdd函數，分別獲得warp內的小於和大於pivot數在數組的位移地址，注意在同一個warp內，這個位移地址是一樣的，然後每個線程將自己的元素放到位移地址，這樣就完成了分割

需要注意的是，這個快排不是in-place的，又涉及到遞迴調用，所以還得處理原數組和緩衝區的調換

由於cuda沒有顯式的鎖，此方法採用了一種特殊的迴圈隊列，本人認為在極端情況下，可能會出現問題

（這裡的代碼有錯，沒有處理原數組和緩衝區的調換，只是協助理解。正確的代碼請參考Samples裡的）

#define QSORT_BLOCKSIZE_SHIFT   9#define QSORT_BLOCKSIZE         (1 << QSORT_BLOCKSIZE_SHIFT)#define BITONICSORT_LEN         1024            // Must be power of 2!#define QSORT_MAXDEPTH          16              // Will force final bitonic stage at depth QSORT_MAXDEPTH+1#define QSORT_STACK_ELEMS   1*1024*1024 // One million stack elements is a HUGE number.typedef struct __align__(128) qsortAtomicData_t{    volatile unsigned int lt_offset;    // Current output offset for <pivot    volatile unsigned int gt_offset;    // Current output offset for >pivot    volatile unsigned int sorted_count; // Total count sorted, for deciding when to launch next wave    volatile unsigned int index;        // Ringbuf tracking index. Can be ignored if not using ringbuf.} qsortAtomicData;////////////////////////////////////////////////////////////////////////////////// A ring-buffer for rapid stack allocation////////////////////////////////////////////////////////////////////////////////typedef struct qsortRingbuf_t{    volatile unsigned int head; //1        // Head pointer - we allocate from here    volatile unsigned int tail; //0        // Tail pointer - indicates last still-in-use element    volatile unsigned int count;//0        // Total count allocated    volatile unsigned int max;  //0        // Max index allocated    unsigned int stacksize;  //           // Wrap-around size of buffer (must be power of 2)    volatile void *stackbase;           // Pointer to the stack we're allocating from} qsortRingbuf;/*for cuda has no lock, so we have to do like this:if alloc , ++headif free , ++tailso [tail, head) contains alloced chunks; head point to the next free chunkcount record the number of chunks had freewe have n chunks, but the index of a chunk is increase when re-allocmax record the maximum index of the free chunksonly if the chunks before max are all free, aka, max == count, we can alter tail value */template<class T>static __device__ void ringbufFree(qsortRingbuf *ringbuf, T *data) {unsigned index = data->index;unsigned count = atomicAdd((unsigned*)&(ringbuf->count), 1) + 1;unsigned max = atomicMax((unsigned*)&(ringbuf->max), index + 1);if (max < (index + 1)) max = index + 1;if (max == count) {atomicMax((unsigned*)&(ringbuf->tail), count);}}template<class T>static __device__ T* ringbufAlloc(qsortRingbuf *ringbuf) {unsigned int loop = 10000;while (((ringbuf->head - ringbuf->tail) >= ringbuf->stacksize) && (loop-- > 0));if (loop == 0) return NULL;unsigned index = atomicAdd((unsigned*)&ringbuf->head, 1);T *ret = (T*)(ringbuf->stackbase) + (index & (ringbuf->stacksize - 1));ret->index = index;return ret;}__global__ void qsort_warp(unsigned *indata,                           unsigned *outdata,                           unsigned int offset,//0                           unsigned int len,//                           qsortAtomicData *atomicData,//stack                           qsortRingbuf *atomicDataStack,//ringbuf                           unsigned int source_is_indata,//true                           unsigned int depth) {//printf("depth = %d", depth);    // Find my data offset, based on warp ID    unsigned int thread_id = threadIdx.x + (blockIdx.x << QSORT_BLOCKSIZE_SHIFT);    //unsigned int warp_id = threadIdx.x >> 5;   // Used for debug only    unsigned int lane_id = threadIdx.x & (warpSize-1);// %32    // Exit if I'm outside the range of sort to be done    if (thread_id >= len)        return;    //    // First part of the algorithm. Each warp counts the number of elements that are    // greater/less than the pivot.    //    // When a warp knows its count, it updates an atomic counter.    //    // Read in the data and the pivot. Arbitrary pivot selection for now.    unsigned pivot = indata[offset + len/2];    unsigned data  = indata[offset + thread_id];    // Count how many are <= and how many are > pivot.    // If all are <= pivot then we adjust the comparison    // because otherwise the sort will move nothing and    // we'll iterate forever.    unsigned int greater = (data > pivot);    unsigned int gt_mask = __ballot(greater);//Evaluate predicate for all active threads of the warp and return an integer whose Nth bit is set if and only if predicate evaluates to non-zero for the Nth thread of the warp and the Nth thread is active.if (gt_mask == 0) {greater = (data >= pivot);gt_mask = __ballot(greater);}unsigned lt_mask = __ballot(!greater);unsigned gt_count = __popc(gt_mask);//count number of 1 in a warp;unsigned lt_count = __popc(lt_mask);//only thread 0 in warp calc//find 2 new positions for this warpunsigned lt_oft, gt_oft;if (lane_id == 0) {if (lt_count > 0)lt_oft = atomicAdd((unsigned*)&atomicData->lt_offset, lt_count);//atomicAdd return old value, not the newer//all the warps will syn call thisif (gt_count > 0)gt_oft = len - (atomicAdd((unsigned*) &atomicData->gt_offset, gt_count) + gt_count);//printf("depth = %d\n", depth);//printf("pivot = %u\n", pivot);//printf("lt_count %u lt_oft %u gt_count %u gt_oft %u  atomicDataGtOffset %u\n", lt_count,lt_oft, gt_count,gt_oft, atomicData->gt_offset);}lt_oft = __shfl((int)lt_oft, 0);gt_oft = __shfl((int)gt_oft, 0);//Everyone pulls the offsets from lane 0__syncthreads(); // Now compute my own personal offset within this. I need to know how many    // threads with a lane ID less than mine are going to write to the same buffer    // as me. We can use popc to implement a single-operation warp scan in this case.    unsigned lane_mask_lt;    asm("mov.u32 %0, %%lanemask_lt;" : "=r"(lane_mask_lt));//bits set in positions less than the thread's lane number the warpunsigned my_mask = greater ? gt_mask : lt_mask;unsigned my_oft = __popc(my_mask & lane_mask_lt);////move datamy_oft += greater ? gt_oft : lt_oft;outdata[offset + my_oft] = data;__syncthreads();//if (lane_id == 0) printf("pivot = %d", pivot);if (lane_id == 0) {/*if (blockIdx.x == 0) {printf("depth = %d\n", depth);for (int i = 0; i < len; ++i)printf("%u ", outdata[offset+i]);printf("\n");}*/unsigned mycount = lt_count + gt_count;//we are the last warp if (atomicAdd((unsigned*)&atomicData->sorted_count, mycount) + mycount == len) {unsigned lt_len = atomicData->lt_offset;unsigned gt_len = atomicData->gt_offset;cudaStream_t lstream, rstream;            cudaStreamCreateWithFlags(&lstream, cudaStreamNonBlocking);            cudaStreamCreateWithFlags(&rstream, cudaStreamNonBlocking);ringbufFree<qsortAtomicData>(atomicDataStack, atomicData);if (lt_len == 0) return;//////////////////////////// leftif (lt_len > BITONICSORT_LEN) {if (depth >= QSORT_MAXDEPTH) {bigBinoticSort<<<1, BITONICSORT_LEN,0, rstream>>>(outdata + offset, lt_len, indata + offset);}else {if ((atomicData = ringbufAlloc<qsortAtomicData>(atomicDataStack)) == NULL)                        printf("Stack-allocation error. Failing left child launch.\n");else {atomicData->lt_offset = atomicData->gt_offset = atomicData->sorted_count = 0;                        unsigned int numblocks = (unsigned int)(lt_len+(QSORT_BLOCKSIZE-1)) >> QSORT_BLOCKSIZE_SHIFT;                        qsort_warp<<< numblocks, QSORT_BLOCKSIZE, 0, lstream >>>(outdata, indata, offset, lt_len, atomicData, atomicDataStack, true, depth+1);}}}else if (lt_len > 1) {unsigned int bitonic_len = 1 << (__btflo(lt_len-1U)+1);                binoticSort<<< 1, bitonic_len, 0, lstream >>>(outdata + offset,lt_len);}////////////////////////// rightif (gt_len > BITONICSORT_LEN) {if (depth >= QSORT_MAXDEPTH)bigBinoticSort<<<1, BITONICSORT_LEN,0, rstream>>>(outdata + offset + lt_len, gt_len, indata + offset + lt_len);else {if ((atomicData = ringbufAlloc<qsortAtomicData>(atomicDataStack)) == NULL)                        printf("Stack allocation error! Failing right-side launch.\n");else {atomicData->lt_offset = atomicData->gt_offset = atomicData->sorted_count = 0;                        unsigned int numblocks = (unsigned int)(gt_len+(QSORT_BLOCKSIZE-1)) >> QSORT_BLOCKSIZE_SHIFT;                        qsort_warp<<< numblocks, QSORT_BLOCKSIZE, 0, rstream >>>(outdata, indata, offset+lt_len, gt_len, atomicData, atomicDataStack, true, depth+1);}}}else if (gt_len > 1) {unsigned int bitonic_len = 1 << (__btflo(gt_len-1U)+1);                binoticSort<<< 1, bitonic_len, 0, rstream >>>(outdata + offset + lt_len,gt_len);}}}}void runqsort(unsigned *gpudata, unsigned *scratchdata, unsigned int count, cudaStream_t stream) {unsigned int stacksize = QSORT_STACK_ELEMS;//1*1024*1024    // This is the stack, for atomic tracking of each sort's status    qsortAtomicData *gpustack;    CUDA_CHECK(cudaMalloc((void **)&gpustack, stacksize * sizeof(qsortAtomicData)));    CUDA_CHECK(cudaMemset(gpustack, 0, sizeof(qsortAtomicData)));     // Only need set first entry to 0    // Create the memory ringbuffer used for handling the stack.    // Initialise everything to where it needs to be.    qsortRingbuf buf;    qsortRingbuf *ringbuf;    CUDA_CHECK(cudaMalloc((void **)&ringbuf, sizeof(qsortRingbuf)));    buf.head = 1;           // We start with one allocation    buf.tail = 0;    buf.count = 0;    buf.max = 0;    buf.stacksize = stacksize;    buf.stackbase = gpustack;    CUDA_CHECK(cudaMemcpy(ringbuf, &buf, sizeof(buf), cudaMemcpyHostToDevice));if (count > BITONICSORT_LEN)//1024    {//QSORT_BLOCKSIZE = 2^9 = 512        unsigned int numblocks = (unsigned int)(count+(QSORT_BLOCKSIZE-1)) >> QSORT_BLOCKSIZE_SHIFT;        qsort_warp<<< numblocks, QSORT_BLOCKSIZE, 0, stream >>>(gpudata, scratchdata, 0U, count, gpustack, ringbuf, true, 0);   }    else    {        binoticSort<<< 1, BITONICSORT_LEN >>>(gpudata, count);CUDA_CHECK(cudaMemcpy(scratchdata, gpudata, sizeof(unsigned)*count, cudaMemcpyDeviceToDevice));    }cudaDeviceSynchronize();}