对于 GTX960M 中大小为 300k 的结构数组，Thrust::sort 慢

Question

我目前正在研究一种 GPU 渲染算法，我需要对这个结构的数组进行排序：

struct RadiosityData {
    vec4 emission;
    vec4 radiosity;
    float nPixLight;
    float nPixCam;
    float __padding[2];
};

我正在使用以下代码对数组进行排序：

thrust::device_ptr<RadiosityData> dev_ptr = thrust::device_pointer_cast(GPUpointer_ssbo);
thrust::sort(dev_ptr, dev_ptr + N);

其中 GPUpointer_ssbo 是来自 cudaOpenGL 互操作的 GPU 指针，N 等于 ~300k。比较是通过：

__host__ __device__ bool operator<(const RadiosityData& lhs, const RadiosityData& rhs) { return (lhs.nPixCam > rhs.nPixCam); };

我的 GTX960M 的排序非常慢：没有排序，我的应用程序每帧大约 10 毫秒，而排序大约需要 35 毫秒。这意味着排序需要大约 25 毫秒。我正在用 VS-NSIGHT 测量执行时间

我知道这个问题可能是 GPU 同步问题，因为我在调用推力之前执行 OpenGL 操作。不过，我对这个说法并不信服，因为如果我用未排序的数组用OpenGL显示数据，总共还是需要10ms，这意味着OpenGL代码本身不存在同步问题。

这种“小”阵列的性能是预期的吗？有没有更好的 GPU 排序算法来解决这类问题？

------------编辑： 我正在使用默认的 VS2019 CUDA 命令编译发布，即：

驱动 API（NVCC 编译类型为 .cubin、.gpu 或 .ptx） 设置 CUDAFE_FLAGS=--sdk_dir "C:\Program Files (x86)\Windows Kits\10\" "C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v10.2\bin\nvcc.exe" --use-local-env -ccbin "C:\Program Files (x86)\Microsoft Visual Studio\2019\Community \VC\Tools\MSVC\14.26.28801\bin\HostX86\x64" -x cu --keep-dir x64\Release -maxrregcount=0 --machine 64 --compile -cudart static -o x64\Release\sortBufferCUDA。 cu.obj "C:\Users\Jose\Desktop\RealTimeDiffuseIlumination\OpenGL-avanzado\sortBufferCUDA.cu"

运行时 API（NVCC 编译类型是混合对象或 .c 文件） 设置 CUDAFE_FLAGS=--sdk_dir "C:\Program Files (x86)\Windows Kits\10\" "C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v10.2\bin\nvcc.exe" --use-local-env -ccbin "C:\Program Files (x86)\Microsoft Visual Studio\2019\Community \VC\Tools\MSVC\14.26.28801\bin\HostX86\x64" -x cu --keep-dir x64\Release -maxrregcount=0 --machine 64 --compile -cudart static -Xcompiler "/EHsc /nologo / fd /FS /Zi " -o x64\Release\sortBufferCUDA.cu.obj "C:\Users\Jose\Desktop\RealTimeDiffuseIlumination\OpenGL-avanzado\sortBufferCUDA.cu"

-------------编辑2：

以下是一个最小的工作示例：

#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>

#include <thrust/device_ptr.h>
#include <thrust/sort.h>
#include <thrust/execution_policy.h>
#include <thrust/extrema.h>
#include <cuda_runtime_api.h>
#include <cuda.h>
#include <thrust/device_vector.h>


struct RadiosityData {
    float emission[4];
    float radiosity[4];
    float nPixLight;
    float nPixCam;
    float __padding[2];
};

extern "C" void CUDAsort();

__host__ __device__ bool operator<(const RadiosityData& lhs, const RadiosityData& rhs) { return (lhs.nPixCam > rhs.nPixCam); };

int pri = 1;

thrust::device_vector<RadiosityData> dev;

void CUDAsort() {
    if (pri == 1) {
        pri = 0;
        dev.resize(300000);

    }
    thrust::sort(dev.begin(), dev.end());

}

int main()
{
    float time;
    cudaEvent_t start, stop;


    while (true) {
        cudaEventCreate(&start);
        cudaEventCreate(&stop);
        cudaEventRecord(start, 0);
        CUDAsort();

        cudaEventRecord(stop, 0);
        cudaEventSynchronize(stop);
        cudaEventElapsedTime(&time, start, stop);
        printf("Time to generate:  %3.1f ms \n", time);
    }
    return 0;
}

Answer 1

要求推力在排序时移动 48 字节结构当然是可能的，但可能不是最有效的方法。

我们可以尝试的是：

1. 将用于排序结构数组 (AoS) 的浮点值拉入 float 数组中
2. 创建一个索引数组来配合这个 0 1 2 3...
3. sort_by_key 浮点值，携带整数索引
4. 使用重新排列的索引数组，从输入到输出执行 AoS 的单个置换副本
5. 将输出数组复制回输入数组，以模拟“就地”排序

看起来工作量很大，但根据我的测试，它实际上要快一点：

$ cat t30.cu
#include <thrust/sort.h>
#include <thrust/device_vector.h>
#include <iostream>
#include <thrust/execution_policy.h>
#include <time.h>
#include <sys/time.h>
#include <cstdlib>
#define USECPSEC 1000000ULL
long long dtime_usec(unsigned long long start){

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

struct RadiosityData {
#ifdef USE_VEC
 float4 emission;
 float4 radiosity;
#else
 float emission[4];
 float radiosity[4];
#endif
        float nPixLight;
 float nPixCam;
        float __padding[2];
};

__global__ void copyKernel(RadiosityData *d, float *f, int *i, int n){
 int idx=threadIdx.x+blockDim.x*blockIdx.x;
 if (idx < n){
  f[idx] = d[idx].nPixCam;
  i[idx] = idx;}
}

__host__ __device__ bool operator<(const RadiosityData &lhs, const RadiosityData  &rhs) { return (lhs.nPixCam > rhs.nPixCam); };
struct my_sort_functor
{
 template <typename T1, typename T2>
 __host__ __device__ bool operator()(T1 lhs, T2 rhs) { return (lhs.nPixCam > rhs.nPixCam); };
};
const int N = 300000;
int main(){
 RadiosityData *GPUpointer_ssbo, *o;
 int sz = N*sizeof(RadiosityData);
 thrust::device_vector<RadiosityData> ii(N);
 GPUpointer_ssbo = thrust::raw_pointer_cast(ii.data());
 thrust::device_ptr<RadiosityData> dev_ptr = thrust::device_pointer_cast(GPUpointer_ssbo);
//method 1: ordinary thrust sort
 long long dt = dtime_usec(0);
 thrust::sort(dev_ptr, dev_ptr+N);
 cudaDeviceSynchronize();
 dt = dtime_usec(dt);
 std::cout << "ordinary sort time: " << dt/(float)USECPSEC << "s" << std::endl;
//method 2: reduced sort-and-copy
 cudaMalloc(&o, sz);
 thrust::device_ptr<RadiosityData> dev_optr = thrust::device_pointer_cast(o);
 for (int i = 0; i < N; i++) {RadiosityData q{0}; q.nPixCam = rand(); ii[i] = q;}
 float *d;
 int *k;
 cudaMalloc(&d, N*sizeof(float));
 cudaMalloc(&k, N*sizeof(int));
 thrust::device_ptr<int> dev_kptr = thrust::device_pointer_cast(k);
 cudaDeviceSynchronize();
 dt = dtime_usec(0);
 copyKernel<<<(N+511)/512, 512>>>(GPUpointer_ssbo, d, k, N);
 thrust::sort_by_key(thrust::device, d, d+N, k);
 thrust::copy(thrust::make_permutation_iterator(dev_ptr, dev_kptr), thrust::make_permutation_iterator(dev_ptr, dev_kptr+N), dev_optr);
 cudaMemcpy(GPUpointer_ssbo, o, sz, cudaMemcpyDeviceToDevice);
 cudaDeviceSynchronize();
 dt = dtime_usec(dt);
 std::cout << "sort+copy time: " << dt/(float)USECPSEC << "s" << std::endl;
}

$ nvcc -o t30 t30.cu -arch=sm_52
$ ./t30
ordinary sort time: 0.009527s
sort+copy time: 0.003143s
$ nvcc -o t30 t30.cu -arch=sm_52 -DUSE_VEC
$ ./t30
ordinary sort time: 0.004409s
sort+copy time: 0.002859s
$

（CUDA 10.1.105，GTX960，fedora 核心 29）

因此，我们观察到改进后的方法大约有 50% 或更好的加速。

如果您只想返回排序的前 M 个元素，使用这种解构复制方法，我们可以通过减少复制操作的大小来做一些进一步的改进。完全排序是在浮点数上完成的，但是在复制 AoS 结果时，只复制前 M 个值：

$ cat t30.cu
#include <thrust/sort.h>
#include <thrust/device_vector.h>
#include <iostream>
#include <thrust/execution_policy.h>
#include <time.h>
#include <sys/time.h>
#include <cstdlib>
#define USECPSEC 1000000ULL
long long dtime_usec(unsigned long long start){

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

struct RadiosityData {
#ifdef USE_VEC
 float4 emission;
 float4 radiosity;
#else
 float emission[4];
 float radiosity[4];
#endif
        float nPixLight;
 float nPixCam;
        float __padding[2];
};

__global__ void copyKernel(RadiosityData *d, float *f, int *i, int n){
 int idx=threadIdx.x+blockDim.x*blockIdx.x;
 if (idx < n){
  f[idx] = d[idx].nPixCam;
  i[idx] = idx;}
}

__host__ __device__ bool operator<(const RadiosityData &lhs, const RadiosityData  &rhs) { return (lhs.nPixCam > rhs.nPixCam); };
struct my_sort_functor
{
 template <typename T1, typename T2>
 __host__ __device__ bool operator()(T1 lhs, T2 rhs) { return (lhs.nPixCam > rhs.nPixCam); };
};
const int N = 300000;
const int M = 1000;  // identifies top-M values to be returned by sort
int main(){
 RadiosityData *GPUpointer_ssbo, *o;
 int sz = N*sizeof(RadiosityData);
 thrust::device_vector<RadiosityData> ii(N);
 GPUpointer_ssbo = thrust::raw_pointer_cast(ii.data());
 thrust::device_ptr<RadiosityData> dev_ptr = thrust::device_pointer_cast(GPUpointer_ssbo);
//method 1: ordinary thrust sort
 long long dt = dtime_usec(0);
 thrust::sort(dev_ptr, dev_ptr+N);
 cudaDeviceSynchronize();
 dt = dtime_usec(dt);
 std::cout << "ordinary sort time: " << dt/(float)USECPSEC << "s" << std::endl;
//method 2: reduced sort-and-copy
 cudaMalloc(&o, sz);
 thrust::device_ptr<RadiosityData> dev_optr = thrust::device_pointer_cast(o);
 for (int i = 0; i < N; i++) {RadiosityData q{0}; q.nPixCam = rand(); ii[i] = q;}
 float *d;
 int *k;
 cudaMalloc(&d, N*sizeof(float));
 cudaMalloc(&k, N*sizeof(int));
 thrust::device_ptr<int> dev_kptr = thrust::device_pointer_cast(k);
        cudaDeviceSynchronize();
 dt = dtime_usec(0);
        copyKernel<<<(N+511)/512, 512>>>(GPUpointer_ssbo, d, k, N);
 thrust::sort_by_key(thrust::device, d, d+N, k);
 thrust::copy_n(thrust::make_permutation_iterator(dev_ptr, dev_kptr), M, dev_optr);
 cudaMemcpy(GPUpointer_ssbo, o, M*sizeof(RadiosityData), cudaMemcpyDeviceToDevice);
        cudaDeviceSynchronize();
 dt = dtime_usec(dt);
 std::cout << "sort+copy time: " << dt/(float)USECPSEC << "s" << std::endl;
}


$ nvcc -o t30 t30.cu -arch=sm_52 -lineinfo -DUSE_VEC
$ ./t30
ordinary sort time: 0.004425s
sort+copy time: 0.001042s
$

其他几点说明：

1. 还观察到，当 4-float 量用矢量类型 (float4) 而非 4 元素数组表示时，推力可以更有效地处理 AoS。我怀疑这可以让编译器识别出更有效的结构复制方法。

2. 还请注意，根据我的测试，为正确的 GPU 架构（在我的情况下为 sm_52）编译似乎会带来一些小的改进。 YMMV。

我不声明此代码或我发布的任何其他代码的正确性。使用我发布的任何代码的任何人都需要自担风险。我只是声称我试图解决原始帖子中的问题，并提供一些解释。我并不是说我的代码没有缺陷，或者它适用于任何特定目的。使用（或不使用）风险自负。