相同的程序，不同的硬件，数值结果略有不同

Question

我有一个程序，使用相同版本的 Julia 在两台不同的电脑上运行。我不想给出程序的具体片段，因为一方面有点复杂，另一方面，我觉得一定有一些东西取决于机器而不是我的程序本身。

以下是我在两台机器上得到的一些结果：

1  1.0000000000582077  15.124999999941792
2  1.9999999999417923  27.536972172616515
3  3.000000000523869  45.722282028989866

---

1  1.0  15.125
2  2.0  27.53697217302397
3  3.0  45.722282028466

问题是舍入误差，我不知道它们是在哪一点引入的，导致结果略有不同。它不可能是格式错误（我在两种情况下使用完全相同的程序），也不可能是机器 epsilon 错误（两种情况下的机器 epsilon 都是>2.220446049250313e-16）。

目前我只能认为我使用的额外软件包，SparseArrays或LinearAlgebra正在FFTW引入此人为错误，但我不确定。

---

编辑：

好吧，这里是一个最小的例子（也许数组不需要像 [a,a] 一样，其中 a 本身就是一个稀疏数组）。我必须修改我的程序才能隔离问题，因此我下面将讨论的结果与我上面写的结果略有不同。然而，两者都显示了我提到的不可忽视的错误。

using LinearAlgebra
using SparseArrays
using FFTW

const Nsites = 10000

function Opbuild1(::Val{true})
    Up = spzeros(Complex{Float64}, Nsites, Nsites)
    N = div(Nsites,2)
    ii = 1
    for p in 0:(N-1)
        Up[ii,ii] = exp(-1.0im*p^2*0.25)
        Up[ii + N, ii + N] = exp(-1.0im*(p - N)^2*0.25)
        ii += 1
    end
    return kron(sparse(Matrix(1.0I, 2, 2)), Up)
end

function Opbuild2(::Val{true})
    pop = spzeros(Float64, Nsites, Nsites)
    N = div(Nsites,2)
    ii = 1
    for p in 0:(N-1)
        pop[ii,ii] = p
        pop[ii+N,ii+N] = (p-N)
        ii += 1
    end
    return kron(sparse(Matrix(1.0I, 2, 2)), pop)
end

function initst1(BoolN, irealiz)
    psi = spzeros(Complex{Float64}, Nsites)
    a = range(680,stop=783,step=1)
    x = a[irealiz]
    psi[x] = 1.0
    psi = kron([1.0,1.0im]/sqrt(2), psi)
    return psi, x
end

function ifftcustom(array)
    array = Array(array)
    array1 = copy(array)
    array1[1:Nsites] = ifft(array[1:Nsites])
    array1[Nsites+1:end] = ifft(array[Nsites+1:end])
    return array1/norm(array1)
end

function fftcustom(array)
    array = Array(array)
    array1 = copy(array)
    array1[1:Nsites] = fft(array[1:Nsites])
    array1[Nsites+1:end] = fft(array[Nsites+1:end])
    return array1/norm(array1)
end

function main()
    boolN = rem(Nsites,2) == 0
    Operator = Opbuild1(Val(boolN))
    Operator_p = Opbuild2(Val(boolN))
    #                                                                                                                                                                                                              
    psi0_inp, pinit = initst1(boolN, 1) #!                                                                                                                                                                         
    psi0_inp = round.(ifftcustom(Operator*psi0_inp), digits=15) # (*)See below
    psi0_inp = round.(fftcustom(psi0_inp), digits=15)
    psi0_inp = round.(psi0_inp, digits=15)
    @show psi0_inp'*Operator_p*psi0_inp
end

main()

所以区别如下。如果(*)我改为运行行 psi0_inp = round.(ifftcustom(psi0_inp), digits=15), 在@show我679.000000000001在两台机器中获得的部分中。另一方面，如果我按照我在代码中编写的方式运行它，在一台机器上我得到，679.0000000000001 + 0.0im但在另一台机器上得到679.0 + 1.2036944776504516e-16im机器上得到。

我不关心 -16im，因为这是双精度中的“零”，但实际上实部的数量级是 -13，这在双精度中不完全是零。

Answer 1

如果机器具有不同的 CPU 架构，则由于不同的 SIMD 指令，特别是不同的并行指令集宽度，结果可能会有所不同。差异是由于以下几点造成的：

1. 与它们近似的纯数学运算不同，浮点运算（如+和 ）*是非关联的，即(x + y) + z可以不同于x + (y + z)。

2. 为了利用 SIMD 并行指令，编译器将重新关联操作作为相当重要的优化。Julia 代码需要显式选择这种重新关联，但像 BLAS（用于稠密线性代数）、SuiteSparse（用于稀疏线性代数）和 FFTW 这样的高性能库肯定会以完美再现性为代价实现此类优化。如果这些库不允许 SIMD 优化，它们的速度会慢 2-8 倍，因此它们都这样做。

3. 在实现求和等操作时，数值运算的关联树取决于 SIMD 指令的宽度。如果“车道宽度”为W，则总和被实现为W个并行总和，它们在最后相加在一起。由于不同CPU架构的宽度W可能不同，因此关联性树可能不同，因此数值结果可能不同。

请注意，这些都不是 Julia 特有的，如果程序是使用这些库（BLAS、SuiteSparse 和 FFTW）以任何语言编写的，也会发生这种情况。