Small discrepancy

[ymtoo]

MWE:

using KernelAbstractions

@kernel function testfunc1_kernel!(As, bs, ys)
    idx = @index(Global)
    @views A = KernelAbstractions.SMatrix{3,4,Float32}(As[:,:,idx])
    @views b = KernelAbstractions.SVector{3,Float32}(bs[:,idx])
    ys[:,idx] .= A[1:3,1:3] * b + A[1:3,4] 
end

@kernel function testfunc2_kernel!(As, bs, ys)
    idx = @index(Global)
    @views R = KernelAbstractions.SMatrix{3,3,Float32}(As[:,1:3,idx])
    @views t = KernelAbstractions.SVector{3,Float32}(As[:,4,idx])
    @views b = KernelAbstractions.SVector{3,Float32}(bs[:,idx])
    ys[:,idx] .= R * b + t
end

function testfunc1(As, bs)
    N =size(As, 3)
    backend = get_backend(As)
    ys = KernelAbstractions.zeros(backend, Float32, 3, N)
    kernel = testfunc1_kernel!(backend)
    kernel(As, bs, ys, ndrange=N)
    return ys
end

function testfunc2(As, bs)
    N =size(As, 3)
    backend = get_backend(As)
    ys = KernelAbstractions.zeros(backend, Float32, 3, N)
    kernel = testfunc2_kernel!(backend)
    kernel(As, bs, ys, ndrange=N)
    return ys
end

N = 100
As = randn(Float32, 3, 4, N)
bs = randn(Float32, 3, N)
ys1 = testfunc1(As, bs)
ys2 = testfunc2(As, bs)

_ys = stack(1:N) do i
    As[1:3,1:3,i] * bs[:,i] + As[1:3,4,i]
end

ys1 == _ys # true
ys2 == _ys # false

Any idea why testfunc2_kernel! introduces small discrepancies?

Julia and package versions:

julia> versioninfo()
Julia Version 1.10.4
Commit 48d4fd48430 (2024-06-04 10:41 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 32 × 13th Gen Intel(R) Core(TM) i9-13980HX
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, goldmont)
Threads: 1 default, 0 interactive, 1 GC (on 32 virtual cores)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 

(julia-kernel-programming) pkg> st
Status `~/Projects/julia-kernel-programming/Project.toml`
  [63c18a36] KernelAbstractions v0.9.20

[bjarthur]

is isapprox(ys2, _ys) true ? floating point arithmetic is not necessarily associative.

[ymtoo]

isapprox(ys2, _ys) returns true.

With the same order of operations for matrix multiplication and vector summation, testfunc1 is able to produce the exact output, i.e., ys1 == _ys is true. The only difference in the behavior of testfunc1 and testfunc2 lies in the way they handle array slicing and the creation of static arrays within the kernels.

[ymtoo]

Closing this issue now as the discrepancy is probably due to floating point arithmetic. Feel free to reopen this issue if necessary.