diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 14bae24d8c020..2220ef250102c 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -677,6 +677,22 @@ eigvals!(A::Hermitian{T,S}, B::Hermitian{T,S}) where {T<:BlasComplex,S<:StridedM eigvecs(A::HermOrSym) = eigvecs(eigen(A)) +function svd(A::RealHermSymComplexHerm, full::Bool=false) + vals, vecs = eigen(A) + I = sortperm(vals; by=abs, rev=true) + permute!(vals, I) + Base.permutecols!!(vecs, I) # left-singular vectors + V = copy(vecs) # right-singular vectors + # shifting -1 from singular values to right-singular vectors + @inbounds for i = 1:length(vals) + if vals[i] < 0 + vals[i] = -vals[i] + for j = 1:size(V,1); V[j,i] = -V[j,i]; end + end + end + return SVD(vecs, vals, V') +end + function svdvals!(A::RealHermSymComplexHerm) vals = eigvals!(A) for i = 1:length(vals) diff --git a/stdlib/LinearAlgebra/test/svd.jl b/stdlib/LinearAlgebra/test/svd.jl index 71f9a41caeef1..7a29d86974a61 100644 --- a/stdlib/LinearAlgebra/test/svd.jl +++ b/stdlib/LinearAlgebra/test/svd.jl @@ -35,25 +35,15 @@ end n = 10 -# Split n into 2 parts for tests needing two matrices -n1 = div(n, 2) -n2 = 2*n1 - Random.seed!(1234321) areal = randn(n,n)/2 aimg = randn(n,n)/2 -a2real = randn(n,n)/2 -a2img = randn(n,n)/2 @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, Int) aa = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal) - aa2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real) asym = aa' + aa # symmetric indefinite - apd = aa' * aa # symmetric positive-definite - for (a, a2) in ((aa, aa2), (view(aa, 1:n, 1:n), view(aa2, 1:n, 1:n))) - ε = εa = eps(abs(float(one(eltya)))) - + for a in (aa, view(aa, 1:n, 1:n)) usv = svd(a) @testset "singular value decomposition" begin @test usv.S === svdvals(usv) @@ -72,28 +62,20 @@ a2img = randn(n,n)/2 @test svdz.Vt ≈ Matrix{eltya}(I, 0, 0) end end - usv = svd(a') - @testset "singular value decomposition of adjoint" begin - @test usv.S === svdvals(usv) - @test usv.U * (Diagonal(usv.S) * usv.Vt) ≈ a' - @test convert(Array, usv) ≈ a' - @test usv.Vt' ≈ usv.V - @test_throws ErrorException usv.Z - b = rand(eltya,n) - @test usv\b ≈ a'\b - end - usv = svd(transpose(a)) - @testset "singular value decomposition of transpose" begin - @test usv.S === svdvals(usv) - @test usv.U * (Diagonal(usv.S) * usv.Vt) ≈ transpose(a) - @test convert(Array, usv) ≈ transpose(a) - @test usv.Vt' ≈ usv.V - @test_throws ErrorException usv.Z - b = rand(eltya,n) - @test usv\b ≈ transpose(a)\b + @testset "singular value decomposition of adjoint/transpose" begin + for transform in (adjoint, transpose) + usv = svd(transform(a)) + @test usv.S === svdvals(usv) + @test usv.U * (Diagonal(usv.S) * usv.Vt) ≈ transform(a) + @test convert(Array, usv) ≈ transform(a) + @test usv.Vt' ≈ usv.V + @test_throws ErrorException usv.Z + b = rand(eltya,n) + @test usv\b ≈ transform(a)\b + end end @testset "Generalized svd" begin - a_svd = a[1:n1, :] + a_svd = a[1:div(n, 2), :] gsvd = svd(a,a_svd) @test gsvd.U*gsvd.D1*gsvd.R*gsvd.Q' ≈ a @test gsvd.V*gsvd.D2*gsvd.R*gsvd.Q' ≈ a_svd @@ -121,6 +103,18 @@ a2img = randn(n,n)/2 @test gsvd.V*gsvd.D2*gsvd.R*gsvd.Q' ≈ c end end + @testset "singular value decomposition of Hermitian/real-Symmetric" begin + for T in (eltya <: Real ? (Symmetric, Hermitian) : (Hermitian,)) + usv = svd(T(asym)) + @test usv.S === svdvals(usv) + @test usv.U * (Diagonal(usv.S) * usv.Vt) ≈ T(asym) + @test convert(Array, usv) ≈ T(asym) + @test usv.Vt' ≈ usv.V + @test_throws ErrorException usv.Z + b = rand(eltya,n) + @test usv\b ≈ T(asym)\b + end + end if eltya <: LinearAlgebra.BlasReal @testset "Number input" begin x, y = randn(eltya, 2)