From 81513300a0465b1fbd96ca08533b35af7df5d5ba Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 20 Nov 2019 16:08:09 +0100 Subject: [PATCH 1/3] add mulstyle feature refactor linear combination multiplication try to get rid of allocations in lincomb mul show coefficients more details in uniform scaling mul opt-out of allocation testing fix mulstyle error in kronecker, improve coverage --- src/LinearMaps.jl | 17 +++++ src/blockmap.jl | 2 + src/linearcombination.jl | 129 ++++++++++++++++---------------------- src/transpose.jl | 2 + src/uniformscalingmap.jl | 23 +++++-- src/wrappedmap.jl | 2 + test/blockmap.jl | 6 +- test/kronecker.jl | 8 ++- test/linearcombination.jl | 22 +++++-- test/runtests.jl | 5 ++ test/transpose.jl | 1 + test/uniformscalingmap.jl | 6 +- 12 files changed, 134 insertions(+), 89 deletions(-) diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index a8adf542..549e773d 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -19,6 +19,23 @@ const MapOrMatrix{T} = Union{LinearMap{T},AbstractMatrix{T}} Base.eltype(::LinearMap{T}) where {T} = T +abstract type MulStyle end + +struct FiveArg <: MulStyle end +struct ThreeArg <: MulStyle end + +mulstyle(::Type{FiveArg}, ::Type{FiveArg}) = FiveArg +mulstyle(::Type{ThreeArg}, ::Type{FiveArg}) = ThreeArg +mulstyle(::Type{FiveArg}, ::Type{ThreeArg}) = ThreeArg +mulstyle(::Type{ThreeArg}, ::Type{ThreeArg}) = ThreeArg +mulstyle(::LinearMap) = ThreeArg # default +if VERSION ≥ v"1.3.0-alpha.115" + mulstyle(::AbstractMatrix) = FiveArg +else + mulstyle(::AbstractMatrix) = ThreeArg +end +mulstyle(A::LinearMap, As::LinearMap...) = mulstyle(mulstyle(A), mulstyle(As...)) + Base.isreal(A::LinearMap) = eltype(A) <: Real LinearAlgebra.issymmetric(::LinearMap) = false # default assumptions LinearAlgebra.ishermitian(A::LinearMap{<:Real}) = issymmetric(A) diff --git a/src/blockmap.jl b/src/blockmap.jl index 04bd4d26..c62152f0 100644 --- a/src/blockmap.jl +++ b/src/blockmap.jl @@ -14,6 +14,8 @@ end BlockMap{T}(maps::As, rows::S) where {T,As<:Tuple{Vararg{LinearMap}},S} = BlockMap{T,As,S}(maps, rows) +mulstyle(A::BlockMap) = mulstyle(A.maps...) + function check_dim(A::LinearMap, dim, n) n == size(A, dim) || throw(DimensionMismatch("Expected $n, got $(size(A, dim))")) return nothing diff --git a/src/linearcombination.jl b/src/linearcombination.jl index 7f7d7340..50f119ea 100644 --- a/src/linearcombination.jl +++ b/src/linearcombination.jl @@ -1,17 +1,20 @@ -struct LinearCombination{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T} +struct LinearCombination{T, MS<:MulStyle, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T} maps::As - function LinearCombination{T, As}(maps::As) where {T, As} + function LinearCombination{T, MS, As}(maps::As) where {T, MS<:MulStyle, As} N = length(maps) sz = size(maps[1]) - for n in 1:N - size(maps[n]) == sz || throw(DimensionMismatch("LinearCombination")) - promote_type(T, eltype(maps[n])) == T || throw(InexactError()) + for Ai in maps + size(Ai) == sz || throw(DimensionMismatch("LinearCombination")) + promote_type(T, eltype(Ai)) == T || throw(InexactError()) end - new{T, As}(maps) + MS === FiveArg && mulstyle(maps...) === ThreeArg && throw("wrong mulstyle in constructor") + new{T, MS, As}(maps) end end -LinearCombination{T}(maps::As) where {T, As} = LinearCombination{T, As}(maps) +LinearCombination{T,MS}(maps::As) where {T, MS<:MulStyle, As} = LinearCombination{T, mulstyle(maps...), As}(maps) + +mulstyle(::LinearCombination{T,MS}) where {T, MS<:MulStyle} = MS # basic methods Base.size(A::LinearCombination) = size(A.maps[1]) @@ -39,18 +42,18 @@ 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₁, A₂)) + return LinearCombination{T, mulstyle(A₁, A₂)}(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...)) + return LinearCombination{T, mulstyle(A₁, A₂)}(tuple(A₁, A₂.maps...)) end Base.:(+)(A₁::LinearCombination, A₂::LinearMap) = +(A₂, A₁) function Base.:(+)(A₁::LinearCombination, A₂::LinearCombination) size(A₁) == size(A₂) || throw(DimensionMismatch("+")) T = promote_type(eltype(A₁), eltype(A₂)) - return LinearCombination{T}(tuple(A₁.maps..., A₂.maps...)) + return LinearCombination{T, mulstyle(A₁, A₂)}(tuple(A₁.maps..., A₂.maps...)) end Base.:(-)(A₁::LinearMap, A₂::LinearMap) = +(A₁, -A₂) @@ -58,87 +61,63 @@ Base.:(-)(A₁::LinearMap, A₂::LinearMap) = +(A₁, -A₂) Base.:(==)(A::LinearCombination, B::LinearCombination) = (eltype(A) == eltype(B) && A.maps == B.maps) # special transposition behavior -LinearAlgebra.transpose(A::LinearCombination) = LinearCombination{eltype(A)}(map(transpose, A.maps)) -LinearAlgebra.adjoint(A::LinearCombination) = LinearCombination{eltype(A)}(map(adjoint, A.maps)) +LinearAlgebra.transpose(A::LinearCombination) = LinearCombination{eltype(A), mulstyle(A)}(map(transpose, A.maps)) +LinearAlgebra.adjoint(A::LinearCombination) = LinearCombination{eltype(A), mulstyle(A)}(map(adjoint, A.maps)) # multiplication with vectors -if VERSION < v"1.3.0-alpha.115" - -function A_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) - # no size checking, will be done by individual maps - A_mul_B!(y, A.maps[1], x) - l = length(A.maps) - if l>1 - z = similar(y) - for n in 2:l - A_mul_B!(z, A.maps[n], x) - y .+= z - end +for Atype in (AbstractVector, AbstractMatrix) + @eval Base.@propagate_inbounds function LinearAlgebra.mul!(y::$Atype, A::LinearCombination, x::$Atype, + α::Number=true, β::Number=false) + @boundscheck check_dim_mul(y, A, x) + return _lincombmul!(y, A, x, α, β) end - return y end -else # 5-arg mul! is available for matrices - -# map types that have an allocation-free 5-arg mul! implementation -const FreeMap = Union{MatrixMap,UniformScalingMap} - -function A_mul_B!(y::AbstractVector, A::LinearCombination{T,As}, x::AbstractVector) where {T, As<:Tuple{Vararg{FreeMap}}} - # no size checking, will be done by individual maps - A_mul_B!(y, A.maps[1], x) - for n in 2:length(A.maps) - mul!(y, A.maps[n], x, true, true) - end - return y -end -function A_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) - # no size checking, will be done by individual maps - A_mul_B!(y, A.maps[1], x) - l = length(A.maps) - if l>1 - z = similar(y) - for n in 2:l - An = A.maps[n] - if An isa FreeMap - mul!(y, An, x, true, true) - else - A_mul_B!(z, A.maps[n], x) - y .+= z - end +@inline function _lincombmul!(y, A::LinearCombination{<:Any,FiveArg}, x, α::Number, β::Number) + if iszero(α) # trivial cases + iszero(β) && (fill!(y, zero(eltype(y))); return y) + isone(β) && return y + # β != 0, 1 + rmul!(y, β) + return y + else + mul!(y, first(A.maps), x, α, β) + @inbounds for An in Base.tail(A.maps) + mul!(y, An, x, α, true) end + return y end - return y end -function LinearAlgebra.mul!(y::AbstractVector, A::LinearCombination{T,As}, x::AbstractVector, α::Number=true, β::Number=false) where {T, As<:Tuple{Vararg{FreeMap}}} - length(y) == size(A, 1) || throw(DimensionMismatch("mul!")) - if isone(α) - iszero(β) && (A_mul_B!(y, A, x); return y) - !isone(β) && rmul!(y, β) - elseif iszero(α) +@inline function _lincombmul!(y, A::LinearCombination{<:Any,ThreeArg}, x, α::Number, β::Number) + if iszero(α) iszero(β) && (fill!(y, zero(eltype(y))); return y) isone(β) && return y # β != 0, 1 rmul!(y, β) return y - else # α != 0, 1 - if iszero(β) - A_mul_B!(y, A, x) - rmul!(y, α) - return y - elseif !isone(β) - rmul!(y, β) - end # β-cases - end # α-cases - - for An in A.maps - mul!(y, An, x, α, true) + else + mul!(y, first(A.maps), x, α, β) + l = length(A.maps) + if l>1 + z = similar(y) + @inbounds for n in 2:l + An = A.maps[n] + muladd!(mulstyle(An), y, An, x, α, z) + end + end + return y end - return y end -end # VERSION +@inline muladd!(::Type{FiveArg}, y, A, x, α, _) = mul!(y, A, x, α, true) +@inline function muladd!(::Type{ThreeArg}, y, A, x, α, z) + A_mul_B!(z, A, x) + y .+= isone(α) ? z : z .* α +end + +A_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = mul!(y, A, x) -At_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = A_mul_B!(y, transpose(A), x) +At_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = mul!(y, transpose(A), x) -Ac_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) +Ac_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = mul!(y, adjoint(A), x) diff --git a/src/transpose.jl b/src/transpose.jl index d367da72..9d158466 100644 --- a/src/transpose.jl +++ b/src/transpose.jl @@ -5,6 +5,8 @@ struct AdjointMap{T, A<:LinearMap{T}} <: LinearMap{T} lmap::A end +mulstyle(A::Union{TransposeMap,AdjointMap}) = mulstyle(A.lmap) + # transposition behavior of LinearMap objects LinearAlgebra.transpose(A::TransposeMap) = A.lmap LinearAlgebra.adjoint(A::AdjointMap) = A.lmap diff --git a/src/uniformscalingmap.jl b/src/uniformscalingmap.jl index 25461858..4bb77f3c 100644 --- a/src/uniformscalingmap.jl +++ b/src/uniformscalingmap.jl @@ -11,6 +11,8 @@ UniformScalingMap(λ::Number, M::Int, N::Int) = UniformScalingMap(λ::T, sz::Dims{2}) where {T} = (sz[1] == sz[2] ? UniformScalingMap(λ, sz[1]) : error("UniformScalingMap needs to be square")) +mulstyle(::UniformScalingMap) = FiveArg + # properties Base.size(A::UniformScalingMap) = (A.M, A.M) Base.isreal(A::UniformScalingMap) = isreal(A.λ) @@ -91,11 +93,22 @@ function _scaling!(y, J::UniformScalingMap, x, α::Number=true, β::Number=false rmul!(y, β) return y else # α != 0, 1 - iszero(β) && (y .= λ .* x .* α; return y) - isone(β) && (y .+= λ .* x .* α; return y) - # β != 0, 1 - y .= y .* β .+ λ .* x .* α - return y + if iszero(β) + iszero(λ) && return fill!(y, zero(eltype(y))) + isone(λ) && return y .= x .* α + y .= λ .* x .* α + return y + elseif isone(β) + iszero(λ) && return y + isone(λ) && return y .+= x .* α + y .+= λ .* x .* α + return y + else # β != 0, 1 + iszero(λ) && (rmul!(y, β); return y) + isone(λ) && (y .= y .* β .+ x .* α; return y) + y .= y .* β .+ λ .* x .* α + return y + end end end diff --git a/src/wrappedmap.jl b/src/wrappedmap.jl index 8ebabc9e..09c94a8b 100644 --- a/src/wrappedmap.jl +++ b/src/wrappedmap.jl @@ -19,6 +19,8 @@ end const MatrixMap{T} = WrappedMap{T,<:AbstractMatrix} +mulstyle(A::WrappedMap) = mulstyle(A.lmap) + LinearAlgebra.transpose(A::MatrixMap{T}) where {T} = WrappedMap{T}(transpose(A.lmap); issymmetric=A._issymmetric, ishermitian=A._ishermitian, isposdef=A._isposdef) LinearAlgebra.adjoint(A::MatrixMap{T}) where {T} = diff --git a/test/blockmap.jl b/test/blockmap.jl index 43947c54..29b0d24d 100644 --- a/test/blockmap.jl +++ b/test/blockmap.jl @@ -6,6 +6,7 @@ using Test, LinearMaps, LinearAlgebra A11 = rand(elty, 10, 10) A12 = rand(elty, 10, n2) L = @inferred hcat(LinearMap(A11), LinearMap(A12)) + @test @inferred(LinearMaps.mulstyle(L)) == matrixstyle @test L isa LinearMaps.BlockMap{elty} A = [A11 A12] x = rand(10+n2) @@ -36,6 +37,7 @@ using Test, LinearMaps, LinearAlgebra A21 = rand(elty, 20, 10) L = @inferred vcat(LinearMap(A11), LinearMap(A21)) @test L isa LinearMaps.BlockMap{elty} + @test @inferred(LinearMaps.mulstyle(L)) == matrixstyle A = [A11; A21] x = rand(10) @test size(L) == size(A) @@ -62,6 +64,7 @@ using Test, LinearMaps, LinearAlgebra A = [A11 A12; A21 A22] @inferred hvcat((2,2), LinearMap(A11), LinearMap(A12), LinearMap(A21), LinearMap(A22)) L = [LinearMap(A11) LinearMap(A12); LinearMap(A21) LinearMap(A22)] + @test @inferred(LinearMaps.mulstyle(L)) == matrixstyle @test @inferred !issymmetric(L) @test @inferred !ishermitian(L) x = rand(30) @@ -102,12 +105,13 @@ using Test, LinearMaps, LinearAlgebra @test Matrix(adjoint(B)) == C' end end - + @testset "adjoint/transpose" begin for elty in (Float32, Float64, ComplexF64), transform in (transpose, adjoint) A12 = rand(elty, 10, 10) A = [I A12; transform(A12) I] L = [I LinearMap(A12); transform(LinearMap(A12)) I] + @test @inferred(LinearMaps.mulstyle(L)) == matrixstyle if elty <: Complex if transform == transpose @test @inferred issymmetric(L) diff --git a/test/kronecker.jl b/test/kronecker.jl index 564c5350..6f8e8605 100644 --- a/test/kronecker.jl +++ b/test/kronecker.jl @@ -9,6 +9,7 @@ using Test, LinearMaps, LinearAlgebra LB = LinearMap(B) LK = @inferred kron(LA, LB) @test @inferred size(LK) == size(K) + @test LinearMaps.mulstyle(LK) == LinearMaps.ThreeArg for i in (1, 2) @test @inferred size(LK, i) == size(K, i) end @@ -31,6 +32,11 @@ using Test, LinearMaps, LinearAlgebra @test @inferred kron(LA, LB)' == @inferred kron(LA', LB') @test (@inferred kron(LA, B)) == (@inferred kron(LA, LB)) == (@inferred kron(A, LB)) @test @inferred ishermitian(kron(LA'LA, LB'LB)) + A = rand(2, 5); B = rand(4, 2) + K = @inferred kron(A, LinearMap(B)) + @test Matrix(K) ≈ kron(A, B) + K = @inferred kron(LinearMap(B), A) + @test Matrix(K) ≈ kron(B, A) A = rand(3, 3); B = rand(2, 2); LA = LinearMap(A); LB = LinearMap(B) @test @inferred issymmetric(kron(LA'LA, LB'LB)) @test @inferred ishermitian(kron(LA'LA, LB'LB)) @@ -59,7 +65,7 @@ using Test, LinearMaps, LinearAlgebra @test Matrix(kronsum(transform(LA), transform(LB))) ≈ transform(KSmat) @test Matrix(transform(LinearMap(kronsum(LA, LB)))) ≈ Matrix(transform(KS)) ≈ transform(KSmat) end - @inferred kronsum(A, A, LB) + @test @inferred(kronsum(A, A, LB)) == @inferred(⊕(A, A, B)) @test Matrix(@inferred LA^⊕(3)) == Matrix(@inferred A^⊕(3)) ≈ Matrix(kronsum(LA, A, A)) @test @inferred(kronsum(LA, LA, LB)) == @inferred(kronsum(LA, kronsum(LA, LB))) == @inferred(kronsum(A, A, B)) @test Matrix(@inferred kronsum(A, B, A, B, A, B)) ≈ Matrix(@inferred kronsum(LA, LB, LA, LB, LA, LB)) diff --git a/test/linearcombination.jl b/test/linearcombination.jl index d3f0aaca..4da621f0 100644 --- a/test/linearcombination.jl +++ b/test/linearcombination.jl @@ -4,7 +4,7 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools CS! = LinearMap{ComplexF64}(cumsum!, (y, x) -> (copyto!(y, x); reverse!(y); cumsum!(y, y)), 10; ismutating=true) - v = rand(10) + v = rand(ComplexF64, 10) u = similar(v) b = @benchmarkable mul!($u, $CS!, $v) @test run(b, samples=3).allocs == 0 @@ -13,12 +13,20 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools @test mul!(u, L, v) ≈ n * cumsum(v) b = @benchmarkable mul!($u, $L, $v) @test run(b, samples=5).allocs <= 1 + for α in (false, true, rand(ComplexF64)), β in (false, true, rand(ComplexF64)) + @test mul!(copy(u), L, v, α, β) ≈ Matrix(L)*v*α + u*β + end A = 2 * rand(ComplexF64, (10, 10)) .- 1 B = rand(ComplexF64, size(A)...) M = @inferred LinearMap(A) N = @inferred LinearMap(B) + @test @inferred(LinearMaps.mulstyle(M)) == matrixstyle + @test @inferred(LinearMaps.mulstyle(N)) == matrixstyle LC = @inferred M + N + @test @inferred(LinearMaps.mulstyle(LC)) == matrixstyle + @test @inferred(LinearMaps.mulstyle(LC + I)) == matrixstyle + @test @inferred(LinearMaps.mulstyle(LC + 2.0*I)) == matrixstyle v = rand(ComplexF64, 10) w = similar(v) b = @benchmarkable mul!($w, $M, $v) @@ -27,10 +35,14 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools b = @benchmarkable mul!($w, $LC, $v) @test run(b, samples=3).allocs == 0 for α in (false, true, rand(ComplexF64)), β in (false, true, rand(ComplexF64)) - b = @benchmarkable mul!($w, $LC, $v, $α, $β) - @test run(b, samples=3).allocs == 0 - b = @benchmarkable mul!($w, $(LC + I), $v, $α, $β) - @test run(b, samples=3).allocs == 0 + if testallocs + b = @benchmarkable mul!($w, $LC, $v, $α, $β) + @test run(b, samples=3).allocs == 0 + b = @benchmarkable mul!($w, $(I + LC), $v, $α, $β) + @test run(b, samples=3).allocs == 0 + b = @benchmarkable mul!($w, $(LC + I), $v, $α, $β) + @test run(b, samples=3).allocs == 0 + end y = rand(ComplexF64, size(v)) @test mul!(copy(y), LC, v, α, β) ≈ Matrix(LC)*v*α + y*β @test mul!(copy(y), LC+I, v, α, β) ≈ Matrix(LC + I)*v*α + y*β diff --git a/test/runtests.jl b/test/runtests.jl index 6851f7b6..31848937 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,9 @@ using Test, LinearMaps +import LinearMaps: FiveArg, ThreeArg + +const matrixstyle = VERSION ≥ v"1.3.0-alpha.115" ? FiveArg : ThreeArg + +const testallocs = false include("linearmaps.jl") diff --git a/test/transpose.jl b/test/transpose.jl index 784fefd7..e1736f14 100644 --- a/test/transpose.jl +++ b/test/transpose.jl @@ -53,6 +53,7 @@ using Test, LinearMaps, LinearAlgebra CS = @inferred LinearMap{ComplexF64}(cumsum, x -> reverse(cumsum(reverse(x))), 10; ismutating=false) for transform in (adjoint, transpose) @test transform(transform(CS)) == CS + @test LinearMaps.mulstyle(transform(CS)) == LinearMaps.mulstyle(CS) end @test transpose(CS) != adjoint(CS) @test adjoint(CS) != transpose(CS) diff --git a/test/uniformscalingmap.jl b/test/uniformscalingmap.jl index ab38ce88..0d09392b 100644 --- a/test/uniformscalingmap.jl +++ b/test/uniformscalingmap.jl @@ -29,8 +29,10 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools Λ = @inferred LinearMap(λ*I, 10) x = rand(10) y = rand(10) - b = @benchmarkable mul!($y, $Λ, $x, $α, $β) - @test run(b, samples=3).allocs == 0 + if testallocs + b = @benchmarkable mul!($y, $Λ, $x, $α, $β) + @test run(b, samples=3).allocs == 0 + end y = deepcopy(x) @inferred mul!(y, Λ, x, α, β) @test y ≈ λ * x * α + β * x From 72abd11b9ed2042af540c32f09736dde27098bd9 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 29 Nov 2019 10:29:00 +0100 Subject: [PATCH 2/3] edits as per review, recursive lincomb multiplication --- src/LinearMaps.jl | 18 ++++---- src/blockmap.jl | 2 +- src/linearcombination.jl | 92 ++++++++++++++++++--------------------- src/transpose.jl | 2 +- src/uniformscalingmap.jl | 2 +- src/wrappedmap.jl | 2 +- test/blockmap.jl | 8 ++-- test/kronecker.jl | 2 +- test/linearcombination.jl | 10 ++--- test/runtests.jl | 4 +- test/transpose.jl | 2 +- 11 files changed, 68 insertions(+), 76 deletions(-) diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index 549e773d..99e2b4fa 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -24,17 +24,17 @@ abstract type MulStyle end struct FiveArg <: MulStyle end struct ThreeArg <: MulStyle end -mulstyle(::Type{FiveArg}, ::Type{FiveArg}) = FiveArg -mulstyle(::Type{ThreeArg}, ::Type{FiveArg}) = ThreeArg -mulstyle(::Type{FiveArg}, ::Type{ThreeArg}) = ThreeArg -mulstyle(::Type{ThreeArg}, ::Type{ThreeArg}) = ThreeArg -mulstyle(::LinearMap) = ThreeArg # default -if VERSION ≥ v"1.3.0-alpha.115" - mulstyle(::AbstractMatrix) = FiveArg +MulStyle(::FiveArg, ::FiveArg) = FiveArg() +MulStyle(::ThreeArg, ::FiveArg) = ThreeArg() +MulStyle(::FiveArg, ::ThreeArg) = ThreeArg() +MulStyle(::ThreeArg, ::ThreeArg) = ThreeArg() +MulStyle(::LinearMap) = ThreeArg() # default +@static if VERSION ≥ v"1.3.0-alpha.115" + MulStyle(::AbstractMatrix) = FiveArg() else - mulstyle(::AbstractMatrix) = ThreeArg + MulStyle(::AbstractMatrix) = ThreeArg() end -mulstyle(A::LinearMap, As::LinearMap...) = mulstyle(mulstyle(A), mulstyle(As...)) +MulStyle(A::LinearMap, As::LinearMap...) = MulStyle(MulStyle(A), MulStyle(As...)) Base.isreal(A::LinearMap) = eltype(A) <: Real LinearAlgebra.issymmetric(::LinearMap) = false # default assumptions diff --git a/src/blockmap.jl b/src/blockmap.jl index c62152f0..d7acf1d9 100644 --- a/src/blockmap.jl +++ b/src/blockmap.jl @@ -14,7 +14,7 @@ end BlockMap{T}(maps::As, rows::S) where {T,As<:Tuple{Vararg{LinearMap}},S} = BlockMap{T,As,S}(maps, rows) -mulstyle(A::BlockMap) = mulstyle(A.maps...) +MulStyle(A::BlockMap) = MulStyle(A.maps...) function check_dim(A::LinearMap, dim, n) n == size(A, dim) || throw(DimensionMismatch("Expected $n, got $(size(A, dim))")) diff --git a/src/linearcombination.jl b/src/linearcombination.jl index 50f119ea..970d5df6 100644 --- a/src/linearcombination.jl +++ b/src/linearcombination.jl @@ -1,20 +1,20 @@ -struct LinearCombination{T, MS<:MulStyle, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T} +struct LinearCombination{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T} maps::As - function LinearCombination{T, MS, As}(maps::As) where {T, MS<:MulStyle, As} + function LinearCombination{T, As}(maps::As) where {T, As} N = length(maps) sz = size(maps[1]) for Ai in maps size(Ai) == sz || throw(DimensionMismatch("LinearCombination")) promote_type(T, eltype(Ai)) == T || throw(InexactError()) end - MS === FiveArg && mulstyle(maps...) === ThreeArg && throw("wrong mulstyle in constructor") - new{T, MS, As}(maps) + new{T, As}(maps) end end -LinearCombination{T,MS}(maps::As) where {T, MS<:MulStyle, As} = LinearCombination{T, mulstyle(maps...), As}(maps) +LinearCombination{T}(maps::As) where {T, As} = + LinearCombination{T, As}(maps) -mulstyle(::LinearCombination{T,MS}) where {T, MS<:MulStyle} = MS +MulStyle(A::LinearCombination) = MulStyle(A.maps...) # basic methods Base.size(A::LinearCombination) = size(A.maps[1]) @@ -42,18 +42,18 @@ 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, mulstyle(A₁, A₂)}(tuple(A₁, A₂)) + 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, mulstyle(A₁, A₂)}(tuple(A₁, A₂.maps...)) + return LinearCombination{T}(tuple(A₁, A₂.maps...)) end Base.:(+)(A₁::LinearCombination, A₂::LinearMap) = +(A₂, A₁) function Base.:(+)(A₁::LinearCombination, A₂::LinearCombination) size(A₁) == size(A₂) || throw(DimensionMismatch("+")) T = promote_type(eltype(A₁), eltype(A₂)) - return LinearCombination{T, mulstyle(A₁, A₂)}(tuple(A₁.maps..., A₂.maps...)) + return LinearCombination{T}(tuple(A₁.maps..., A₂.maps...)) end Base.:(-)(A₁::LinearMap, A₂::LinearMap) = +(A₁, -A₂) @@ -61,59 +61,51 @@ Base.:(-)(A₁::LinearMap, A₂::LinearMap) = +(A₁, -A₂) Base.:(==)(A::LinearCombination, B::LinearCombination) = (eltype(A) == eltype(B) && A.maps == B.maps) # special transposition behavior -LinearAlgebra.transpose(A::LinearCombination) = LinearCombination{eltype(A), mulstyle(A)}(map(transpose, A.maps)) -LinearAlgebra.adjoint(A::LinearCombination) = LinearCombination{eltype(A), mulstyle(A)}(map(adjoint, A.maps)) +LinearAlgebra.transpose(A::LinearCombination) = + LinearCombination{eltype(A)}(map(transpose, A.maps)) +LinearAlgebra.adjoint(A::LinearCombination) = + LinearCombination{eltype(A)}(map(adjoint, A.maps)) -# multiplication with vectors +# multiplication with vectors & matrices for Atype in (AbstractVector, AbstractMatrix) @eval Base.@propagate_inbounds function LinearAlgebra.mul!(y::$Atype, A::LinearCombination, x::$Atype, α::Number=true, β::Number=false) @boundscheck check_dim_mul(y, A, x) - return _lincombmul!(y, A, x, α, β) - end -end - -@inline function _lincombmul!(y, A::LinearCombination{<:Any,FiveArg}, x, α::Number, β::Number) - if iszero(α) # trivial cases - iszero(β) && (fill!(y, zero(eltype(y))); return y) - isone(β) && return y - # β != 0, 1 - rmul!(y, β) - return y - else - mul!(y, first(A.maps), x, α, β) - @inbounds for An in Base.tail(A.maps) - mul!(y, An, x, α, true) + if iszero(α) # trivial cases + iszero(β) && (fill!(y, zero(eltype(y))); return y) + isone(β) && return y + # β != 0, 1 + rmul!(y, β) + return y + else + mul!(y, first(A.maps), x, α, β) + return _mul!(MulStyle(A), y, A, x, α, β) end - return y end end -@inline function _lincombmul!(y, A::LinearCombination{<:Any,ThreeArg}, x, α::Number, β::Number) - if iszero(α) - iszero(β) && (fill!(y, zero(eltype(y))); return y) - isone(β) && return y - # β != 0, 1 - rmul!(y, β) - return y - else - mul!(y, first(A.maps), x, α, β) - l = length(A.maps) - if l>1 - z = similar(y) - @inbounds for n in 2:l - An = A.maps[n] - muladd!(mulstyle(An), y, An, x, α, z) - end - end - return y - end +@inline _mul!(::FiveArg, y, A::LinearCombination, x, α::Number, β::Number) = + __mul!(y, Base.tail(A.maps), x, α, nothing) +@inline function _mul!(::ThreeArg, y, A::LinearCombination, x, α::Number, β::Number) + z = similar(y) + __mul!(y, Base.tail(A.maps), x, α, z) end -@inline muladd!(::Type{FiveArg}, y, A, x, α, _) = mul!(y, A, x, α, true) -@inline function muladd!(::Type{ThreeArg}, y, A, x, α, z) +@inline __mul!(y, As::Tuple{Vararg{LinearMap}}, x, α, z) = + __mul!(__mul!(y, first(As), x, α, z), Base.tail(As), x, α, z) +@inline __mul!(y, A::Tuple{LinearMap}, x, α, z) = __mul!(y, first(A), x, α, z) +@inline __mul!(y, A::LinearMap, x, α, z) = muladd!(MulStyle(A), y, A, x, α, z) + +@inline muladd!(::FiveArg, y, A, x, α, _) = mul!(y, A, x, α, true) +@inline function muladd!(::ThreeArg, y, A, x, α, z) + # TODO: replace by mul!(z, A, x) A_mul_B!(z, A, x) - y .+= isone(α) ? z : z .* α + if isone(α) + y .+= z + else + y .+= z .* α + end + return y end A_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = mul!(y, A, x) diff --git a/src/transpose.jl b/src/transpose.jl index 9d158466..60caf3c2 100644 --- a/src/transpose.jl +++ b/src/transpose.jl @@ -5,7 +5,7 @@ struct AdjointMap{T, A<:LinearMap{T}} <: LinearMap{T} lmap::A end -mulstyle(A::Union{TransposeMap,AdjointMap}) = mulstyle(A.lmap) +MulStyle(A::Union{TransposeMap,AdjointMap}) = MulStyle(A.lmap) # transposition behavior of LinearMap objects LinearAlgebra.transpose(A::TransposeMap) = A.lmap diff --git a/src/uniformscalingmap.jl b/src/uniformscalingmap.jl index 4bb77f3c..3a8cb0c6 100644 --- a/src/uniformscalingmap.jl +++ b/src/uniformscalingmap.jl @@ -11,7 +11,7 @@ UniformScalingMap(λ::Number, M::Int, N::Int) = UniformScalingMap(λ::T, sz::Dims{2}) where {T} = (sz[1] == sz[2] ? UniformScalingMap(λ, sz[1]) : error("UniformScalingMap needs to be square")) -mulstyle(::UniformScalingMap) = FiveArg +MulStyle(::UniformScalingMap) = FiveArg() # properties Base.size(A::UniformScalingMap) = (A.M, A.M) diff --git a/src/wrappedmap.jl b/src/wrappedmap.jl index 09c94a8b..212c0de7 100644 --- a/src/wrappedmap.jl +++ b/src/wrappedmap.jl @@ -19,7 +19,7 @@ end const MatrixMap{T} = WrappedMap{T,<:AbstractMatrix} -mulstyle(A::WrappedMap) = mulstyle(A.lmap) +MulStyle(A::WrappedMap) = MulStyle(A.lmap) LinearAlgebra.transpose(A::MatrixMap{T}) where {T} = WrappedMap{T}(transpose(A.lmap); issymmetric=A._issymmetric, ishermitian=A._ishermitian, isposdef=A._isposdef) diff --git a/test/blockmap.jl b/test/blockmap.jl index 29b0d24d..734a9edc 100644 --- a/test/blockmap.jl +++ b/test/blockmap.jl @@ -6,7 +6,7 @@ using Test, LinearMaps, LinearAlgebra A11 = rand(elty, 10, 10) A12 = rand(elty, 10, n2) L = @inferred hcat(LinearMap(A11), LinearMap(A12)) - @test @inferred(LinearMaps.mulstyle(L)) == matrixstyle + @test @inferred(LinearMaps.MulStyle(L)) === matrixstyle @test L isa LinearMaps.BlockMap{elty} A = [A11 A12] x = rand(10+n2) @@ -37,7 +37,7 @@ using Test, LinearMaps, LinearAlgebra A21 = rand(elty, 20, 10) L = @inferred vcat(LinearMap(A11), LinearMap(A21)) @test L isa LinearMaps.BlockMap{elty} - @test @inferred(LinearMaps.mulstyle(L)) == matrixstyle + @test @inferred(LinearMaps.MulStyle(L)) === matrixstyle A = [A11; A21] x = rand(10) @test size(L) == size(A) @@ -64,7 +64,7 @@ using Test, LinearMaps, LinearAlgebra A = [A11 A12; A21 A22] @inferred hvcat((2,2), LinearMap(A11), LinearMap(A12), LinearMap(A21), LinearMap(A22)) L = [LinearMap(A11) LinearMap(A12); LinearMap(A21) LinearMap(A22)] - @test @inferred(LinearMaps.mulstyle(L)) == matrixstyle + @test @inferred(LinearMaps.MulStyle(L)) === matrixstyle @test @inferred !issymmetric(L) @test @inferred !ishermitian(L) x = rand(30) @@ -111,7 +111,7 @@ using Test, LinearMaps, LinearAlgebra A12 = rand(elty, 10, 10) A = [I A12; transform(A12) I] L = [I LinearMap(A12); transform(LinearMap(A12)) I] - @test @inferred(LinearMaps.mulstyle(L)) == matrixstyle + @test @inferred(LinearMaps.MulStyle(L)) === matrixstyle if elty <: Complex if transform == transpose @test @inferred issymmetric(L) diff --git a/test/kronecker.jl b/test/kronecker.jl index 6f8e8605..69e0367e 100644 --- a/test/kronecker.jl +++ b/test/kronecker.jl @@ -9,7 +9,7 @@ using Test, LinearMaps, LinearAlgebra LB = LinearMap(B) LK = @inferred kron(LA, LB) @test @inferred size(LK) == size(K) - @test LinearMaps.mulstyle(LK) == LinearMaps.ThreeArg + @test LinearMaps.MulStyle(LK) === LinearMaps.ThreeArg() for i in (1, 2) @test @inferred size(LK, i) == size(K, i) end diff --git a/test/linearcombination.jl b/test/linearcombination.jl index 4da621f0..fd755010 100644 --- a/test/linearcombination.jl +++ b/test/linearcombination.jl @@ -21,12 +21,12 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools B = rand(ComplexF64, size(A)...) M = @inferred LinearMap(A) N = @inferred LinearMap(B) - @test @inferred(LinearMaps.mulstyle(M)) == matrixstyle - @test @inferred(LinearMaps.mulstyle(N)) == matrixstyle + @test @inferred(LinearMaps.MulStyle(M)) === matrixstyle + @test @inferred(LinearMaps.MulStyle(N)) === matrixstyle LC = @inferred M + N - @test @inferred(LinearMaps.mulstyle(LC)) == matrixstyle - @test @inferred(LinearMaps.mulstyle(LC + I)) == matrixstyle - @test @inferred(LinearMaps.mulstyle(LC + 2.0*I)) == matrixstyle + @test @inferred(LinearMaps.MulStyle(LC)) === matrixstyle + @test @inferred(LinearMaps.MulStyle(LC + I)) === matrixstyle + @test @inferred(LinearMaps.MulStyle(LC + 2.0*I)) === matrixstyle v = rand(ComplexF64, 10) w = similar(v) b = @benchmarkable mul!($w, $M, $v) diff --git a/test/runtests.jl b/test/runtests.jl index 31848937..b5444f0c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,9 +1,9 @@ using Test, LinearMaps import LinearMaps: FiveArg, ThreeArg -const matrixstyle = VERSION ≥ v"1.3.0-alpha.115" ? FiveArg : ThreeArg +const matrixstyle = VERSION ≥ v"1.3.0-alpha.115" ? FiveArg() : ThreeArg() -const testallocs = false +const testallocs = true include("linearmaps.jl") diff --git a/test/transpose.jl b/test/transpose.jl index e1736f14..ad71cccf 100644 --- a/test/transpose.jl +++ b/test/transpose.jl @@ -53,7 +53,7 @@ using Test, LinearMaps, LinearAlgebra CS = @inferred LinearMap{ComplexF64}(cumsum, x -> reverse(cumsum(reverse(x))), 10; ismutating=false) for transform in (adjoint, transpose) @test transform(transform(CS)) == CS - @test LinearMaps.mulstyle(transform(CS)) == LinearMaps.mulstyle(CS) + @test LinearMaps.MulStyle(transform(CS)) === LinearMaps.MulStyle(CS) end @test transpose(CS) != adjoint(CS) @test adjoint(CS) != transpose(CS) From 7cd8c9810d570634050cb9d39e95d1af200fc8e5 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 29 Nov 2019 10:47:10 +0100 Subject: [PATCH 3/3] revert some minor layout modifications --- src/linearcombination.jl | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/src/linearcombination.jl b/src/linearcombination.jl index 970d5df6..e9e7696b 100644 --- a/src/linearcombination.jl +++ b/src/linearcombination.jl @@ -3,16 +3,15 @@ struct LinearCombination{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T} function LinearCombination{T, As}(maps::As) where {T, As} N = length(maps) sz = size(maps[1]) - for Ai in maps - size(Ai) == sz || throw(DimensionMismatch("LinearCombination")) - promote_type(T, eltype(Ai)) == T || throw(InexactError()) + for n in 1:N + size(maps[n]) == sz || throw(DimensionMismatch("LinearCombination")) + promote_type(T, eltype(maps[n])) == T || throw(InexactError()) end new{T, As}(maps) end end -LinearCombination{T}(maps::As) where {T, As} = - LinearCombination{T, As}(maps) +LinearCombination{T}(maps::As) where {T, As} = LinearCombination{T, As}(maps) MulStyle(A::LinearCombination) = MulStyle(A.maps...) @@ -61,10 +60,8 @@ Base.:(-)(A₁::LinearMap, A₂::LinearMap) = +(A₁, -A₂) Base.:(==)(A::LinearCombination, B::LinearCombination) = (eltype(A) == eltype(B) && A.maps == B.maps) # special transposition behavior -LinearAlgebra.transpose(A::LinearCombination) = - LinearCombination{eltype(A)}(map(transpose, A.maps)) -LinearAlgebra.adjoint(A::LinearCombination) = - LinearCombination{eltype(A)}(map(adjoint, A.maps)) +LinearAlgebra.transpose(A::LinearCombination) = LinearCombination{eltype(A)}(map(transpose, A.maps)) +LinearAlgebra.adjoint(A::LinearCombination) = LinearCombination{eltype(A)}(map(adjoint, A.maps)) # multiplication with vectors & matrices for Atype in (AbstractVector, AbstractMatrix)