Non-symmetric inverse of Cholesky decomposition

[astrozot]

Because of rounding errors, in some cases explicitly symmetric matrices might produce non-symmetric inverses:

using StaticArrays, LinearAlgebra

c = cholesky(@SMatrix[π 3.1; 3.1 12.0])
d = inv(c)
d[1,2] - d[2,1] # returns -1.3877787807814457e-17

This can be a serious issue and it can break various codes: for example, cholesky(d) would fail. Note that the same code works without problems if one uses standard arrays instead of static ones:

c1 = cholesky([π 3.1; 3.1 12.0])
d1 = inv(c1)
d1[1,2] - d1[2,1] # returns 0.0

A fix is to force, in the definition of inv(::Cholesky)

https://github.com/JuliaArrays/StaticArrays.jl/blob/c4092a1a8866313d0cc0d93f0bb870e5d82a2070/src/cholesky.jl#L72-L74

the symmetry of the result:

function inv(A::Cholesky{T,<:StaticMatrix{N,N,T}}) where {N,T}
    B = A.U \ (A.U' \ SDiagonal{N}(I))
    return (B .+ B') ./ T(2)
end