CUDA GPU与CPU浮点运算结果不一致的排查及统一方法问询

编写的易并行内核 #

function kernel!(u, r, β::Number)
    u .= u .* β .+ r
end

CPU与GPU的一致性测试代码 #

using Test
using CUDA
using LinearAlgebra: norm

@testset "What is my GPU doing?" begin
    N = 1024
    u = rand(N)
    u_cu = CuVector(u)
    @test Vector(u_cu) == u  # 初始数据一致
    r = rand(N)
    r_cu = CuVector(r)
    @test Vector(r_cu) == r  # 初始数据一致
    β = 3.0
    u .= u .* β .+ r
    u_cu .= u_cu .* β .+ r_cu
    @test_broken u == u_cu   # 结果不一致
    @test norm(u - Vector(u_cu)) < 1.0e-14  # 误差处于机器epsilon级别
    @test eltype(u_cu) === eltype(u)  # 元素类型一致
end

测试结果显示CPU与GPU的运算结果存在微小误差，推测是浮点运算不满足结合律导致某一方重排了运算逻辑。以下是针对问题的解答：

问题解答 #

1. 确认差异原因：运算重排还是其他因素？ #

• 验证FMA的影响：从你补充的FMA测试结果可以直接得出结论：GPU自动将u_cu .* β .+ r_cu融合成了FMA（融合乘加）指令，而CPU上分开执行乘、加操作的结果和显式调用fma的结果不同，说明CPU未自动做FMA融合。FMA会一次性完成a*b+c，仅进行一次舍入；而分开计算(a*b)+c会经历两次舍入，这正是浮点结合律不满足导致的差异核心原因。
• 查看编译后的指令：
  • 对于GPU：把广播操作包装成显式内核，再用@device_code_ptx查看指令：

    function broadcast_kernel!(u, r, β)
        idx = threadIdx().x + (blockIdx().x - 1) * blockDim().x
        if idx <= length(u)
            u[idx] = u[idx] * β + r[idx]
        end
        return
    end
    
    N = 1024
    u_cu = CuVector(rand(N))
    r_cu = CuVector(rand(N))
    β = 3.0
    @device_code_ptx @cuda threads=256 blocks=cld(N,256) broadcast_kernel!(u_cu, r_cu, β)
    

    可以在输出中查找fma.rn之类的指令，确认GPU是否使用了FMA。
  • 对于CPU：用@code_native查看生成的汇编代码，判断是否是独立的mul和add指令，而非FMA。
• 排除其他因素：你的测试已经验证了初始数据拷贝一致、元素类型相同，因此可以排除数据传输、类型转换等误差来源。

2. 让CPU与GPU运算结果完全一致 #

• 强制CPU使用FMA：在支持FMA的CPU上，通过编译器选项开启FMA融合，让CPU的乘加逻辑和GPU对齐。在Julia中可以：
  • 启动时添加参数：julia --cpu-fma=on
  • 运行时设置环境变量：ENV["JULIA_CPU_FMA"] = true
    开启后，CPU上的u .* β .+ r会被编译成FMA指令，与GPU结果一致。
• 强制GPU禁用FMA：如果不想修改CPU行为，可在GPU端禁用自动FMA融合：

  @fastmath false begin
      u_cu .= u_cu .* β .+ r_cu
  end
  

  注意需要通过查看PTX代码确认GPU编译器是否真的禁用了FMA，部分GPU架构可能仍会保留部分优化。
• 使用高精度浮点类型：如果对精度要求极高，可以使用DoubleFloats等扩展精度类型，大幅降低舍入误差的影响，但会牺牲运算性能。

辅助问题：查看广播操作的PTX/SASS #

自动广播操作是由CUDA后端隐式生成内核的，无法直接用@device_code_ptx捕获。解决方法是手动实现等价的显式内核（如上面的broadcast_kernel!），然后对该显式内核使用@device_code_ptx或@device_code_sass工具，即可查看对应的PTX或SASS代码，其逻辑与自动广播完全一致。

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