在 CUDA 中检查矩阵稳定性的有效方法

Question

许多算法迭代直到达到某个收敛标准（例如特定矩阵的稳定性）。在许多情况下，每次迭代都必须启动一个 CUDA 内核。我的问题是：那么如何有效且准确地确定矩阵在最后一次内核调用过程中是否发生了变化？以下是三种似乎同样不令人满意的可能性：

• 每次在内核中修改矩阵时写入一个全局标志。这可行，但效率极低，并且在技术上不是线程安全的。
• 使用原子操作执行与上述相同的操作。同样，这似乎效率低下，因为在最坏的情况下，每个线程都会发生一次全局写入。
• 使用归约内核来计算矩阵的某些参数（例如总和、均值、方差）。在某些情况下，这可能会更快，但似乎仍然是矫枉过正。此外，还可以设想矩阵已更改但总和/均值/方差未更改的情况（例如交换了两个元素）。

是否存在上述三个选项中的任何一个或替代方案，被认为是最佳实践和/或通常更有效？

Answer 1

如果浏览器崩溃，我也会回到我在 2012 年发布的答案。

基本思想是，您可以使用 warp 投票指令执行简单、廉价的缩减，然后使用每个块的零个或一个原子操作来更新主机在每次内核启动后可以读取的固定映射标志。使用映射标志消除了在每次内核启动后显式设备到主机传输的需要。

这需要内核中每个 warp 一个字的共享内存，这是一个很小的开销，如果您提供每个块的 warp 数量作为模板参数，一些模板技巧可以允许循环展开。

一个完整的工作示例（使用 C++ 主机代码，我目前无法访问工作的 PyCUDA 安装）如下所示：

#include <cstdlib>
#include <vector>
#include <algorithm>
#include <assert.h>

__device__ unsigned int process(int & val)
{
    return (++val < 10);
}

template<int nwarps>
__global__ void kernel(int *inout, unsigned int *kchanged)
{
    __shared__ int wchanged[nwarps];
    unsigned int laneid = threadIdx.x % warpSize;
    unsigned int warpid = threadIdx.x / warpSize;

    // Do calculations then check for change/convergence 
    // and set tchanged to be !=0 if required
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    unsigned int tchanged = process(inout[idx]);

    // Simple blockwise reduction using voting primitives
    // increments kchanged is any thread in the block 
    // returned tchanged != 0
    tchanged = __any(tchanged != 0);
    if (laneid == 0) {
        wchanged[warpid] = tchanged;
    }
    __syncthreads();

    if (threadIdx.x == 0) {
        int bchanged = 0;
#pragma unroll
        for(int i=0; i<nwarps; i++) {
            bchanged |= wchanged[i];
        }
        if (bchanged) {
            atomicAdd(kchanged, 1);
        }
    }
}

int main(void)
{
    const int N = 2048;
    const int min = 5, max = 15;
    std::vector<int> data(N);
    for(int i=0; i<N; i++) {
        data[i] = min + (std::rand() % (int)(max - min + 1));
    }

    int* _data;
    size_t datasz = sizeof(int) * (size_t)N;
    cudaMalloc<int>(&_data, datasz);
    cudaMemcpy(_data, &data[0], datasz, cudaMemcpyHostToDevice);

    unsigned int *kchanged, *_kchanged;
    cudaHostAlloc((void **)&kchanged, sizeof(unsigned int), cudaHostAllocMapped);
    cudaHostGetDevicePointer((void **)&_kchanged, kchanged, 0);

    const int nwarps = 4;
    dim3 blcksz(32*nwarps), grdsz(16);

    // Loop while the kernel signals it needs to run again
    do {
        *kchanged = 0;
        kernel<nwarps><<<grdsz, blcksz>>>(_data, _kchanged);
        cudaDeviceSynchronize(); 
    } while (*kchanged != 0); 

    cudaMemcpy(&data[0], _data, datasz, cudaMemcpyDeviceToHost);
    cudaDeviceReset();

    int minval = *std::min_element(data.begin(), data.end());
    assert(minval == 10);

    return 0;
}

这里，kchanged 是内核用来向主机发出信号它需要再次运行的标志。内核一直运行，直到输入中的每个条目都增加到阈值以上。在每个线程处理结束时，它参与一次warp 投票，之后每个warp 中的一个线程将投票结果加载到共享内存中。一个线程减少扭曲结果，然后自动更新kchanged 值。主机线程等待设备完成，然后可以直接从映射的主机变量中读取结果。

您应该能够根据您的应用程序要求调整它

Answer 2

我会回到我原来的建议。我已经用我自己的答案更新了the related question，我认为这是正确的。

在全局内存中创建一个标志：

__device__ int flag;

在每次迭代中，

1. 将标志初始化为零（在主机代码中）：

   int init_val = 0;
   cudaMemcpyToSymbol(flag, &init_val, sizeof(int));
   

2. 在您的内核设备代码中，如果对矩阵进行了更改，请将标志修改为 1：

   __global void iter_kernel(float *matrix){
   
   ...
     if (new_val[i] != matrix[i]){
       matrix[i] = new_val[i];
       flag = 1;}
   ...
   }
   

3. 调用内核后，在迭代结束时（在宿主代码中），测试修改：

   int modified = 0;
   cudaMemcpyFromSymbol(&modified, flag, sizeof(int));
   if (modified){
     ...
     }
   

即使在单独的块甚至单独的网格中的多个线程正在写入flag 值，只要它们所做的唯一事情是写入相同的值（即在这种情况下为 1），就没有危险。写入不会“丢失”，flag 变量中不会出现虚假值。

以这种方式测试 float 或 double 数量是否相等是有问题的，但这似乎不是您问题的重点。如果您有首选方法来声明“修改”，请改用该方法（例如在容差范围内测试相等性，也许）。

此方法的一些明显改进是为每个线程创建一个（本地）标志变量，并让每个线程在每个内核更新一次全局标志变量，而不是在每次修改时更新。这将导致每个内核每个线程最多一次全局写入。另一种方法是在共享内存中为每个块保留一个标志变量，并让所有线程简单地更新该变量。在块完成时，对全局内存进行一次写入（如果需要）以更新全局标志。在这种情况下，我们不需要求助于复杂的归约，因为整个内核只有一个布尔结果，并且我们可以容忍多个线程写入共享或全局变量，只要所有线程都写入相同价值。

我看不出使用原子的任何理由，或者它对任何事情都有什么好处。

至少与其中一种优化方法（例如，每个块共享标志）相比，缩减内核似乎有点过头了。它会有你提到的缺点，例如任何小于 CRC 或类似复杂计算的东西都可能将两个不同的矩阵结果别名为“相同”。