为什么 Thrust Library 的交集返回意外结果？

Question

我正在使用库 Thrust 来获取两个较大整数集的交集。在使用 2 个小输入的测试中，我得到了正确的结果，但是当我使用具有 10^8 和 65535*1024 元素的两组时，我得到了一个空集。谁能解释一下这个问题？将前两个变量更改为较小的值，推力返回预期的交集。我的代码如下。

#include <thrust/set_operations.h>
#include <thrust/device_vector.h>
#include <thrust/device_ptr.h>
#include <iostream>
#include <stdio.h>


int main() {
    int sizeArrayLonger = 100*1000*1000;
    int sizeArraySmaller = 65535*1024;
    int length_result = sizeArraySmaller;    
    int* list = (int*) malloc(4*sizeArrayLonger);
    int* list_smaller = (int*) malloc(4*sizeArraySmaller);
    int* result = (int*) malloc(4*length_result);

    int* list_gpu;
    int* list_smaller_gpu;
    int* result_gpu;

    // THE NEXT TWO FORS TRANSFORMS THE SMALLER ARRAY IN A SUBSET OF THE LARGER ARRAY
    for (int i=0; i < sizeArraySmaller; i++) {
        list_smaller[i] = i+1;
        list[i] = i+1;
    }
    for (int i=sizeArraySmaller; i < sizeArrayLonger; i++) {
        list[i] = i+1;
    }

    cudaMalloc(&list_gpu, sizeof(int) * sizeArrayLonger);
    cudaMalloc(&list_smaller_gpu, sizeof(int) * sizeArraySmaller);
    cudaMalloc(&result_gpu, sizeof(int) * length_result);

    cudaMemcpy(list_gpu, list, sizeof(int) * sizeArrayLonger, cudaMemcpyHostToDevice);
    cudaMemcpy(list_smaller_gpu, list_smaller, sizeof(int) * sizeArraySmaller, cudaMemcpyHostToDevice);
    cudaMemset(result_gpu, 0, sizeof(int) * length_result);

    typedef thrust::device_ptr<int> device_ptr;

    thrust::set_intersection(device_ptr(list_gpu), device_ptr(list_gpu + sizeArrayLonger), device_ptr(list_smaller_gpu),
        device_ptr(list_smaller_gpu + sizeArraySmaller), device_ptr(result_gpu), thrust::less<int>() );

    // MOVING TO CPU THE MARKER ARRAY OF ELEMENTS OF INTERSECTION SET
    cudaMemcpy(result, result_gpu, sizeof(int)*length_result, cudaMemcpyDeviceToHost);

    cudaDeviceSynchronize();

    // THIS LOOP ITERATES ALL ARRAY NAMED "result" WHERE THE POSITION ARE MARKED WITH 1
    int counter = 0;
    for (int i=0; i < length_result; i++)
        if (result[i]) {
            printf("\n-> %d", result[i]);
            counter++;
        }

    printf("\nTHRUST -> Total of elements: %d\n", counter);

    cudaDeviceReset();

    return 0;
}

Answer 1

似乎 OP 最近没有访问过，所以我将扩展我的 cmets 以供其他读者使用。 （我希望得到一些确认，即在编译期间指定正在使用的设备的计算目标也可以修复 OP 的观察结果。）

根据我的测试，OP的代码会：

• 如果为 cc2.0 设备编译并在 cc2.0 设备上运行，则通过。
• 如果为 cc3.0 设备编译并在 cc3.x 设备上运行，则通过。
• 如果为 cc2.0 设备编译并在 cc3.x 设备上运行，则会失败。

最后一个结果有点不直观。通常，由于the runtime JIT mechanism，我们喜欢将使用 PTX（例如nvcc -arch=sm_20 ... 或类似）编译的 CUDA 代码视为与未来架构的前向兼容。

但是，存在一个陷阱（以及与推力相关的问题。）CUDA 代码查询它们实际运行的设备（例如通过cudaGetDeviceProperties）并做出决策（例如内核）并不少见配置决策）基于使用的设备。具体来说，在这种情况下，thrust 是在底层启动一个内核，并根据实际使用的设备决定为该内核选择的网格 x 维度的大小。 CC 2.x 设备的此参数限制为 65535，但 CC 3.x 和更高版本的设备have a much higher limit。所以，在这种情况下，对于足够大的数据集，如果thrust检测到它在cc3.0设备上运行，它会为这个特定的内核配置一个大于65535的网格x维度。（对于足够小的数据集，它不会这样做，因此这个可能的错误不会浮出水面。因此，问题与问题大小有松散的联系。）

如果我们在二进制文件中嵌入了 cc 2.x 和 cc 3.x PTX（或适当的 SASS），那么仍然不会有问题。但是，如果我们只在二进制文件中嵌入了 cc2.x PTX，那么 JIT 进程将使用它来创建适合在 cc 3.x 设备上运行的机器代码，如果那是正在使用的设备。 但是这个前向 JIT 编译的 SASS 仍然受到 CC 2.x 的限制，包括 65535 的网格 X 尺寸限制。但是 cudaGetDeviceProperties 返回设备是 cc3.x 设备的事实，因此如果将此信息用于此特定决策（可接受的网格 X 尺寸），则会产生误导。

由于此序列，内核配置不正确，内核启动失败并出现特定类型的非粘性 CUDA 运行时 API 错误。这种类型的非粘性错误不会破坏 CUDA 上下文，因此仍然允许进一步的 CUDA 操作，并且未来的 CUDA API 调用将不会返回此错误。为了在 CUDA 内核启动后捕获此类错误，有必要在内核启动后发出 cudaGetLastError() 或 cudaPeekAtLastError() 调用，如 proper cuda error checking 所建议的那样。未能执行此操作意味着错误“丢失”并且无法从未来的 CUDA API 调用中发现（cudaGetLastError() 或 cudaPeekAtLastError() 除外），因为它们不会在状态返回值中指示存在此错误或内核启动失败。

以上大部分内容都可以通过仔细使用 cuda 分析工具来发现，例如nvprof，在通过和失败的情况下，以及cuda-memcheck。在过去的案例中，cuda-memcheck 没有报告任何错误，分析器显示对cudaLaunch 的 8 次调用以及在 GPU 上实际执行的 8 个内核。在失败的情况下，cuda-memcheck 报告了 2 个上述类型的内核启动失败，并且分析器显示了 8 次 cudaLaunch 调用，但实际上只有 6 个内核在 GPU 上执行。在 cc2.x GPU 上运行时，失败的内核配置为 65535 的网格 X 维度，而在 cc3.x GPU 上运行时配置为更大的数字。

因此，通过适当的 cuda 错误检查，上述序列虽然不一定可取，但至少会因显式错误而失败。但是 OP 的代码会静默失败 - 在失败的情况下它会返回错误的结果，但推力不会引发任何类型的错误。

事实证明，在引擎盖下，对从集合操作启动的内核（至少是这个，特别是这个）进行的推力错误检查具有这种特殊的错误检查差距。

通过仔细研究分析器的输出，我们可以发现哪些文件包含在这种情况下推力用于启动关闭的代码（即内核启动的实际来源）。 （您也可以通过仔细跟踪模板序列来解决这个问题。）在特定的失败情况下，我相信内核启动是由here 引起的。如果我们查看其中一个内核启动，我们会看到类似this：

#ifndef __CUDA_ARCH__ 
  kernel<<<(unsigned int) num_blocks, (unsigned int) block_size, (unsigned int) smem_size, stream(thrust::detail::derived_cast(exec))>>>(f); 
#else 
  ...
#endif // __CUDA_ARCH__ 
  synchronize_if_enabled("launch_closure_by_value"); 

synchronize_if_enabled（在这个特定的代码路径中）将在内核启动后立即被调用。那个函数可以找到here:

inline __host__ __device__ 
void synchronize_if_enabled(const char *message) 
{ 
// XXX this could potentially be a runtime decision 
//     note we always have to synchronize in __device__ code 
#if __THRUST_SYNCHRONOUS || defined(__CUDA_ARCH__) 
  synchronize(message); 
#else 
  // WAR "unused parameter" warning 
  (void) message; 
#endif

调用synchronize()：

inline __host__ __device__ 
void synchronize(const char *message) 
{ 
  throw_on_error(cudaDeviceSynchronize(), message); 
} // end synchronize() 

我们在synchronize() 中看到throw_on_error 调用cudaDeviceSynchronize()，这消除了先前表示错误配置的内核启动尝试的非粘性错误11，实际上返回cudaSuccess（因为cudaDeviceSynchronize() 操作本身是，其实是成功的。）

所以总结是存在两个问题：

1. Thrust（在这种情况下）做出关于内核启动配置的运行时决定，如果执行设备是 cc3.0 或更高版本并且代码是为 cc2.x 编译的（仅），这将是不正确的。

2. 对该特定 set_intersection 调用的推力错误检查存在缺陷，因为它没有适当的机制来捕获与错误配置的内核启动相关的非粘性 CUDA 运行时 API 错误（错误 11）。

因此，如果您打算在 cc3.0 或更高版本的设备上运行，建议始终编译指定 cc3.0 或更高版本目标的推力代码（至少）。 （当然，您可以同时指定一个 cc2.x 和一个 cc3.x 目标，并选择适当的 nvcc 命令行开关。）Thrust 在后台使用各种启动机制，并且并非所有（也许大多数）都不受这种特殊缺陷的影响（#2），但在我看来（对我来说）这个特殊的 set_intersection 调用 在这个时候受到这种缺陷的影响（推力v1.8)。

（对我来说）不清楚是否有办法系统地解决上述第一个问题（#1）。我已将上述第二个问题（#2）提请推力开发人员注意（通过 RFE 或错误报告。）

作为一种解决方法，推力开发人员可以在他们的推力应用程序中插入对cudaGetLastError() 的调用（可能在最后），以防止这种类型的错误“静默”。