Add FillMaps

docs/src/types.md

@@ -84,6 +84,10 @@ Base.cat
 SparseArrays.blockdiag
 ```
 
+### `FillMap``
+
+Type for lazily representing constantly filled matrices.
+
 ## Methods
 
 ### Multiplication methods

src/LinearMaps.jl

@@ -239,18 +239,22 @@ include("composition.jl") # composition of linear maps
 include("functionmap.jl") # using a function as linear map
 include("blockmap.jl") # block linear maps
 include("kronecker.jl") # Kronecker product of linear maps
+include("fillmap.jl") # linear maps representing constantly filled matrices
 include("conversion.jl") # conversion of linear maps to matrices
 include("show.jl") # show methods for LinearMap objects
 
 """
     LinearMap(A::LinearMap; kwargs...)::WrappedMap
     LinearMap(A::AbstractMatrix; kwargs...)::WrappedMap
     LinearMap(J::UniformScaling, M::Int)::UniformScalingMap
+    LinearMap(λ::Number, M::Int, N::Int) = FillMap(λ, (M, N))::FillMap
+    LinearMap(λ::Number, dims::Dims{2}) = FillMap(λ, dims)::FillMap
     LinearMap{T=Float64}(f, [fc,], M::Int, N::Int = M; kwargs...)::FunctionMap
 
 Construct a linear map object, either from an existing `LinearMap` or `AbstractMatrix` `A`,
 with the purpose of redefining its properties via the keyword arguments `kwargs`;
-a `UniformScaling` object `J` with specified (square) dimension `M`; or
+a `UniformScaling` object `J` with specified (square) dimension `M`; from a `Number`
+object to lazily represent filled matrices; or
 from a function or callable object `f`. In the latter case, one also needs to specify
 the size of the equivalent matrix representation `(M, N)`, i.e., for functions `f` acting
 on length `N` vectors and producing length `M` vectors (with default value `N=M`).
@@ -281,6 +285,8 @@ LinearMap(f, M::Int; kwargs...) = LinearMap{Float64}(f, M; kwargs...)
 LinearMap(f, M::Int, N::Int; kwargs...) = LinearMap{Float64}(f, M, N; kwargs...)
 LinearMap(f, fc, M::Int; kwargs...) = LinearMap{Float64}(f, fc, M; kwargs...)
 LinearMap(f, fc, M::Int, N::Int; kwargs...) = LinearMap{Float64}(f, fc, M, N; kwargs...)
+LinearMap(λ::Number, M::Int, N::Int) = FillMap(λ, (M, N))
+LinearMap(λ::Number, dims::Dims{2}) = FillMap(λ, dims)
 
 LinearMap{T}(A::MapOrMatrix; kwargs...) where {T} = WrappedMap{T}(A; kwargs...)
 LinearMap{T}(f, args...; kwargs...) where {T} = FunctionMap{T}(f, args...; kwargs...)

src/conversion.jl

@@ -171,3 +171,6 @@ function SparseArrays.sparse(L::KroneckerSumMap)
     IB = sparse(Diagonal(ones(Bool, size(B, 1))))
     return kron(convert(AbstractMatrix, A), IB) + kron(IA, convert(AbstractMatrix, B))
 end
+
+# FillMap
+Base.Matrix{T}(A::FillMap) where {T} = fill(T(A.λ), size(A))

src/fillmap.jl

@@ -0,0 +1,56 @@
+struct FillMap{T} <: LinearMap{T}
+    λ::T
+    size::Dims{2}
+    function FillMap(λ::T, dims::Dims{2}) where {T}
+        (dims[1]>=0 && dims[2]>=0) || throw(ArgumentError("dims of FillMap must be non-negative"))
+        return new{T}(λ, dims)
+    end
+end
+
+# properties
+Base.size(A::FillMap) = A.size
+MulStyle(A::FillMap) = FiveArg()
+LinearAlgebra.issymmetric(A::FillMap) = A.size[1] == A.size[2]
+LinearAlgebra.ishermitian(A::FillMap) = isreal(A.λ) && A.size[1] == A.size[2]
+LinearAlgebra.isposdef(A::FillMap) = (size(A, 1) == size(A, 2) == 1 && isposdef(A.λ))
+Base.:(==)(A::FillMap, B::FillMap) = A.λ == B.λ && A.size == B.size
+
+LinearAlgebra.adjoint(A::FillMap) = FillMap(adjoint(A.λ), reverse(A.size))
+LinearAlgebra.transpose(A::FillMap) = FillMap(transpose(A.λ), reverse(A.size))
+
+function Base.:(*)(A::FillMap, x::AbstractVector)
+    T = typeof(oneunit(eltype(A)) * (zero(eltype(x)) + zero(eltype(x))))
+    return fill(iszero(A.λ) ? zero(T) : A.λ*sum(x), A.size[1])
+end
+
+function _unsafe_mul!(y::AbstractVecOrMat, A::FillMap, x::AbstractVector)
+    return fill!(y, iszero(A.λ) ? zero(eltype(y)) : A.λ*sum(x))
+end
+
+function _unsafe_mul!(y::AbstractVecOrMat, A::FillMap, x::AbstractVector, α::Number, β::Number)
+    if iszero(α)
+        !isone(β) && rmul!(y, β)
+    else
+        temp = A.λ * sum(x) * α
+        if iszero(β)
+            y .= temp
+        elseif isone(β)
+            y .+= temp
+        else
+            y .= y .* β .+ temp
+        end
+    end
+    return y
+end
+
+Base.:(+)(A::FillMap, B::FillMap) = A.size == B.size ? FillMap(A.λ + B.λ, A.size) : throw(DimensionMismatch())
+Base.:(-)(A::FillMap) = FillMap(-A.λ, A.size)
+Base.:(*)(λ::Number, A::FillMap) = FillMap(λ * A.λ, size(A))
+Base.:(*)(A::FillMap, λ::Number) = FillMap(A.λ * λ, size(A))
+Base.:(*)(λ::RealOrComplex, A::FillMap) = FillMap(λ * A.λ, size(A))
+Base.:(*)(A::FillMap, λ::RealOrComplex) = FillMap(A.λ * λ, size(A))
+
+function Base.:(*)(A::FillMap, B::FillMap)
+    check_dim_mul(A, B)
+    return FillMap(A.λ*B.λ*size(A, 2), (size(A, 1), size(B, 2)))
+end

src/show.jl

@@ -39,6 +39,9 @@ end
 function _show(io::IO, A::ScaledMap{T}, i) where {T}
     " with scale: $(A.λ) of\n" * map_show(io, A.lmap, i+2)
 end
+function _show(io::IO, A::FillMap{T}, _) where {T}
+    " with fill value: $(A.λ)"
+end
 
 # helper functions
 function _show_typeof(A::LinearMap{T}) where {T}

test/fillmap.jl

@@ -0,0 +1,40 @@
+using LinearMaps, LinearAlgebra, Test
+
+@testset "filled maps" begin
+    M, N = 2, 3
+    μ = rand()
+    for λ in (true, false, 3, μ, μ + 2im)
+        L = LinearMap(λ, (M, N))
+        @test L == LinearMap(λ, M, N)
+        @test occursin("$M×$N LinearMaps.FillMap{$(typeof(λ))} with fill value: $λ", sprint((t, s) -> show(t, "text/plain", s), L))
+        @test LinearMaps.MulStyle(L) === LinearMaps.FiveArg()
+        A = fill(λ, (M, N))
+        x = rand(typeof(λ) <: Real ? Float64 : ComplexF64, 3)
+        X = rand(typeof(λ) <: Real ? Float64 : ComplexF64, 3, 4)
+        w = similar(x, 2)
+        W = similar(X, 2, 4)
+        @test size(L) == (M, N)
+        @test adjoint(L) == LinearMap(adjoint(λ), (3,2))
+        @test transpose(L) == LinearMap(λ, (3,2))
+        @test Matrix(L) == A
+        @test L * x ≈ A * x
+        @test mul!(w, L, x) ≈ A * x
+        @test mul!(W, L, X) ≈ A * X
+        for α in (true, false, 1, 0, randn()), β in (true, false, 1, 0, randn())
+            @test mul!(copy(w), L, x, α, β) ≈ fill(λ * sum(x) * α, M) + w * β
+            @test mul!(copy(W), L, X, α, β) ≈ λ * reduce(vcat, sum(X, dims=1) for _ in 1:2) * α + W * β
+        end
+    end
+    @test issymmetric(LinearMap(μ + 1im, (3, 3)))
+    @test ishermitian(LinearMap(μ + 0im, (3, 3)))
+    @test isposdef(LinearMap(μ, (1,1))) == isposdef(μ)
+    @test !isposdef(LinearMap(μ, (3,3)))
+    α = rand()
+    β = rand()
+    @test LinearMap(μ, (M, N)) + LinearMap(α, (M, N)) == LinearMap(μ + α, (M, N))
+    @test LinearMap(μ, (M, N)) - LinearMap(α, (M, N)) == LinearMap(μ - α, (M, N))
+    @test α*LinearMap(μ, (M, N)) == LinearMap(α * μ, (M, N))
+    @test LinearMap(μ, (M, N))*α == LinearMap(μ * α, (M, N))
+    @test LinearMap(μ, (M, N))*LinearMap(μ, (N, M)) == LinearMap(μ^2*N, (M, M))
+    @test Matrix(LinearMap(μ, (M, N))*LinearMap(μ, (N, M))) ≈ fill(μ, (M, N))*fill(μ, (N, M))
+end

test/numbertypes.jl

@@ -45,6 +45,7 @@ using Test, LinearMaps, LinearAlgebra, Quaternions
     @test M == (L * L * γ) * β == (L * L * α) * β == (L * L * α) * β.λ
     @test length(M.maps) == 3
     @test M.maps[1].λ == γ*β.λ
+    @test γ*LinearMap(γ, (3, 4)) == LinearMap(γ^2, (3, 4)) == LinearMap(γ, (3, 4))*γ
 
     # exercise non-RealOrComplex scalar operations
     @test Array(γ * (L'*L)) ≈ γ * (A'*A) # CompositeMap

test/runtests.jl

@@ -30,3 +30,5 @@ include("kronecker.jl")
 include("conversion.jl")
 
 include("left.jl")
+
+include("fillmap.jl")