miscellaneous updates (#66)

* require one-based indexing in indexing muls

* add inference test

* make alpha and beta Bool by default in 5-arg mul!

* prevent UniformScalingMaps with negative size

* add some docstrings

* remove unnecessary type specification

* use isone/iszero

Author: dkarrasch

src/LinearMaps.jl

@@ -5,6 +5,13 @@ export LinearMap
 using LinearAlgebra
 using SparseArrays
 
+if VERSION < v"1.2-"
+    import Base: has_offset_axes
+    require_one_based_indexing(A...) = !has_offset_axes(A...) || throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))
+else
+    import Base: require_one_based_indexing
+end
+
 abstract type LinearMap{T} end
 
 Base.eltype(::LinearMap{T}) where {T} = T
@@ -20,32 +27,32 @@ Base.size(A::LinearMap, n) = (n==1 || n==2 ? size(A)[n] : error("LinearMap objec
 Base.length(A::LinearMap) = size(A)[1] * size(A)[2]
 
 Base.:(*)(A::LinearMap, x::AbstractVector) = mul!(similar(x, promote_type(eltype(A), eltype(x)), size(A, 1)), A, x)
-function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap{T}, x::AbstractVector, α::Number=one(T), β::Number=zero(T)) where {T}
+function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap, x::AbstractVector, α::Number=true, β::Number=false)
     length(y) == size(A, 1) || throw(DimensionMismatch("mul!"))
-    if α == 1
-        β == 0 && (A_mul_B!(y, A, x); return y)
-        β == 1 && (y .+= A * x; return y)
+    if isone(α)
+        iszero(β) && (A_mul_B!(y, A, x); return y)
+        isone(β) && (y .+= A * x; return y)
         # β != 0, 1
         rmul!(y, β)
         y .+= A * x
         return y
-    elseif α == 0
-        β == 0 && (fill!(y, zero(eltype(y))); return y)
-        β == 1 && return y
+    elseif iszero(α)
+        iszero(β) && (fill!(y, zero(eltype(y))); return y)
+        isone(β) && return y
         # β != 0, 1
         rmul!(y, β)
         return y
     else # α != 0, 1
-        β == 0 && (A_mul_B!(y, A, x); rmul!(y, α); return y)
-        β == 1 && (y .+= rmul!(A * x, α); return y)
+        iszero(β) && (A_mul_B!(y, A, x); rmul!(y, α); return y)
+        isone(β) && (y .+= rmul!(A * x, α); return y)
         # β != 0, 1
         rmul!(y, β)
         y .+= rmul!(A * x, α)
         return y
     end
 end
 # the following is of interest in, e.g., subspace-iteration methods
-function LinearAlgebra.mul!(Y::AbstractMatrix, A::LinearMap{T}, X::AbstractMatrix, α::Number=one(T), β::Number=zero(T)) where {T}
+function LinearAlgebra.mul!(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix, α::Number=true, β::Number=false)
     (size(Y, 1) == size(A, 1) && size(X, 1) == size(A, 2) && size(Y, 2) == size(X, 2)) || throw(DimensionMismatch("mul!"))
     @inbounds @views for i = 1:size(X, 2)
         mul!(Y[:, i], A, X[:, i], α, β)

src/blockmap.jl

@@ -50,7 +50,7 @@ end
 Base.size(A::BlockMap) = (last(A.rowranges[end]), last(A.colranges[end]))
 
 ############
-# hcat
+# concatenation
 ############
 
 for k in 1:8 # is 8 sufficient?
@@ -63,6 +63,32 @@ for k in 1:8 # is 8 sufficient?
     @eval Base.hvcat(rows::Tuple{Vararg{Int}}, $(Is...), $L, As::Union{LinearMap,UniformScaling}...) = _hvcat(rows, $(args...), As...)
 end
 
+############
+# hcat
+############
+"""
+    hcat(As::Union{LinearMap,UniformScaling}...)
+
+Construct a `BlockMap <: LinearMap` object, a (lazy) representation of the
+horizontal concatenation of the arguments. `UniformScaling` objects are promoted
+to `LinearMap` automatically. To avoid fallback to the generic [`Base.hcat`](@ref),
+there must be a `LinearMap` object among the first 8 arguments.
+
+# Examples
+```jldoctest; setup=(using LinearMaps)
+julia> CS = LinearMap{Int}(cumsum, 3)::LinearMaps.FunctionMap;
+
+julia> L = [CS LinearMap(ones(Int, 3, 3))]::LinearMaps.BlockMap;
+
+julia> L * ones(Int, 6)
+3-element Array{Int64,1}:
+ 4
+ 5
+ 6
+```
+"""
+Base.hcat
+
 function _hcat(As::Union{LinearMap,UniformScaling}...)
     T = promote_type(map(eltype, As)...)
     nbc = length(As)
@@ -82,6 +108,31 @@ end
 ############
 # vcat
 ############
+"""
+    vcat(As::Union{LinearMap,UniformScaling}...)
+
+Construct a `BlockMap <: LinearMap` object, a (lazy) representation of the
+vertical concatenation of the arguments. `UniformScaling` objects are promoted
+to `LinearMap` automatically. To avoid fallback to the generic [`Base.vcat`](@ref),
+there must be a `LinearMap` object among the first 8 arguments.
+
+# Examples
+```jldoctest; setup=(using LinearMaps)
+julia> CS = LinearMap{Int}(cumsum, 3)::LinearMaps.FunctionMap;
+
+julia> L = [CS; LinearMap(ones(Int, 3, 3))]::LinearMaps.BlockMap;
+
+julia> L * ones(Int, 3)
+6-element Array{Int64,1}:
+ 1
+ 2
+ 3
+ 3
+ 3
+ 3
+```
+"""
+Base.vcat
 
 function _vcat(As::Union{LinearMap,UniformScaling}...)
     T = promote_type(map(eltype, As)...)
@@ -103,6 +154,35 @@ end
 ############
 # hvcat
 ############
+"""
+    hvcat(rows::Tuple{Vararg{Int}}, As::Union{LinearMap,UniformScaling}...)
+
+Construct a `BlockMap <: LinearMap` object, a (lazy) representation of the
+horizontal-vertical concatenation of the arguments. The first argument specifies
+the number of arguments to concatenate in each block row. `UniformScaling` objects
+are promoted to `LinearMap` automatically. To avoid fallback to the generic
+[`Base.hvcat`](@ref), there must be a `LinearMap` object among the first 8 arguments.
+
+# Examples
+```jldoctest; setup=(using LinearMaps)
+julia> CS = LinearMap{Int}(cumsum, 3)::LinearMaps.FunctionMap;
+
+julia> L = [CS CS; CS CS]::LinearMaps.BlockMap;
+
+julia> L.rows
+(2, 2)
+
+julia> L * ones(Int, 6)
+6-element Array{Int64,1}:
+ 2
+ 4
+ 6
+ 2
+ 4
+ 6
+```
+"""
+Base.hvcat
 
 function _hvcat(rows::Tuple{Vararg{Int}}, As::Union{LinearMap,UniformScaling}...)
     nr = length(rows)
@@ -221,6 +301,7 @@ LinearAlgebra.adjoint(A::BlockMap)  = AdjointMap(A)
 ############
 
 function A_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector)
+    require_one_based_indexing(y, x)
     m, n = size(A)
     @boundscheck (m == length(y) && n == length(x)) || throw(DimensionMismatch("A_mul_B!"))
     maps, rows, yinds, xinds = A.maps, A.rows, A.rowranges, A.colranges
@@ -238,6 +319,7 @@ function A_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector)
 end
 
 function At_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector)
+    require_one_based_indexing(y, x)
     m, n = size(A)
     @boundscheck (n == length(y) && m == length(x)) || throw(DimensionMismatch("At_mul_B!"))
     maps, rows, xinds, yinds = A.maps, A.rows, A.rowranges, A.colranges
@@ -262,6 +344,7 @@ function At_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector)
 end
 
 function Ac_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector)
+    require_one_based_indexing(y, x)
     m, n = size(A)
     @boundscheck (n == length(y) && m == length(x)) || throw(DimensionMismatch("At_mul_B!"))
     maps, rows, xinds, yinds = A.maps, A.rows, A.rowranges, A.colranges

src/composition.jl

@@ -70,10 +70,26 @@ Base.:(/)(A::LinearMap, α::Number) = A * inv(α)
 Base.:(-)(A::LinearMap) = -1 * A
 
 # composition of linear maps
-function Base.:(*)(A₁::CompositeMap, A₂::CompositeMap)
+"""
+    A::LinearMap * B::LinearMap
+
+Construct a `CompositeMap <: LinearMap`, a (lazy) representation of the product
+of the two operators. Products of `LinearMap`/`CompositeMap` objects and
+`LinearMap`/`CompositeMap` objects are reduced to a single `CompositeMap`.
+In products of `LinearMap`s and `AbstractMatrix`/`UniformScaling` objects, the
+latter get promoted to `LinearMap`s automatically.
+
+# Examples
+```jldoctest; setup=(using LinearAlgebra, LinearMaps)
+julia> CS = LinearMap{Int}(cumsum, 3)::LinearMaps.FunctionMap;
+
+julia> LinearMap(ones(Int, 3, 3)) * CS * I * rand(3, 3);
+```
+"""
+function Base.:(*)(A₁::LinearMap, A₂::LinearMap)
     size(A₁, 2) == size(A₂, 1) || throw(DimensionMismatch("*"))
     T = promote_type(eltype(A₁), eltype(A₂))
-    return CompositeMap{T}(tuple(A₂.maps..., A₁.maps...))
+    return CompositeMap{T}(tuple(A₂, A₁))
 end
 function Base.:(*)(A₁::LinearMap, A₂::CompositeMap)
     size(A₁, 2) == size(A₂, 1) || throw(DimensionMismatch("*"))
@@ -85,10 +101,10 @@ function Base.:(*)(A₁::CompositeMap, A₂::LinearMap)
     T = promote_type(eltype(A₁), eltype(A₂))
     return CompositeMap{T}(tuple(A₂, A₁.maps...))
 end
-function Base.:(*)(A₁::LinearMap, A₂::LinearMap)
+function Base.:(*)(A₁::CompositeMap, A₂::CompositeMap)
     size(A₁, 2) == size(A₂, 1) || throw(DimensionMismatch("*"))
     T = promote_type(eltype(A₁), eltype(A₂))
-    return CompositeMap{T}(tuple(A₂, A₁))
+    return CompositeMap{T}(tuple(A₂.maps..., A₁.maps...))
 end
 Base.:(*)(A₁::LinearMap, A₂::UniformScaling) = A₁ * A₂.λ
 Base.:(*)(A₁::UniformScaling, A₂::LinearMap) = A₁.λ * A₂

src/functionmap.jl

@@ -50,19 +50,19 @@ function Base.:(*)(A::FunctionMap, x::AbstractVector)
 end
 
 function A_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
-    (length(x) == A.N && length(y) == A.M) || throw(DimensionMismatch())
+    (length(x) == A.N && length(y) == A.M) || throw(DimensionMismatch("A_mul_B!"))
     ismutating(A) ? A.f(y, x) : copyto!(y, A.f(x))
     return y
 end
 
 function At_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
     issymmetric(A) && return A_mul_B!(y, A, x)
-    (length(x) == A.M && length(y) == A.N) || throw(DimensionMismatch())
+    (length(x) == A.M && length(y) == A.N) || throw(DimensionMismatch("At_mul_B!"))
     if A.fc !== nothing
         if !isreal(A)
             x = conj(x)
         end
-        (ismutating(A) ? A.fc(y, x) : copyto!(y, A.fc(x)))
+        ismutating(A) ? A.fc(y, x) : copyto!(y, A.fc(x))
         if !isreal(A)
             conj!(y)
         end
@@ -74,9 +74,10 @@ end
 
 function Ac_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
     ishermitian(A) && return A_mul_B!(y, A, x)
-    (length(x) == A.M && length(y) == A.N) || throw(DimensionMismatch())
+    (length(x) == A.M && length(y) == A.N) || throw(DimensionMismatch("Ac_mul_B!"))
     if A.fc !== nothing
-        return (ismutating(A) ? A.fc(y, x) : copyto!(y, A.fc(x)))
+        ismutating(A) ? A.fc(y, x) : copyto!(y, A.fc(x))
+        return y
     else
         error("adjoint not implemented for $A")
     end

src/linearcombination.jl

@@ -20,21 +20,37 @@ LinearAlgebra.ishermitian(A::LinearCombination) = all(ishermitian, A.maps) # suf
 LinearAlgebra.isposdef(A::LinearCombination) = all(isposdef, A.maps) # sufficient but not necessary
 
 # adding linear maps
-function Base.:(+)(A₁::LinearCombination, A₂::LinearCombination)
+"""
+    A::LinearMap + B::LinearMap
+
+Construct a `LinearCombination <: LinearMap`, a (lazy) representation of the sum
+of the two operators. Sums of `LinearMap`/`LinearCombination` objects and
+`LinearMap`/`LinearCombination` objects are reduced to a single `LinearCombination`.
+In sums of `LinearMap`s and `AbstractMatrix`/`UniformScaling` objects, the latter
+get promoted to `LinearMap`s automatically.
+
+# Examples
+```jldoctest; setup=(using LinearAlgebra, LinearMaps)
+julia> CS = LinearMap{Int}(cumsum, 3)::LinearMaps.FunctionMap;
+
+julia> LinearMap(ones(Int, 3, 3)) + CS + I + rand(3, 3);
+```
+"""
+function Base.:(+)(A₁::LinearMap, A₂::LinearMap)
     size(A₁) == size(A₂) || throw(DimensionMismatch("+"))
     T = promote_type(eltype(A₁), eltype(A₂))
-    return LinearCombination{T}(tuple(A₁.maps..., A₂.maps...))
+    return LinearCombination{T}(tuple(A₁, A₂))
 end
 function Base.:(+)(A₁::LinearMap, A₂::LinearCombination)
     size(A₁) == size(A₂) || throw(DimensionMismatch("+"))
     T = promote_type(eltype(A₁), eltype(A₂))
     return LinearCombination{T}(tuple(A₁, A₂.maps...))
 end
 Base.:(+)(A₁::LinearCombination, A₂::LinearMap) = +(A₂, A₁)
-function Base.:(+)(A₁::LinearMap, A₂::LinearMap)
+function Base.:(+)(A₁::LinearCombination, A₂::LinearCombination)
     size(A₁) == size(A₂) || throw(DimensionMismatch("+"))
     T = promote_type(eltype(A₁), eltype(A₂))
-    return LinearCombination{T}(tuple(A₁, A₂))
+    return LinearCombination{T}(tuple(A₁.maps..., A₂.maps...))
 end
 Base.:(-)(A₁::LinearMap, A₂::LinearMap) = +(A₁, -A₂)
 

src/transpose.jl

@@ -14,7 +14,7 @@ LinearAlgebra.adjoint(A::LinearMap{<:Real}) = transpose(A)
 LinearAlgebra.adjoint(A::LinearMap) = AdjointMap(A)
 
 # properties
-Base.size(A::Union{TransposeMap, AdjointMap}) = (size(A.lmap, 2), size(A.lmap, 1))
+Base.size(A::Union{TransposeMap, AdjointMap}) = reverse(size(A.lmap))
 LinearAlgebra.issymmetric(A::Union{TransposeMap, AdjointMap}) = issymmetric(A.lmap)
 LinearAlgebra.ishermitian(A::Union{TransposeMap, AdjointMap}) = ishermitian(A.lmap)
 LinearAlgebra.isposdef(A::Union{TransposeMap, AdjointMap}) = isposdef(A.lmap)

src/uniformscalingmap.jl

@@ -1,10 +1,14 @@
 struct UniformScalingMap{T} <: LinearMap{T}
     λ::T
     M::Int
+    function UniformScalingMap(λ::Number, M::Int)
+        M < 0 && throw(ArgumentError("size of UniformScalingMap must be non-negative, got $M"))
+        return new{typeof(λ)}(λ, M)
+    end
 end
-UniformScalingMap(λ::T, M::Int, N::Int) where {T} =
+UniformScalingMap(λ::Number, M::Int, N::Int) =
     (M == N ? UniformScalingMap(λ, M) : error("UniformScalingMap needs to be square"))
-UniformScalingMap(λ::T, sz::Tuple{Int, Int}) where {T} =
+UniformScalingMap(λ::Number, sz::Tuple{Int, Int}) =
     (sz[1] == sz[2] ? UniformScalingMap(λ, sz[1]) : error("UniformScalingMap needs to be square"))
 
 # properties
@@ -49,7 +53,7 @@ function A_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector)
 end
 end # VERSION
 
-function LinearAlgebra.mul!(y::AbstractVector, J::UniformScalingMap{T}, x::AbstractVector, α::Number=one(T), β::Number=zero(T)) where {T}
+function LinearAlgebra.mul!(y::AbstractVector, J::UniformScalingMap, x::AbstractVector, α::Number=true, β::Number=false)
     @boundscheck (length(x) == length(y) == J.M || throw(DimensionMismatch("mul!")))
     λ = J.λ
     @inbounds if isone(α)

test/composition.jl

@@ -15,6 +15,7 @@ LinearAlgebra.mul!(y::Vector, A::Union{SimpleFunctionMap,SimpleComplexFunctionMa
 
 @testset "composition" begin
     F = @inferred LinearMap(cumsum, y -> reverse(cumsum(reverse(x))), 10; ismutating=false)
+    FC = @inferred LinearMap{ComplexF64}(cumsum, y -> reverse(cumsum(reverse(x))), 10; ismutating=false)
     A = 2 * rand(ComplexF64, (10, 10)) .- 1
     B = rand(size(A)...)
     M = @inferred 1 * LinearMap(A)
@@ -32,6 +33,10 @@ LinearAlgebra.mul!(y::Vector, A::Union{SimpleFunctionMap,SimpleComplexFunctionMa
     @test @inferred ishermitian(N * N')
     @test @inferred !issymmetric(M' * M)
     @test @inferred ishermitian(M' * M)
+    @test @inferred issymmetric(F'F)
+    @test @inferred ishermitian(F'F)
+    @test @inferred !issymmetric(FC'FC)
+    @test @inferred ishermitian(FC'FC)
     @test @inferred isposdef(transpose(F) * F * 3)
     @test @inferred isposdef(transpose(F) * 3 * F)
     @test @inferred !isposdef(-5*transpose(F) * F)

test/functionmap.jl

@@ -57,6 +57,7 @@ using Test, LinearMaps, LinearAlgebra
     @test_throws ErrorException CS!'v
     @test_throws ErrorException adjoint(CS!) * v
     CS! = LinearMap{ComplexF64}(cumsum!, (y, x) -> (copyto!(y, x); reverse!(y); cumsum!(y, y)), 10; ismutating=true)
+    @inferred adjoint(CS!)
     @test @inferred LinearMaps.ismutating(CS!)
     @test @inferred CS! * v == cv
     @test @inferred *(CS!, v) == cv

test/uniformscalingmap.jl

@@ -1,6 +1,7 @@
 using Test, LinearMaps, LinearAlgebra, BenchmarkTools
 
 @testset "identity/scaling map" begin
+    @test_throws ArgumentError LinearMaps.UniformScalingMap(true, -1)
     A = 2 * rand(ComplexF64, (10, 10)) .- 1
     B = rand(size(A)...)
     M = @inferred 1 * LinearMap(A)