From 9ca631f3266f875ceee08448e003965563cb727f Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sun, 22 Dec 2024 13:33:58 +0100 Subject: [PATCH 1/4] Make unitary matrices in `svd`/`eigen` of `Diagonal` be unitless --- src/diagonal.jl | 15 +++++++++------ test/diagonal.jl | 10 +++++++--- 2 files changed, 16 insertions(+), 9 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index 280caec1..d81a3dfb 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -931,15 +931,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)) + # 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)) @@ -961,7 +964,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) @@ -1019,13 +1022,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 diff --git a/test/diagonal.jl b/test/diagonal.jl index 1cd2e9df..9549e46d 100644 --- a/test/diagonal.jl +++ b/test/diagonal.jl @@ -478,14 +478,14 @@ Random.seed!(1) U, s, V = F @test map(x -> x.val, Matrix(F)) ≈ map(x -> x.val, Du) @test svdvals(Du) == s - @test U isa AbstractMatrix{<:Furlong{0}} - @test V isa AbstractMatrix{<:Furlong{0}} + @test U isa AbstractMatrix{<:AbstractFloat} + @test V isa AbstractMatrix{<:AbstractFloat} @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{<:AbstractFloat} end end @@ -858,6 +858,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 From 85a4a2a08e9142817d050227ad79e6f5d686a0c3 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sun, 22 Dec 2024 18:29:50 +0100 Subject: [PATCH 2/4] Update test/diagonal.jl --- test/diagonal.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/diagonal.jl b/test/diagonal.jl index 9549e46d..fcaf0bb7 100644 --- a/test/diagonal.jl +++ b/test/diagonal.jl @@ -478,7 +478,7 @@ Random.seed!(1) U, s, V = F @test map(x -> x.val, Matrix(F)) ≈ map(x -> x.val, Du) @test svdvals(Du) == s - @test U isa AbstractMatrix{<:AbstractFloat} + @test U isa AbstractMatrix{<:Union{Real,Complex}} @test V isa AbstractMatrix{<:AbstractFloat} @test s isa AbstractVector{<:Furlong{1}} E = eigen(Du) From b642425e20e64f80a4f6ddea96a5e57a5897f32e Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 24 Jan 2025 20:41:31 +0100 Subject: [PATCH 3/4] return tests --- test/unitful.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/unitful.jl b/test/unitful.jl index cd4ade52..8d6c4260 100644 --- a/test/unitful.jl +++ b/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 From 530fcee4fd7016506eb27ce7c03d828857ca2492 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 27 Jan 2025 11:35:53 +0100 Subject: [PATCH 4/4] Update src/diagonal.jl --- src/diagonal.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index d1a9ea00..8a70d103 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -911,7 +911,7 @@ function pinv(D::Diagonal{T}, tol::Real) where T end _ortho_eltype(T) = Base.promote_op(/, T, T) -_ortho_eltype(T::Type{<:Number}) = typeof(one(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?