Special case for Zeros(5)'*b and dot(::AbstractFill,::AbstractFill) (#160)

* Special case for Zeros(5)'*b and dot(::AbstractFill,::AbstractFill)

* more overloads

* add tests

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "FillArrays"
 uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "0.12.5"
+version = "0.12.6"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/FillArrays.jl

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

src/fillalgebra.jl

@@ -83,6 +83,8 @@ end
 *(a::Zeros{<:Any,2}, b::AbstractVector) = mult_zeros(a, b)
 *(a::AbstractVector, b::Zeros{<:Any,2}) = mult_zeros(a, b)
 
+*(a::Zeros{<:Any,1}, b::AdjOrTransAbsVec) = mult_zeros(a, b)
+
 *(a::Zeros{<:Any,1}, b::Diagonal) = mult_zeros(a, b)
 *(a::Zeros{<:Any,2}, b::Diagonal) = mult_zeros(a, b)
 *(a::Diagonal, b::Zeros{<:Any,1}) = mult_zeros(a, b)
@@ -117,16 +119,29 @@ function *(a::StridedMatrix{T}, b::Fill{T, 2}) where T
     fill!(fB, b.value)
     return a*fB
 end
-function _adjvec_mul_zeros(a::Adjoint{T}, b::Zeros{S, 1}) where {T, S}
+function _adjvec_mul_zeros(a, b)
     la, lb = length(a), length(b)
     if la ≠ lb
         throw(DimensionMismatch("dot product arguments have lengths $la and $lb"))
     end
-    return zero(Base.promote_op(*, T, S))
+    return zero(Base.promote_op(*, eltype(a), eltype(b)))
 end
 
+*(a::AdjointAbsVec{<:Any,<:Zeros{<:Any,1}}, b::AbstractMatrix) = (b' * a')'
+*(a::AdjointAbsVec{<:Any,<:Zeros{<:Any,1}}, b::Zeros{<:Any,2}) = (b' * a')'
+*(a::TransposeAbsVec{<:Any,<:Zeros{<:Any,1}}, b::AbstractMatrix) = transpose(transpose(b) * transpose(a))
+*(a::TransposeAbsVec{<:Any,<:Zeros{<:Any,1}}, b::Zeros{<:Any,2}) = transpose(transpose(b) * transpose(a))
+
+*(a::AbstractVector, b::AdjOrTransAbsVec{<:Any,<:Zeros{<:Any,1}}) = a * permutedims(parent(b))
+*(a::AbstractMatrix, b::AdjOrTransAbsVec{<:Any,<:Zeros{<:Any,1}}) = a * permutedims(parent(b))
+*(a::Zeros{<:Any,1}, b::AdjOrTransAbsVec{<:Any,<:Zeros{<:Any,1}}) = a * permutedims(parent(b))
+*(a::Zeros{<:Any,2}, b::AdjOrTransAbsVec{<:Any,<:Zeros{<:Any,1}}) = a * permutedims(parent(b))
+
 *(a::AdjointAbsVec, b::Zeros{<:Any, 1}) = _adjvec_mul_zeros(a, b)
 *(a::AdjointAbsVec{<:Number}, b::Zeros{<:Number, 1}) = _adjvec_mul_zeros(a, b)
+*(a::TransposeAbsVec, b::Zeros{<:Any, 1}) = _adjvec_mul_zeros(a, b)
+*(a::TransposeAbsVec{<: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
@@ -138,6 +153,39 @@ function *(a::Transpose{T, <:AbstractVector{T}}, b::Zeros{T, 1}) where T<:Real
 end
 *(a::Transpose{T, <:AbstractMatrix{T}}, b::Zeros{T, 1}) where T<:Real = mult_zeros(a, b)
 
+# treat zero separately to support ∞-vectors
+function _zero_dot(a, b)
+    axes(a) == axes(b) || throw(DimensionMismatch("dot product arguments have lengths $(length(a)) and $(length(b))"))
+    zero(promote_type(eltype(a),eltype(b)))
+end
+
+_fill_dot(a::Zeros, b::Zeros) = _zero_dot(a, b)
+_fill_dot(a::Zeros, b) = _zero_dot(a, b)
+_fill_dot(a, b::Zeros) = _zero_dot(a, b)
+_fill_dot(a::Zeros, b::AbstractFill) = _zero_dot(a, b)
+_fill_dot(a::AbstractFill, b::Zeros) = _zero_dot(a, b)
+
+function _fill_dot(a::AbstractFill, b::AbstractFill)
+    axes(a) == axes(b) || throw(DimensionMismatch("dot product arguments have lengths $(length(a)) and $(length(b))"))
+    getindex_value(a)getindex_value(b)*length(b)
+end
+
+# support types with fast sum
+function _fill_dot(a::AbstractFill, b)
+    axes(a) == axes(b) || throw(DimensionMismatch("dot product arguments have lengths $(length(a)) and $(length(b))"))
+    getindex_value(a)sum(b)
+end
+
+function _fill_dot(a, b::AbstractFill)
+    axes(a) == axes(b) || throw(DimensionMismatch("dot product arguments have lengths $(length(a)) and $(length(b))"))
+    sum(a)getindex_value(b)
+end
+
+
+dot(a::AbstractFill{<:Any,1}, b::AbstractFill{<:Any,1}) = _fill_dot(a, b)
+dot(a::AbstractFill{<:Any,1}, b::AbstractVector) = _fill_dot(a, b)
+dot(a::AbstractVector, b::AbstractFill{<:Any,1}) = _fill_dot(a, b)
+
 function dot(u::AbstractVector, E::Eye, v::AbstractVector)
     length(u) == size(E,1) && length(v) == size(E,2) ||
         throw(DimensionMismatch("dot product arguments have dimensions $(length(u))×$(size(E))×$(length(v))"))

test/runtests.jl

@@ -451,29 +451,49 @@ end
     @test [1,2,3]*Zeros(1,3) ≡ Zeros(3,3)
     @test_throws MethodError [1,2,3]*Zeros(3) # Not defined for [1,2,3]*[0,0,0] either
 
-    # Check multiplication by Adjoint vectors works as expected.
-    @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.
-    @test transpose(randn(4, 3)) * Zeros(4) === Zeros(3)
-    @test transpose(randn(4)) * Zeros(4) === zero(Float64)
-    @test transpose([1, 2, 3]) * Zeros{Int}(3) === zero(Int)
-    @test_throws DimensionMismatch transpose(randn(4)) * Zeros(3)
-
-    @test +(Zeros{Float64}(3, 5)) === Zeros{Float64}(3, 5)
-    @test -(Zeros{Float32}(5, 2)) === Zeros{Float32}(5, 2)
-
-    # `Zeros` are closed under addition and subtraction (both unary and binary).
+    @testset "Check multiplication by Adjoint vectors works as expected." begin
+        @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)
+        @test Zeros(5)' * randn(5,3) ≡ Zeros(5)'*Zeros(5,3) ≡ Zeros(5)'*Ones(5,3) ≡ Zeros(3)'
+        @test Zeros(5)' * randn(5) ≡ Zeros(5)' * Zeros(5) ≡ Zeros(5)' * Ones(5) ≡ 0.0 
+        @test Zeros(5) * Zeros(6)' ≡ Zeros(5,1) * Zeros(6)' ≡ Zeros(5,6)
+        @test randn(5) * Zeros(6)' ≡ randn(5,1) * Zeros(6)' ≡ Zeros(5,6)
+        @test Zeros(5) * randn(6)' ≡ Zeros(5,6)
+
+        @test ([[1,2]])' * Zeros{SVector{2,Int}}(1) ≡ 0
+        @test_broken ([[1,2,3]])' * Zeros{SVector{2,Int}}(1)
+    end
+
+    @testset "Check multiplication by Transpose-d vectors works as expected." begin
+        @test transpose(randn(4, 3)) * Zeros(4) === Zeros(3)
+        @test transpose(randn(4)) * Zeros(4) === zero(Float64)
+        @test transpose([1, 2, 3]) * Zeros{Int}(3) === zero(Int)
+        @test_throws DimensionMismatch transpose(randn(4)) * Zeros(3)
+        @test transpose(Zeros(5)) * randn(5,3) ≡ transpose(Zeros(5))*Zeros(5,3) ≡ transpose(Zeros(5))*Ones(5,3) ≡ transpose(Zeros(3))
+        @test transpose(Zeros(5)) * randn(5) ≡ transpose(Zeros(5)) * Zeros(5) ≡ transpose(Zeros(5)) * Ones(5) ≡ 0.0 
+        @test randn(5) * transpose(Zeros(6)) ≡ randn(5,1) * transpose(Zeros(6)) ≡ Zeros(5,6)
+        @test Zeros(5) * transpose(randn(6)) ≡ Zeros(5,6)
+        @test transpose(randn(5)) * Zeros(5) ≡ 0.0
+
+        @test transpose([[1,2]]) * Zeros{SVector{2,Int}}(1) ≡ 0
+        @test_broken transpose([[1,2,3]]) * Zeros{SVector{2,Int}}(1)
+    end
+
     z1, z2 = Zeros{Float64}(4), Zeros{Int}(4)
-    @test +(z1) === z1
-    @test -(z1) === z1
 
-    test_addition_and_subtraction([z1, z2], [z1, z2], Zeros)
-    test_addition_and_subtraction_dim_mismatch(z1, Zeros{Float64}(4, 2))
+    @testset "`Zeros` are closed under addition and subtraction (both unary and binary)." begin
+        @test +(Zeros{Float64}(3, 5)) === Zeros{Float64}(3, 5)
+        @test -(Zeros{Float32}(5, 2)) === Zeros{Float32}(5, 2)
+
+        @test +(z1) === z1
+        @test -(z1) === z1
+
+        test_addition_and_subtraction([z1, z2], [z1, z2], Zeros)
+        test_addition_and_subtraction_dim_mismatch(z1, Zeros{Float64}(4, 2))
+    end
 
     # `Zeros` +/- `Fill`s should yield `Fills`.
     fill1, fill2 = Fill(5.0, 4), Fill(5, 4)
@@ -502,36 +522,41 @@ end
         @test op(Zeros{Float64}(4, 5), Zeros{Int}(4, 5)) === Zeros{Float64}(4, 5)
     end
 
-    # Zeros +/- dense where + / - have different results.
-    @test +(Zeros(3, 5), X) == X && +(X, Zeros(3, 5)) == X
-    @test !(Zeros(3, 5) + X === X) && !(X + Zeros(3, 5) === X)
-    @test -(Zeros(3, 5), X) == -X
-
-    # Addition with different eltypes.
-    @test +(Zeros{Float32}(3, 5), X) isa Matrix{Float64}
-    @test !(+(Zeros{Float32}(3, 5), X) === X)
-    @test +(Zeros{Float32}(3, 5), X) == X
-    @test !(+(Zeros{ComplexF64}(3, 5), X) === X)
-    @test +(Zeros{ComplexF64}(3, 5), X) == X
-
-    # Subtraction with different eltypes.
-    @test -(Zeros{Float32}(3, 5), X) isa Matrix{Float64}
-    @test -(Zeros{Float32}(3, 5), X) == -X
-    @test -(Zeros{ComplexF64}(3, 5), X) == -X
-
-    # Tests for ranges.
-    X = randn(5)
-    @test !(Zeros(5) + X === X)
-    @test Zeros{Int}(5) + (1:5) === (1:5) && (1:5) + Zeros{Int}(5) === (1:5)
-    @test Zeros(5) + (1:5) === (1.0:1.0:5.0) && (1:5) + Zeros(5) === (1.0:1.0:5.0)
-    @test (1:5) - Zeros{Int}(5) === (1:5)
-    @test Zeros{Int}(5) - (1:5) === -1:-1:-5
-    @test Zeros(5) - (1:5) === -1.0:-1.0:-5.0
-
-    # test Base.zero
-    @test zero(Zeros(10)) == Zeros(10)
-    @test zero(Ones(10,10)) == Zeros(10,10)
-    @test zero(Fill(0.5, 10, 10)) == Zeros(10,10)
+    @testset "Zeros +/- dense where + / - have different results." begin
+        @test +(Zeros(3, 5), X) == X && +(X, Zeros(3, 5)) == X
+        @test !(Zeros(3, 5) + X === X) && !(X + Zeros(3, 5) === X)
+        @test -(Zeros(3, 5), X) == -X
+    end
+
+    @testset "Addition with different eltypes." begin
+        @test +(Zeros{Float32}(3, 5), X) isa Matrix{Float64}
+        @test !(+(Zeros{Float32}(3, 5), X) === X)
+        @test +(Zeros{Float32}(3, 5), X) == X
+        @test !(+(Zeros{ComplexF64}(3, 5), X) === X)
+        @test +(Zeros{ComplexF64}(3, 5), X) == X
+    end
+
+    @testset "Subtraction with different eltypes." begin
+        @test -(Zeros{Float32}(3, 5), X) isa Matrix{Float64}
+        @test -(Zeros{Float32}(3, 5), X) == -X
+        @test -(Zeros{ComplexF64}(3, 5), X) == -X
+    end
+
+    @testset "Tests for ranges." begin
+        X = randn(5)
+        @test !(Zeros(5) + X === X)
+        @test Zeros{Int}(5) + (1:5) === (1:5) && (1:5) + Zeros{Int}(5) === (1:5)
+        @test Zeros(5) + (1:5) === (1.0:1.0:5.0) && (1:5) + Zeros(5) === (1.0:1.0:5.0)
+        @test (1:5) - Zeros{Int}(5) === (1:5)
+        @test Zeros{Int}(5) - (1:5) === -1:-1:-5
+        @test Zeros(5) - (1:5) === -1.0:-1.0:-5.0
+    end
+
+    @testset "test Base.zero" begin
+        @test zero(Zeros(10)) == Zeros(10)
+        @test zero(Ones(10,10)) == Zeros(10,10)
+        @test zero(Fill(0.5, 10, 10)) == Zeros(10,10)
+    end
 end
 
 @testset "maximum/minimum/svd/sort" begin
@@ -1135,25 +1160,7 @@ end
     @test E*(1:5) ≡ 1.0:5.0
     @test (1:5)'E == (1.0:5)'
     @test E*E ≡ E
-end
-
-@testset "count" begin
-    @test count(Ones{Bool}(10)) == count(Fill(true,10)) == 10
-    @test count(Zeros{Bool}(10)) == count(Fill(false,10)) == 0
-    @test count(x -> 1 ≤ x < 2, Fill(1.3,10)) == 10
-    @test count(x -> 1 ≤ x < 2, Fill(2.0,10)) == 0
-end
-
-@testset "norm" begin
-    for a in (Zeros{Int}(5), Zeros(5,3), Zeros(2,3,3),
-                Ones{Int}(5), Ones(5,3), Ones(2,3,3),
-                Fill(2.3,5), Fill([2.3,4.2],5), Fill(4)),
-        p in (-Inf, 0, 0.1, 1, 2, 3, Inf)
-        @test norm(a,p) ≈ norm(Array(a),p)
-    end
-end
 
-@testset "multiplication" begin
     for T in (Float64, ComplexF64)
         fv = T == Float64 ? Float64(1.6) : ComplexF64(1.6, 1.3)
         n  = 10
@@ -1172,6 +1179,22 @@ end
     end
 end
 
+@testset "count" begin
+    @test count(Ones{Bool}(10)) == count(Fill(true,10)) == 10
+    @test count(Zeros{Bool}(10)) == count(Fill(false,10)) == 0
+    @test count(x -> 1 ≤ x < 2, Fill(1.3,10)) == 10
+    @test count(x -> 1 ≤ x < 2, Fill(2.0,10)) == 0
+end
+
+@testset "norm" begin
+    for a in (Zeros{Int}(5), Zeros(5,3), Zeros(2,3,3),
+                Ones{Int}(5), Ones(5,3), Ones(2,3,3),
+                Fill(2.3,5), Fill([2.3,4.2],5), Fill(4)),
+        p in (-Inf, 0, 0.1, 1, 2, 3, Inf)
+        @test norm(a,p) ≈ norm(Array(a),p)
+    end
+end
+
 @testset "dot products" begin
     n = 15
     o = Ones(1:n)
@@ -1187,6 +1210,17 @@ end
     @test dot(u, 2D, v) == 2dot(u, v)
     @test dot(u, Z, v) == 0
 
+    @test dot(Zeros(5), Zeros{ComplexF16}(5)) ≡ zero(ComplexF64)
+    @test dot(Zeros(5), Ones{ComplexF16}(5)) ≡ zero(ComplexF64)
+    @test dot(Ones{ComplexF16}(5), Zeros(5)) ≡ zero(ComplexF64)
+    @test dot(randn(5), Zeros{ComplexF16}(5)) ≡ dot(Zeros{ComplexF16}(5), randn(5)) ≡ zero(ComplexF64)
+
+    @test dot(Fill(1,5), Fill(2.0,5)) ≡ 10.0
+
+    let N = 2^big(1000) # fast dot for fast sum
+        @test dot(Fill(2,N),1:N) == dot(Fill(2,N),1:N) == dot(1:N,Fill(2,N)) == 2*sum(1:N)
+    end
+
     @test_throws DimensionMismatch dot(u[1:end-1], D, v)
     @test_throws DimensionMismatch dot(u[1:end-1], D, v[1:end-1])
 
@@ -1195,6 +1229,9 @@ end
 
     @test_throws DimensionMismatch dot(u, Z, v[1:end-1])
     @test_throws DimensionMismatch dot(u, Z, v[1:end-1])
+
+    @test_throws DimensionMismatch dot(Zeros(5), Zeros(6))
+    @test_throws DimensionMismatch dot(Zeros(5), randn(6))
 end
 
 @testset "print" begin