Float16 input now gives Float 16 factors for svd,eigen and cholesky. (JuliaLang#41352)

Co-authored-by: Daniel Karrasch <[email protected]>
Co-authored-by: Fredrik Bagge Carlson <[email protected]>

Author: 3 people

stdlib/LinearAlgebra/src/cholesky.jl

@@ -388,6 +388,11 @@ true
 cholesky(A::Union{StridedMatrix,RealHermSymComplexHerm{<:Real,<:StridedMatrix}},
     ::Val{false}=Val(false); check::Bool = true) = cholesky!(cholcopy(A); check = check)
 
+function cholesky(A::Union{StridedMatrix{Float16},RealHermSymComplexHerm{Float16,<:StridedMatrix}}, ::Val{false}=Val(false); check::Bool = true)
+    X = cholesky!(cholcopy(A); check = check)
+    return Cholesky{Float16}(X)
+end
+
 
 ## With pivoting
 """

stdlib/LinearAlgebra/src/eigen.jl

@@ -275,6 +275,17 @@ function eigen(A::AbstractMatrix{T}; permute::Bool=true, scale::Bool=true, sortb
     isdiag(AA) && return eigen(Diagonal(AA); permute=permute, scale=scale, sortby=sortby)
     return eigen!(AA; permute=permute, scale=scale, sortby=sortby, jvl=jvl, jvr=jvr, jce=jce, jcv=jcv)
 end
+function eigen(A::AbstractMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where {T <: Union{Float16,Complex{Float16}}}
+    AA = copy_oftype(A, eigtype(T))
+    isdiag(AA) && return eigen(Diagonal(AA); permute=permute, scale=scale, sortby=sortby)
+    A = eigen!(AA; permute, scale, sortby)
+    values = convert(AbstractVector{isreal(A.values) ? Float16 : Complex{Float16}}, A.values)
+    vectors = convert(AbstractMatrix{isreal(A.vectors) ? Float16 : Complex{Float16}}, A.vectors)
+    vectorsl = convert(AbstractMatrix{isreal(A.vectors) ? Float16 : Complex{Float16}}, A.vectorsl)
+    rconde = convert(Vector{Float16}, A.rconde)
+    rcondv = convert(Vector{Float16}, A.rcondv)
+    return Eigen(values, vectors, vectorsl, A.unitary, rconde, rcondv)
+end
 eigen(x::Number) = Eigen([x], fill(one(x), 1, 1))
 
 """

stdlib/LinearAlgebra/src/svd.jl

@@ -175,6 +175,10 @@ true
 function svd(A::StridedVecOrMat{T}; full::Bool = false, alg::Algorithm = default_svd_alg(A)) where {T}
     svd!(copy_oftype(A, eigtype(T)), full = full, alg = alg)
 end
+function svd(A::StridedVecOrMat{T}; full::Bool = false, alg::Algorithm = default_svd_alg(A)) where {T <: Union{Float16,Complex{Float16}}}
+    A = svd!(copy_oftype(A, eigtype(T)), full = full, alg = alg)
+    return SVD{T}(A)
+end
 function svd(x::Number; full::Bool = false, alg::Algorithm = default_svd_alg(x))
     SVD(x == 0 ? fill(one(x), 1, 1) : fill(x/abs(x), 1, 1), [abs(x)], fill(one(x), 1, 1))
 end

stdlib/LinearAlgebra/test/cholesky.jl

@@ -483,4 +483,17 @@ end
     @test F\b == F'\b
 end
 
+@testset "Float16" begin
+    A = Float16[4. 12. -16.; 12. 37. -43.; -16. -43. 98.]
+    B = cholesky(A)
+    B32 = cholesky(Float32.(A))
+    @test B isa Cholesky{Float16, Matrix{Float16}}
+    @test B.U isa UpperTriangular{Float16, Matrix{Float16}}
+    @test B.L isa LowerTriangular{Float16, Matrix{Float16}}
+    @test B.UL isa UpperTriangular{Float16, Matrix{Float16}}
+    @test B.U ≈ B32.U
+    @test B.L ≈ B32.L
+    @test B.UL ≈ B32.UL
+end
+
 end # module TestCholesky

stdlib/LinearAlgebra/test/eigen.jl

@@ -177,4 +177,36 @@ end
     @test isequal(eigen(A), eigen(A))
 end
 
+@testset "Float16" begin
+    A = Float16[4. 12. -16.; 12. 37. -43.; -16. -43. 98.]
+    B = eigen(A)
+    B32 = eigen(Float32.(A))
+    C = Float16[3 -2; 4 -1]
+    D = eigen(C)
+    D32 = eigen(Float32.(C))
+    F = eigen(complex(C))
+    F32 = eigen(complex(Float32.(C)))
+    @test B isa Eigen{Float16, Float16, Matrix{Float16}, Vector{Float16}, Vector{Float16}}
+    @test B.values isa Vector{Float16}
+    @test B.vectors isa Matrix{Float16}
+    @test B.vectorsl isa Matrix{Float16}
+    @test B.values ≈ B32.values
+    @test B.vectors ≈ B32.vectors
+    @test B.vectorsl ≈ B32.vectorsl
+    @test D isa Eigen{ComplexF16, ComplexF16, Matrix{ComplexF16}, Vector{ComplexF16}, Vector{Float16}}
+    @test D.values isa Vector{ComplexF16}
+    @test D.vectors isa Matrix{ComplexF16}
+    @test D.vectorsl isa Matrix{ComplexF16}
+    @test D.values ≈ D32.values
+    @test D.vectors ≈ D32.vectors
+    @test D.vectorsl ≈ D32.vectorsl
+    @test F isa Eigen{ComplexF16, ComplexF16, Matrix{ComplexF16}, Vector{ComplexF16}, Vector{Float16}}
+    @test F.values isa Vector{ComplexF16}
+    @test F.vectors isa Matrix{ComplexF16}
+    @test F.vectorsl isa Matrix{ComplexF16}
+    @test F.values ≈ F32.values
+    @test F.vectors ≈ F32.vectors
+    @test F.vectorsl ≈ F32.vectorsl
+end
+
 end # module TestEigen

stdlib/LinearAlgebra/test/svd.jl

@@ -248,4 +248,17 @@ end
     end
 end
 
+@testset "Float16" begin
+    A = Float16[4. 12. -16.; 12. 37. -43.; -16. -43. 98.]
+    B = svd(A)
+    B32 = svd(Float32.(A))
+    @test B isa SVD{Float16, Float16, Matrix{Float16}}
+    @test B.U isa Matrix{Float16}
+    @test B.Vt isa Matrix{Float16}
+    @test B.S isa Vector{Float16}
+    @test B.U ≈ B32.U
+    @test B.Vt ≈ B32.Vt
+    @test B.S ≈ B32.S
+end
+
 end # module TestSVD