Various fixes for complex SVD

Author: c42f

src/StaticArrays.jl

@@ -9,7 +9,7 @@ import Base: getindex, setindex!, size, similar, vec, show,
              mapreduce, broadcast, broadcast!, conj, transpose, ctranspose,
              hcat, vcat, ones, zeros, eye, one, cross, vecdot, reshape, fill,
              fill!, det, inv, eig, eigvals, expm, sqrtm, trace, vecnorm, norm, dot, diagm, diag,
-             lu, svd, svdvals!, svdfact,
+             lu, svd, svdvals, svdfact,
              sum, diff, prod, count, any, all, minimum,
              maximum, extrema, mean, copy, rand, randn, randexp, rand!, randn!,
              randexp!, normalize, normalize!, read, read!, write

src/svd.jl

@@ -13,14 +13,18 @@ getindex(::SVD, ::Symbol) = error("In order to avoid type instability, StaticArr
 
 function svdvals(A::StaticMatrix)
     sv = svdvals(Matrix(A))
-    similar_type(A, eltype(sv), Size(diagsize(A)))(sv)
+    # We should be using `T2=eltype(sv)`, but it's not inferrable for complex
+    # eltypes.  See https://github.com/JuliaLang/julia/pull/22443
+    T = eltype(A)
+    T2 = promote_type(Float32, real(typeof(one(T)/norm(one(T)))))
+    similar_type(A, T2, Size(diagsize(A)))(sv)
 end
 
 function svdfact(A::StaticMatrix)
     # "Thin" SVD only for now.
     f = svdfact(Matrix(A))
     U = similar_type(A, eltype(f.U), Size(Size(A)[1], diagsize(A)))(f.U)
-    S = similar_type(A, Size(diagsize(A)))(f.S)
+    S = similar_type(A, eltype(f.S), Size(diagsize(A)))(f.S)
     Vt = similar_type(A, eltype(f.Vt), Size(diagsize(A), Size(A)[2]))(f.Vt)
     SVD(U,S,Vt)
 end

test/svd.jl

@@ -1,36 +1,47 @@
 using StaticArrays, Base.Test
 
 @testset "SVD factorization" begin
+    m3 = @SMatrix Float64[3 9 4; 6 6 2; 3 7 9]
+    m3c = Complex128.(m3)
+    m23 = @SMatrix Float64[3 9 4; 6 6 2]
+
     @testset "svd" begin
         @testinf svdvals(@SMatrix [2 0; 0 0])::StaticVector ≊ [2, 0]
         @testinf svdvals((@SMatrix [2 -2; 1 1]) / sqrt(2)) ≊ [2, 1]
 
-        m3 = @SMatrix Float64[3 9 4; 6 6 2; 3 7 9]
-        @testinf svdvals(m3) ≈ svdvals(Matrix(m3))
+        @testinf svdvals(m3) ≊ svdvals(Matrix(m3))
+        @testinf svdvals(m3c) isa SVector{3,Float64}
 
-        @testinf svd(m3)[1]::StaticMatrix ≈ svd(Matrix(m3))[1]
-        @testinf svd(m3)[2]::StaticVector ≈ svd(Matrix(m3))[2]
-        @testinf svd(m3)[3]::StaticMatrix ≈ svd(Matrix(m3))[3]
+        @testinf svd(m3)[1]::StaticMatrix ≊ svd(Matrix(m3))[1]
+        @testinf svd(m3)[2]::StaticVector ≊ svd(Matrix(m3))[2]
+        @testinf svd(m3)[3]::StaticMatrix ≊ svd(Matrix(m3))[3]
     end
 
     @testset "svdfact" begin
         @test_throws ErrorException svdfact(@SMatrix [1 0; 0 1])[:U]
 
         @testinf svdfact(@SMatrix [2 0; 0 0]).U === eye(SMatrix{2,2})
-        @testinf svdfact(@SMatrix [2 0; 0 0]).S === SVector(2, 0)
+        @testinf svdfact(@SMatrix [2 0; 0 0]).S === SVector(2.0, 0.0)
         @testinf svdfact(@SMatrix [2 0; 0 0]).Vt === eye(SMatrix{2,2})
 
         @testinf svdfact((@SMatrix [2 -2; 1 1]) / sqrt(2)).U ≊ [-1 0; 0 1]
         @testinf svdfact((@SMatrix [2 -2; 1 1]) / sqrt(2)).S ≊ [2, 1]
         @testinf svdfact((@SMatrix [2 -2; 1 1]) / sqrt(2)).Vt ≊ [-1 1; 1 1]/sqrt(2)
 
-        m23 = @SMatrix Float64[3 9 4; 6 6 2]
         @testinf svdfact(m23).U  ≊ svdfact(Matrix(m23))[:U]
         @testinf svdfact(m23).S  ≊ svdfact(Matrix(m23))[:S]
         @testinf svdfact(m23).Vt ≊ svdfact(Matrix(m23))[:Vt]
 
         @testinf svdfact(m23').U  ≊ svdfact(Matrix(m23'))[:U]
         @testinf svdfact(m23').S  ≊ svdfact(Matrix(m23'))[:S]
         @testinf svdfact(m23').Vt ≊ svdfact(Matrix(m23'))[:Vt]
+
+        @testinf svdfact(m3c).U  ≊ svdfact(Matrix(m3c))[:U]
+        @testinf svdfact(m3c).S  ≊ svdfact(Matrix(m3c))[:S]
+        @testinf svdfact(m3c).Vt ≊ svdfact(Matrix(m3c))[:Vt]
+
+        @testinf svdfact(m3c).U  isa SMatrix{3,3,Complex128}
+        @testinf svdfact(m3c).S  isa SVector{3,Float64}
+        @testinf svdfact(m3c).Vt isa SMatrix{3,3,Complex128}
     end
 end