如何使用 Thrust 对矩阵的行进行排序？

Posted 2023-03-23

【中文标题】如何使用 Thrust 对矩阵的行进行排序？

【英文标题】：How to use Thrust to sort the rows of a matrix?

【发布时间】：2015-03-24 20:45:08

【问题描述】：

我有一个 5000x500 的矩阵，我想用 cuda 分别对每一行进行排序。我可以使用arrayfire，但这只是thrust::sort的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>());
                    
                
            
        

有没有办法在推力中融合操作以使排序并行运行？事实上，我正在寻找一种将 for 循环迭代融合到其中的通用方法。

【问题讨论】：

你能适应How to normalize matrix columns in CUDA with max performance?中的方法吗？

我会尝试在对thrust::for_each的调用中嵌套对thrust::sort的调用。

我正在尝试理解这两种方法...谢谢。

好的！我放弃。打算用简单的方法来做。

【参考方案1】：

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

使用矢量化排序。如果要由嵌套 for 循环排序的区域不重叠，则可以使用 2 个背靠背稳定排序操作进行矢量化排序，如 here 所述。

推力 v1.8（可用于 CUDA 7 RC，或通过从 thrust github repository 直接下载包括 support for nesting thrust algorithms，通过在传递给另一个推力算法的自定义函子中包含推力算法调用。如果您使用thrust::for_each 操作来选择您需要执行的单个排序，您可以通过在传递给 thrust::for_each 的函子中包含 thrust::sort 操作来使用单个推力算法调用来运行这些排序。

李>

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

原始的 sort-in-a-loop 方法 矢量化/批量排序 嵌套排序

在每种情况下，我们都对相同的 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
$

注意事项：

这些结果会因 GPU 而异。 “嵌套”时间/方法在支持动态并行的 GPU 上可能会有很大差异，因为这会影响推力如何运行嵌套排序函数。要使用动态并行进行测试，请将编译开关从 -arch=sm_20 更改为 -arch=sm_35 -rdc=true -lcudadevrt 此代码需要 CUDA 7 RC。我用的是 Fedora 20。 嵌套排序方法也会从设备端分配，因此我们必须使用cudaDeviceSetLimit大幅增加设备分配堆。 如果您使用动态并行，并且根据您所运行的 GPU 类型，使用cudaDeviceSetLimit 保留的内存量可能需要增加 8 倍。

【讨论】：

非常感谢！您不仅回答了我的问题，还向我展示了如何正确地做很多事情。我必须补充一下，我看到了这个：***.com/questions/9037906/… 并意识到正常的推力排序会比其他任何东西都快（除了手动编码的并行排序方法）。我的数据是无符号的，我数据的 16 个最高有效位全为零。所以我只是将行号放在 16 msb 中并排序它会推进排序。

嗨，你知道如何沿着矩阵行进行 argsort 吗？这意味着找出矩阵中行元素的顺序而不是排序矩阵。