diff --git a/Project.toml b/Project.toml index b5c826ea..80b21214 100644 --- a/Project.toml +++ b/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" diff --git a/src/FillArrays.jl b/src/FillArrays.jl index c5d2c3b1..34eb8d73 100644 --- a/src/FillArrays.jl +++ b/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 diff --git a/src/fillalgebra.jl b/src/fillalgebra.jl index 148db495..3df1fd21 100644 --- a/src/fillalgebra.jl +++ b/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))")) diff --git a/test/runtests.jl b/test/runtests.jl index 953ab1d4..4f415485 100644 --- a/test/runtests.jl +++ b/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