Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Use copyto! in converting Diagonal/Bidiagonal/Tridiagonal to Matrix #53912

Merged
merged 2 commits into from
Apr 4, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 7 additions & 12 deletions stdlib/LinearAlgebra/src/bidiag.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
11 changes: 6 additions & 5 deletions stdlib/LinearAlgebra/src/diagonal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
33 changes: 18 additions & 15 deletions stdlib/LinearAlgebra/src/tridiag.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
13 changes: 13 additions & 0 deletions stdlib/LinearAlgebra/test/bidiag.jl
Original file line number Diff line number Diff line change
Expand Up @@ -884,4 +884,17 @@ end
@test mul!(C1, B, sv, 1, 2) == mul!(C2, B, v, 1 ,2)
end

@testset "Matrix conversion for non-numeric and undef" begin
B = Bidiagonal(Vector{BigInt}(undef, 4), fill(big(3), 3), :U)
M = Matrix(B)
B[diagind(B)] .= 4
M[diagind(M)] .= 4
@test diag(B) == diag(M)

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
13 changes: 13 additions & 0 deletions stdlib/LinearAlgebra/test/diagonal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1292,4 +1292,17 @@ end
@test yadj == x'
end

@testset "Matrix conversion for non-numeric and undef" begin
D = Diagonal(Vector{BigInt}(undef, 4))
M = Matrix(D)
D[diagind(D)] .= 4
M[diagind(M)] .= 4
@test diag(D) == diag(M)

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
2 changes: 2 additions & 0 deletions stdlib/LinearAlgebra/test/special.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
8 changes: 8 additions & 0 deletions stdlib/LinearAlgebra/test/tridiag.jl
Original file line number Diff line number Diff line change
Expand Up @@ -830,4 +830,12 @@ end
@test axes(B) === (ax, ax)
end

@testset "Matrix conversion for non-numeric and undef" begin
T = Tridiagonal(fill(big(3), 3), Vector{BigInt}(undef, 4), fill(big(3), 3))
M = Matrix(T)
T[diagind(T)] .= 4
M[diagind(M)] .= 4
@test diag(T) == diag(M)
end

end # module TestTridiagonal