Thrust::sort GTX960M 中大小为 300k 的结构数组速度慢

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 interop 的 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 同步问题，因为我在调用 thrust 之前执行 OpenGL 操作。尽管如此，我并不相信这个论点，因为如果我用未排序的数组用OpenGL显示数据，它仍然总共需要10ms，这意味着OpenGL代码本身没有同步问题。

这种 "small" 阵列的性能是否符合预期？有没有更好的GPU排序算法可以解决这类问题？

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

驱动程序API（NVCC 编译类型为 .cubin、.gpu 或 .ptx）设置 CUDAFE_FLAGS=--sdk_dir "C:\Program Files (x86)\Windows Kits\" "C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v10.2\bin\nvcc.exe" --use-local-env -ccbin "C:\Program Files (x86)\Microsoft Visual Studio19\Community\VC\Tools\MSVC.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\" "C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v10.2\bin\nvcc.exe" --use-local-env -ccbin "C:\Program Files (x86)\Microsoft Visual Studio19\Community\VC\Tools\MSVC.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

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

我们可以尝试的是：

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

看起来工作量很大，但根据我的测试，实际上速度更快了一点：

$ 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 core 29)

因此我们观察到改进方法的速度提高了大约 50% 或更好。

如果您只想return排序前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
$

其他一些注意事项：

还观察到，当 4-float 量用矢量类型 (float4) 而不是 4- 表示时，推力可以更有效地处理 AoS元素数组。我怀疑这允许编译器识别更有效的结构复制方法。
另请注意，根据我的测试，为正确的 GPU 架构（sm_52 在我的例子中）编译似乎有一点改进。 YMMV.

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

Thrust::sort GTX960M 中大小为 300k 的结构数组速度慢

Thrust::sort slow for array of structs of size 300k in GTX960M

opengl

cuda

thrust