Diagonal * AbstractFill matrix special case (#108)

* Diagonal * AbstractFill matrix special case

* Fix Fill * Fill

* Test on 1.5

* Fixes for Julia v1.5

* Adjoint of StaticArray test

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "FillArrays"
 uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "0.9"
+version = "0.9.1"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -12,7 +12,8 @@ julia = "1"
 
 [extras]
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Base64"]
+test = ["Test", "Base64", "StaticArrays"]

src/FillArrays.jl

@@ -9,7 +9,7 @@ import Base: size, getindex, setindex!, IndexStyle, checkbounds, convert,
     show
 
 import LinearAlgebra: rank, svdvals!, tril, triu, tril!, triu!, diag, transpose, adjoint, fill!,
-    norm2, norm1, normInf, normMinusInf, normp, lmul!, rmul!, diagzero, AbstractTriangular
+    norm2, norm1, normInf, normMinusInf, normp, lmul!, rmul!, diagzero, AbstractTriangular, AdjointAbsVec
 
 import Base.Broadcast: broadcasted, DefaultArrayStyle, broadcast_shape
 

src/fillalgebra.jl

@@ -36,25 +36,20 @@ end
 function mult_fill(a::AbstractFill, b::AbstractFill{<:Any,2})
     axes(a, 2) ≠ axes(b, 1) &&
         throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
-    return Fill(getindex_value(a)*getindex_value(b), (axes(a, 1), axes(b, 2)))
+    return Fill(getindex_value(a)*getindex_value(b)*size(a,2), (axes(a, 1), axes(b, 2)))
 end
 
 function mult_fill(a::AbstractFill, b::AbstractFill{<:Any,1})
     axes(a, 2) ≠ axes(b, 1) &&
         throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
-    return Fill(getindex_value(a)*getindex_value(b), (axes(a, 1),))
+    return Fill(getindex_value(a)*getindex_value(b)*size(a,2), (axes(a, 1),))
 end
 
-function mult_ones(a, b::AbstractMatrix)
+function mult_ones(a::AbstractVector, b::AbstractMatrix)
     axes(a, 2) ≠ axes(b, 1) &&
         throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
     return Ones{promote_type(eltype(a), eltype(b))}((axes(a, 1), axes(b, 2)))
 end
-function mult_ones(a, b::AbstractVector)
-    axes(a, 2) ≠ axes(b, 1) &&
-        throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
-    return Ones{promote_type(eltype(a), eltype(b))}((axes(a, 1),))
-end
 
 function mult_zeros(a, b::AbstractMatrix)
     axes(a, 2) ≠ axes(b, 1) &&
@@ -72,8 +67,6 @@ end
 *(a::AbstractFill{<:Any,2}, b::AbstractFill{<:Any,1}) = mult_fill(a,b)
 
 *(a::Ones{<:Any,1}, b::Ones{<:Any,2}) = mult_ones(a, b)
-*(a::Ones{<:Any,2}, b::Ones{<:Any,2}) = mult_ones(a, b)
-*(a::Ones{<:Any,2}, b::Ones{<:Any,1}) = mult_ones(a, b)
 
 *(a::Zeros{<:Any,1}, b::Zeros{<:Any,2}) = mult_zeros(a, b)
 *(a::Zeros{<:Any,2}, b::Zeros{<:Any,2}) = mult_zeros(a, b)
@@ -98,6 +91,14 @@ end
 *(a::Zeros{<:Any,2}, b::Diagonal) = mult_zeros(a, b)
 *(a::Diagonal, b::Zeros{<:Any,1}) = mult_zeros(a, b)
 *(a::Diagonal, b::Zeros{<:Any,2}) = mult_zeros(a, b)
+function *(a::Diagonal, b::AbstractFill{<:Any,2}) 
+    size(a,2) == size(b,1) || throw(DimensionMismatch("A has dimensions $(size(a)) but B has dimensions $(size(b))"))
+    a.diag .* b # use special broadcast
+end
+function *(a::AbstractFill{<:Any,2}, b::Diagonal) 
+    size(a,2) == size(b,1) || throw(DimensionMismatch("A has dimensions $(size(a)) but B has dimensions $(size(b))"))
+    a .* permutedims(b.diag) # use special broadcast
+end
 
 *(a::Adjoint{T, <:StridedMatrix{T}},   b::Fill{T, 1}) where T = reshape(sum(conj.(parent(a)); dims=1) .* b.value, size(parent(a), 2))
 *(a::Transpose{T, <:StridedMatrix{T}}, b::Fill{T, 1}) where T = reshape(sum(parent(a); dims=1) .* b.value, size(parent(a), 2))
@@ -120,13 +121,16 @@ function *(a::StridedMatrix{T}, b::Fill{T, 2}) where T
     fill!(fB, b.value)
     return a*fB
 end
-function *(a::Adjoint{T, <:AbstractVector{T}}, b::Zeros{S, 1}) where {T, S}
+function _adjvec_mul_zeros(a::Adjoint{T}, b::Zeros{S, 1}) where {T, S}
     la, lb = length(a), length(b)
     if la ≠ lb
         throw(DimensionMismatch("dot product arguments have lengths $la and $lb"))
     end
-    return zero(promote_type(T, S))
+    return zero(Base.promote_op(*, T, S))
 end
+
+*(a::AdjointAbsVec, b::Zeros{<:Any, 1}) = _adjvec_mul_zeros(a, b)
+*(a::AdjointAbsVec{<:Number}, b::Zeros{<:Number, 1}) = _adjvec_mul_zeros(a, b)
 *(a::Adjoint{T, <:AbstractMatrix{T}} where T, b::Zeros{<:Any, 1}) = mult_zeros(a, b)
 
 function *(a::Transpose{T, <:AbstractVector{T}}, b::Zeros{T, 1}) where T<:Real

src/fillbroadcast.jl

@@ -100,4 +100,4 @@ end
 broadcasted(::DefaultArrayStyle{N}, op, r::AbstractFill{T,N}, x::Number) where {T,N} = Fill(op(getindex_value(r),x), size(r))
 broadcasted(::DefaultArrayStyle{N}, op, x::Number, r::AbstractFill{T,N}) where {T,N} = Fill(op(x, getindex_value(r)), size(r))
 broadcasted(::DefaultArrayStyle{N}, op, r::AbstractFill{T,N}, x::Ref) where {T,N} = Fill(op(getindex_value(r),x[]), size(r))
-broadcasted(::DefaultArrayStyle{N}, op, x::Ref, r::AbstractFill{T,N}) where {T,N} = Fill(op(x[], getindex_value(r)), size(r))
+broadcasted(::DefaultArrayStyle{N}, op, x::Ref, r::AbstractFill{T,N}) where {T,N} = Fill(op(x[], getindex_value(r)), size(r))

test/runtests.jl

@@ -1,4 +1,4 @@
-using FillArrays, LinearAlgebra, SparseArrays, Random, Base64, Test
+using FillArrays, LinearAlgebra, SparseArrays, StaticArrays, Random, Base64, Test
 import FillArrays: AbstractFill, RectDiagonal, SquareEye
 
 @testset "fill array constructors and convert" begin
@@ -354,6 +354,7 @@ end
     @test randn(4, 3)' * Zeros(4) === Zeros(3)
     @test randn(4)' * Zeros(4) === zero(Float64)
     @test [1, 2, 3]' * Zeros{Int}(3) === zero(Int)
+    @test [SVector(1,2)', SVector(2,3)', SVector(3,4)']' * Zeros{Int}(3) === SVector(0,0)
     @test_throws DimensionMismatch randn(4)' * Zeros(3)
 
     # Check multiplication by Transpose-d vectors works as expected.
@@ -866,22 +867,22 @@ end
 @testset "mult" begin
     @test Fill(2,10)*Fill(3,1,12) == Vector(Fill(2,10))*Matrix(Fill(3,1,12))
     @test Fill(2,10)*Fill(3,1,12) ≡ Fill(6,10,12)
-    @test Fill(2,3,10)*Fill(3,10,12) ≡ Fill(6,3,12)
-    @test Fill(2,3,10)*Fill(3,10) ≡ Fill(6,3)
+    @test Fill(2,3,10)*Fill(3,10,12) ≡ Fill(60,3,12)
+    @test Fill(2,3,10)*Fill(3,10) ≡ Fill(60,3)
     @test_throws DimensionMismatch Fill(2,10)*Fill(3,2,12)
     @test_throws DimensionMismatch Fill(2,3,10)*Fill(3,2,12)
 
     @test Ones(10)*Fill(3,1,12) ≡ Fill(3.0,10,12)
-    @test Ones(10,3)*Fill(3,3,12) ≡ Fill(3.0,10,12)
-    @test Ones(10,3)*Fill(3,3) ≡ Fill(3.0,10)
+    @test Ones(10,3)*Fill(3,3,12) ≡ Fill(9.0,10,12)
+    @test Ones(10,3)*Fill(3,3) ≡ Fill(9.0,10)
 
     @test Fill(2,10)*Ones(1,12) ≡ Fill(2.0,10,12)
-    @test Fill(2,3,10)*Ones(10,12) ≡ Fill(2.0,3,12)
-    @test Fill(2,3,10)*Ones(10) ≡ Fill(2.0,3)
+    @test Fill(2,3,10)*Ones(10,12) ≡ Fill(20.0,3,12)
+    @test Fill(2,3,10)*Ones(10) ≡ Fill(20.0,3)
 
     @test Ones(10)*Ones(1,12) ≡ Ones(10,12)
-    @test Ones(3,10)*Ones(10,12) ≡ Ones(3,12)
-    @test Ones(3,10)*Ones(10) ≡ Ones(3)
+    @test Ones(3,10)*Ones(10,12) ≡ Fill(10.0,3,12)
+    @test Ones(3,10)*Ones(10) ≡ Fill(10.0,3)
 
     @test Zeros(10)*Fill(3,1,12) ≡   Zeros(10,12)
     @test Zeros(10,3)*Fill(3,3,12) ≡ Zeros(10,12)
@@ -924,6 +925,16 @@ end
     @test D*Zeros(1,1) ≡ Zeros(1,1)
     @test D*Zeros(1) ≡ Zeros(1)
 
+    D = Diagonal(Fill(2,10))
+    @test D * Ones(10) ≡ Fill(2.0,10)
+    @test D * Ones(10,5) ≡ Fill(2.0,10,5)
+    @test Ones(5,10) * D ≡ Fill(2.0,5,10)
+
+    # following test is broken in Base as of Julia v1.5
+    @test_skip @test_throws DimensionMismatch Diagonal(Fill(1,1)) * Ones(10)
+    @test_throws DimensionMismatch Diagonal(Fill(1,1)) * Ones(10,5)
+    @test_throws DimensionMismatch Ones(5,10) * Diagonal(Fill(1,1))
+
     E = Eye(5)
     @test E*(1:5) ≡ 1.0:5.0
     @test (1:5)'E == (1.0:5)'