Merge branch 'master' into permutedims

Author: jishnub

Project.toml

@@ -22,9 +22,11 @@ FillArraysStatisticsExt = "Statistics"
 [compat]
 Aqua = "0.8"
 Base64 = "1.6"
+Documenter = "1"
 Infinities = "0.1"
 LinearAlgebra = "1.6"
 PDMats = "0.11.17"
+Quaternions = "0.7"
 Random = "1.6"
 ReverseDiff = "1"
 SparseArrays = "1.6"
@@ -36,13 +38,15 @@ julia = "1.6"
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Infinities = "e1ba4f0e-776d-440f-acd9-e1d2e9742647"
 PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
+Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Test", "Base64", "Infinities", "PDMats", "ReverseDiff", "SparseArrays", "StaticArrays", "Statistics"]
+test = ["Aqua", "Test", "Base64", "Infinities", "PDMats", "ReverseDiff", "SparseArrays", "StaticArrays", "Statistics", "Quaternions", "Documenter"]

docs/Project.toml

@@ -6,5 +6,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
-Documenter = "0.27"
+Documenter = "1"
+Random = "1"
+SparseArrays = "1"
 StaticArrays = "1"

docs/src/index.md

@@ -69,7 +69,7 @@ The elements of a `Fill` array don't need to be restricted to numbers, and these
 
 ```jldoctest
 julia> f = Fill("hello", 2,5)
-2×5 Fill{String}, with entries equal to hello
+2×5 Fill{String}, with entries equal to "hello"
 
 julia> Array(f)
 2×5 Matrix{String}:

src/FillArrays.jl

@@ -4,14 +4,14 @@ module FillArrays
 using LinearAlgebra
 import Base: size, getindex, setindex!, IndexStyle, checkbounds, convert,
     +, -, *, /, \, diff, sum, cumsum, maximum, minimum, sort, sort!,
-    any, all, axes, isone, iterate, unique, allunique, permutedims, inv,
+    any, all, axes, isone, iszero, iterate, unique, allunique, permutedims, inv,
     copy, vec, setindex!, count, ==, reshape, map, zero,
     show, view, in, mapreduce, one, reverse, promote_op, promote_rule, repeat,
     parent, similar, issorted
 
 import LinearAlgebra: rank, svdvals!, tril, triu, tril!, triu!, diag, transpose, adjoint, fill!,
     dot, norm2, norm1, normInf, normMinusInf, normp, lmul!, rmul!, diagzero, AdjointAbsVec, TransposeAbsVec,
-    issymmetric, ishermitian, AdjOrTransAbsVec, checksquare, mul!, kron
+    issymmetric, ishermitian, AdjOrTransAbsVec, checksquare, mul!, kron, AbstractTriangular
 
 
 import Base.Broadcast: broadcasted, DefaultArrayStyle, broadcast_shape, BroadcastStyle, Broadcasted
@@ -94,14 +94,10 @@ Typically created by `Fill` or `Zeros` or `Ones`
 
 ```jldoctest
 julia> Fill(7, (2,3))
-2×3 Fill{Int64,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}}:
- 7  7  7
- 7  7  7
-
-julia> Fill{Float64, 1, Tuple{UnitRange{Int64}}}(7., (1:2,))
-2-element Fill{Float64,1,Tuple{UnitRange{Int64}}} with indices 1:2:
- 7.0
- 7.0
+2×3 Fill{Int64}, with entries equal to 7
+
+julia> Fill{Float64, 1, Tuple{UnitRange{Int64}}}(7.0, (1:2,))
+2-element Fill{Float64, 1, Tuple{UnitRange{Int64}}} with indices 1:2, with entries equal to 7.0
 ```
 """
 struct Fill{T, N, Axes} <: AbstractFill{T, N, Axes}
@@ -621,9 +617,10 @@ end
 #########
 
 function isone(AF::AbstractFillMatrix)
-    isone(getindex_value(AF)) || return false
     (n,m) = size(AF)
     n != m && return false
+    (n == 0 || m == 0) && return true
+    isone(getindex_value(AF)) || return false
     n == 1 && return true
     return false
 end
@@ -715,12 +712,9 @@ function Base.show(io::IO, ::MIME"text/plain", x::Union{Eye, AbstractFill})
         return show(io, x)
     end
     summary(io, x)
-    if x isa Union{AbstractZeros, AbstractOnes, Eye}
-        # then no need to print entries
-    elseif length(x) > 1
-        print(io, ", with entries equal to ", getindex_value(x))
-    else
-        print(io, ", with entry equal to ", getindex_value(x))
+    if !(x isa Union{AbstractZeros, AbstractOnes, Eye})
+        print(io, ", with ", length(x) > 1 ? "entries" : "entry", " equal to ")
+        show(io, getindex_value(x))
     end
 end
 
@@ -742,7 +736,13 @@ function Base.show(io::IO, x::AbstractFill)  # for example (Fill(π,3),)
     join(io, size(x), ", ")
     print(io, ")")
 end
-Base.show(io::IO, x::Eye) = print(io, "Eye(", size(x,1), ")")
+function Base.show(io::IO, x::Eye)
+    print(io, "Eye(", size(x,1))
+    if size(x,1) != size(x,2)
+        print(io, ",", size(x,2))
+    end
+    print(io, ")")
+end
 
 Base.array_summary(io::IO, ::Zeros{T}, inds::Tuple{Vararg{Base.OneTo}}) where T =
     print(io, Base.dims2string(length.(inds)), " Zeros{$T}")

src/fillalgebra.jl

@@ -6,10 +6,17 @@ vec(a::AbstractFill) = fillsimilar(a, length(a))
 # cannot do this for vectors since that would destroy scalar dot product
 
 
-transpose(a::Union{AbstractOnesMatrix, AbstractZerosMatrix}) = fillsimilar(a, reverse(axes(a)))
-adjoint(a::Union{AbstractOnesMatrix, AbstractZerosMatrix}) = fillsimilar(a, reverse(axes(a)))
-transpose(a::FillMatrix{T}) where T = Fill{T}(transpose(a.value), reverse(a.axes))
-adjoint(a::FillMatrix{T}) where T = Fill{T}(adjoint(a.value), reverse(a.axes))
+for OP in (:transpose, :adjoint)
+    @eval begin
+        function $OP(a::AbstractZerosMatrix)
+            v = getindex_value(a)
+            T = typeof($OP(v))
+            Zeros{T}(reverse(axes(a)))
+        end
+        $OP(a::AbstractOnesMatrix) = fillsimilar(a, reverse(axes(a)))
+        $OP(a::FillMatrix) = Fill($OP(a.value), reverse(a.axes))
+    end
+end
 
 permutedims(a::AbstractFillMatrix) = fillsimilar(a, reverse(a.axes))
 
@@ -91,9 +98,18 @@ mult_ones(a, b) = mult_ones(a, b, mult_axes(a, b))
 *(a::AbstractFillMatrix, b::AbstractZerosMatrix) = mult_zeros(a, b)
 *(a::AbstractFillMatrix, b::AbstractZerosVector) = mult_zeros(a, b)
 
-*(a::AbstractZerosMatrix, b::AbstractMatrix) = mult_zeros(a, b)
-*(a::AbstractMatrix, b::AbstractZerosVector) = mult_zeros(a, b)
-*(a::AbstractMatrix, b::AbstractZerosMatrix) = mult_zeros(a, b)
+for MT in (:AbstractMatrix, :AbstractTriangular)
+    @eval *(a::AbstractZerosMatrix, b::$MT) = mult_zeros(a, b)
+    @eval *(a::$MT, b::AbstractZerosMatrix) = mult_zeros(a, b)
+end
+# Odd way to deal with the type-parameters to avoid ambiguities
+for MT in (:(AbstractMatrix{T}), :(Transpose{<:Any, <:AbstractMatrix{T}}), :(Adjoint{<:Any, <:AbstractMatrix{T}}),
+            :(AbstractTriangular{T}))
+    @eval *(a::$MT, b::AbstractZerosVector) where {T} = mult_zeros(a, b)
+end
+for MT in (:(Transpose{<:Any, <:AbstractVector}), :(Adjoint{<:Any, <:AbstractVector}))
+    @eval *(a::$MT, b::AbstractZerosMatrix) = mult_zeros(a, b)
+end
 *(a::AbstractZerosMatrix, b::AbstractVector) = mult_zeros(a, b)
 
 function lmul_diag(a::Diagonal, b)
@@ -139,60 +155,56 @@ end
 function mul!(y::AbstractVector, A::AbstractFillMatrix, b::AbstractFillVector, alpha::Number, beta::Number)
     check_matmul_sizes(y, A, b)
 
-    αAb = alpha * getindex_value(A) * getindex_value(b) * length(b)
+    Abα = Ref(getindex_value(A) * getindex_value(b) * alpha * length(b))
 
     if iszero(beta)
-        y .= αAb
+        y .= Abα
     else
-        y .= αAb .+ beta .* y
+        y .= Abα .+ y .* beta
     end
     y
 end
 
 function mul!(y::StridedVector, A::StridedMatrix, b::AbstractFillVector, alpha::Number, beta::Number)
     check_matmul_sizes(y, A, b)
 
-    αb = alpha * getindex_value(b)
+    bα = Ref(getindex_value(b) * alpha)
 
     if iszero(beta)
-        y .= zero(eltype(y))
-        for col in eachcol(A)
-            y .+= αb .* col
-        end
+        y .= Ref(zero(eltype(y)))
     else
-        lmul!(beta, y)
-        for col in eachcol(A)
-            y .+= αb .* col
-        end
+        rmul!(y, beta)
+    end
+    for Acol in eachcol(A)
+        @. y += Acol * bα
     end
     y
 end
 
 function mul!(y::StridedVector, A::AbstractFillMatrix, b::StridedVector, alpha::Number, beta::Number)
     check_matmul_sizes(y, A, b)
 
-    αA = alpha * getindex_value(A)
+    Abα = Ref(getindex_value(A) * sum(b) * alpha)
 
     if iszero(beta)
-        y .= αA .* sum(b)
+        y .= Abα
     else
-        y .= αA .* sum(b) .+ beta .* y
+        y .= Abα .+ y .* beta
     end
     y
 end
 
-function _mul_adjtrans!(y::AbstractVector, A::AbstractMatrix, b::AbstractVector, alpha, beta, f)
-    α = alpha * getindex_value(b)
-
+function _mul_adjtrans!(y::AbstractVector, A::AbstractMatrix, b::AbstractFillVector, alpha, beta, f)
+    bα = getindex_value(b) * alpha
     At = f(A)
 
     if iszero(beta)
-        for (ind, col) in zip(eachindex(y), eachcol(At))
-            y[ind] = α .* f(sum(col))
+        for (ind, Atcol) in zip(eachindex(y), eachcol(At))
+            y[ind] = f(sum(Atcol)) * bα
         end
     else
-        for (ind, col) in zip(eachindex(y), eachcol(At))
-            y[ind] = α .* f(sum(col)) .+ beta .* y[ind]
+        for (ind, Atcol) in zip(eachindex(y), eachcol(At))
+            y[ind] = f(sum(Atcol)) * bα .+ y[ind] .* beta
         end
     end
     y
@@ -207,11 +219,11 @@ end
 
 function mul!(C::AbstractMatrix, A::AbstractFillMatrix, B::AbstractFillMatrix, alpha::Number, beta::Number)
     check_matmul_sizes(C, A, B)
-    αAB = alpha * getindex_value(A) * getindex_value(B) * size(B,1)
+    ABα = getindex_value(A) * getindex_value(B) * alpha * size(B,1)
     if iszero(beta)
-        C .= αAB
+        C .= ABα
     else
-        C .= αAB .+ beta .* C
+        C .= ABα .+ C .* beta
     end
     C
 end
@@ -237,7 +249,7 @@ _firstcol(C::Union{Adjoint, Transpose}) = view(parent(C), 1, :)
 function _mulfill!(C, A, B::AbstractFillMatrix, alpha, beta)
     check_matmul_sizes(C, A, B)
     if iszero(size(B,2))
-        return lmul!(beta, C)
+        return rmul!(C, beta)
     end
     mul!(_firstcol(C), A, view(B, :, 1), alpha, beta)
     copyfirstcol!(C)
@@ -288,13 +300,25 @@ function _adjvec_mul_zeros(a, b)
     return a1 * b[1]
 end
 
-*(a::AdjointAbsVec{<:Any,<:AbstractZerosVector}, b::AbstractMatrix) = (b' * a')'
+for MT in (:AbstractMatrix, :AbstractTriangular, :(Adjoint{<:Any,<:TransposeAbsVec}))
+    @eval *(a::AdjointAbsVec{<:Any,<:AbstractZerosVector}, b::$MT) = (b' * a')'
+end
+# ambiguity
+function *(a::AdjointAbsVec{<:Any,<:AbstractZerosVector}, b::TransposeAbsVec{<:Any,<:AdjointAbsVec})
+    # change from Transpose ∘ Adjoint to Adjoint ∘ Transpose
+    b2 = adjoint(transpose(adjoint(transpose(b))))
+    a * b2
+end
 *(a::AdjointAbsVec{<:Any,<:AbstractZerosVector}, b::AbstractZerosMatrix) = (b' * a')'
-*(a::TransposeAbsVec{<:Any,<:AbstractZerosVector}, b::AbstractMatrix) = transpose(transpose(b) * transpose(a))
+for MT in (:AbstractMatrix, :AbstractTriangular, :(Transpose{<:Any,<:AdjointAbsVec}))
+    @eval *(a::TransposeAbsVec{<:Any,<:AbstractZerosVector}, b::$MT) = transpose(transpose(b) * transpose(a))
+end
 *(a::TransposeAbsVec{<:Any,<:AbstractZerosVector}, b::AbstractZerosMatrix) = transpose(transpose(b) * transpose(a))
 
 *(a::AbstractVector, b::AdjOrTransAbsVec{<:Any,<:AbstractZerosVector}) = a * permutedims(parent(b))
-*(a::AbstractMatrix, b::AdjOrTransAbsVec{<:Any,<:AbstractZerosVector}) = a * permutedims(parent(b))
+for MT in (:AbstractMatrix, :AbstractTriangular)
+    @eval *(a::$MT, b::AdjOrTransAbsVec{<:Any,<:AbstractZerosVector}) = a * permutedims(parent(b))
+end
 *(a::AbstractZerosVector, b::AdjOrTransAbsVec{<:Any,<:AbstractZerosVector}) = a * permutedims(parent(b))
 *(a::AbstractZerosMatrix, b::AdjOrTransAbsVec{<:Any,<:AbstractZerosVector}) = a * permutedims(parent(b))
 
@@ -305,7 +329,8 @@ end
 
 *(a::Adjoint{T, <:AbstractMatrix{T}} where T, b::AbstractZeros{<:Any, 1}) = mult_zeros(a, b)
 
-*(D::Diagonal, a::AdjointAbsVec{<:Any,<:AbstractZerosVector}) = (a' * D')'
+*(D::Diagonal, a::Adjoint{<:Any,<:AbstractZerosVector}) = (a' * D')'
+*(D::Diagonal, a::Transpose{<:Any,<:AbstractZerosVector}) = transpose(transpose(a) * transpose(D))
 *(a::AdjointAbsVec{<:Any,<:AbstractZerosVector}, D::Diagonal) = (D' * a')'
 *(a::TransposeAbsVec{<:Any,<:AbstractZerosVector}, D::Diagonal) = transpose(D*transpose(a))
 function _triple_zeromul(x, D::Diagonal, y)
@@ -323,7 +348,7 @@ end
 *(x::TransposeAbsVec{<:Any,<:AbstractZerosVector}, D::Diagonal, y::AbstractZerosVector) = _triple_zeromul(x, D, y)
 
 
-function *(a::Transpose{T, <:AbstractVector{T}}, b::AbstractZerosVector{T}) where T<:Real
+function *(a::Transpose{T, <:AbstractVector}, b::AbstractZerosVector{T}) where T<:Real
     la, lb = length(a), length(b)
     if la ≠ lb
         throw(DimensionMismatch("dot product arguments have lengths $la and $lb"))
@@ -496,3 +521,11 @@ function kron(f::AbstractFillVecOrMat, g::AbstractFillVecOrMat)
     sz = _kronsize(f, g)
     return _kron(f, g, sz)
 end
+
+# bandedness
+function LinearAlgebra.istriu(A::AbstractFillMatrix, k::Integer = 0)
+    iszero(A) || k <= -(size(A,1)-1)
+end
+function LinearAlgebra.istril(A::AbstractFillMatrix, k::Integer = 0)
+    iszero(A) || k >= size(A,2)-1
+end

src/trues.jl

@@ -7,17 +7,24 @@ Typically created using `Trues(dims)` or `Trues(dims...)`
 
 # Example
 ```jldoctest
-julia> Trues(1,3)
-1×3 Ones{Bool,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}} = true
+julia> T = Trues(1,3)
+1×3 Ones{Bool}
 
-julia> Trues((2,3))
-2×3 Ones{Bool,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}} = true
+julia> Array(T)
+1×3 Matrix{Bool}:
+ 1  1  1
 ```
 """
 const Trues = Ones{Bool, N, Axes} where {N, Axes}
 
 
-""" `Falses = Zeros{Bool, N, Axes}` (see `Trues`) """
+"""
+    Falses = Zeros{Bool, N, Axes}
+
+Lazy version of `falses` with axes.
+
+See also: [`Trues`](@ref)
+"""
 const Falses = Zeros{Bool, N, Axes} where {N, Axes}