From 28aaafc46349b5e49c33d22775ea07dae39dbd03 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Wed, 15 May 2024 23:32:33 +0530 Subject: [PATCH] Reland " Use copyto! in converting Diagonal/Bidiagonal/Tridiagonal to Matrix #53912 " (#54459) --- stdlib/LinearAlgebra/src/bidiag.jl | 19 ++++++--------- stdlib/LinearAlgebra/src/diagonal.jl | 11 +++++---- stdlib/LinearAlgebra/src/tridiag.jl | 33 +++++++++++++++------------ stdlib/LinearAlgebra/test/bidiag.jl | 7 ++++++ stdlib/LinearAlgebra/test/diagonal.jl | 7 ++++++ stdlib/LinearAlgebra/test/special.jl | 2 ++ 6 files changed, 47 insertions(+), 32 deletions(-) diff --git a/stdlib/LinearAlgebra/src/bidiag.jl b/stdlib/LinearAlgebra/src/bidiag.jl index 7236488563ca7..53fbff5839630 100644 --- a/stdlib/LinearAlgebra/src/bidiag.jl +++ b/stdlib/LinearAlgebra/src/bidiag.jl @@ -195,19 +195,14 @@ end #Converting from Bidiagonal to dense Matrix function Matrix{T}(A::Bidiagonal) where T - n = size(A, 1) - B = Matrix{T}(undef, n, n) - n == 0 && return B - n > 1 && fill!(B, zero(T)) - @inbounds for i = 1:n - 1 - B[i,i] = A.dv[i] - if A.uplo == 'U' - B[i,i+1] = A.ev[i] - else - B[i+1,i] = A.ev[i] - end + B = Matrix{T}(undef, size(A)) + if haszero(T) # optimized path for types with zero(T) defined + size(B,1) > 1 && fill!(B, zero(T)) + copyto!(view(B, diagind(B)), A.dv) + copyto!(view(B, diagind(B, A.uplo == 'U' ? 1 : -1)), A.ev) + else + copyto!(B, A) end - B[n,n] = A.dv[n] return B end Matrix(A::Bidiagonal{T}) where {T} = Matrix{promote_type(T, typeof(zero(T)))}(A) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index d68a691f8a149..bb722d20bb05b 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -116,11 +116,12 @@ AbstractMatrix{T}(D::Diagonal{T}) where {T} = copy(D) Matrix(D::Diagonal{T}) where {T} = Matrix{promote_type(T, typeof(zero(T)))}(D) Array(D::Diagonal{T}) where {T} = Matrix(D) function Matrix{T}(D::Diagonal) where {T} - n = size(D, 1) - B = Matrix{T}(undef, n, n) - n > 1 && fill!(B, zero(T)) - @inbounds for i in 1:n - B[i,i] = D.diag[i] + B = Matrix{T}(undef, size(D)) + if haszero(T) # optimized path for types with zero(T) defined + size(B,1) > 1 && fill!(B, zero(T)) + copyto!(view(B, diagind(B)), D.diag) + else + copyto!(B, D) end return B end diff --git a/stdlib/LinearAlgebra/src/tridiag.jl b/stdlib/LinearAlgebra/src/tridiag.jl index 081232b08a46a..ffdb56ba2dd77 100644 --- a/stdlib/LinearAlgebra/src/tridiag.jl +++ b/stdlib/LinearAlgebra/src/tridiag.jl @@ -135,13 +135,17 @@ function Matrix{T}(M::SymTridiagonal) where T n = size(M, 1) Mf = Matrix{T}(undef, n, n) n == 0 && return Mf - n > 2 && fill!(Mf, zero(T)) - @inbounds for i = 1:n-1 - Mf[i,i] = symmetric(M.dv[i], :U) - Mf[i+1,i] = transpose(M.ev[i]) - Mf[i,i+1] = M.ev[i] + if haszero(T) # optimized path for types with zero(T) defined + n > 2 && fill!(Mf, zero(T)) + @inbounds for i = 1:n-1 + Mf[i,i] = symmetric(M.dv[i], :U) + Mf[i+1,i] = transpose(M.ev[i]) + Mf[i,i+1] = M.ev[i] + end + Mf[n,n] = symmetric(M.dv[n], :U) + else + copyto!(Mf, M) end - Mf[n,n] = symmetric(M.dv[n], :U) return Mf end Matrix(M::SymTridiagonal{T}) where {T} = Matrix{promote_type(T, typeof(zero(T)))}(M) @@ -586,15 +590,14 @@ axes(M::Tridiagonal) = (ax = axes(M.d,1); (ax, ax)) function Matrix{T}(M::Tridiagonal) where {T} A = Matrix{T}(undef, size(M)) - n = length(M.d) - n == 0 && return A - n > 2 && fill!(A, zero(T)) - for i in 1:n-1 - A[i,i] = M.d[i] - A[i+1,i] = M.dl[i] - A[i,i+1] = M.du[i] - end - A[n,n] = M.d[n] + if haszero(T) # optimized path for types with zero(T) defined + size(A,1) > 2 && fill!(A, zero(T)) + copyto!(view(A, diagind(A)), M.d) + copyto!(view(A, diagind(A,1)), M.du) + copyto!(view(A, diagind(A,-1)), M.dl) + else + copyto!(A, M) + end A end Matrix(M::Tridiagonal{T}) where {T} = Matrix{promote_type(T, typeof(zero(T)))}(M) diff --git a/stdlib/LinearAlgebra/test/bidiag.jl b/stdlib/LinearAlgebra/test/bidiag.jl index 80a1b27ead829..cf6c0845194d2 100644 --- a/stdlib/LinearAlgebra/test/bidiag.jl +++ b/stdlib/LinearAlgebra/test/bidiag.jl @@ -904,4 +904,11 @@ end @test mul!(C1, B, sv, 1, 2) == mul!(C2, B, v, 1 ,2) end +@testset "Matrix conversion for non-numeric" begin + B = Bidiagonal(fill(Diagonal([1,3]), 3), fill(Diagonal([1,3]), 2), :U) + M = Matrix{eltype(B)}(B) + @test M isa Matrix{eltype(B)} + @test M == B +end + end # module TestBidiagonal diff --git a/stdlib/LinearAlgebra/test/diagonal.jl b/stdlib/LinearAlgebra/test/diagonal.jl index d45aa5e7fad98..366ecdecb9b24 100644 --- a/stdlib/LinearAlgebra/test/diagonal.jl +++ b/stdlib/LinearAlgebra/test/diagonal.jl @@ -1298,4 +1298,11 @@ end @test yadj == x' end +@testset "Matrix conversion for non-numeric" begin + D = Diagonal(fill(Diagonal([1,3]), 2)) + M = Matrix{eltype(D)}(D) + @test M isa Matrix{eltype(D)} + @test M == D +end + end # module TestDiagonal diff --git a/stdlib/LinearAlgebra/test/special.jl b/stdlib/LinearAlgebra/test/special.jl index c6fb187add3fb..6df27e00cc81a 100644 --- a/stdlib/LinearAlgebra/test/special.jl +++ b/stdlib/LinearAlgebra/test/special.jl @@ -113,6 +113,8 @@ Random.seed!(1) Base.convert(::Type{TypeWithZero}, ::TypeWithoutZero) = TypeWithZero() Base.zero(::Type{<:Union{TypeWithoutZero, TypeWithZero}}) = TypeWithZero() LinearAlgebra.symmetric(::TypeWithoutZero, ::Symbol) = TypeWithoutZero() + LinearAlgebra.symmetric_type(::Type{TypeWithoutZero}) = TypeWithoutZero + Base.copy(A::TypeWithoutZero) = A Base.transpose(::TypeWithoutZero) = TypeWithoutZero() d = fill(TypeWithoutZero(), 3) du = fill(TypeWithoutZero(), 2)