diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index a76e2c2a..bd994b33 100644 --- a/src/LinearMaps.jl +++ b/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,24 +27,24 @@ 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, α) @@ -45,7 +52,7 @@ function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap{T}, x::AbstractVecto 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], α, β) diff --git a/src/blockmap.jl b/src/blockmap.jl index 42e71443..04bd4d26 100644 --- a/src/blockmap.jl +++ b/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 diff --git a/src/composition.jl b/src/composition.jl index f7df3a14..0fefc57c 100644 --- a/src/composition.jl +++ b/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₂ diff --git a/src/functionmap.jl b/src/functionmap.jl index 6f1dbd84..097af290 100644 --- a/src/functionmap.jl +++ b/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 diff --git a/src/linearcombination.jl b/src/linearcombination.jl index 55e2e6d3..7f7d7340 100644 --- a/src/linearcombination.jl +++ b/src/linearcombination.jl @@ -20,10 +20,26 @@ 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("+")) @@ -31,10 +47,10 @@ function Base.:(+)(A₁::LinearMap, A₂::LinearCombination) 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₂) diff --git a/src/transpose.jl b/src/transpose.jl index 8de084c5..d367da72 100644 --- a/src/transpose.jl +++ b/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) diff --git a/src/uniformscalingmap.jl b/src/uniformscalingmap.jl index 01c96470..7c2de56a 100644 --- a/src/uniformscalingmap.jl +++ b/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(α) diff --git a/test/composition.jl b/test/composition.jl index 16ec06e3..749ad10c 100644 --- a/test/composition.jl +++ b/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) diff --git a/test/functionmap.jl b/test/functionmap.jl index 04336915..168965ef 100644 --- a/test/functionmap.jl +++ b/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 diff --git a/test/uniformscalingmap.jl b/test/uniformscalingmap.jl index 437eb3af..140592c9 100644 --- a/test/uniformscalingmap.jl +++ b/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)