CUDA内核中的Heisenbug，全局内存访问答案

【问题标题】：Heisenbug in CUDA kernel, global memory accessCUDA内核中的Heisenbug，全局内存访问
【发布时间】：2015-01-28 22:36:58
【问题描述】：

大约两年前，我编写了一个内核，用于同时处理多个数值网格。出现了一些非常奇怪的行为，导致了错误的结果。当在内核中利用 printf()-statements 寻找错误时，错误消失了。

由于截止日期的限制，我一直保持这种方式，尽管最近我发现这不是合适的编码风格。所以我重新审视了我的内核并将其归结为您在下面看到的内容。

__launch_bounds__(672, 2)
__global__ void heisenkernel(float *d_u, float *d_r, float *d_du, int radius,
        int numNodesPerGrid, int numBlocksPerSM, int numGridsPerSM, int numGrids)
{
    __syncthreads();
    int id_sm           = blockIdx.x /   numBlocksPerSM;                                    // (arbitrary) ID of Streaming Multiprocessor (SM) this thread works upon           - (constant over lifetime of thread)
    int id_blockOnSM    = blockIdx.x % numBlocksPerSM;                                      // Block number on this specific SM                                                 - (constant over lifetime of thread)
    int id_r            = id_blockOnSM  * (blockDim.x - 2*radius) + threadIdx.x - radius;   // Grid point number this thread is to work upon                                    - (constant over lifetime of thread)
    int id_grid         = id_sm         * numGridsPerSM;                                    // Grid ID this thread is to work upon                                              - (not constant over lifetime of thread)

    while(id_grid < numGridsPerSM * (id_sm + 1))    // this loops over numGridsPerSM grids
    {
        __syncthreads();
        int id_numInArray       = id_grid * numNodesPerGrid + id_r;     // Entry in array this thread is responsible for (read and possibly write)  - (not constant over lifetime of thread)
        float uchange           = 0.0f;
        //uchange                   = 1.0f;                                 // if this line is uncommented, results will be computed correctly ("Solution 1")
        float du                = 0.0f;

        if((threadIdx.x > radius-1) && (threadIdx.x < blockDim.x - radius) && (id_r < numNodesPerGrid) && (id_grid < numGrids))
        {
            if (id_r == 0)  // FO-forward difference
                du = (d_u[id_numInArray+1] - d_u[id_numInArray])/(d_r[id_numInArray+1] - d_r[id_numInArray]);
            else if (id_r == numNodesPerGrid - 1)  // FO-rearward difference
                du = (d_u[id_numInArray] - d_u[id_numInArray-1])/(d_r[id_numInArray] - d_r[id_numInArray-1]);
            else if (id_r == 1 || id_r == numNodesPerGrid - 2) //SO-central difference
                du = (d_u[id_numInArray+1] - d_u[id_numInArray-1])/(d_r[id_numInArray+1] - d_r[id_numInArray-1]);
            else if(id_r > 1 && id_r < numNodesPerGrid - 2)
                du = d_fourpoint_constant * ((d_u[id_numInArray+1] - d_u[id_numInArray-1])/(d_r[id_numInArray+1] - d_r[id_numInArray-1])) + (1-d_fourpoint_constant) * ((d_u[id_numInArray+2] - d_u[id_numInArray-2])/(d_r[id_numInArray+2] - d_r[id_numInArray-2]));
            else
                du = 0;
        }

        __syncthreads();
        if((threadIdx.x > radius-1 && threadIdx.x < blockDim.x - radius) && (id_r < numNodesPerGrid) && (id_grid < numGrids))
        {
            d_u[    id_numInArray] = d_u[id_numInArray] * uchange;          // if this line is commented out, results will be computed correctly ("Solution 2")
            d_du[   id_numInArray] = du;
        }

    __syncthreads();
    ++id_grid;
}

这个内核计算一些数值在所有网格点的导数，用于许多数字一维网格。

需要考虑的事项：（请参阅底部的完整代码库）

一个网格由 1300 个网格点组成
每个网格必须由两个块处理（由于内存/寄存器限制）
每个块在 37 个网格上连续工作（或更好：网格的一半，while 循环负责处理）
每个线程负责每个网格中的相同网格点
为了计算导数，线程需要访问来自四个下一个网格点的数据
为了保持块相互独立，在网格上引入了一个小的重叠（每个网格的网格点 666、667、668、669 由来自不同块的两个线程读取，尽管只有一个线程是写信给他们，正是这种重叠出现了问题）
由于沸腾过程，块两边的两个线程不做任何计算，原来它们负责将相应的网格值写入共享内存

网格的值存储在u_arr、du_arr 和r_arr（及其对应的设备数组d_u、d_du 和d_r）中。每个网格在每个数组中占用 1300 个连续值。内核中的 while 循环为每个块迭代超过 37 个网格。

为了评估内核的工作情况，每个网格都使用完全相同的值进行初始化，因此确定性程序将为每个网格生成相同的结果。我的代码不会发生这种情况。

海森堡的怪异：

我将网格 0 的计算值与其他每个网格进行了比较，在重叠处（网格点 666-669）存在差异，但不一致。有些网格具有正确的值，有些则没有。连续两次运行会将不同的网格标记为错误。首先想到的是，在这个重叠处的两个线程尝试同时写入内存，尽管情况似乎并非如此（我检查了......并重新检查了）。

注释或取消注释行或使用printf() 进行调试将改变程序的结果也是：当“询问”负责相关网格点的线程时，他们告诉我一切都很好，而且他们实际上是正确的。一旦我强制线程打印出它的变量，它们就会被正确计算（更重要的是：存储）。使用 Nsight Eclipse 进行调试也是如此。

Memcheck / Racecheck：

cuda-memcheck（memcheck 和 racecheck）报告没有内存/竞争条件问题，尽管即使使用其中一种工具也有能力影响结果的正确性。 Valgrind 给出了一些警告，尽管我认为它们与我无法影响的 CUDA API 有关，并且似乎与我的问题无关。

（更新） 正如所指出的，cuda-memcheck --tool racecheck 仅适用于共享内存竞争条件，而手头的问题在 d_u 上存在竞争条件，即全局内存。

测试环境：

原始内核已在不同的 CUDA 设备和不同的计算能力（2.0、3.0 和 3.5）上进行了测试，每个配置中都会出现错误（以某种形式）。

我的（主要）测试系统如下：

2 x GTX 460，在运行 X-server 的 GPU 以及另一个
驱动版本：340.46
Cuda 工具包 6.5
Linux 内核 3.11.0-12-generic (Linux Mint 16 - Xfce)

解决方案状态：

到目前为止，我很确定某些内存访问是罪魁祸首，可能是编译器的一些优化或使用未初始化的值，而且我显然不了解一些基本的 CUDA 范例。内核中的printf() 语句（通过一些黑魔法也必须利用设备和主机内存）和 memcheck 算法（cuda-memcheck 和 valgrind）影响的事实行为指向同一方向。

我很抱歉这个有点复杂的内核，但我尽可能地简化了原始内核和调用，这就是我所能得到的。到目前为止，我已经学会欣赏这个问题，我期待着了解这里发生了什么。

在代码中标记了两个强制内核按预期工作的“解决方案”。

（更新） 正如下面正确答案中提到的，我的代码的问题是线程块边界的竞争条件。由于每个网格上有两个块在工作，并且不能保证哪个块首先工作，因此导致下面概述的行为。它还解释了使用代码中提到的“解决方案1”时的正确结果，因为uchange = 1.0 时输入/输出值d_u 不会改变。

简单的解决方案是将这个内核分成两个内核，一个计算d_u，另一个计算导数d_du。最好只有一个内核调用而不是两个，尽管我不知道如何使用-arch=sm_20 来实现这一点。使用-arch=sm_35 可能会使用动态并行来实现这一点，尽管第二次内核调用的开销可以忽略不计。

heisenbug.cu：

#include <cuda.h>
#include <cuda_runtime.h>
#include <stdio.h>

const float r_sol = 6.955E8f;
__constant__ float d_fourpoint_constant = 0.2f;

__launch_bounds__(672, 2)
__global__ void heisenkernel(float *d_u, float *d_r, float *d_du, int radius,
        int numNodesPerGrid, int numBlocksPerSM, int numGridsPerSM, int numGrids)
{
    __syncthreads();
    int id_sm           = blockIdx.x / numBlocksPerSM;                                      // (arbitrary) ID of Streaming Multiprocessor (SM) this thread works upon           - (constant over lifetime of thread)
    int id_blockOnSM    = blockIdx.x % numBlocksPerSM;                                      // Block number on this specific SM                                                 - (constant over lifetime of thread)
    int id_r            = id_blockOnSM  * (blockDim.x - 2*radius) + threadIdx.x - radius;   // Grid point number this thread is to work upon                                    - (constant over lifetime of thread)
    int id_grid         = id_sm         * numGridsPerSM;                                    // Grid ID this thread is to work upon                                              - (not constant over lifetime of thread)

    while(id_grid < numGridsPerSM * (id_sm + 1))    // this loops over numGridsPerSM grids
    {
        __syncthreads();
        int id_numInArray       = id_grid * numNodesPerGrid + id_r;     // Entry in array this thread is responsible for (read and possibly write)  - (not constant over lifetime of thread)
        float uchange           = 0.0f;
        //uchange                   = 1.0f;                                 // if this line is uncommented, results will be computed correctly ("Solution 1")
        float du                = 0.0f;

        if((threadIdx.x > radius-1) && (threadIdx.x < blockDim.x - radius) && (id_r < numNodesPerGrid) && (id_grid < numGrids))
        {
            if (id_r == 0)  // FO-forward difference
                du = (d_u[id_numInArray+1] - d_u[id_numInArray])/(d_r[id_numInArray+1] - d_r[id_numInArray]);
            else if (id_r == numNodesPerGrid - 1)  // FO-rearward difference
                du = (d_u[id_numInArray] - d_u[id_numInArray-1])/(d_r[id_numInArray] - d_r[id_numInArray-1]);
            else if (id_r == 1 || id_r == numNodesPerGrid - 2) //SO-central difference
                du = (d_u[id_numInArray+1] - d_u[id_numInArray-1])/(d_r[id_numInArray+1] - d_r[id_numInArray-1]);
            else if(id_r > 1 && id_r < numNodesPerGrid - 2)
                du = d_fourpoint_constant * ((d_u[id_numInArray+1] - d_u[id_numInArray-1])/(d_r[id_numInArray+1] - d_r[id_numInArray-1])) + (1-d_fourpoint_constant) * ((d_u[id_numInArray+2] - d_u[id_numInArray-2])/(d_r[id_numInArray+2] - d_r[id_numInArray-2]));
            else
                du = 0;
        }

        __syncthreads();
        if((threadIdx.x > radius-1 && threadIdx.x < blockDim.x - radius) && (id_r < numNodesPerGrid) && (id_grid < numGrids))
        {
            d_u[    id_numInArray] = d_u[id_numInArray] * uchange;          // if this line is commented out, results will be computed correctly ("Solution 2")
            d_du[   id_numInArray] = du;
        }

        __syncthreads();
        ++id_grid;
    }
}

bool gridValuesEqual(float *matarray, uint id0, uint id1, const char *label, int numNodesPerGrid){

    bool retval = true;
    for(uint i=0; i<numNodesPerGrid; ++i)
        if(matarray[id0 * numNodesPerGrid + i] != matarray[id1 * numNodesPerGrid + i])
        {
            printf("value %s at position %u of grid %u not equal that of grid %u: %E != %E, diff: %E\n",
                    label, i, id0, id1, matarray[id0 * numNodesPerGrid + i], matarray[id1 * numNodesPerGrid + i],
                    matarray[id0 * numNodesPerGrid + i] - matarray[id1 * numNodesPerGrid + i]);
            retval = false;
        }
    return retval;
}

int main(int argc, const char* argv[])
{
    float *d_u;
    float *d_du;
    float *d_r;

    float *u_arr;
    float *du_arr;
    float *r_arr;

    int numNodesPerGrid = 1300;
    int numBlocksPerSM  = 2;
    int numGridsPerSM   = 37;
    int numSM           = 7;
    int TPB             = 672;
    int radius          = 2;
    int numGrids        = 259;
    int memsize_grid    = sizeof(float) * numNodesPerGrid;

    int numBlocksPerGrid    = numNodesPerGrid / (TPB - 2 * radius) + (numNodesPerGrid%(TPB - 2 * radius) == 0 ? 0 : 1);

    printf("---------------------------------------------------------------------------\n");
    printf("--- Heisenbug Extermination Tracker ---------------------------------------\n");
    printf("---------------------------------------------------------------------------\n\n");

    cudaSetDevice(0);
    cudaDeviceReset();

    cudaMalloc((void **) &d_u,      memsize_grid * numGrids);
    cudaMalloc((void **) &d_du,     memsize_grid * numGrids);
    cudaMalloc((void **) &d_r,      memsize_grid * numGrids);

    u_arr   = new float[numGrids * numNodesPerGrid];
    du_arr  = new float[numGrids * numNodesPerGrid];
    r_arr   = new float[numGrids * numNodesPerGrid];

    for(uint k=0; k<numGrids; ++k)
        for(uint i=0; i<numNodesPerGrid; ++i)
        {
            uint index  = k * numNodesPerGrid + i;

            if (i < 585)
                r_arr[index] = i * (6000.0f);
            else
            {
                if (i == 585)
                    r_arr[index] = r_arr[index - 1] + 8.576E-6f * r_sol;
                else
                    r_arr[index] = r_arr[index - 1] + 1.02102f  * ( r_arr[index - 1] - r_arr[index - 2] );
            }

            u_arr[index]    = 1E-10f * (i+1);
            du_arr[index]   = 0.0f;
        }

    /*
    printf("\n\nbefore kernel start\n\n");
    for(uint k=0; k<numGrids; ++k)
        printf("matrix->du_arr[k*paramH.numNodes + 668]:\t%E\n", du_arr[k*numNodesPerGrid + 668]);//*/

    bool equal = true;
    for(int k=1; k<numGrids; ++k)
    {
        equal &= gridValuesEqual(u_arr, 0, k, "u", numNodesPerGrid);
        equal &= gridValuesEqual(du_arr, 0, k, "du", numNodesPerGrid);
        equal &= gridValuesEqual(r_arr, 0, k, "r", numNodesPerGrid);
    }

    if(!equal)
        printf("Input values are not identical for different grids!\n\n");
    else
        printf("All grids contain the same values at same grid points.!\n\n");

    cudaMemcpy(d_u, u_arr,      memsize_grid * numGrids, cudaMemcpyHostToDevice);
    cudaMemcpy(d_du, du_arr,    memsize_grid * numGrids, cudaMemcpyHostToDevice);
    cudaMemcpy(d_r, r_arr,      memsize_grid * numGrids, cudaMemcpyHostToDevice);

    printf("Configuration:\n\n");
    printf("numNodesPerGrid:\t%i\nnumBlocksPerSM:\t\t%i\nnumGridsPerSM:\t\t%i\n", numNodesPerGrid, numBlocksPerSM, numGridsPerSM);
    printf("numSM:\t\t\t\t%i\nTPB:\t\t\t\t%i\nradius:\t\t\t\t%i\nnumGrids:\t\t\t%i\nmemsize_grid:\t\t%i\n", numSM, TPB, radius, numGrids, memsize_grid);
    printf("numBlocksPerGrid:\t%i\n\n", numBlocksPerGrid);
    printf("Kernel launch parameters:\n\n");
    printf("moduleA2_3<<<%i, %i, %i>>>(...)\n\n", numBlocksPerSM * numSM, TPB, 0);
    printf("Launching Kernel...\n\n");

    heisenkernel<<<numBlocksPerSM * numSM, TPB, 0>>>(d_u, d_r, d_du, radius, numNodesPerGrid, numBlocksPerSM, numGridsPerSM, numGrids);
    cudaDeviceSynchronize();

    cudaMemcpy(u_arr, d_u,      memsize_grid * numGrids, cudaMemcpyDeviceToHost);
    cudaMemcpy(du_arr, d_du,    memsize_grid * numGrids, cudaMemcpyDeviceToHost);
    cudaMemcpy(r_arr, d_r,      memsize_grid * numGrids, cudaMemcpyDeviceToHost);

    /*
    printf("\n\nafter kernel finished\n\n");
    for(uint k=0; k<numGrids; ++k)
        printf("matrix->du_arr[k*paramH.numNodes + 668]:\t%E\n", du_arr[k*numNodesPerGrid + 668]);//*/

    equal = true;
    for(int k=1; k<numGrids; ++k)
    {
        equal &= gridValuesEqual(u_arr, 0, k, "u", numNodesPerGrid);
        equal &= gridValuesEqual(du_arr, 0, k, "du", numNodesPerGrid);
        equal &= gridValuesEqual(r_arr, 0, k, "r", numNodesPerGrid);
    }

    if(!equal)
        printf("Results are wrong!!\n");
    else
        printf("All went well!\n");

    cudaFree(d_u);
    cudaFree(d_du);
    cudaFree(d_r);

    delete [] u_arr;
    delete [] du_arr;
    delete [] r_arr;

    return 0;
}

生成文件：

CUDA            = 1
DEFINES         = 

ifeq ($(CUDA), 1)
    DEFINES     += -DCUDA
    CUDAPATH    = /usr/local/cuda-6.5
    CUDAINCPATH = -I$(CUDAPATH)/include
    CUDAARCH    = -arch=sm_20
endif

CXX             = g++
CXXFLAGS        = -pipe -g -std=c++0x -fPIE -O0 $(DEFINES)
VALGRIND        = valgrind
VALGRIND_FLAGS  = -v --leak-check=yes --log-file=out.memcheck
CUDAMEMCHECK    = cuda-memcheck
CUDAMC_FLAGS    = --tool memcheck
RACECHECK       = $(CUDAMEMCHECK)
RACECHECK_FLAGS = --tool racecheck  
INCPATH         = -I. $(CUDAINCPATH)
LINK            = g++
LFLAGS          = -O0
LIBS            = 

ifeq ($(CUDA), 1)
    NVCC        = $(CUDAPATH)/bin/nvcc
    LIBS        += -L$(CUDAPATH)/lib64/ 
    LIBS        += -lcuda -lcudart -lcudadevrt
    NVCCFLAGS   = -g -G -O0 --ptxas-options=-v
    NVCCFLAGS   += -lcuda -lcudart -lcudadevrt -lineinfo --machine 64 -x cu $(CUDAARCH) $(DEFINES)
endif 

all: 
    $(NVCC) $(NVCCFLAGS) $(INCPATH) -c -o $(DST_DIR)heisenbug.o $(SRC_DIR)heisenbug.cu
    $(LINK) $(LFLAGS) -o heisenbug heisenbug.o $(LIBS)

clean:
    rm heisenbug.o
    rm heisenbug

memrace: all
    ./heisenbug > out
    $(VALGRIND) $(VALGRIND_FLAGS) ./heisenbug > out.memcheck.log
    $(CUDAMEMCHECK) $(CUDAMC_FLAGS) ./heisenbug > out.cudamemcheck
    $(RACECHECK) $(RACECHECK_FLAGS) ./heisenbug > out.racecheck

【问题讨论】：

标签： cuda nvcc

【解决方案1】：

请注意，在您的整个文章中，我没有看到明确提出的问题，因此我正在回复：

我期待了解这里发生了什么。

您在d_u 上存在竞争条件。

由您自己声明：

•为了保持块相互独立，在网格上引入了一个小的重叠（每个网格的网格点 666、667、668、669 由来自不同块的两个线程读取，尽管只有一个线程正在写信给他们，正是这个重叠出现了问题）

另外，如果你注释掉写到d_u，根据你在代码中的说法，问题就消失了。

CUDA 线程块可以按任何顺序执行。您至少有 2 个不同的块正在从网格点 666、667、668、669 读取数据。根据实际发生的情况，结果会有所不同：

两个块都在任何写入发生之前读取值。
一个块读取值，然后发生写入，然后另一个块读取值。

如果一个块正在读取另一个块可以写入的值，则这些块不是相互独立的（与您的陈述相反）。这种情况下block执行的顺序会决定结果，CUDA没有指定block执行的顺序。

注意cuda-memcheck 与-tool racecheck 选项only captures race conditions related to __shared__ memory usage。您发布的内核不使用 __shared__ 内存，因此我不希望 cuda-memcheck 报告任何内容。

cuda-memcheck，为了收集它的数据，确实影响了块执行的顺序，所以影响行为也就不足为奇了。

内核中printf 表示代价高昂的函数调用，写入全局内存缓冲区。所以它也会影响执行行为/模式。而且，如果您要打印大量数据，超出输出的缓冲区行数，则在缓冲区溢出的情况下，效果（就执行时间而言）是极其昂贵的。

顺便说一句，据我所知，Linux Mint 是not a supported distro for CUDA。但是我认为这与您的问题无关；我可以在受支持的配置上重现该行为。

【讨论】：

到底发生了什么，感谢您在另一个系统上进行测试。