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

Question

我有一个5000x500的矩阵，想要使用cuda将每一行分别排序。我可以使用arrayfire，但这只是一个对thrust::sort进行循环的过程，效率不高。

参考链接：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>());
                    }
                }
            }
        }

有没有一种方法可以在thrust中融合操作，以便使排序并行运行？实际上，我正在寻找一种将for循环迭代融合成的通用方法。

- amir

---

你能否适应如何在CUDA中以最大性能规范化矩阵列？的方法？ - Vitality

我会尝试在thrust::for_each中嵌套调用thrust::sort。 - Jared Hoberock

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

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

Answer 1

我可以想到两种可能性，其中一种已经被@JaredHoberock建议过了。我不知道在thrust中融合for循环迭代的通用方法，但第二种方法是更通用的方法。我猜在这种情况下第一种方法会更快。

1. 使用向量化排序。如果要按嵌套的for循环排序的区域没有重叠部分，您可以使用2个背靠背的稳定排序操作进行向量化排序，如此处讨论的那样。
2. Thrust v1.8（可通过CUDA 7 RC或直接从thrust github存储库下载）包括通过在传递给另一个thrust算法的自定义函数对象中包含thrust算法调用来支持嵌套thrust算法的支持。如果使用thrust::for_each操作选择需要执行的单个排序，您可以通过将thrust::sort操作包含在传递给thrust::for_each的函数对象中，使用单个thrust算法调用运行这些排序。

这里是三种方法的完全比较：

1. 原始的循环排序方法
2. 向量化/批处理排序
3. 嵌套排序

在每种情况下，我们都要对相同的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不同而显著不同。
2. 如果GPU支持动态并行处理，则“嵌套”时间/方法可能会在很大程度上变化，因为这将影响Thrust运行嵌套排序函数的方式。要使用动态并行处理进行测试，请将编译开关从-arch=sm_20更改为-arch=sm_35 -rdc=true -lcudadevrt
3. 此代码需要CUDA 7 RC。我使用的是Fedora 20。
4. 嵌套排序方法也将从设备端分配内存，因此我们必须使用cudaDeviceSetLimit大幅增加设备分配堆。
5. 如果您正在使用动态并行处理，并且取决于您所运行的GPU类型，则可能需要将cudaDeviceSetLimit保留的内存量增加8倍。

- Robert Crovella

---

非常感谢！你不仅回答了我的问题，还向我展示了如何正确地做很多事情。我必须补充一点，我看到了这个链接：https://dev59.com/C1_Va4cB1Zd3GeqPWsV9，并意识到普通的thrust排序会比其他任何方法都要快（除了手写的并行排序方法）。我的数据是无符号的，而且我的数据的16个最高有效位都是零。所以我只需将行号放在16个最高有效位中，并使用thrust排序即可。 - amir

你好，你知道如何在矩阵的行上执行argsort吗？这意味着要找出矩阵行中元素的顺序，而不是排序后的矩阵。 - coin cheung