为何CUDA复制远线程索引数组值时结果不一致？

问题现象 #

当offset超过约20500时，以下CUDA代码的GPU结果与CPU模拟结果不一致（阈值因平台/GPU架构而异）：

测试代码 #

#ifndef __CUDACC__ 
#define __CUDACC__
#endif
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "device_functions.h"

#include <memory>
#include <stdio.h>
#include <cmath>

constexpr unsigned int N = 100000; // array size
constexpr unsigned int BLOCK_DIM = 256;
constexpr unsigned int GRID_DIM = (N + BLOCK_DIM - 1) / BLOCK_DIM;
constexpr unsigned int offset = 30000; // please set the offset here

__global__ void foo(float* a, float* b) {
    unsigned int tid = threadIdx.x + blockDim.x * blockIdx.x;
    if (tid < N) {
        b[tid] *= a[tid];
    }
    __syncthreads();
    if (tid < 10) {
        b[tid] = b[tid + offset];
    }
}

int main(void) {
    float* a, * b, *resCpu, * dev_a, * dev_b, *resGpu;

    // Allocate on host
    a = (float*)malloc(sizeof(float) * N);
    b = (float*)malloc(sizeof(float) * N);
    resCpu = (float*)malloc(sizeof(float) * N);
    resGpu = (float*)malloc(sizeof(float) * N);

    // Allocate on device
    cudaMalloc(&dev_a, sizeof(float) * N);
    cudaMalloc(&dev_b, sizeof(float) * N);

    // Initialize arrays
    float phase = 0.2f;
    for (unsigned int i = 0; i < N; i++) {
        a[i] = std::sin((float)i);
        b[i] = std::sin((float)i + phase);
    }

    // Replicate GPU algorithm on CPU 
    for (unsigned int i = 0; i < N; i++) {
        resCpu[i] = a[i] * b[i];
    }
    for (unsigned int j = 0; j < N / 2; j++) {
        resCpu[j] = resCpu[j + offset];
    }

    // Transfer data to GPU
    cudaMemcpy(dev_a, a, sizeof(float) * N, cudaMemcpyHostToDevice);
    cudaMemcpy(dev_b, b, sizeof(float) * N, cudaMemcpyHostToDevice);

    // Kernel call
    foo <<< GRID_DIM, BLOCK_DIM >>> (dev_a, dev_b);

    // Copy result to host
    cudaMemcpy(resGpu, dev_b, sizeof(float) * N, cudaMemcpyDeviceToHost);

    // Print and compare CPU and GPU results
    for (int i = 0; i < 10; i++) {
        printf("cpu[%d] = %.4f\n", i, resCpu[i]);
    }
    printf("\n");
    for (int i = 0; i < 10; i++) {
        printf("gpu[%d] = %.4f\n", i, resGpu[i]);
    }

    // Free memory 
    free(a);
    free(b);
    free(resCpu);
    cudaFree(dev_a);
    cudaFree(dev_b);
    cudaFree(resGpu);
}

输出对比 #

• offset=30000时（未定义行为）：

cpu[0] = 0.7263
cpu[1] = 0.7926
cpu[2] = 0.0022
cpu[3] = 0.5937
cpu[4] = 0.8918
cpu[5] = 0.0522
cpu[6] = 0.4529
cpu[7] = 0.9590
cpu[8] = 0.1371
cpu[9] = 0.3151

gpu[0] = -0.9048
gpu[1] = -0.8472
gpu[2] = -0.0106
gpu[3] = 0.8357
gpu[4] = 0.9137
gpu[5] = 0.1516
gpu[6] = -0.7498
gpu[7] = -0.9619
gpu[8] = -0.2896
gpu[9] = 0.6489

• offset=100时（结果一致）：

cpu[0] = 0.1645
cpu[1] = 0.2804
cpu[2] = 0.9900
cpu[3] = 0.2836
cpu[4] = 0.1619
cpu[5] = 0.9696
cpu[6] = 0.4190
cpu[7] = 0.0695
cpu[8] = 0.9110
cpu[9] = 0.5601

gpu[0] = 0.1645
gpu[1] = 0.2804
gpu[2] = 0.9900
gpu[3] = 0.2836
gpu[4] = 0.1619
gpu[5] = 0.9696
gpu[6] = 0.4190
gpu[7] = 0.0695
gpu[8] = 0.9110
gpu[9] = 0.5601

问题核心原因 #

1. __syncthreads()的局限性 #

__syncthreads()仅能同步当前线程块内的所有线程，无法跨线程块保证内存操作的完成顺序。当offset较小时，tid+offset对应的线程和tid=0-9的线程属于同一个线程块，__syncthreads()能确保前者完成b[tid] *= a[tid]的写操作后，后者才读取该值。但当offset足够大时，tid+offset对应的线程属于完全不同的线程块，__syncthreads()对这些线程没有约束。

2. CUDA弱一致性内存模型 #

CUDA的全局内存采用弱一致性模型：线程块对全局内存的写操作会先存入自身SM（流式多处理器）的L1缓存，不会立即同步到全局内存或其他SM的缓存。当tid=0-9的线程读取b[tid+offset]时，对应的线程块可能还在执行，或者其写操作仍停留在自身SM的缓存中，导致读取到的是内存中的旧值（初始化时的b值，而非乘以a后的新值）。

3. 硬件层面的缓存隔离 #

每个SM拥有独立的L1/共享缓存资源，线程块之间的缓存数据不会自动同步。当offset超过单个线程块覆盖的内存范围（比如256线程对应1024字节float数组，30000远大于这个值），tid+offset的线程在另一个SM上运行，其对b的修改无法被当前SM的线程感知，除非有显式的同步操作强制缓存写回全局内存。

解决方法 #

• 将操作拆分为两个独立内核：第一个内核完成所有b[tid] *= a[tid]的计算，调用cudaDeviceSynchronize()确保所有线程块执行完毕后，再启动第二个内核完成b[tid] = b[tid+offset]的赋值。
• 若要在单个内核中实现，需在__syncthreads()后添加__threadfence()，强制线程块将缓存中的写操作同步到全局内存，确保其他线程块能读取到最新值。

内容的提问来源于stack exchange，提问作者chckx592