Improve performance of svd and eigen of Diagonals (#43856)

stdlib/LinearAlgebra/src/diagonal.jl

@@ -675,20 +675,38 @@ 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
- Eigen(sorteig!(eigvals(D), eigvecs(D), sortby)...)
+ Td = Base.promote_op(/, eltype(D), eltype(D))
+ λ = eigvals(D)
+ if !isnothing(sortby)
+ p = sortperm(λ; alg=QuickSort, by=sortby)
+ λ = λ[p] # make a copy, otherwise this permutes D.diag
+ evecs = zeros(Td, size(D))
+ @inbounds for i in eachindex(p)
+ evecs[p[i],i] = one(Td)
+ end
+ else
+ evecs = Matrix{Td}(I, size(D))
+ end
+ Eigen(λ, evecs)
 end
 
 #Singular system
 svdvals(D::Diagonal{<:Number}) = sort!(abs.(D.diag), rev = true)
 svdvals(D::Diagonal) = [svdvals(v) for v in D.diag]
-function svd(D::Diagonal{T}) where T<:Number
- S = abs.(D.diag)
- piv = sortperm(S, rev = true)
- U = Diagonal(D.diag ./ S)
- Up = hcat([U[:,i] for i = 1:length(D.diag)][piv]...)
- V = Diagonal(fill!(similar(D.diag), one(T)))
- Vp = hcat([V[:,i] for i = 1:length(D.diag)][piv]...)
- return SVD(Up, S[piv], copy(Vp'))
+function svd(D::Diagonal{T}) where {T<:Number}
+ d = D.diag
+ s = abs.(d)
+ piv = sortperm(s, rev = true)
+ S = s[piv]
+ Td = typeof(oneunit(T)/oneunit(T))
+ U = zeros(Td, size(D))
+ Vt = copy(U)
+ for i in 1:length(d)
+ j = piv[i]
+ U[j,i] = d[j] / S[i]
+ Vt[i,j] = one(Td)
+ end
+ return SVD(U, S, Vt)
 end
 
 # disambiguation methods: * and / of Diagonal and Adj/Trans AbsVec

stdlib/LinearAlgebra/test/diagonal.jl

@@ -5,6 +5,10 @@ module TestDiagonal
 using Test, LinearAlgebra, Random
 using LinearAlgebra: BlasFloat, BlasComplex
 
+const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
+isdefined(Main, :Furlongs) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "Furlongs.jl"))
+using .Main.Furlongs
+
 n=12 #Size of matrix problem to test
 Random.seed!(1)
 
@@ -344,8 +348,12 @@ Random.seed!(1)
 
  @testset "Eigensystem" begin
  eigD = eigen(D)
- @test Diagonal(eigD.values) ≈ D
+ @test Diagonal(eigD.values) == D
  @test eigD.vectors == Matrix(I, size(D))
+ eigsortD = eigen(D, sortby=LinearAlgebra.eigsortby)
+ @test eigsortD.values !== D.diag
+ @test eigsortD.values == sort(D.diag, by=LinearAlgebra.eigsortby)
+ @test Matrix(eigsortD) == D
  end
 
  @testset "ldiv" begin
@@ -411,6 +419,22 @@ Random.seed!(1)
  @test svd(D).V == V
  end
 
+ @testset "svd/eigen with Diagonal{Furlong}" begin
+ Du = Furlong.(D)
+ @test Du isa Diagonal{<:Furlong{1}}
+ F = svd(Du)
+ 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 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}}
+ end
 end
 
 @testset "rdiv! (#40887)" begin