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

Question

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

• 每次在内核中修改矩阵时写入一个全局标志。这可行，但效率非常低，并且在技术上不是线程安全的。
• 使用原子操作来完成与上面相同的操作。同样，这似乎效率低下，因为在最坏的情况下，每个线程都会发生一次全局写入。
• 使用缩减内核计算矩阵的某些参数(例如总和、均值、方差)。在某些情况下这可能会更快，但仍然看起来有点矫枉过正。此外，可以设想矩阵发生变化但总和/均值/方差没有变化的情况(例如，交换两个元素)。

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

Answer 1

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

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

这需要内核中每个 warp 共享一个字的内存，这是一个很小的开销，如果您提供每个 block 的 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 值。主机线程等待设备完成，然后可以直接从映射的主机变量中读取结果。

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