Make unitary matrices in svd/eigen of Diagonal be unitless (#1155)

Author: dkarrasch

src/diagonal.jl

@@ -911,15 +911,18 @@ function pinv(D::Diagonal{T}, tol::Real) where T
     Diagonal(Di)
 end
 
+_ortho_eltype(T) = Base.promote_op(/, T, T)
+_ortho_eltype(T::Type{<:Number}) = typeof(one(T)/one(T))
+
 # TODO Docstrings for eigvals, eigvecs, eigen all mention permute, scale, sortby as keyword args
 # but not all of them below provide them. Do we need to fix that?
 #Eigensystem
 eigvals(D::Diagonal{<:Number}; permute::Bool=true, scale::Bool=true) = copy(D.diag)
 eigvals(D::Diagonal; permute::Bool=true, scale::Bool=true) =
     reduce(vcat, eigvals(x) for x in D.diag) #For block matrices, etc.
-function eigvecs(D::Diagonal{T}) where T<:AbstractMatrix
+function eigvecs(D::Diagonal{T}) where {T<:AbstractMatrix}
     diag_vecs = [ eigvecs(x) for x in D.diag ]
-    matT = reduce((a,b) -> promote_type(typeof(a),typeof(b)), diag_vecs)
+    matT = promote_type(map(typeof, diag_vecs)...)
     ncols_diag = [ size(x, 2) for x in D.diag ]
     nrows = size(D, 1)
     vecs = Matrix{Vector{eltype(matT)}}(undef, nrows, sum(ncols_diag))
@@ -941,7 +944,7 @@ function eigen(D::Diagonal; permute::Bool=true, scale::Bool=true, sortby::Union{
     if any(!isfinite, D.diag)
         throw(ArgumentError("matrix contains Infs or NaNs"))
     end
-    Td = Base.promote_op(/, eltype(D), eltype(D))
+    Td = _ortho_eltype(eltype(D))
     λ = eigvals(D)
     if !isnothing(sortby)
         p = sortperm(λ; alg=QuickSort, by=sortby)
@@ -999,13 +1002,13 @@ function svd(D::Diagonal{T}) where {T<:Number}
     s = abs.(d)
     piv = sortperm(s, rev = true)
     S = s[piv]
-    Td  = typeof(oneunit(T)/oneunit(T))
+    Td  = _ortho_eltype(T)
     U = zeros(Td, size(D))
     Vt = copy(U)
     for i in 1:length(d)
         j = piv[i]
-        U[j,i] = iszero(d[j]) ? oneunit(Td) : d[j] / S[i]
-        Vt[i,j] = oneunit(Td)
+        U[j,i] = iszero(d[j]) ? one(Td) : d[j] / S[i]
+        Vt[i,j] = one(Td)
     end
     return SVD(U, S, Vt)
 end

test/diagonal.jl

@@ -839,6 +839,10 @@ end
     @test eigD.values == evals
     @test eigD.vectors ≈ evecs
     @test D * eigD.vectors ≈ eigD.vectors * Diagonal(eigD.values)
+
+    # test concrete types
+    D = Diagonal([I2 for _ in 1:4])
+    @test eigen(D) isa Eigen{Vector{Float64}, Float64, Matrix{Vector{Float64}}, Vector{Float64}}
 end
 
 @testset "linear solve for block diagonal matrices" begin

test/unitful.jl

@@ -79,14 +79,14 @@ end
         U, s, V = F
         @test map(getval, Matrix(F)) ≈ map(getval, Du)
         @test svdvals(Du) == s
-        @test U isa AbstractMatrix{<:Furlong{0}}
-        @test V isa AbstractMatrix{<:Furlong{0}}
+        @test U isa AbstractMatrix{<:Union{Real,Complex}}
+        @test V isa AbstractMatrix{<:Union{Real,Complex}}
         @test s isa AbstractVector{<:Furlong{1}}
         E = eigen(Du)
         vals, vecs = E
         @test Matrix(E) == Du
         @test vals isa AbstractVector{<:Furlong{1}}
-        @test vecs isa AbstractMatrix{<:Furlong{0}}
+        @test vecs isa AbstractMatrix{<:Union{Real,Complex}}
     end
 end