From 4660ee7c422eedd9970ee5af68282184be119830 Mon Sep 17 00:00:00 2001 From: MasonProtter Date: Wed, 22 Apr 2020 12:55:09 -0600 Subject: [PATCH 1/9] special case empty covec-diagonal-vec product --- stdlib/LinearAlgebra/src/diagonal.jl | 18 ++++++++++++++---- stdlib/LinearAlgebra/test/diagonal.jl | 5 +++++ 2 files changed, 19 insertions(+), 4 deletions(-) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index d67c0caab0c0f..40ae83c2596be 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -659,10 +659,20 @@ end # disambiguation methods: * of Diagonal and Adj/Trans AbsVec *(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal) = Adjoint(map((t,s) -> t'*s, D.diag, parent(x))) *(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal) = Transpose(map((t,s) -> transpose(t)*s, D.diag, parent(x))) -*(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) = - mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y)) -*(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) = - mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y)) +function *(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) + if all(isempty.((x, D, y))) + return zero(promote_type(eltype.((x, D, y))...)) + else + return mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y)) + end +end +function *(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) + if all(isempty.((x, D, y))) + return zero(promote_type(eltype.((x, D, y))...)) + else + return mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y)) + end +end function dot(x::AbstractVector, D::Diagonal, y::AbstractVector) mapreduce(t -> dot(t[1], t[2], t[3]), +, zip(x, D.diag, y)) end diff --git a/stdlib/LinearAlgebra/test/diagonal.jl b/stdlib/LinearAlgebra/test/diagonal.jl index fdb070dd70aab..b11e6354ab2bc 100644 --- a/stdlib/LinearAlgebra/test/diagonal.jl +++ b/stdlib/LinearAlgebra/test/diagonal.jl @@ -711,4 +711,9 @@ end @test s1 == prod(sign, d) end +@testset "Empty (#35424)" begin + @test zeros(0)'*Diagonal(zeros(0))*zeros(0) === 0.0 + @test zeros(0)'*Diagonal(zeros(Complex{Int}, 0))*zeros(0) === 0.0 + 0.0im +end + end # module TestDiagonal From 7910d08bfd7736ce65fbdec265d25e664097ddfc Mon Sep 17 00:00:00 2001 From: Mason Protter Date: Wed, 22 Apr 2020 17:50:10 -0600 Subject: [PATCH 2/9] Update stdlib/LinearAlgebra/src/diagonal.jl Co-Authored-By: Takafumi Arakaki --- stdlib/LinearAlgebra/src/diagonal.jl | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index 40ae83c2596be..ab8dee2d71e8f 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -667,11 +667,8 @@ function *(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) end end function *(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) - if all(isempty.((x, D, y))) - return zero(promote_type(eltype.((x, D, y))...)) - else - return mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y)) - end + return mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y); + init = zero(promote_type(map(eltype, (x, D, y))...))) end function dot(x::AbstractVector, D::Diagonal, y::AbstractVector) mapreduce(t -> dot(t[1], t[2], t[3]), +, zip(x, D.diag, y)) From df36bf9f6b4632f8b22139b6b704c66f52ca39e3 Mon Sep 17 00:00:00 2001 From: Mason Protter Date: Wed, 22 Apr 2020 17:50:27 -0600 Subject: [PATCH 3/9] Update stdlib/LinearAlgebra/src/diagonal.jl Co-Authored-By: Takafumi Arakaki --- stdlib/LinearAlgebra/src/diagonal.jl | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index ab8dee2d71e8f..46f3f3eff1659 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -660,11 +660,8 @@ end *(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal) = Adjoint(map((t,s) -> t'*s, D.diag, parent(x))) *(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal) = Transpose(map((t,s) -> transpose(t)*s, D.diag, parent(x))) function *(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) - if all(isempty.((x, D, y))) - return zero(promote_type(eltype.((x, D, y))...)) - else - return mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y)) - end + return mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y); + init = zero(promote_type(map(eltype, (x, D, y))...))) end function *(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) return mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y); From 310abcdde6f2f7677373e2e5f34d72f5edcb09b1 Mon Sep 17 00:00:00 2001 From: MasonProtter Date: Wed, 22 Apr 2020 19:10:54 -0600 Subject: [PATCH 4/9] fix --- stdlib/LinearAlgebra/src/diagonal.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index 46f3f3eff1659..7f3585a42e19a 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -664,11 +664,12 @@ function *(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) init = zero(promote_type(map(eltype, (x, D, y))...))) end function *(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) - return mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y); - init = zero(promote_type(map(eltype, (x, D, y))...))) + mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y), + init=zero(promote_type(eltype(x), eltype(D), eltype(y)))) end function dot(x::AbstractVector, D::Diagonal, y::AbstractVector) - mapreduce(t -> dot(t[1], t[2], t[3]), +, zip(x, D.diag, y)) + mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y), + init=zero(promote_type(eltype(x), eltype(D), eltype(y)))) end function cholesky!(A::Diagonal, ::Val{false} = Val(false); check::Bool = true) From ea09f49684b6bcdb1e4bd4ed427f2bc962a55a07 Mon Sep 17 00:00:00 2001 From: MasonProtter Date: Wed, 22 Apr 2020 22:04:12 -0600 Subject: [PATCH 5/9] revert and don't use broadcast --- stdlib/LinearAlgebra/src/diagonal.jl | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index 7f3585a42e19a..35548023e1f16 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -660,16 +660,21 @@ end *(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal) = Adjoint(map((t,s) -> t'*s, D.diag, parent(x))) *(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal) = Transpose(map((t,s) -> transpose(t)*s, D.diag, parent(x))) function *(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) - return mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y); - init = zero(promote_type(map(eltype, (x, D, y))...))) + if all(isempty(x), isempty(D), isempty(y)) + return zero(promote_type(eltype(x), eltype(D), eltype(y))) + else + return mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y)) + end end function *(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) - mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y), - init=zero(promote_type(eltype(x), eltype(D), eltype(y)))) + if all(isempty(x), isempty(D), isempty(y)) + return zero(promote_type(eltype(x), eltype(D), eltype(y))) + else + return mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y)) + end end function dot(x::AbstractVector, D::Diagonal, y::AbstractVector) - mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y), - init=zero(promote_type(eltype(x), eltype(D), eltype(y)))) + mapreduce(t -> dot(t[1], t[2], t[3]), +, zip(x, D.diag, y)) end function cholesky!(A::Diagonal, ::Val{false} = Val(false); check::Bool = true) From a6c3a609784b5905abe8d3dc4902aa10cdef61c1 Mon Sep 17 00:00:00 2001 From: MasonProtter Date: Wed, 22 Apr 2020 23:13:40 -0600 Subject: [PATCH 6/9] fix all method and reuse function body across 3 methods --- stdlib/LinearAlgebra/src/diagonal.jl | 19 +++++++------------ stdlib/LinearAlgebra/test/diagonal.jl | 3 ++- 2 files changed, 9 insertions(+), 13 deletions(-) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index 35548023e1f16..d91e95ff92a2b 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -659,23 +659,18 @@ end # disambiguation methods: * of Diagonal and Adj/Trans AbsVec *(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal) = Adjoint(map((t,s) -> t'*s, D.diag, parent(x))) *(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal) = Transpose(map((t,s) -> transpose(t)*s, D.diag, parent(x))) -function *(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) - if all(isempty(x), isempty(D), isempty(y)) - return zero(promote_type(eltype(x), eltype(D), eltype(y))) - else - return mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y)) - end -end -function *(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) - if all(isempty(x), isempty(D), isempty(y)) +*(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) = _dot(x, D, y) +*(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) = _dot(x, D, y) +dot(x::AbstractVector, D::Diagonal, y::AbstractVector) = _dot(x, D, y) + +function _dot(x, D::Diagonal, y) + if all((isempty(x), isempty(D), isempty(y))) return zero(promote_type(eltype(x), eltype(D), eltype(y))) else return mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y)) end end -function dot(x::AbstractVector, D::Diagonal, y::AbstractVector) - mapreduce(t -> dot(t[1], t[2], t[3]), +, zip(x, D.diag, y)) -end + function cholesky!(A::Diagonal, ::Val{false} = Val(false); check::Bool = true) info = 0 diff --git a/stdlib/LinearAlgebra/test/diagonal.jl b/stdlib/LinearAlgebra/test/diagonal.jl index b11e6354ab2bc..532b556847835 100644 --- a/stdlib/LinearAlgebra/test/diagonal.jl +++ b/stdlib/LinearAlgebra/test/diagonal.jl @@ -713,7 +713,8 @@ end @testset "Empty (#35424)" begin @test zeros(0)'*Diagonal(zeros(0))*zeros(0) === 0.0 - @test zeros(0)'*Diagonal(zeros(Complex{Int}, 0))*zeros(0) === 0.0 + 0.0im + @test transpose(zeros(0))*Diagonal(zeros(Complex{Int}, 0))*zeros(0) === 0.0 + 0.0im + @test dot(zeros(Int32, 0), Digonal(zeros(Int, 0)), zeros(Int16, 0)) === 0 end end # module TestDiagonal From bb51ff68ea5d352c882445be2bb88d4f54493100 Mon Sep 17 00:00:00 2001 From: MasonProtter Date: Thu, 23 Apr 2020 00:23:42 -0600 Subject: [PATCH 7/9] fix typo --- stdlib/LinearAlgebra/test/diagonal.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/test/diagonal.jl b/stdlib/LinearAlgebra/test/diagonal.jl index 532b556847835..98b2fca354fd6 100644 --- a/stdlib/LinearAlgebra/test/diagonal.jl +++ b/stdlib/LinearAlgebra/test/diagonal.jl @@ -714,7 +714,7 @@ end @testset "Empty (#35424)" begin @test zeros(0)'*Diagonal(zeros(0))*zeros(0) === 0.0 @test transpose(zeros(0))*Diagonal(zeros(Complex{Int}, 0))*zeros(0) === 0.0 + 0.0im - @test dot(zeros(Int32, 0), Digonal(zeros(Int, 0)), zeros(Int16, 0)) === 0 + @test dot(zeros(Int32, 0), Diagonal(zeros(Int, 0)), zeros(Int16, 0)) === 0 end end # module TestDiagonal From c9074228d52fc0caab23a3cd60580502883f1907 Mon Sep 17 00:00:00 2001 From: MasonProtter Date: Thu, 23 Apr 2020 08:45:10 -0600 Subject: [PATCH 8/9] fix change of behaviour in dot --- stdlib/LinearAlgebra/src/diagonal.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index d91e95ff92a2b..ce9801d367563 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -659,15 +659,15 @@ end # disambiguation methods: * of Diagonal and Adj/Trans AbsVec *(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal) = Adjoint(map((t,s) -> t'*s, D.diag, parent(x))) *(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal) = Transpose(map((t,s) -> transpose(t)*s, D.diag, parent(x))) -*(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) = _dot(x, D, y) -*(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) = _dot(x, D, y) -dot(x::AbstractVector, D::Diagonal, y::AbstractVector) = _dot(x, D, y) +*(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) = _mapreduce_prod(*, x, D, y) +*(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) = _mapreduce_prod(*, x, D, y) +dot(x::AbstractVector, D::Diagonal, y::AbstractVector) = _mapreduce_prod(dot, x, D, y) -function _dot(x, D::Diagonal, y) - if all((isempty(x), isempty(D), isempty(y))) +function _mapreduce_prod(f, x, D::Diagonal, y) + if isempty(x) && isempty(D) && isempty(y) return zero(promote_type(eltype(x), eltype(D), eltype(y))) else - return mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y)) + return mapreduce(t -> f(t[1], t[2], t[3]), +, zip(x, D.diag, y)) end end From b53e804d2a82efe7c31a60beffdd4f71166e21bb Mon Sep 17 00:00:00 2001 From: Mason Protter Date: Thu, 23 Apr 2020 11:32:11 -0600 Subject: [PATCH 9/9] Update stdlib/LinearAlgebra/src/diagonal.jl Co-Authored-By: Daniel Karrasch --- stdlib/LinearAlgebra/src/diagonal.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index ce9801d367563..f3b4ac17eec78 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -665,7 +665,7 @@ dot(x::AbstractVector, D::Diagonal, y::AbstractVector) = _mapreduce_prod(dot, x, function _mapreduce_prod(f, x, D::Diagonal, y) if isempty(x) && isempty(D) && isempty(y) - return zero(promote_type(eltype(x), eltype(D), eltype(y))) + return zero(Base.promote_op(f, eltype(x), eltype(D), eltype(y))) else return mapreduce(t -> f(t[1], t[2], t[3]), +, zip(x, D.diag, y)) end