CUDA推力排序或CUB :: DeviceRadixSort

Question

我有一個由float4數組表示的粒子池，其中w分量是粒子當前的生命週期[0，1]。

我需要根據粒子的壽命以降序對這個數組進行排序，這樣我就可以爲當前有多少粒子「活躍」（壽命大於0）保留一個精確計數器。我需要這個計數器，因爲當我需要激活更多的粒子（隨機發生）時，它將允許我索引到陣列中的正確位置。

我的微粒數組存儲在設備內存中，似乎我應該能夠對數組進行排序而不必將數組傳輸到主機內存。

我在網上找到例子並沒有太多的運氣，這說明我可以用Thrust或CUB做到這一點。另外，我對使用Thrust猶豫不決，因爲我不知道如何防止它退化爲合併排序（這比基數排序慢得多），因爲我需要根據w組件進行排序。至於CUB，我根本沒有找到任何資源來說明我如何做到這一點。

我還希望保存在w組件中的生命週期，因爲這使我的生活在我的代碼的其他部分更容易。

有沒有簡單的方法來做到這一點？謝謝你的幫助。

Answer 1

在任一幼崽或推力，我們可以排序的.w「鑰匙」只是，那種做一個鍵值，其中值只是一個線性遞增指數：

0, 1, 2, 3, ... 

然後我們可以使用合成重新排列索引序列以便在一個步驟中重新排列原始的float4陣列（有效地按.w排序）。這將允許您保留基數分選速度（無論是幼體還是推力），並且也可能相當有效，因爲float4數量只需移動/重新排列一次，而不是在分類操作過程中連續移動。

這裏是一個完整的推力示例，在32M元素上，演示了一個「普通」推力排序，使用函數來指定排序.w元素（sort_f4_w），接着是上述方法。在這種情況下，我的特定的設置（Fedora的20，CUDA 7，Quadro5000），第二個方法似乎是約5倍快：

$ cat t686.cu 
#include <iostream> 
#include <vector_types.h> 
#include <stdlib.h> 
#include <thrust/host_vector.h> 
#include <thrust/device_vector.h> 
#include <thrust/sort.h> 
#include <thrust/iterator/transform_iterator.h> 
#include <thrust/iterator/permutation_iterator.h> 
#include <thrust/sequence.h> 
#include <thrust/copy.h> 
#include <thrust/equal.h> 

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

unsigned long long dtime_usec(unsigned long long start){ 

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

#define DSIZE (32*1048576) 

struct sort_f4_w 
{ 
    __host__ __device__ 
    bool operator()(const float4 &a, const float4 &b) const { 
    return (a.w < b.w);} 
}; 
// functor to extract the .w element from a float4 
struct f4_to_fw : public thrust::unary_function<float4, float> 
{ 
    __host__ __device__ 
    float operator()(const float4 &a) const { 
    return a.w;} 
}; 
// functor to extract the .x element from a float4 
struct f4_to_fx : public thrust::unary_function<float4, float> 
{ 
    __host__ __device__ 
    float operator()(const float4 &a) const { 
    return a.x;} 
}; 


bool validate(thrust::device_vector<float4> &d1, thrust::device_vector<float4> &d2){ 
    return thrust::equal(thrust::make_transform_iterator(d1.begin(), f4_to_fx()), thrust::make_transform_iterator(d1.end(), f4_to_fx()), thrust::make_transform_iterator(d2.begin(), f4_to_fx())); 
} 


int main(){ 
    unsigned long long t1_time, t2_time; 
    float4 *mydata = new float4[DSIZE]; 
    for (int i = 0; i < DSIZE; i++){ 
    mydata[i].x = i; 
    mydata[i].y = i; 
    mydata[i].z = i; 
    mydata[i].w = rand()/(float)RAND_MAX;} 

    thrust::host_vector<float4> h_data(mydata, mydata+DSIZE); 
    // do once as a warm-up run, then report timings on second run 
    for (int i = 0; i < 2; i++){ 
    thrust::device_vector<float4> d_data1 = h_data; 
    thrust::device_vector<float4> d_data2 = h_data; 

    // first time sort using typical thrust approach 
    t1_time = dtime_usec(0); 
    thrust::sort(d_data1.begin(), d_data1.end(), sort_f4_w()); 
    cudaDeviceSynchronize(); 
    t1_time = dtime_usec(t1_time); 
    // now extract keys and create index values, sort, then rearrange 
    t2_time = dtime_usec(0); 
    thrust::device_vector<float> keys(DSIZE); 
    thrust::device_vector<int> vals(DSIZE); 
    thrust::copy(thrust::make_transform_iterator(d_data2.begin(), f4_to_fw()), thrust::make_transform_iterator(d_data2.end(), f4_to_fw()), keys.begin()); 
    thrust::sequence(vals.begin(), vals.end()); 
    thrust::sort_by_key(keys.begin(), keys.end(), vals.begin()); 
    thrust::device_vector<float4> result(DSIZE); 
    thrust::copy(thrust::make_permutation_iterator(d_data2.begin(), vals.begin()), thrust::make_permutation_iterator(d_data2.begin(), vals.end()), result.begin()); 
    cudaDeviceSynchronize(); 
    t2_time = dtime_usec(t2_time); 
    if (!validate(d_data1, result)){ 
     std::cout << "Validation failure " << std::endl; 
     } 
    } 
    std::cout << "thrust t1 time: " << t1_time/(float)USECPSEC << "s, t2 time: " << t2_time/(float)USECPSEC << std::endl; 
} 


$ nvcc -o t686 t686.cu 
$ ./t686 
thrust t1 time: 0.731456s, t2 time: 0.149959 
$