eigen and svd with Symmetric{Float16, ...}

[btmit]

eigen and svd of Symmetric{Float16, ...} returns values with element type Float32. This is not the case when the input is not wrapped in Symmetric. cholesky decompositions don't have this issue.

a = rand(Float16, 3, 3);
a = a*a'

svd(a)
# SVD{Float16, Float16, Matrix{Float16}, Vector{Float16}}

eigen(a)
# Eigen{Float16, Float16, Matrix{Float16}, Vector{Float16}}

cholesky(a)
# Cholesky{Float16, Matrix{Float16}}

svd(Symmetric(a))
# SVD{Float32, Float32, Matrix{Float32}, Vector{Float32}}

eigen(Symmetric(a))
# Eigen{Float32, Float32, Matrix{Float32}, Vector{Float32}}

cholesky(Symmetric(a))
# Cholesky{Float16, Matrix{Float16}}

[aravindh-krishnamoorthy]

It seems that for svd, in

julia/stdlib/LinearAlgebra/src/symmetric.jl

Line 758 in 10db278

function svd(A::RealHermSymComplexHerm; full::Bool=false)

and for eigen in

julia/stdlib/LinearAlgebra/src/symmetriceigen.jl

Line 11 in 10db278

function eigen(A::RealHermSymComplexHerm; sortby::Union{Function,Nothing}=nothing)

the matrix is promoted to Float32 for LAPACK but the result is not cast back to Float16.

Working cases: For plain vanilla svd(a), the result is cast back correctly in svd.jl and eigen(a) uses a specialised Float16 routine in eigen.jl. Furthermore, for cholesky, there's a special method for AbstractMatrix{Float16} here:

julia/stdlib/LinearAlgebra/src/cholesky.jl

Line 405 in 10db278

function cholesky(A::AbstractMatrix{Float16}, ::NoPivot=NoPivot(); check::Bool = true)

so the cholesky(Symmetric(a)) function is Ok.

[jishnub]

A PR to add an analogous method is welcome