如何使用Thrust来sortingmatrix的行？

我有一个5000x500matrix，我想用cuda分别对每一行进行sorting。 我可以使用arrayfire，但是这仅仅是一个for循环，它不应该是有效的。

https://github.com/arrayfire/arrayfire/blob/devel/src/backend/cuda/kernel/sort.hpp

for(dim_type w = 0; w < val.dims[3]; w++) { dim_type valW = w * val.strides[3]; for(dim_type z = 0; z < val.dims[2]; z++) { dim_type valWZ = valW + z * val.strides[2]; for(dim_type y = 0; y < val.dims[1]; y++) { dim_type valOffset = valWZ + y * val.strides[1]; if(isAscending) { thrust::sort(val_ptr + valOffset, val_ptr + valOffset + val.dims[0]); } else { thrust::sort(val_ptr + valOffset, val_ptr + valOffset + val.dims[0], thrust::greater<T>()); } } } } 

有没有一种方法可以将推力操作融合起来，从而实现并行运行？ 事实上，我正在寻找的是通用的方式来融合循环迭代。

Answer 1

我可以想到两种可能性，其中之一是@JaredHoberock已经提出的。 我不知道一种通用方法来融合推力中的循环迭代，但第二种方法是更一般的方法。 我的猜测是，在这种情况下，第一种方法是两种方法中较快的方法。

1. 使用vector化sorting。 如果要通过嵌套for循环进行sorting的区域不重叠，则可以使用2个背对背稳定sorting操作进行vector化sorting，如此处所述。

2. Thrust v1.8（通过CUDA 7 RC提供，或者通过直接从github库中下载，包括支持嵌套推力algorithm ，通过在自定义函子中包含一个推力algorithm调用传递给另一个推力algorithm。 thrust::for_each操作来select需要执行的各种sorting，可以通过在您传递给thrust::for_each的函子中包含thrust::sort操作，使用单个推式algorithm调用来运行这些sorting。

以下是3种方法之间的比较：

1. 原始的分类方法
2. vector化/批量sorting
3. 嵌套sorting

在每种情况下，我们正在对相同的16000套1000单位进行分类。

 $ cat t617.cu #include <thrust/device_vector.h> #include <thrust/device_ptr.h> #include <thrust/host_vector.h> #include <thrust/sort.h> #include <thrust/execution_policy.h> #include <thrust/generate.h> #include <thrust/equal.h> #include <thrust/sequence.h> #include <thrust/for_each.h> #include <iostream> #include <stdlib.h> #define NSORTS 16000 #define DSIZE 1000 int my_mod_start = 0; int my_mod(){ return (my_mod_start++)/DSIZE; } bool validate(thrust::device_vector<int> &d1, thrust::device_vector<int> &d2){ return thrust::equal(d1.begin(), d1.end(), d2.begin()); } struct sort_functor { thrust::device_ptr<int> data; int dsize; __host__ __device__ void operator()(int start_idx) { thrust::sort(thrust::device, data+(dsize*start_idx), data+(dsize*(start_idx+1))); } }; #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; } int main(){ cudaDeviceSetLimit(cudaLimitMallocHeapSize, (16*DSIZE*NSORTS)); thrust::host_vector<int> h_data(DSIZE*NSORTS); thrust::generate(h_data.begin(), h_data.end(), rand); thrust::device_vector<int> d_data = h_data; // first time a loop thrust::device_vector<int> d_result1 = d_data; thrust::device_ptr<int> r1ptr = thrust::device_pointer_cast<int>(d_result1.data()); unsigned long long mytime = dtime_usec(0); for (int i = 0; i < NSORTS; i++) thrust::sort(r1ptr+(i*DSIZE), r1ptr+((i+1)*DSIZE)); cudaDeviceSynchronize(); mytime = dtime_usec(mytime); std::cout << "loop time: " << mytime/(float)USECPSEC << "s" << std::endl; //vectorized sort thrust::device_vector<int> d_result2 = d_data; thrust::host_vector<int> h_segments(DSIZE*NSORTS); thrust::generate(h_segments.begin(), h_segments.end(), my_mod); thrust::device_vector<int> d_segments = h_segments; mytime = dtime_usec(0); thrust::stable_sort_by_key(d_result2.begin(), d_result2.end(), d_segments.begin()); thrust::stable_sort_by_key(d_segments.begin(), d_segments.end(), d_result2.begin()); cudaDeviceSynchronize(); mytime = dtime_usec(mytime); std::cout << "vectorized time: " << mytime/(float)USECPSEC << "s" << std::endl; if (!validate(d_result1, d_result2)) std::cout << "mismatch 1!" << std::endl; //nested sort thrust::device_vector<int> d_result3 = d_data; sort_functor f = {d_result3.data(), DSIZE}; thrust::device_vector<int> idxs(NSORTS); thrust::sequence(idxs.begin(), idxs.end()); mytime = dtime_usec(0); thrust::for_each(idxs.begin(), idxs.end(), f); cudaDeviceSynchronize(); mytime = dtime_usec(mytime); std::cout << "nested time: " << mytime/(float)USECPSEC << "s" << std::endl; if (!validate(d_result1, d_result3)) std::cout << "mismatch 2!" << std::endl; return 0; } $ nvcc -arch=sm_20 -std=c++11 -o t617 t617.cu $ ./t617 loop time: 8.51577s vectorized time: 0.068802s nested time: 0.567959s $ 

笔记：

1. 从GPU到GPU，这些结果会有很大的不同。
2. 在可以支持dynamic并行性的GPU上，“嵌套”时间/方法可能会有很大差异，因为这会影响嵌套sorting函数的推力。 要使用dynamic并行性进行testing，请将编译开关从-arch=sm_20 -arch=sm_35 -rdc=true -lcudadevrt
3. 此代码需要CUDA 7 RC。 我使用了Fedora 20。
4. 嵌套sorting方法也将从设备端分配，因此我们必须使用cudaDeviceSetLimit大幅增加设备分配堆。
5. 如果您正在使用dynamic并行机制，并且根据您正在运行的GPU的types，则可能需要增加一个额外的因子8来增加用cudaDeviceSetLimit保留的内存量。