Add svd for real Symmetric and Hermitian matrices (#32017)

* add svd for realsym/hermitian, and tests

* shift improved code from svd.jl to symmetric.jl

* minimize allocations

* combine tests

* move test out of loop

* clean up svd tests

Author: dkarrasch

stdlib/LinearAlgebra/src/symmetric.jl

@@ -689,6 +689,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)

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)