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

Question

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

• 每次在内核中修改矩阵时写入全局标志。这种方法虽然可行，但效率极低，而且不是线程安全的。
• 使用原子操作来执行与上述相同的操作。同样，这似乎效率低下，因为在最坏情况下，每个线程都要进行一次全局写入。
• 使用约简内核计算矩阵的某个参数（例如总和、平均值、方差）。在某些情况下，这可能更快，但仍然似乎有点过度。此外，可以想象出矩阵已更改但总和/平均值/方差未更改的情况（例如，两个元素被交换）。

以上三种方法中是否有最佳实践或通常更有效的替代方法？

- user2398029

---

@talonmies 最后一次迭代中矩阵未被修改。 - user2398029

talonmies的建议很好。如果不小心处理，追求这种方法可能并不值得，因此我关于在全局内存中使用单个布尔标志的评论可能是不明智的。我对你的代码结构尤其是算法迭代方面了解不够。 - Robert Crovella

@louism：我几天前写了一篇回答，包括一些演示代码，但是浏览器崩溃了，很抱歉。我会做一个sum-reduce，但是你可以使用warp投票原语（例如__any()warp vote）来高效地完成它。然后，您只需要对每个块内的每个warp的结果进行非常简单的约简，并对每个块进行单个原子添加以更新全局标志。如果标志位于零拷贝内存中，则无需显式复制即可检查主机上的结果。 - talonmies

@talonmies 非常感谢，这肯定为我指明了方向。您能否更详细地解释一下“如果标志位在零拷贝内存中，则无需显式复制即可在主机上检查结果”的含义？ - user2398029

零拷贝内存（又称“固定映射”）是主机内存通过PCI Express总线映射到GPU地址空间中的一种技术。您无需使用复制即可将其读回主机，GPU会为您完成此操作。 - talonmies

Answer 1

我将回到2012年本应发布的答案，但因为浏览器崩溃而未能发表。

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

此操作需要每个内核中每个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中的一个线程将投票结果加载到共享内存中。一个线程减少了warp结果，然后原子更新了kchanged值。主机线程等待设备完成，然后可以直接从映射的主机变量中读取结果。

您应该能够根据您的应用程序要求进行适应。

- talonmies

---

我不清楚在没有拓扑结构的情况下如何评估矩阵所经历的变化。要生成拓扑结构，需要一种范数，而矩阵范数通常需要缩减操作。因此，我认为回答这个问题的唯一方法是通过缩减计算当前和上一迭代步骤中矩阵之间差异的范数，然后设置一个标志，如果差异的范数超过某个预定的精度，则返回一个可靠的理论基础。 - Vitality

这里的隐含假设是每个线程至少在输出中发出一个条目。因此，每个线程都能够确定其输出是否在一个内核运行过程中发生了变化（这是评论中给定的要求）。确定如何检测更改的标准是应用程序特定的，大多数情况下并不重要。上面的所有代码所做的就是以最便宜的方式减少线程范围状态标志，以便主机知道是否需要进一步迭代。 - talonmies

Answer 2

我会回到我的原始建议。我已经更新了相关问题，并提供了我自己的答案，我相信是正确的。

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

__device__ int flag;

在每次迭代中，

1. initialize the flag to zero (in host code):

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

2. In your kernel device code, modify the flag to 1 if a change is made to the matrix:

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

3. after calling the kernel, at the end of the iteration (in host code), test for modification:

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

即使多个线程在不同的块或甚至不同的网格中写入flag值，只要它们所做的唯一操作是写入相同的值（例如，在这种情况下为1），则没有任何危险。写入不会“丢失”，并且在flag变量中不会出现虚假值。

以这种方式测试float或double类型的相等性是有问题的，但这似乎并不是您的问题所在。如果您有首选的声明“修改”的方法，请使用该方法（例如，在容忍度范围内测试相等性）。

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

我看不到使用原子操作的原因，也不知道它如何有益于任何事情。

与优化方法（例如，每个块的共享flag）之一相比，约简内核似乎过于复杂。而且，它会有您提到的缺点，例如，如果计算不到CRC或类似的复杂度，则任何小于这种复杂度的计算可能会将两个不同的矩阵结果视为“相同”。

- Robert Crovella