oneMKL: test failures with rotate and reflect

[kballeda]

Hi @maleadt , I am facing following errors when I tested latest master in my local setup.
Looks like couple of rotate/reflect primitives are failing for complexf64 type at level-1.

~/Kali/2023/2002_fp16_rebase/oneAPI.jl$ $JULIA --project -L test/setup.jl test/onemkl.jl 
rotate: Test Failed at /home/kali/Kali/2023/2002_fp16_rebase/oneAPI.jl/test/onemkl.jl:33
  Expression: testf(rotate!, rand(T, m), rand(T, m), rand(real(T)), rand(T))
Stacktrace:
 [1] macro expansion
   @ ~/Kali/oneAPI.jl/julia-1.8.5/share/julia/stdlib/v1.8/Test/src/Test.jl:464 [inlined]
 [2] macro expansion
   @ ~/Kali/2023/2002_fp16_rebase/oneAPI.jl/test/onemkl.jl:33 [inlined]
 [3] macro expansion
   @ ~/Kali/oneAPI.jl/julia-1.8.5/share/julia/stdlib/v1.8/Test/src/Test.jl:1363 [inlined]
 [4] macro expansion
   @ ~/Kali/2023/2002_fp16_rebase/oneAPI.jl/test/onemkl.jl:32 [inlined]
 [5] macro expansion
   @ ~/Kali/oneAPI.jl/julia-1.8.5/share/julia/stdlib/v1.8/Test/src/Test.jl:1439 [inlined]
 [6] macro expansion
   @ ~/Kali/2023/2002_fp16_rebase/oneAPI.jl/test/onemkl.jl:12 [inlined]
 [7] macro expansion
   @ ~/Kali/oneAPI.jl/julia-1.8.5/share/julia/stdlib/v1.8/Test/src/Test.jl:1363 [inlined]
 [8] top-level scope
   @ ~/Kali/2023/2002_fp16_rebase/oneAPI.jl/test/onemkl.jl:12
reflect: Test Failed at /home/kali/Kali/2023/2002_fp16_rebase/oneAPI.jl/test/onemkl.jl:38
  Expression: testf(reflect!, rand(T, m), rand(T, m), rand(real(T)), rand(T))
Stacktrace:
 [1] macro expansion
   @ ~/Kali/oneAPI.jl/julia-1.8.5/share/julia/stdlib/v1.8/Test/src/Test.jl:464 [inlined]
 [2] macro expansion
   @ ~/Kali/2023/2002_fp16_rebase/oneAPI.jl/test/onemkl.jl:38 [inlined]
 [3] macro expansion
   @ ~/Kali/oneAPI.jl/julia-1.8.5/share/julia/stdlib/v1.8/Test/src/Test.jl:1363 [inlined]
 [4] macro expansion
   @ ~/Kali/2023/2002_fp16_rebase/oneAPI.jl/test/onemkl.jl:37 [inlined]
 [5] macro expansion
   @ ~/Kali/oneAPI.jl/julia-1.8.5/share/julia/stdlib/v1.8/Test/src/Test.jl:1439 [inlined]
 [6] macro expansion
   @ ~/Kali/2023/2002_fp16_rebase/oneAPI.jl/test/onemkl.jl:12 [inlined]
 [7] macro expansion
   @ ~/Kali/oneAPI.jl/julia-1.8.5/share/julia/stdlib/v1.8/Test/src/Test.jl:1363 [inlined]
 [8] top-level scope
   @ ~/Kali/2023/2002_fp16_rebase/oneAPI.jl/test/onemkl.jl:12

Test Environment:
OS: "Ubuntu 20.04.5 LTS"
CPU: i9-9900K
GPU: Intel(R) UHD Graphics 630 [0x3e98]

[maleadt]

Looks like a bug in oneMKL? MWE:

using oneAPI, LinearAlgebra

function main(; n = 1, T = ComplexF64)
    x = ones(T, n)
    y = ones(T, n)
    c = one(real(T))
    s = one(T)

    d_x = oneArray(x)
    d_y = oneArray(y)

    rotate!(x, y, c, s)

    #rotate!(d_x, d_y, c, s)
    queue = global_queue(context(d_x), device(d_x))
    incx = stride(x, 1)
    incy = stride(y, 1)
    @assert length(d_x) >= 1 + (n - 1)*abs(incx)
    @assert length(d_y) >= 1 + (n - 1)*abs(incy)
    if T == ComplexF64
        oneMKL.onemklZrot(sycl_queue(queue), n, d_x, incx, d_y, incy, c, s)
    elseif T == ComplexF32
        oneMKL.onemklCrot(sycl_queue(queue), n, d_x, incx, d_y, incy, c, s)
    end

    @show x
    @show Array(d_x)
    @assert x == Array(d_x)
end

The C++ wrappers seem fine... Even with ComplexF32, it occasionally fails:

julia> main(; T=ComplexF32)
x = ComplexF32[2.0f0 + 0.0f0im]
Array(d_x) = ComplexF32[1.0f0 + 0.0f0im]
ERROR: AssertionError: x == Array(d_x)
Stacktrace:
 [1] main(; n::Int64, T::Type)
   @ Main ~/Julia/pkg/oneAPI/wip.jl:30
 [2] top-level scope
   @ REPL[10]:1

julia> main(; T=ComplexF32)
x = ComplexF32[2.0f0 + 0.0f0im]
Array(d_x) = ComplexF32[2.0f0 + 0.0f0im]

julia> main(; T=ComplexF32)
x = ComplexF32[2.0f0 + 0.0f0im]
Array(d_x) = ComplexF32[2.0f0 + 0.0f0im]

julia> main(; T=ComplexF32)
x = ComplexF32[2.0f0 + 0.0f0im]
Array(d_x) = ComplexF32[2.0f0 + 0.0f0im]

julia> main(; T=ComplexF32)
x = ComplexF32[2.0f0 + 0.0f0im]
Array(d_x) = ComplexF32[2.0f0 + 0.0f0im]

julia> main(; T=ComplexF32)
x = ComplexF32[2.0f0 + 0.0f0im]
Array(d_x) = ComplexF32[1.0f0 + 0.0f0im]
ERROR: AssertionError: x == Array(d_x)
Stacktrace:
 [1] main(; n::Int64, T::Type)
   @ Main ~/Julia/pkg/oneAPI/wip.jl:30
 [2] top-level scope
   @ REPL[10]:1

cc @amontoison, in case you're already relying on these.

[amontoison]

We use the rot routine with c and s when we have the relation |c|^2 + |s|^2 = 1 because they represent cosines and sines of a 2x2 Givens rotation.
Maybe oneMKL relies on this relation to reduce numerical errors and we should update our tests to satisfy it.

[maleadt]

Maybe oneMKL relies on this relation to reduce numerical errors and we should update our tests to satisfy it.

Would that explain it switching from ComplexF32[1.0f0 + 0.0f0im] to ComplexF32[2.0f0 + 0.0f0im] randomly? I'm not using random inputs in the MWE above.

[amontoison]

I fixed a typo with #284 but the error is not related.

using oneAPI, LinearAlgebra, MKL

function main(; n = 1, T = ComplexF64)
    x = ones(T, n)
    y = ones(T, n)
    c = rand(real(T))
    s = T(sqrt(1 - c^2)) 

    @test c*c + s*conj(s) ≈ 1
    d_x = oneArray(x)
    d_y = oneArray(y)
    
    incx = stride(x, 1)
    incy = stride(y, 1)
    BLAS.rot!(n, x, incx, y, incy, c, s)

    queue = global_queue(context(d_x), device(d_x))
    incx = stride(d_x, 1)
    incy = stride(d_y, 1)
    @assert length(d_x) >= 1 + (n - 1)*abs(incx)
    @assert length(d_y) >= 1 + (n - 1)*abs(incy)

    (T == Float32)    && oneMKL.onemklSrot(sycl_queue(queue), n, d_x, incx, d_y, incy, c, s)
    (T == Float64)    && oneMKL.onemklDrot(sycl_queue(queue), n, d_x, incx, d_y, incy, c, s)
    (T == ComplexF32) && oneMKL.onemklCrot(sycl_queue(queue), n, d_x, incx, d_y, incy, c, s)
    (T == ComplexF64) && oneMKL.onemklZrot(sycl_queue(queue), n, d_x, incx, d_y, incy, c, s)

    @show x
    @show Array(d_x)
    @assert x == Array(d_x)
end

What I don't understand is that we have a different result if we use MKL on CPU.
Is it not the same source code that is used by both versions (CPU / GPU)?

I can't explain that we have random switch from ComplexF32[1.0f0 + 0.0f0im] to ComplexF32[2.0f0 + 0.0f0im] without random inputs in your example above. 🤔
Do we have a way to contact Intel like Nvidia?

[maleadt]

@kballeda is from Intel so probably can elaborate. Otherwise @pengtu.

[pengtu]

@maleadt @amontoison We will take this to the oneMKL team. Thank you for the investigation and the very good test cases.

[maleadt]

IIRC this has been fixed by an update of the support library.