CUDA并行Jacobi算法不收敛问题求助（附可复现代码）

安全阿岳

2026-6-29

CUDA并行Jacobi算法不收敛问题排查思路

我实现的基于CUDA的并行Jacobi迭代算法无法收敛，但串行版本运行完全正常，检查代码后未发现明显错误。以下是可复现代码：

constexpr int N{ 1024 };
constexpr int block_size = 256;
constexpr float epsilon = 1e-4;


void fillMatrix(float* M, const int size) {
    for (int i = 0; i < size; ++i) {
        M[i * size + i] = size + 1;
        for (int j = 0; j<size; ++j) {
            if (j != i) {
                M[j + size * i] = 1;
            }
        }
    }
}

void fillVector(float* M, const int size) {
    for (int i = 0; i < size; ++i) {
        M[i] = i;
    }
}

void solveJacobi(const float* A, const float* b, float* xNew, float* xOld, const int size) {
    for (int i = 0; i < size; ++i) {
        xOld[i] = xNew[i];
        xNew[i] = 0.0;
        for (int j = 0; j < size; ++j) {
            if (j != i) {
                xNew[i] -= A[i * size + j] * xOld[j];
            }
        }
        xNew[i] += b[i];
        xNew[i] /= A[i * size + i];
    }
}

void computeNorm(const float* xNew, const float* xOld, const int size, float* norm) {
    *norm = 0;
    for (int i = 0; i < size; ++i) {
        *norm += (xOld[i] - xNew[i])*(xOld[i] - xNew[i]);
    }
    *norm = sqrt(*norm);
}

__global__ void solveJacobiCuda(const float* A, const float* b, float* xNew, float* xOld, const int size){
    int idx = threadIdx.x+blockDim.x*blockIdx.x; // create thread x index
    if ((idx < size)){
        xOld[idx] = xNew[idx];
        xNew[idx] = 0;
        for(int i=0; i < size; ++i){
            if(i != idx){
                xNew[idx] -= A[idx * size + i] * xOld[i];
            }
        }
        xNew[idx] += b[idx];
        xNew[idx] /= A[idx * size + idx];
    }
}

int main() {
    float* h_A = new float[N * N];
    float* h_b = new float[N];
    float* h_bStar = new float[N];
    float* h_xOld = new float[N];
    float* h_xNew = new float[N];
    float* h_norm = new float{std::numeric_limits<float>::max()};

    fillMatrix(h_A, N);
    fillVector(h_b, N);
    fillVector(h_xOld,N);
    fillVector(h_xNew,N);
    
    float* d_A;
    float* d_b;
    float* d_xOld;
    float* d_xNew;

    cudaMalloc(&d_A, N*N*sizeof(float));
    cudaMalloc(&d_b, N*sizeof(float));
    cudaMalloc(&d_xOld, N*sizeof(float));
    cudaMalloc(&d_xNew, N*sizeof(float));

    cudaMemcpy(d_A, h_A, N*N*sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, h_b, N*sizeof(float), cudaMemcpyHostToDevice);

    
    int iterationCounter{ 0 };
    while (*h_norm > epsilon && iterationCounter < 10000) {
        ++iterationCounter;

        //solveJacobi(h_A, h_b, h_xNew, h_xOld, N);
        cudaMemcpy(d_xOld, h_xOld, N*sizeof(float), cudaMemcpyHostToDevice);
        cudaMemcpy(d_xNew, h_xNew, N*sizeof(float), cudaMemcpyHostToDevice);
        // Launch kernel
        solveJacobiCuda<<<(N+block_size-1)/block_size, block_size>>>(d_A, d_b, d_xNew, d_xOld, N);
        cudaCheckErrors("kernel launch failure");
        cudaMemcpy(h_xNew, d_xNew, N*sizeof(float), cudaMemcpyDeviceToHost);
        cudaMemcpy(h_xOld, d_xOld, N*sizeof(float), cudaMemcpyDeviceToHost);

        computeNorm(h_xNew, h_xOld, N, h_norm);
    }

    for(int i=0; i<N; ++i){
        h_bStar[i] = 0;
        for(int j=0; j<N; ++j){
            h_bStar[i] += h_A[i*N + j]*h_xNew[j];
        }
    }
        
    cout << "Jacobi finished" << endl;
    cout << "N_iterations = " << iterationCounter << endl;
    cout << "Final norm = " << *h_norm << endl;
    cout << "b*[0] = " << h_bStar[0] << endl;
    cout << "b*[N] = " << h_bStar[N-1] << endl;

    return 0;
}

排查思路

核心逻辑错误：xOld与xNew的更新顺序
Jacobi迭代的核心要求是：当前轮所有x的计算必须完全依赖上一轮的全量结果。串行代码中，先把xNew的所有值复制到xOld，再逐个计算新的xNew，全程xOld不会被修改。但并行kernel里，每个线程直接执行xOld[idx] = xNew[idx]，这会导致不同线程的xOld被实时覆盖——比如线程A刚更新了xOld[0]，线程B计算时可能用到的是刚被修改的xOld[0]，而非上一轮的原始值，直接破坏迭代逻辑。
修正方案：迭代开始前，在device上完成xNew到xOld的整体复制（用cudaMemcpy(d_xOld, d_xNew, N*sizeof(float), cudaMemcpyDeviceToDevice)），kernel内只读取xOld的旧值，计算并写入xNew，全程不修改xOld。
冗余数据传输与同步问题
现有代码每次迭代都将host端的xOld、xNew复制到device，不仅性能低下，还可能因数据传输顺序出错导致迭代异常。正确流程应为：
1. 初始化时将xOld、xNew一次性传到device
2. 迭代步骤：先在device内部完成xOld = xNew的复制
3. 启动Jacobi kernel，仅读取xOld，写入xNew
4. 仅在需要计算收敛norm时，将xNew回传到host（或直接在device上计算norm，避免频繁数据传输）
数组索引正确性验证
检查矩阵索引：串行代码中A[i * size + j]对应行优先存储的A[i][j]，并行kernel中A[idx * size + i]对应A[idx][i]，这部分是正确的，因为每个线程负责计算第idx行的xNew值。
收敛判断的一致性
确保norm计算时，xOld是上一轮的结果，xNew是当前轮的计算结果，不要搞混两者的传输顺序。比如kernel执行后，应先回传xNew到host的h_xNew，而h_xOld保留上一轮的值，再调用computeNorm。
完善CUDA错误检查
现有代码仅检查了kernel启动错误，需在cudaMalloc、cudaMemcpy等所有CUDA操作后添加错误检查，避免内存分配/传输失败导致的隐性问题。示例代码：
```
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
    fprintf(stderr, "CUDA error: %s\n", cudaGetErrorString(err));
    exit(EXIT_FAILURE);
}
```

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

CUDA并行Jacobi算法不收敛问题求助（附可复现代码）

排查思路 #

排查思路