为什么矩阵求逆计算中的中间值在 C# 和 CUDA C++ 之间略有不同？

Question

我编写了两个程序来使用 Gaussian elimination 计算矩阵的逆，第一个程序是用 C# 编写的，第二个程序是用 CUDA C++ 编写的。这两个程序遵循完全相同的程序并给出相同的最终结果。但是，当我检查中间步骤中的值时，我发现值略有不同，小于 1e-5 相对误差。

这是两个程序的每个代码的一部分。

C#

int i, j, i1, n, y, z;
double[,] M = new double[n, n];
double[,] inv = new double[n, n];
for (i = 0; i < n; i++)
    inv[i, i] = 1;
for (i = 0; i < n; i++)
{
    for (j = i + 1; j < n; j++)
        M[i, j] /= M[i, i];
    for (j = 0; j < n; j++)
        inv[i, j] /= M[i, i];
    if (i != n - 1)
    {
        for (i1 = i + 1; i1 < n; i1++)
            if (Math.Abs(M[i1, i]) >= 1e-9)
            {
                for (j = i + 1; j < n; j++)
                    M[i1, j] -= M[i1, i] * M[i, j];
                for (j = 0; j < n; j++)
                    inv[i1, j] -= M[i1, i] * inv[i, j];
            }
        f = new StreamWriter("CPU.csv");
        for (y = 0; y < n; y++)
        {
            for (z = 0; z < n; z++)
                f.Write(M[y, z].ToString() + ",");
            for (z = 0; z < n; z++)
                f.Write(ans[y, z].ToString() + ",");
            f.WriteLine();
        }
        f.Close();
    }
}
for (i = n - 1; i > 0; i--)
{
    for (i1 = 0; i1 < i; i1++)
        if (Math.Abs(M[i1, i]) >= 1e-9)
            for (j = 0; j < n; j++)
                inv[i1, j] -= M[i1, i] * inv[i, j];
}

CUDA C++

int i, j;
double v;
double* d_A, * d_B, * d_v, * Z;
size = n * n * sizeof(double);
cudaMalloc(&d_A, size);
cudaMemcpy(d_A, A, size, cudaMemcpyHostToDevice);
cudaMalloc(&d_B, size);
cudaMalloc(&d_v, sizeof(double));
Z = new double[n * n];
Unity <<<1, n>>> (d_B, n);
cudaDeviceSynchronize();
for (i = 0; i < n; i++)
{
    GetVal <<<1, 1>>> (d_A, i * (n + 1), d_v);
    cudaMemcpy(&v, d_v, sizeof(double), cudaMemcpyDeviceToHost);
    if (i != n - 1)
        DivideRow <<<1, n - i - 1>>> (d_A, i * (n + 1) + 1, n - i - 1, v);
    DivideRow <<<1, n>>> (d_B, i * n, n, v);
    cudaDeviceSynchronize();
    cudaMemcpy(Z, d_A, size, cudaMemcpyDeviceToHost);
    cudaMemcpy(B, d_B, size, cudaMemcpyDeviceToHost);
    if (i != n - 1)
    {
        dim3 GridA(1, 1);
        dim3 BlockA(n - i - 1, n - i - 1);
        dim3 GridB(1, 1);
        dim3 BlockB(n - i - 1, n);
        ModifyRow <<<GridA, BlockA>>> (d_A, i, i, i + 1, n - i - 1, n - i - 1);
        ModifyRow <<<GridB, BlockB>>> (d_A, n, i, i, d_B, i + 1, 0, n - i - 1, n);
        cudaDeviceSynchronize();
        cudaMemcpy(Z, d_A, size, cudaMemcpyDeviceToHost);
        cudaMemcpy(B, d_B, size, cudaMemcpyDeviceToHost);
        myfile.open("GPU.csv");
        for (x = 0; x < n; x++)
        {
            for (y = 0; y < n; y++)
                myfile << Z[x * n + y] << ",";
            for (y = 0; y < n; y++)
                myfile << B[x * n + y] << ",";
            myfile << "\n";
        }
        myfile.close();
    }
}
cudaFree(d_v);
for (i = n - 1; i > 0; i--)
{
    dim3 GridB(1, 1);
    dim3 BlockB(i, n);
    ModifyRow <<<GridB, BlockB>>> (d_A, n, i, i, d_B, 0, 0, i, n);
    cudaDeviceSynchronize();
    cudaMemcpy(Z, d_A, size, cudaMemcpyDeviceToHost);
    cudaMemcpy(B, d_B, size, cudaMemcpyDeviceToHost);
}
cudaMemcpy(B, d_B, size, cudaMemcpyDeviceToHost);
cudaFree(d_A);
cudaFree(d_B);

我比较了 CPU.csv 和 GPU.csv 文件中的值，发现了这些差异。

这可能是什么原因？ CUDA C++ 中的计算精度是否低于 C#？

Answer 1

从 NVIDIA documentation 开始（大约下降 2/3）：

The consequence [of rounding] is that different math libraries cannot be expected to compute exactly the same result for a given input. This applies to GPU programming as well. Functions compiled for the GPU will use the NVIDIA CUDA math library implementation while functions compiled for the CPU will use the host compiler math library implementation (e.g., glibc on Linux). Because these implementations are independent and neither is guaranteed to be correctly rounded, the results will often differ slightly.

告诉你所有你需要知道的，真的。

为什么矩阵求逆计算中的中间值在 C# 和 CUDA C++ 之间略有不同？

Why the intermediate values in matrix inverse calculations are slightly different between C# and CUDA C++?

c#

c++

cuda

precision

matrix-inverse