JuliaLinearAlgebra / MKL.jl

Intel MKL linear algebra backend for Julia
Other
205 stars 32 forks source link

Discrepancy between OpenBLAS and MKL result #141

Closed andreasvarga closed 10 months ago

andreasvarga commented 10 months ago

Following the discussions here, I started to make my code in PeriodicSystems compliant with MKL. During tests, I came accross to the following issue:

julia> using LinearAlgebra

julia> B = [0.6950637459946359 0.4709911648119093]
1×2 Matrix{Float64}:
 0.695064  0.470991

julia> norm(B*B' - B*Matrix(B'),Inf)
0.0

julia> using MKL

julia> norm(B*B' - B*Matrix(B'),Inf)
1.1102230246251565e-16

While B*Matrix(B') is always a symmetric matrix, when computed with OpenBLAS, this is not anymore true with MKL. Is any explanation for this? This is unexpected for me, being simply a dot product computed probably in two different ways.

In my tests I relied on the fact that B*Matrix(B') is symmetric. For MKL compliance I needed to change to something like (B*Matrix(B')+B*Matrix(B'))/2` to ensure exact symmetry (not a big trouble, but still unexpected somehow).

Note that issymmetric(B*B') is always true for both OpenBLAS and MKL,

andreasvarga commented 10 months ago

I think this issue is not worth to spend time on it.