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.

Sign up for GitHub

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

Jump to bottom

Fix multiplication of Triangular and Diagonal matrix #33334

Merged
merged 5 commits into from
Oct 10, 2019
Merged

Fix multiplication of Triangular and Diagonal matrix #33334

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
5d395c2
Fix multiplication of Triangular and Diagonal matrix
goggle Sep 20, 2019
bc367b8
Adopt more multiplication methods
goggle Sep 20, 2019
2032b6f
Add `similar` methods for `Adjoint` and `Transpose`
goggle Sep 20, 2019
cfc0a23
Enable type tests
goggle Sep 20, 2019
074607e
Write out `similar` methods explicitly
goggle Sep 20, 2019
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
Prev Previous commit
Next Next commit
Adopt more multiplication methods
  • Loading branch information
@goggle
goggle committed Sep 20, 2019
commit bc367b813913b0c7bbe18a4fbc1619d209a7ab1b
12 changes: 8 additions & 4 deletions stdlib/LinearAlgebra/src/diagonal.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -216,7 +216,8 @@ function lmul!(D::Diagonal, B::UnitUpperTriangular)
end

*(D::Adjoint{<:Any,<:Diagonal}, B::Diagonal) = Diagonal(adjoint.(D.parent.diag) .* B.diag)
*(A::Adjoint{<:Any,<:AbstractTriangular}, D::Diagonal) = rmul!(copy(A), D)
*(A::Adjoint{<:Any,<:AbstractTriangular}, D::Diagonal) =
rmul!(copyto!(similar(A, promote_op(*, eltype(A), eltype(D.diag))), A), D)
function *(adjA::Adjoint{<:Any,<:AbstractMatrix}, D::Diagonal)
A = adjA.parent
Ac = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
Expand All @@ -225,7 +226,8 @@ function *(adjA::Adjoint{<:Any,<:AbstractMatrix}, D::Diagonal)
end

*(D::Transpose{<:Any,<:Diagonal}, B::Diagonal) = Diagonal(transpose.(D.parent.diag) .* B.diag)
*(A::Transpose{<:Any,<:AbstractTriangular}, D::Diagonal) = rmul!(copy(A), D)
*(A::Transpose{<:Any,<:AbstractTriangular}, D::Diagonal) =
rmul!(copyto!(similar(A, promote_op(*, eltype(A), eltype(D.diag))), A), D)
function *(transA::Transpose{<:Any,<:AbstractMatrix}, D::Diagonal)
A = transA.parent
At = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
Expand All @@ -234,7 +236,8 @@ function *(transA::Transpose{<:Any,<:AbstractMatrix}, D::Diagonal)
end

*(D::Diagonal, B::Adjoint{<:Any,<:Diagonal}) = Diagonal(D.diag .* adjoint.(B.parent.diag))
*(D::Diagonal, B::Adjoint{<:Any,<:AbstractTriangular}) = lmul!(D, collect(B))
*(D::Diagonal, B::Adjoint{<:Any,<:AbstractTriangular}) =
lmul!(D, copyto!(similar(B, promote_op(*, eltype(B), eltype(D.diag))), B))
*(D::Diagonal, adjQ::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) = (Q = adjQ.parent; rmul!(Array(D), adjoint(Q)))
function *(D::Diagonal, adjA::Adjoint{<:Any,<:AbstractMatrix})
A = adjA.parent
Expand All @@ -244,7 +247,8 @@ function *(D::Diagonal, adjA::Adjoint{<:Any,<:AbstractMatrix})
end

*(D::Diagonal, B::Transpose{<:Any,<:Diagonal}) = Diagonal(D.diag .* transpose.(B.parent.diag))
*(D::Diagonal, B::Transpose{<:Any,<:AbstractTriangular}) = lmul!(D, copy(B))
*(D::Diagonal, B::Transpose{<:Any,<:AbstractTriangular}) =
lmul!(D, copyto!(similar(B, promote_op(*, eltype(B), eltype(D.diag))), B))
function *(D::Diagonal, transA::Transpose{<:Any,<:AbstractMatrix})
A = transA.parent
At = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
Expand Down
33 changes: 23 additions & 10 deletions stdlib/LinearAlgebra/test/diagonal.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -444,30 +444,23 @@ end
@test Matrix(1.0I, 5, 5) \ Diagonal(fill(1.,5)) == Matrix(I, 5, 5)

@testset "Triangular and Diagonal" begin
for T in (LowerTriangular(randn(5,5)), LinearAlgebra.UnitLowerTriangular(randn(5,5)))
D = Diagonal(randn(5))
@test T*D == Array(T)*Array(D)
@test T'D == Array(T)'*Array(D)
@test transpose(T)*D == transpose(Array(T))*Array(D)
@test D*T' == Array(D)*Array(T)'
@test D*transpose(T) == Array(D)*transpose(Array(T))
@test D*T == Array(D)*Array(T)
end
function _test_matrix(type)
if type == Int
return rand(1:9, 5, 5)
else
return randn(type, 5, 5)
end
end
types = (Float32, Float64, Int, ComplexF64)
types = (Float64, Int, ComplexF64)
for ta in types
D = Diagonal(_test_matrix(ta))
for tb in types
B = _test_matrix(tb)
Tmats = (LowerTriangular(B), UnitLowerTriangular(B), UpperTriangular(B), UnitUpperTriangular(B))
restypes = (LowerTriangular, LowerTriangular, UpperTriangular, UpperTriangular)
for (T, rtype) in zip(Tmats, restypes)
adjtype = (rtype == LowerTriangular) ? UpperTriangular : LowerTriangular

# Triangular * Diagonal
R = T * D
@test R == Array(T) * Array(D)
Expand All @@ -477,6 +470,26 @@ end
R = D * T
@test R == Array(D) * Array(T)
@test isa(R, rtype)

# Adjoint of Triangular * Diagonal
R = T' * D
@test R == Array(T)' * Array(D)
# @test isa(R, adjtype)

# Diagonal * Adjoint of Triangular
R = D * T'
@test R == Array(D) * Array(T)'
# @test isa(R, adjtype)

# Transpose of Triangular * Diagonal
R = transpose(T) * D
@test R == transpose(Array(T)) * Array(D)
# @test isa(R, adjtype)

# Diagonal * Transpose of Triangular
R = D * transpose(T)
@test R == Array(D) * transpose(Array(T))
# @test isa(R, adjtype)
end
end
end
Expand Down