From 3feab935be0222d1514a5134af40cff26680cde1 Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Wed, 29 Jul 2020 14:02:02 -0400 Subject: [PATCH 01/18] Left vector multiply --- src/wrappedmap.jl | 7 +++++++ test/wrappedmap.jl | 30 ++++++++++++++++++++++++++++++ 2 files changed, 37 insertions(+) diff --git a/src/wrappedmap.jl b/src/wrappedmap.jl index 574827ea..aa29a7c7 100644 --- a/src/wrappedmap.jl +++ b/src/wrappedmap.jl @@ -1,3 +1,5 @@ +import LinearAlgebra: AdjointAbsVec, TransposeAbsVec + struct WrappedMap{T, A<:MapOrMatrix} <: LinearMap{T} lmap::A _issymmetric::Bool @@ -71,3 +73,8 @@ Base.:(-)(A₁::AbstractMatrix, A₂::LinearMap) = -(WrappedMap(A₁), A₂) Base.:(*)(A₁::LinearMap, A₂::AbstractMatrix) = *(A₁, WrappedMap(A₂)) Base.:(*)(A₁::AbstractMatrix, A₂::LinearMap) = *(WrappedMap(A₁), A₂) + +# An AdjointAbsVec isa AbstractMatrix but we handle it like a vector (Issue#99) +# which is an exception to the left multiplication rule that makes a WrappedMap +Base.:(*)(y::AdjointAbsVec, A::LinearMap) = adjoint(*(A', y')) +Base.:(*)(y::TransposeAbsVec, A::LinearMap) = transpose(transpose(A) * transpose(y)) diff --git a/test/wrappedmap.jl b/test/wrappedmap.jl index 57171b25..d07b13eb 100644 --- a/test/wrappedmap.jl +++ b/test/wrappedmap.jl @@ -15,3 +15,33 @@ using Test, LinearMaps, LinearAlgebra @test @inferred isposdef(MA) @test @inferred isposdef(MB) end + +# y'*L is an exception to the left multiplication rule that makes a WrappedMap +@testset "left mul vec" begin + T = ComplexF32 + A = rand(T, 4, 5) + x = rand(T, 5) + y = rand(T, 4) + L = LinearMap{T}(A) + @test y'A ≈ y'L + @test (y' * A) * x ≈ y' * (L * x) # associative + @test transpose(y) * A ≈ transpose(y) * L + + LL = L * [L' L'] # stress test with CompositeMap & BlockMap + AA = A * [A' A'] + x = rand(T, 8) + @test y'AA ≈ y'LL + @test (y' * AA) * x ≈ y' * (LL * x) # associative + @test transpose(y)*AA ≈ transpose(y)*LL + + # mul! versions + b1 = y'*L + b2 = similar(b1) + mul!(b2, y', A) + @test b1 ≈ b2 + + b1 = transpose(y)*L + b2 = similar(b1) + mul!(b2, transpose(y), A) + @test b1 ≈ b2 +end From e45289d60c2bf14e8b47c2a94ff05a3a708e9c7b Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Wed, 29 Jul 2020 14:02:43 -0400 Subject: [PATCH 02/18] Update README and bump version --- Project.toml | 2 +- README.md | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 7ed11620..f1af504e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "LinearMaps" uuid = "7a12625a-238d-50fd-b39a-03d52299707e" -version = "2.7.0" +version = "2.8.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/README.md b/README.md index 155b33bc..24f0cc5e 100644 --- a/README.md +++ b/README.md @@ -8,6 +8,10 @@ A Julia package for defining and working with linear maps, also known as linear transformations or linear operators acting on vectors. The only requirement for a LinearMap is that it can act on a vector (by multiplication) efficiently. +## What's new in v2.8 +* Left multiplying by a transpose or adjoint vector (e.g., `y'*A`) + produces a transpose or adjoint vector output, rather than a `LinearMap`. + ## What's new in v2.7 * Potential reduction of memory allocations in multiplication of `LinearCombination`s, `BlockMap`s, and real- or complex-scaled `LinearMap`s. From 57d78c4c230d34229a138899e57d66fe66b4ac40 Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Wed, 29 Jul 2020 23:49:31 -0400 Subject: [PATCH 03/18] Support mul! for left vectors --- src/LinearMaps.jl | 25 ++++++++++++----- src/blockmap.jl | 16 +++++------ src/composition.jl | 14 +++++----- src/functionmap.jl | 6 ++--- src/kronecker.jl | 18 ++++++------- src/left.jl | 35 ++++++++++++++++++++++++ src/linearcombination.jl | 6 ++--- src/scaledmap.jl | 8 +++--- src/transpose.jl | 12 ++++----- src/wrappedmap.jl | 13 +++------ test/left.jl | 58 ++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 2 ++ test/wrappedmap.jl | 30 --------------------- 13 files changed, 157 insertions(+), 86 deletions(-) create mode 100644 src/left.jl create mode 100644 test/left.jl diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index 30419ca5..8be3546f 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -4,6 +4,7 @@ export LinearMap export ⊗, kronsum, ⊕ using LinearAlgebra +import LinearAlgebra: AdjointAbsVec, TransposeAbsVec using SparseArrays if VERSION < v"1.2-" @@ -18,6 +19,10 @@ abstract type LinearMap{T} end const MapOrMatrix{T} = Union{LinearMap{T},AbstractMatrix{T}} const RealOrComplex = Union{Real,Complex} +# valid types for left vector multiplication: +const VecIn = Union{AdjointAbsVec, TransposeAbsVec} +const VecOut = Union{AbstractVector, AdjointAbsVec, TransposeAbsVec} + Base.eltype(::LinearMap{T}) where {T} = T abstract type MulStyle end @@ -48,7 +53,7 @@ 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] # check dimension consistency for y = A*x and Y = A*X -function check_dim_mul(y::AbstractVector, A::LinearMap, x::AbstractVector) +function check_dim_mul(y::VecOut, A::LinearMap, x::AbstractVector) # @info "checked vector dimensions" # uncomment for testing m, n = size(A) (m == length(y) && n == length(x)) || throw(DimensionMismatch("mul!")) @@ -61,6 +66,13 @@ function check_dim_mul(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix) return nothing end +# check dimension consistency for left multiplication x = y'*A +function check_dim_mul(x::V, y::V, A::LinearMap) where {V <: VecIn} + m, n = size(A) + ((1,m) == size(y) && (1,n) == size(x)) || throw(DimensionMismatch("left mul!")) + return nothing +end + # conversion of AbstractMatrix to LinearMap convert_to_lmaps_(A::AbstractMatrix) = LinearMap(A) convert_to_lmaps_(A::LinearMap) = A @@ -73,11 +85,11 @@ function Base.:(*)(A::LinearMap, x::AbstractVector) size(A, 2) == length(x) || throw(DimensionMismatch("mul!")) return @inbounds A_mul_B!(similar(x, promote_type(eltype(A), eltype(x)), size(A, 1)), A, x) end -function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap, x::AbstractVector) +function LinearAlgebra.mul!(y::VecOut, A::LinearMap, x::AbstractVector) @boundscheck check_dim_mul(y, A, x) return @inbounds A_mul_B!(y, A, x) end -function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap, x::AbstractVector, α::Number, β::Number) +function LinearAlgebra.mul!(y::VecOut, A::LinearMap, x::AbstractVector, α::Number, β::Number) @boundscheck check_dim_mul(y, A, x) if isone(α) iszero(β) && (A_mul_B!(y, A, x); return y) @@ -114,10 +126,11 @@ Base.@propagate_inbounds function LinearAlgebra.mul!(Y::AbstractMatrix, A::Linea return Y end -A_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, A, x) -At_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, transpose(A), x) -Ac_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, adjoint(A), x) +A_mul_B!(y::VecOut, A::AbstractMatrix, x::AbstractVector) = mul!(y, A, x) +At_mul_B!(y::VecOut, A::AbstractMatrix, x::AbstractVector) = mul!(y, transpose(A), x) +Ac_mul_B!(y::VecOut, A::AbstractMatrix, x::AbstractVector) = mul!(y, adjoint(A), x) +include("left.jl") # left multiplication by a transpose or adjoint vector include("wrappedmap.jl") # wrap a matrix of linear map in a new type, thereby allowing to alter its properties include("uniformscalingmap.jl") # the uniform scaling map, to be able to make linear combinations of LinearMap objects and multiples of I include("transpose.jl") # transposing linear maps diff --git a/src/blockmap.jl b/src/blockmap.jl index becc2de4..ff88e066 100644 --- a/src/blockmap.jl +++ b/src/blockmap.jl @@ -382,19 +382,19 @@ end # multiplication with vectors & matrices ############ -Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) = +Base.@propagate_inbounds A_mul_B!(y::VecOut, A::BlockMap, x::AbstractVector) = mul!(y, A, x, true, false) -Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::TransposeMap{<:Any,<:BlockMap}, x::AbstractVector) = +Base.@propagate_inbounds A_mul_B!(y::VecOut, A::TransposeMap{<:Any,<:BlockMap}, x::AbstractVector) = mul!(y, A, x, true, false) -Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) = +Base.@propagate_inbounds At_mul_B!(y::VecOut, A::BlockMap, x::AbstractVector) = mul!(y, transpose(A), x, true, false) -Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::AdjointMap{<:Any,<:BlockMap}, x::AbstractVector) = +Base.@propagate_inbounds A_mul_B!(y::VecOut, A::AdjointMap{<:Any,<:BlockMap}, x::AbstractVector) = mul!(y, A, x, true, false) -Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) = +Base.@propagate_inbounds Ac_mul_B!(y::VecOut, A::BlockMap, x::AbstractVector) = mul!(y, adjoint(A), x, true, false) for Atype in (AbstractVector, AbstractMatrix) @@ -495,13 +495,13 @@ LinearAlgebra.transpose(A::BlockDiagonalMap{T}) where {T} = BlockDiagonalMap{T}( Base.:(==)(A::BlockDiagonalMap, B::BlockDiagonalMap) = (eltype(A) == eltype(B) && A.maps == B.maps) -Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::BlockDiagonalMap, x::AbstractVector) = +Base.@propagate_inbounds A_mul_B!(y::VecOut, A::BlockDiagonalMap, x::AbstractVector) = mul!(y, A, x, true, false) -Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::BlockDiagonalMap, x::AbstractVector) = +Base.@propagate_inbounds At_mul_B!(y::VecOut, A::BlockDiagonalMap, x::AbstractVector) = mul!(y, transpose(A), x, true, false) -Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::BlockDiagonalMap, x::AbstractVector) = +Base.@propagate_inbounds Ac_mul_B!(y::VecOut, A::BlockDiagonalMap, x::AbstractVector) = mul!(y, adjoint(A), x, true, false) for Atype in (AbstractVector, AbstractMatrix) diff --git a/src/composition.jl b/src/composition.jl index b3ce93fc..9d91edee 100644 --- a/src/composition.jl +++ b/src/composition.jl @@ -129,23 +129,23 @@ LinearAlgebra.adjoint(A::CompositeMap{T}) where {T} = CompositeMap{T}(map(adjo Base.:(==)(A::CompositeMap, B::CompositeMap) = (eltype(A) == eltype(B) && A.maps == B.maps) # multiplication with vectors -function A_mul_B!(y::AbstractVector, A::CompositeMap{T,<:Tuple{LinearMap}}, x::AbstractVector) where {T} +function A_mul_B!(y::VecOut, A::CompositeMap{T,<:Tuple{LinearMap}}, x::AbstractVector) where {T} return A_mul_B!(y, A.maps[1], x) end -function A_mul_B!(y::AbstractVector, A::CompositeMap{T,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector) where {T} +function A_mul_B!(y::VecOut, A::CompositeMap{T,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector) where {T} _compositemul!(y, A, x, similar(y, size(A.maps[1], 1))) end -function A_mul_B!(y::AbstractVector, A::CompositeMap{T,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector) where {T} +function A_mul_B!(y::VecOut, A::CompositeMap{T,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector) where {T} _compositemul!(y, A, x, similar(y, size(A.maps[1], 1)), similar(y, size(A.maps[2], 1))) end -function _compositemul!(y::AbstractVector, A::CompositeMap{T,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector, z::AbstractVector) where {T} +function _compositemul!(y::VecOut, A::CompositeMap{T,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector, z::AbstractVector) where {T} # no size checking, will be done by individual maps A_mul_B!(z, A.maps[1], x) A_mul_B!(y, A.maps[2], z) return y end -function _compositemul!(y::AbstractVector, A::CompositeMap{T,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector, source::AbstractVector, dest::AbstractVector) where {T} +function _compositemul!(y::VecOut, A::CompositeMap{T,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector, source::AbstractVector, dest::AbstractVector) where {T} # no size checking, will be done by individual maps N = length(A.maps) A_mul_B!(source, A.maps[1], x) @@ -166,6 +166,6 @@ function _compositemul!(y::AbstractVector, A::CompositeMap{T,<:Tuple{Vararg{Line return y end -At_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) +At_mul_B!(y::VecOut, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) -Ac_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) +Ac_mul_B!(y::VecOut, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) diff --git a/src/functionmap.jl b/src/functionmap.jl index 4e1b086f..937dc86c 100644 --- a/src/functionmap.jl +++ b/src/functionmap.jl @@ -107,13 +107,13 @@ function Base.:(*)(A::TransposeMap{<:Any,<:FunctionMap}, x::AbstractVector) end end -function A_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector) +function A_mul_B!(y::VecOut, A::FunctionMap, x::AbstractVector) (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) +function At_mul_B!(y::VecOut, A::FunctionMap, x::AbstractVector) (issymmetric(A) || (isreal(A) && ishermitian(A))) && return A_mul_B!(y, A, x) (length(x) == A.M && length(y) == A.N) || throw(DimensionMismatch("At_mul_B!")) if A.fc !== nothing @@ -135,7 +135,7 @@ function At_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector) end end -function Ac_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector) +function Ac_mul_B!(y::VecOut, A::FunctionMap, x::AbstractVector) ishermitian(A) && return A_mul_B!(y, A, x) (length(x) == A.M && length(y) == A.N) || throw(DimensionMismatch("Ac_mul_B!")) if A.fc !== nothing diff --git a/src/kronecker.jl b/src/kronecker.jl index 80fd502a..be4db862 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -131,7 +131,7 @@ end # multiplication with vectors ################# -Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T} +Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T} require_one_based_indexing(y) @boundscheck check_dim_mul(y, L, x) A, B = L.maps @@ -139,7 +139,7 @@ Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T, _kronmul!(y, B, X, transpose(A), T) return y end -Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T}, x::AbstractVector) where {T} +Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::KroneckerMap{T}, x::AbstractVector) where {T} require_one_based_indexing(y) @boundscheck check_dim_mul(y, L, x) A = first(L.maps) @@ -150,7 +150,7 @@ Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T} end # mixed-product rule, prefer the right if possible # (A₁ ⊗ A₂ ⊗ ... ⊗ Aᵣ) * (B₁ ⊗ B₂ ⊗ ... ⊗ Bᵣ) = (A₁B₁) ⊗ (A₂B₂) ⊗ ... ⊗ (AᵣBᵣ) -Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{<:Any,<:Tuple{KroneckerMap,KroneckerMap}}, x::AbstractVector) +Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::CompositeMap{<:Any,<:Tuple{KroneckerMap,KroneckerMap}}, x::AbstractVector) B, A = L.maps if length(A.maps) == length(B.maps) && all(M -> check_dim_mul(M[1], M[2]), zip(A.maps, B.maps)) A_mul_B!(y, kron(map(*, A.maps, B.maps)...), x) @@ -160,7 +160,7 @@ Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{<: end # mixed-product rule, prefer the right if possible # (A₁ ⊗ B₁)*(A₂⊗B₂)*...*(Aᵣ⊗Bᵣ) = (A₁*A₂*...*Aᵣ) ⊗ (B₁*B₂*...*Bᵣ) -Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{T,<:Tuple{Vararg{KroneckerMap{<:Any,<:Tuple{LinearMap,LinearMap}}}}}, x::AbstractVector) where {T} +Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::CompositeMap{T,<:Tuple{Vararg{KroneckerMap{<:Any,<:Tuple{LinearMap,LinearMap}}}}}, x::AbstractVector) where {T} As = map(AB -> AB.maps[1], L.maps) Bs = map(AB -> AB.maps[2], L.maps) As1, As2 = Base.front(As), Base.tail(As) @@ -173,10 +173,10 @@ Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{T, end end -Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = +Base.@propagate_inbounds At_mul_B!(y::VecOut, A::KroneckerMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) -Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = +Base.@propagate_inbounds Ac_mul_B!(y::VecOut, A::KroneckerMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) ############### @@ -261,7 +261,7 @@ LinearAlgebra.transpose(A::KroneckerSumMap{T}) where {T} = KroneckerSumMap{T}(ma Base.:(==)(A::KroneckerSumMap, B::KroneckerSumMap) = (eltype(A) == eltype(B) && A.maps == B.maps) -Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerSumMap, x::AbstractVector) +Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::KroneckerSumMap, x::AbstractVector) @boundscheck check_dim_mul(y, L, x) A, B = L.maps ma, na = size(A) @@ -273,8 +273,8 @@ Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerSumMap return y end -Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = +Base.@propagate_inbounds At_mul_B!(y::VecOut, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) -Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = +Base.@propagate_inbounds Ac_mul_B!(y::VecOut, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) diff --git a/src/left.jl b/src/left.jl new file mode 100644 index 00000000..75c14378 --- /dev/null +++ b/src/left.jl @@ -0,0 +1,35 @@ +# left.jl +# "special cases" related left vector multiplication like x = y'*A +# The subtlety is that y' is a AdjointAbsVec or a TransposeAbsVec +# which is a subtype of AbstractMatrix but it is really "vector like" +# so we want handle it like a vector (Issue#99) +# So this is an exception to the left multiplication rule by a AbstractMatrix +# that usually makes a WrappedMap. +# The "transpose" versions may be of dubious use, but the adjoint versions +# are useful for ensuring that (y'A)*x ≈ y'*(A*x) are both scalars. + +import LinearAlgebra: AdjointAbsVec, TransposeAbsVec + +Base.:(*)(y::AdjointAbsVec, A::LinearMap) = adjoint(*(A', y')) +Base.:(*)(y::TransposeAbsVec, A::LinearMap) = transpose(transpose(A) * transpose(y)) + +# mul!(x, y', A) +LinearAlgebra.mul!(x::AdjointAbsVec, y::AdjointAbsVec, A::LinearMap) = + mul!(x, y, A, true, false) + +# not sure if we need bounds checks and propagate inbounds stuff here + +# mul!(x, y', A, α, β) +function LinearAlgebra.mul!(x::AdjointAbsVec, y::AdjointAbsVec, A::LinearMap, + α::Number, β::Number) + check_dim_mul(x, y, A) + mul!(x, A', y', α, β) + return conj!(x) +end + +# mul!(x, transpose(y), A, α, β) +function LinearAlgebra.mul!(x::TransposeAbsVec, y::TransposeAbsVec, A::LinearMap, + α::Number, β::Number) + check_dim_mul(x, y, A) + return mul!(x, transpose(A), transpose(y), α, β) +end diff --git a/src/linearcombination.jl b/src/linearcombination.jl index 153c861f..705c8bbe 100644 --- a/src/linearcombination.jl +++ b/src/linearcombination.jl @@ -117,8 +117,8 @@ end return y end -A_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = mul!(y, A, x, true, false) +A_mul_B!(y::VecOut, A::LinearCombination, x::AbstractVector) = mul!(y, A, x, true, false) -At_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = mul!(y, transpose(A), x, true, false) +At_mul_B!(y::VecOut, A::LinearCombination, x::AbstractVector) = mul!(y, transpose(A), x, true, false) -Ac_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = mul!(y, adjoint(A), x, true, false) +Ac_mul_B!(y::VecOut, A::LinearCombination, x::AbstractVector) = mul!(y, adjoint(A), x, true, false) diff --git a/src/scaledmap.jl b/src/scaledmap.jl index 14b914c1..dca651c3 100644 --- a/src/scaledmap.jl +++ b/src/scaledmap.jl @@ -54,15 +54,15 @@ Base.:(*)(A::ScaledMap, B::LinearMap) = A.λ * (A.lmap * B) Base.:(*)(A::LinearMap, B::ScaledMap) = (A * B.lmap) * B.λ # multiplication with vectors -function A_mul_B!(y::AbstractVector, A::ScaledMap, x::AbstractVector) +function A_mul_B!(y::VecOut, A::ScaledMap, x::AbstractVector) # no size checking, will be done by map mul!(y, A.lmap, x, A.λ, false) end -function LinearAlgebra.mul!(y::AbstractVector, A::ScaledMap, x::AbstractVector, α::Number, β::Number) +function LinearAlgebra.mul!(y::VecOut, A::ScaledMap, x::AbstractVector, α::Number, β::Number) # no size checking, will be done by map mul!(y, A.lmap, x, A.λ * α, β) end -At_mul_B!(y::AbstractVector, A::ScaledMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) -Ac_mul_B!(y::AbstractVector, A::ScaledMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) +At_mul_B!(y::VecOut, A::ScaledMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) +Ac_mul_B!(y::VecOut, A::ScaledMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) diff --git a/src/transpose.jl b/src/transpose.jl index 60caf3c2..0a60d6bf 100644 --- a/src/transpose.jl +++ b/src/transpose.jl @@ -32,18 +32,18 @@ Base.:(==)(A::LinearMap, B::TransposeMap) = issymmetric(A) && B.lmap == A Base.:(==)(A::LinearMap, B::AdjointMap) = ishermitian(A) && B.lmap == A # multiplication with vector -A_mul_B!(y::AbstractVector, A::TransposeMap, x::AbstractVector) = +A_mul_B!(y::VecOut, A::TransposeMap, x::AbstractVector) = issymmetric(A.lmap) ? A_mul_B!(y, A.lmap, x) : At_mul_B!(y, A.lmap, x) -At_mul_B!(y::AbstractVector, A::TransposeMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x) +At_mul_B!(y::VecOut, A::TransposeMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x) -Ac_mul_B!(y::AbstractVector, A::TransposeMap, x::AbstractVector) = +Ac_mul_B!(y::VecOut, A::TransposeMap, x::AbstractVector) = isreal(A.lmap) ? A_mul_B!(y, A.lmap, x) : (A_mul_B!(y, A.lmap, conj(x)); conj!(y)) -A_mul_B!(y::AbstractVector, A::AdjointMap, x::AbstractVector) = +A_mul_B!(y::VecOut, A::AdjointMap, x::AbstractVector) = ishermitian(A.lmap) ? A_mul_B!(y, A.lmap, x) : Ac_mul_B!(y, A.lmap, x) -At_mul_B!(y::AbstractVector, A::AdjointMap, x::AbstractVector) = +At_mul_B!(y::VecOut, A::AdjointMap, x::AbstractVector) = isreal(A.lmap) ? A_mul_B!(y, A.lmap, x) : (A_mul_B!(y, A.lmap, conj(x)); conj!(y)) -Ac_mul_B!(y::AbstractVector, A::AdjointMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x) +Ac_mul_B!(y::VecOut, A::AdjointMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x) diff --git a/src/wrappedmap.jl b/src/wrappedmap.jl index aa29a7c7..321380dc 100644 --- a/src/wrappedmap.jl +++ b/src/wrappedmap.jl @@ -1,5 +1,3 @@ -import LinearAlgebra: AdjointAbsVec, TransposeAbsVec - struct WrappedMap{T, A<:MapOrMatrix} <: LinearMap{T} lmap::A _issymmetric::Bool @@ -39,13 +37,13 @@ LinearAlgebra.ishermitian(A::WrappedMap) = A._ishermitian LinearAlgebra.isposdef(A::WrappedMap) = A._isposdef # multiplication with vectors & matrices -A_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x) +A_mul_B!(y::VecOut, A::WrappedMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x) Base.:(*)(A::WrappedMap, x::AbstractVector) = *(A.lmap, x) -At_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) = +At_mul_B!(y::VecOut, A::WrappedMap, x::AbstractVector) = (issymmetric(A) || (isreal(A) && ishermitian(A))) ? A_mul_B!(y, A.lmap, x) : At_mul_B!(y, A.lmap, x) -Ac_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) = +Ac_mul_B!(y::VecOut, A::WrappedMap, x::AbstractVector) = ishermitian(A) ? A_mul_B!(y, A.lmap, x) : Ac_mul_B!(y, A.lmap, x) if VERSION ≥ v"1.3.0-alpha.115" @@ -73,8 +71,3 @@ Base.:(-)(A₁::AbstractMatrix, A₂::LinearMap) = -(WrappedMap(A₁), A₂) Base.:(*)(A₁::LinearMap, A₂::AbstractMatrix) = *(A₁, WrappedMap(A₂)) Base.:(*)(A₁::AbstractMatrix, A₂::LinearMap) = *(WrappedMap(A₁), A₂) - -# An AdjointAbsVec isa AbstractMatrix but we handle it like a vector (Issue#99) -# which is an exception to the left multiplication rule that makes a WrappedMap -Base.:(*)(y::AdjointAbsVec, A::LinearMap) = adjoint(*(A', y')) -Base.:(*)(y::TransposeAbsVec, A::LinearMap) = transpose(transpose(A) * transpose(y)) diff --git a/test/left.jl b/test/left.jl new file mode 100644 index 00000000..9630400d --- /dev/null +++ b/test/left.jl @@ -0,0 +1,58 @@ +using Test, LinearMaps +#, LinearAlgebra + +# y'*L is an exception to the left multiplication rule that makes a WrappedMap + +function left_tester(L::LinearMap{T}) where {T} + M,N = size(L) + A = Matrix(L) + + x = rand(T, N) + y = rand(T, M) + + # * + @test y' * A ≈ y' * L + @test (y' * A) * x ≈ y' * (L * x) # associative + @test transpose(y) * A ≈ transpose(y) * L + + # mul! + b1 = y' * A + b2 = similar(b1) + mul!(b2, y', L) # 3-arg + @test b1 ≈ b2 + mul!(b2, y', L, true, false) # 5-arg + @test b1 ≈ b2 + + b1 = transpose(y) * A + b2 = similar(b1) + mul!(b2, transpose(y), L) # 3-arg + @test b1 ≈ b2 + mul!(b2, transpose(y), L, true, false) # 5-arg + @test b1 ≈ b2 + + true +end + + +@testset "left mul vec" begin + T = ComplexF32 + M,N = 5,5 + L = LinearMap{T}(cumsum, reverse ∘ cumsum ∘ reverse, N) + + @test left_tester(L) # FunctionMap + @test left_tester(L'*L) # CompositeMap + @test left_tester(2L) # ScaledMap + @test left_tester(kron(L,L')) # KroneckerMap + +#= +todo: fails with DimensionMismatch + W = LinearMap(randn(T,5,4)) # + @test left_tester(W) # WrappedMap +=# + +#= +todo: fails with stack overflow + @test left_tester([L L]) # BlockMap + @test left_tester(2L+3L') # LinearCombination +=# +end diff --git a/test/runtests.jl b/test/runtests.jl index d09d4bb2..ad807aa9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -28,3 +28,5 @@ include("blockmap.jl") include("kronecker.jl") include("conversion.jl") + +include("left.jl") diff --git a/test/wrappedmap.jl b/test/wrappedmap.jl index d07b13eb..57171b25 100644 --- a/test/wrappedmap.jl +++ b/test/wrappedmap.jl @@ -15,33 +15,3 @@ using Test, LinearMaps, LinearAlgebra @test @inferred isposdef(MA) @test @inferred isposdef(MB) end - -# y'*L is an exception to the left multiplication rule that makes a WrappedMap -@testset "left mul vec" begin - T = ComplexF32 - A = rand(T, 4, 5) - x = rand(T, 5) - y = rand(T, 4) - L = LinearMap{T}(A) - @test y'A ≈ y'L - @test (y' * A) * x ≈ y' * (L * x) # associative - @test transpose(y) * A ≈ transpose(y) * L - - LL = L * [L' L'] # stress test with CompositeMap & BlockMap - AA = A * [A' A'] - x = rand(T, 8) - @test y'AA ≈ y'LL - @test (y' * AA) * x ≈ y' * (LL * x) # associative - @test transpose(y)*AA ≈ transpose(y)*LL - - # mul! versions - b1 = y'*L - b2 = similar(b1) - mul!(b2, y', A) - @test b1 ≈ b2 - - b1 = transpose(y)*L - b2 = similar(b1) - mul!(b2, transpose(y), A) - @test b1 ≈ b2 -end From c97d467d5783a479cd52bfe9ec87ebeb33386dd8 Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Thu, 30 Jul 2020 00:02:43 -0400 Subject: [PATCH 04/18] Suspend testing 1.0 for now --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 67f187f1..685b8152 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -14,7 +14,7 @@ jobs: strategy: matrix: version: - - '1.0' +# - '1.0' - '1' - 'nightly' os: From f87cdcdfe9fa67ba6fedfad4ae6000eeaabacdc3 Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Thu, 30 Jul 2020 06:41:02 -0400 Subject: [PATCH 05/18] Add mul! to test --- test/left.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/left.jl b/test/left.jl index 9630400d..1ffb7d05 100644 --- a/test/left.jl +++ b/test/left.jl @@ -1,5 +1,5 @@ using Test, LinearMaps -#, LinearAlgebra +using LinearAlgebra: mul! # y'*L is an exception to the left multiplication rule that makes a WrappedMap From 16caf902e7eaff0e7e2a36ce44904c0b0ef23d4c Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Thu, 30 Jul 2020 06:42:01 -0400 Subject: [PATCH 06/18] Restore 1.0 test --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 685b8152..67f187f1 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -14,7 +14,7 @@ jobs: strategy: matrix: version: -# - '1.0' + - '1.0' - '1' - 'nightly' os: From 0775fa39fc6845a81355424193f817a2970416b5 Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Thu, 30 Jul 2020 07:03:50 -0400 Subject: [PATCH 07/18] Abandon pre 1.4 --- test/left.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/left.jl b/test/left.jl index 1ffb7d05..79086d1c 100644 --- a/test/left.jl +++ b/test/left.jl @@ -3,6 +3,8 @@ using LinearAlgebra: mul! # y'*L is an exception to the left multiplication rule that makes a WrappedMap +@static if VERSION ≥ v"1.4.0" + function left_tester(L::LinearMap{T}) where {T} M,N = size(L) A = Matrix(L) @@ -56,3 +58,5 @@ todo: fails with stack overflow @test left_tester(2L+3L') # LinearCombination =# end + +end # version From 702f4bb0e876df6bcb7064bc3b2ea387b02f414c Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Thu, 30 Jul 2020 21:59:15 -0400 Subject: [PATCH 08/18] Working version (!) with VecOut = AbstractVecOrMat --- src/LinearMaps.jl | 24 +++++++++--------------- src/left.jl | 19 +++++++++++-------- test/left.jl | 13 +++---------- 3 files changed, 23 insertions(+), 33 deletions(-) diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index 8be3546f..41d16847 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -4,7 +4,6 @@ export LinearMap export ⊗, kronsum, ⊕ using LinearAlgebra -import LinearAlgebra: AdjointAbsVec, TransposeAbsVec using SparseArrays if VERSION < v"1.2-" @@ -20,8 +19,7 @@ const MapOrMatrix{T} = Union{LinearMap{T},AbstractMatrix{T}} const RealOrComplex = Union{Real,Complex} # valid types for left vector multiplication: -const VecIn = Union{AdjointAbsVec, TransposeAbsVec} -const VecOut = Union{AbstractVector, AdjointAbsVec, TransposeAbsVec} +const VecOut = AbstractVecOrMat # todo - temporary Base.eltype(::LinearMap{T}) where {T} = T @@ -52,24 +50,20 @@ Base.ndims(::LinearMap) = 2 Base.size(A::LinearMap, n) = (n==1 || n==2 ? size(A)[n] : error("LinearMap objects have only 2 dimensions")) Base.length(A::LinearMap) = size(A)[1] * size(A)[2] -# check dimension consistency for y = A*x and Y = A*X -function check_dim_mul(y::VecOut, A::LinearMap, x::AbstractVector) +# check dimension consistency for right multiply: y = A*x and Y = A*X +function check_dim_mul(Y::AbstractVecOrMat, A::LinearMap, X::AbstractVecOrMat) # @info "checked vector dimensions" # uncomment for testing m, n = size(A) - (m == length(y) && n == length(x)) || throw(DimensionMismatch("mul!")) - return nothing -end -function check_dim_mul(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix) - # @info "checked matrix dimensions" # uncomment for testing - m, n = size(A) - (m == size(Y, 1) && n == size(X, 1) && size(Y, 2) == size(X, 2)) || throw(DimensionMismatch("mul!")) + (m == size(Y, 1) && n == size(X, 1) && size(Y, 2) == size(X, 2)) || + throw(DimensionMismatch("mul! $(size(X)) $(size(Y)) $(size(A))")) return nothing end -# check dimension consistency for left multiplication x = y'*A -function check_dim_mul(x::V, y::V, A::LinearMap) where {V <: VecIn} +# check dimension consistency for left multiply: X = Y*A (e.g., Y=y') +function check_dim_mul(X::AbstractVecOrMat, Y::AbstractVecOrMat, A::LinearMap) m, n = size(A) - ((1,m) == size(y) && (1,n) == size(x)) || throw(DimensionMismatch("left mul!")) + (n == size(X, 2) && m == size(Y, 2) && size(Y, 1) == size(X, 1)) || + throw(DimensionMismatch("left mul!")) return nothing end diff --git a/src/left.jl b/src/left.jl index 75c14378..3e9ff35c 100644 --- a/src/left.jl +++ b/src/left.jl @@ -10,6 +10,7 @@ import LinearAlgebra: AdjointAbsVec, TransposeAbsVec +# x = y'*A ⇐⇒ x' = (A'*y) Base.:(*)(y::AdjointAbsVec, A::LinearMap) = adjoint(*(A', y')) Base.:(*)(y::TransposeAbsVec, A::LinearMap) = transpose(transpose(A) * transpose(y)) @@ -17,19 +18,21 @@ Base.:(*)(y::TransposeAbsVec, A::LinearMap) = transpose(transpose(A) * transpose LinearAlgebra.mul!(x::AdjointAbsVec, y::AdjointAbsVec, A::LinearMap) = mul!(x, y, A, true, false) -# not sure if we need bounds checks and propagate inbounds stuff here +# todo: not sure if we need bounds checks and propagate inbounds stuff here # mul!(x, y', A, α, β) -function LinearAlgebra.mul!(x::AdjointAbsVec, y::AdjointAbsVec, A::LinearMap, - α::Number, β::Number) +# the key here is that "adjoint" and "transpose" are lazy "views" +function LinearAlgebra.mul!(x::AdjointAbsVec, y::AdjointAbsVec, + A::LinearMap, α::Number, β::Number) check_dim_mul(x, y, A) - mul!(x, A', y', α, β) - return conj!(x) + mul!(adjoint(x), A', y', α, β) + return adjoint(x) end # mul!(x, transpose(y), A, α, β) -function LinearAlgebra.mul!(x::TransposeAbsVec, y::TransposeAbsVec, A::LinearMap, - α::Number, β::Number) +function LinearAlgebra.mul!(x::TransposeAbsVec, y::TransposeAbsVec, + A::LinearMap, α::Number, β::Number) check_dim_mul(x, y, A) - return mul!(x, transpose(A), transpose(y), α, β) + mul!(transpose(x), transpose(A), transpose(y), α, β) + return transpose(x) end diff --git a/test/left.jl b/test/left.jl index 79086d1c..4478c48f 100644 --- a/test/left.jl +++ b/test/left.jl @@ -45,18 +45,11 @@ end @test left_tester(L'*L) # CompositeMap @test left_tester(2L) # ScaledMap @test left_tester(kron(L,L')) # KroneckerMap + @test left_tester(2L+3L') # LinearCombination + @test left_tester([L L]) # BlockMap -#= -todo: fails with DimensionMismatch - W = LinearMap(randn(T,5,4)) # + W = LinearMap(randn(T,5,4)) @test left_tester(W) # WrappedMap -=# - -#= -todo: fails with stack overflow - @test left_tester([L L]) # BlockMap - @test left_tester(2L+3L') # LinearCombination -=# end end # version From afc6b8f9a5a6ab7db72b517699c281095bcb1f9f Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Thu, 30 Jul 2020 22:18:00 -0400 Subject: [PATCH 09/18] Nix useless VecOut --- src/LinearMaps.jl | 13 +++++-------- src/blockmap.jl | 16 ++++++++-------- src/composition.jl | 14 +++++++------- src/functionmap.jl | 6 +++--- src/kronecker.jl | 18 +++++++++--------- src/linearcombination.jl | 6 +++--- src/scaledmap.jl | 8 ++++---- src/transpose.jl | 12 ++++++------ src/wrappedmap.jl | 6 +++--- 9 files changed, 48 insertions(+), 51 deletions(-) diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index 41d16847..27b95bd9 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -18,9 +18,6 @@ abstract type LinearMap{T} end const MapOrMatrix{T} = Union{LinearMap{T},AbstractMatrix{T}} const RealOrComplex = Union{Real,Complex} -# valid types for left vector multiplication: -const VecOut = AbstractVecOrMat # todo - temporary - Base.eltype(::LinearMap{T}) where {T} = T abstract type MulStyle end @@ -79,11 +76,11 @@ function Base.:(*)(A::LinearMap, x::AbstractVector) size(A, 2) == length(x) || throw(DimensionMismatch("mul!")) return @inbounds A_mul_B!(similar(x, promote_type(eltype(A), eltype(x)), size(A, 1)), A, x) end -function LinearAlgebra.mul!(y::VecOut, A::LinearMap, x::AbstractVector) +function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap, x::AbstractVector) @boundscheck check_dim_mul(y, A, x) return @inbounds A_mul_B!(y, A, x) end -function LinearAlgebra.mul!(y::VecOut, A::LinearMap, x::AbstractVector, α::Number, β::Number) +function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap, x::AbstractVector, α::Number, β::Number) @boundscheck check_dim_mul(y, A, x) if isone(α) iszero(β) && (A_mul_B!(y, A, x); return y) @@ -120,9 +117,9 @@ Base.@propagate_inbounds function LinearAlgebra.mul!(Y::AbstractMatrix, A::Linea return Y end -A_mul_B!(y::VecOut, A::AbstractMatrix, x::AbstractVector) = mul!(y, A, x) -At_mul_B!(y::VecOut, A::AbstractMatrix, x::AbstractVector) = mul!(y, transpose(A), x) -Ac_mul_B!(y::VecOut, A::AbstractMatrix, x::AbstractVector) = mul!(y, adjoint(A), x) +A_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, A, x) +At_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, transpose(A), x) +Ac_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, adjoint(A), x) include("left.jl") # left multiplication by a transpose or adjoint vector include("wrappedmap.jl") # wrap a matrix of linear map in a new type, thereby allowing to alter its properties diff --git a/src/blockmap.jl b/src/blockmap.jl index ff88e066..becc2de4 100644 --- a/src/blockmap.jl +++ b/src/blockmap.jl @@ -382,19 +382,19 @@ end # multiplication with vectors & matrices ############ -Base.@propagate_inbounds A_mul_B!(y::VecOut, A::BlockMap, x::AbstractVector) = +Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) = mul!(y, A, x, true, false) -Base.@propagate_inbounds A_mul_B!(y::VecOut, A::TransposeMap{<:Any,<:BlockMap}, x::AbstractVector) = +Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::TransposeMap{<:Any,<:BlockMap}, x::AbstractVector) = mul!(y, A, x, true, false) -Base.@propagate_inbounds At_mul_B!(y::VecOut, A::BlockMap, x::AbstractVector) = +Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) = mul!(y, transpose(A), x, true, false) -Base.@propagate_inbounds A_mul_B!(y::VecOut, A::AdjointMap{<:Any,<:BlockMap}, x::AbstractVector) = +Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::AdjointMap{<:Any,<:BlockMap}, x::AbstractVector) = mul!(y, A, x, true, false) -Base.@propagate_inbounds Ac_mul_B!(y::VecOut, A::BlockMap, x::AbstractVector) = +Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) = mul!(y, adjoint(A), x, true, false) for Atype in (AbstractVector, AbstractMatrix) @@ -495,13 +495,13 @@ LinearAlgebra.transpose(A::BlockDiagonalMap{T}) where {T} = BlockDiagonalMap{T}( Base.:(==)(A::BlockDiagonalMap, B::BlockDiagonalMap) = (eltype(A) == eltype(B) && A.maps == B.maps) -Base.@propagate_inbounds A_mul_B!(y::VecOut, A::BlockDiagonalMap, x::AbstractVector) = +Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::BlockDiagonalMap, x::AbstractVector) = mul!(y, A, x, true, false) -Base.@propagate_inbounds At_mul_B!(y::VecOut, A::BlockDiagonalMap, x::AbstractVector) = +Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::BlockDiagonalMap, x::AbstractVector) = mul!(y, transpose(A), x, true, false) -Base.@propagate_inbounds Ac_mul_B!(y::VecOut, A::BlockDiagonalMap, x::AbstractVector) = +Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::BlockDiagonalMap, x::AbstractVector) = mul!(y, adjoint(A), x, true, false) for Atype in (AbstractVector, AbstractMatrix) diff --git a/src/composition.jl b/src/composition.jl index 9d91edee..b3ce93fc 100644 --- a/src/composition.jl +++ b/src/composition.jl @@ -129,23 +129,23 @@ LinearAlgebra.adjoint(A::CompositeMap{T}) where {T} = CompositeMap{T}(map(adjo Base.:(==)(A::CompositeMap, B::CompositeMap) = (eltype(A) == eltype(B) && A.maps == B.maps) # multiplication with vectors -function A_mul_B!(y::VecOut, A::CompositeMap{T,<:Tuple{LinearMap}}, x::AbstractVector) where {T} +function A_mul_B!(y::AbstractVector, A::CompositeMap{T,<:Tuple{LinearMap}}, x::AbstractVector) where {T} return A_mul_B!(y, A.maps[1], x) end -function A_mul_B!(y::VecOut, A::CompositeMap{T,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector) where {T} +function A_mul_B!(y::AbstractVector, A::CompositeMap{T,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector) where {T} _compositemul!(y, A, x, similar(y, size(A.maps[1], 1))) end -function A_mul_B!(y::VecOut, A::CompositeMap{T,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector) where {T} +function A_mul_B!(y::AbstractVector, A::CompositeMap{T,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector) where {T} _compositemul!(y, A, x, similar(y, size(A.maps[1], 1)), similar(y, size(A.maps[2], 1))) end -function _compositemul!(y::VecOut, A::CompositeMap{T,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector, z::AbstractVector) where {T} +function _compositemul!(y::AbstractVector, A::CompositeMap{T,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector, z::AbstractVector) where {T} # no size checking, will be done by individual maps A_mul_B!(z, A.maps[1], x) A_mul_B!(y, A.maps[2], z) return y end -function _compositemul!(y::VecOut, A::CompositeMap{T,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector, source::AbstractVector, dest::AbstractVector) where {T} +function _compositemul!(y::AbstractVector, A::CompositeMap{T,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector, source::AbstractVector, dest::AbstractVector) where {T} # no size checking, will be done by individual maps N = length(A.maps) A_mul_B!(source, A.maps[1], x) @@ -166,6 +166,6 @@ function _compositemul!(y::VecOut, A::CompositeMap{T,<:Tuple{Vararg{LinearMap}}} return y end -At_mul_B!(y::VecOut, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) +At_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) -Ac_mul_B!(y::VecOut, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) +Ac_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) diff --git a/src/functionmap.jl b/src/functionmap.jl index 937dc86c..4e1b086f 100644 --- a/src/functionmap.jl +++ b/src/functionmap.jl @@ -107,13 +107,13 @@ function Base.:(*)(A::TransposeMap{<:Any,<:FunctionMap}, x::AbstractVector) end end -function A_mul_B!(y::VecOut, A::FunctionMap, x::AbstractVector) +function A_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector) (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::VecOut, A::FunctionMap, x::AbstractVector) +function At_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector) (issymmetric(A) || (isreal(A) && ishermitian(A))) && return A_mul_B!(y, A, x) (length(x) == A.M && length(y) == A.N) || throw(DimensionMismatch("At_mul_B!")) if A.fc !== nothing @@ -135,7 +135,7 @@ function At_mul_B!(y::VecOut, A::FunctionMap, x::AbstractVector) end end -function Ac_mul_B!(y::VecOut, A::FunctionMap, x::AbstractVector) +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("Ac_mul_B!")) if A.fc !== nothing diff --git a/src/kronecker.jl b/src/kronecker.jl index be4db862..80fd502a 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -131,7 +131,7 @@ end # multiplication with vectors ################# -Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T} +Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T} require_one_based_indexing(y) @boundscheck check_dim_mul(y, L, x) A, B = L.maps @@ -139,7 +139,7 @@ Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::KroneckerMap{T,<:NTuple _kronmul!(y, B, X, transpose(A), T) return y end -Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::KroneckerMap{T}, x::AbstractVector) where {T} +Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T}, x::AbstractVector) where {T} require_one_based_indexing(y) @boundscheck check_dim_mul(y, L, x) A = first(L.maps) @@ -150,7 +150,7 @@ Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::KroneckerMap{T}, x::Abs end # mixed-product rule, prefer the right if possible # (A₁ ⊗ A₂ ⊗ ... ⊗ Aᵣ) * (B₁ ⊗ B₂ ⊗ ... ⊗ Bᵣ) = (A₁B₁) ⊗ (A₂B₂) ⊗ ... ⊗ (AᵣBᵣ) -Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::CompositeMap{<:Any,<:Tuple{KroneckerMap,KroneckerMap}}, x::AbstractVector) +Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{<:Any,<:Tuple{KroneckerMap,KroneckerMap}}, x::AbstractVector) B, A = L.maps if length(A.maps) == length(B.maps) && all(M -> check_dim_mul(M[1], M[2]), zip(A.maps, B.maps)) A_mul_B!(y, kron(map(*, A.maps, B.maps)...), x) @@ -160,7 +160,7 @@ Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::CompositeMap{<:Any,<:Tu end # mixed-product rule, prefer the right if possible # (A₁ ⊗ B₁)*(A₂⊗B₂)*...*(Aᵣ⊗Bᵣ) = (A₁*A₂*...*Aᵣ) ⊗ (B₁*B₂*...*Bᵣ) -Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::CompositeMap{T,<:Tuple{Vararg{KroneckerMap{<:Any,<:Tuple{LinearMap,LinearMap}}}}}, x::AbstractVector) where {T} +Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{T,<:Tuple{Vararg{KroneckerMap{<:Any,<:Tuple{LinearMap,LinearMap}}}}}, x::AbstractVector) where {T} As = map(AB -> AB.maps[1], L.maps) Bs = map(AB -> AB.maps[2], L.maps) As1, As2 = Base.front(As), Base.tail(As) @@ -173,10 +173,10 @@ Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::CompositeMap{T,<:Tuple{ end end -Base.@propagate_inbounds At_mul_B!(y::VecOut, A::KroneckerMap, x::AbstractVector) = +Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) -Base.@propagate_inbounds Ac_mul_B!(y::VecOut, A::KroneckerMap, x::AbstractVector) = +Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) ############### @@ -261,7 +261,7 @@ LinearAlgebra.transpose(A::KroneckerSumMap{T}) where {T} = KroneckerSumMap{T}(ma Base.:(==)(A::KroneckerSumMap, B::KroneckerSumMap) = (eltype(A) == eltype(B) && A.maps == B.maps) -Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::KroneckerSumMap, x::AbstractVector) +Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerSumMap, x::AbstractVector) @boundscheck check_dim_mul(y, L, x) A, B = L.maps ma, na = size(A) @@ -273,8 +273,8 @@ Base.@propagate_inbounds function A_mul_B!(y::VecOut, L::KroneckerSumMap, x::Abs return y end -Base.@propagate_inbounds At_mul_B!(y::VecOut, A::KroneckerSumMap, x::AbstractVector) = +Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) -Base.@propagate_inbounds Ac_mul_B!(y::VecOut, A::KroneckerSumMap, x::AbstractVector) = +Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) diff --git a/src/linearcombination.jl b/src/linearcombination.jl index 705c8bbe..153c861f 100644 --- a/src/linearcombination.jl +++ b/src/linearcombination.jl @@ -117,8 +117,8 @@ end return y end -A_mul_B!(y::VecOut, A::LinearCombination, x::AbstractVector) = mul!(y, A, x, true, false) +A_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = mul!(y, A, x, true, false) -At_mul_B!(y::VecOut, A::LinearCombination, x::AbstractVector) = mul!(y, transpose(A), x, true, false) +At_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = mul!(y, transpose(A), x, true, false) -Ac_mul_B!(y::VecOut, A::LinearCombination, x::AbstractVector) = mul!(y, adjoint(A), x, true, false) +Ac_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = mul!(y, adjoint(A), x, true, false) diff --git a/src/scaledmap.jl b/src/scaledmap.jl index dca651c3..14b914c1 100644 --- a/src/scaledmap.jl +++ b/src/scaledmap.jl @@ -54,15 +54,15 @@ Base.:(*)(A::ScaledMap, B::LinearMap) = A.λ * (A.lmap * B) Base.:(*)(A::LinearMap, B::ScaledMap) = (A * B.lmap) * B.λ # multiplication with vectors -function A_mul_B!(y::VecOut, A::ScaledMap, x::AbstractVector) +function A_mul_B!(y::AbstractVector, A::ScaledMap, x::AbstractVector) # no size checking, will be done by map mul!(y, A.lmap, x, A.λ, false) end -function LinearAlgebra.mul!(y::VecOut, A::ScaledMap, x::AbstractVector, α::Number, β::Number) +function LinearAlgebra.mul!(y::AbstractVector, A::ScaledMap, x::AbstractVector, α::Number, β::Number) # no size checking, will be done by map mul!(y, A.lmap, x, A.λ * α, β) end -At_mul_B!(y::VecOut, A::ScaledMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) -Ac_mul_B!(y::VecOut, A::ScaledMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) +At_mul_B!(y::AbstractVector, A::ScaledMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) +Ac_mul_B!(y::AbstractVector, A::ScaledMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) diff --git a/src/transpose.jl b/src/transpose.jl index 0a60d6bf..60caf3c2 100644 --- a/src/transpose.jl +++ b/src/transpose.jl @@ -32,18 +32,18 @@ Base.:(==)(A::LinearMap, B::TransposeMap) = issymmetric(A) && B.lmap == A Base.:(==)(A::LinearMap, B::AdjointMap) = ishermitian(A) && B.lmap == A # multiplication with vector -A_mul_B!(y::VecOut, A::TransposeMap, x::AbstractVector) = +A_mul_B!(y::AbstractVector, A::TransposeMap, x::AbstractVector) = issymmetric(A.lmap) ? A_mul_B!(y, A.lmap, x) : At_mul_B!(y, A.lmap, x) -At_mul_B!(y::VecOut, A::TransposeMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x) +At_mul_B!(y::AbstractVector, A::TransposeMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x) -Ac_mul_B!(y::VecOut, A::TransposeMap, x::AbstractVector) = +Ac_mul_B!(y::AbstractVector, A::TransposeMap, x::AbstractVector) = isreal(A.lmap) ? A_mul_B!(y, A.lmap, x) : (A_mul_B!(y, A.lmap, conj(x)); conj!(y)) -A_mul_B!(y::VecOut, A::AdjointMap, x::AbstractVector) = +A_mul_B!(y::AbstractVector, A::AdjointMap, x::AbstractVector) = ishermitian(A.lmap) ? A_mul_B!(y, A.lmap, x) : Ac_mul_B!(y, A.lmap, x) -At_mul_B!(y::VecOut, A::AdjointMap, x::AbstractVector) = +At_mul_B!(y::AbstractVector, A::AdjointMap, x::AbstractVector) = isreal(A.lmap) ? A_mul_B!(y, A.lmap, x) : (A_mul_B!(y, A.lmap, conj(x)); conj!(y)) -Ac_mul_B!(y::VecOut, A::AdjointMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x) +Ac_mul_B!(y::AbstractVector, A::AdjointMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x) diff --git a/src/wrappedmap.jl b/src/wrappedmap.jl index 321380dc..574827ea 100644 --- a/src/wrappedmap.jl +++ b/src/wrappedmap.jl @@ -37,13 +37,13 @@ LinearAlgebra.ishermitian(A::WrappedMap) = A._ishermitian LinearAlgebra.isposdef(A::WrappedMap) = A._isposdef # multiplication with vectors & matrices -A_mul_B!(y::VecOut, A::WrappedMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x) +A_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x) Base.:(*)(A::WrappedMap, x::AbstractVector) = *(A.lmap, x) -At_mul_B!(y::VecOut, A::WrappedMap, x::AbstractVector) = +At_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) = (issymmetric(A) || (isreal(A) && ishermitian(A))) ? A_mul_B!(y, A.lmap, x) : At_mul_B!(y, A.lmap, x) -Ac_mul_B!(y::VecOut, A::WrappedMap, x::AbstractVector) = +Ac_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) = ishermitian(A) ? A_mul_B!(y, A.lmap, x) : Ac_mul_B!(y, A.lmap, x) if VERSION ≥ v"1.3.0-alpha.115" From 4c770c458d22db97433274f3fb5b1f2241d0562f Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Thu, 30 Jul 2020 22:48:42 -0400 Subject: [PATCH 10/18] Try 1.0 again --- test/left.jl | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/test/left.jl b/test/left.jl index 4478c48f..45674407 100644 --- a/test/left.jl +++ b/test/left.jl @@ -3,16 +3,16 @@ using LinearAlgebra: mul! # y'*L is an exception to the left multiplication rule that makes a WrappedMap -@static if VERSION ≥ v"1.4.0" +@static if true # VERSION ≥ v"1.4.0" function left_tester(L::LinearMap{T}) where {T} - M,N = size(L) + M,N = size(L) A = Matrix(L) x = rand(T, N) y = rand(T, M) - # * + # * @test y' * A ≈ y' * L @test (y' * A) * x ≈ y' * (L * x) # associative @test transpose(y) * A ≈ transpose(y) * L @@ -32,24 +32,24 @@ function left_tester(L::LinearMap{T}) where {T} mul!(b2, transpose(y), L, true, false) # 5-arg @test b1 ≈ b2 - true + true end @testset "left mul vec" begin T = ComplexF32 - M,N = 5,5 - L = LinearMap{T}(cumsum, reverse ∘ cumsum ∘ reverse, N) - - @test left_tester(L) # FunctionMap - @test left_tester(L'*L) # CompositeMap - @test left_tester(2L) # ScaledMap - @test left_tester(kron(L,L')) # KroneckerMap - @test left_tester(2L+3L') # LinearCombination - @test left_tester([L L]) # BlockMap - - W = LinearMap(randn(T,5,4)) - @test left_tester(W) # WrappedMap + M,N = 5,5 + L = LinearMap{T}(cumsum, reverse ∘ cumsum ∘ reverse, N) + + @test left_tester(L) # FunctionMap + @test left_tester(L'*L) # CompositeMap + @test left_tester(2L) # ScaledMap + @test left_tester(kron(L,L')) # KroneckerMap + @test left_tester(2L+3L') # LinearCombination + @test left_tester([L L]) # BlockMap + + W = LinearMap(randn(T,5,4)) + @test left_tester(W) # WrappedMap end end # version From d006c82dee113385c0513940526805d1693c3716 Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Thu, 30 Jul 2020 23:06:09 -0400 Subject: [PATCH 11/18] Narrow 1.0 exclusion --- test/left.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/test/left.jl b/test/left.jl index 45674407..f38c3191 100644 --- a/test/left.jl +++ b/test/left.jl @@ -3,7 +3,6 @@ using LinearAlgebra: mul! # y'*L is an exception to the left multiplication rule that makes a WrappedMap -@static if true # VERSION ≥ v"1.4.0" function left_tester(L::LinearMap{T}) where {T} M,N = size(L) @@ -27,10 +26,12 @@ function left_tester(L::LinearMap{T}) where {T} b1 = transpose(y) * A b2 = similar(b1) - mul!(b2, transpose(y), L) # 3-arg +@static if VERSION ≥ v"1.4.0" + mul!(b2, transpose(y), L) # 3-arg # todo: fails in Julia 1.0 @test b1 ≈ b2 mul!(b2, transpose(y), L, true, false) # 5-arg @test b1 ≈ b2 +end # version true end @@ -51,5 +52,3 @@ end W = LinearMap(randn(T,5,4)) @test left_tester(W) # WrappedMap end - -end # version From e8f57bb01f94b87c4971b10b89fee9538525cd5a Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 31 Jul 2020 17:29:44 +0200 Subject: [PATCH 12/18] minor fixes --- src/LinearMaps.jl | 23 ++++++++++------------- src/left.jl | 30 +++++++++++++++++------------- test/left.jl | 6 ++---- 3 files changed, 29 insertions(+), 30 deletions(-) diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index 27b95bd9..2b94cad7 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -47,20 +47,17 @@ Base.ndims(::LinearMap) = 2 Base.size(A::LinearMap, n) = (n==1 || n==2 ? size(A)[n] : error("LinearMap objects have only 2 dimensions")) Base.length(A::LinearMap) = size(A)[1] * size(A)[2] -# check dimension consistency for right multiply: y = A*x and Y = A*X -function check_dim_mul(Y::AbstractVecOrMat, A::LinearMap, X::AbstractVecOrMat) +# check dimension consistency for multiplication C = A*B +function check_dim_mul(C, A, B) # @info "checked vector dimensions" # uncomment for testing - m, n = size(A) - (m == size(Y, 1) && n == size(X, 1) && size(Y, 2) == size(X, 2)) || - throw(DimensionMismatch("mul! $(size(X)) $(size(Y)) $(size(A))")) - return nothing -end - -# check dimension consistency for left multiply: X = Y*A (e.g., Y=y') -function check_dim_mul(X::AbstractVecOrMat, Y::AbstractVecOrMat, A::LinearMap) - m, n = size(A) - (n == size(X, 2) && m == size(Y, 2) && size(Y, 1) == size(X, 1)) || - throw(DimensionMismatch("left mul!")) + mA, nA = size(A) # A always has two dimensions + mB, nB = size(B, 1), size(B, 2) + if mB != nA + throw(DimensionMismatch("left factor has dimensions ($mA,$nA), right factor has dimensions ($mB,$nB)")) + end + if size(C, 1) != mA || size(C, 2) != nB + throw(DimensionMismatch("result has dimensions $(size(C)), needs ($mA,$nB)")) + end return nothing end diff --git a/src/left.jl b/src/left.jl index 3e9ff35c..84f2cbcf 100644 --- a/src/left.jl +++ b/src/left.jl @@ -15,24 +15,28 @@ Base.:(*)(y::AdjointAbsVec, A::LinearMap) = adjoint(*(A', y')) Base.:(*)(y::TransposeAbsVec, A::LinearMap) = transpose(transpose(A) * transpose(y)) # mul!(x, y', A) -LinearAlgebra.mul!(x::AdjointAbsVec, y::AdjointAbsVec, A::LinearMap) = - mul!(x, y, A, true, false) - -# todo: not sure if we need bounds checks and propagate inbounds stuff here +Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::AdjointAbsVec, A::LinearMap) + @boundscheck check_dim_mul(x, y, A) + @inbounds mul!(adjoint(x), A', y') + return adjoint(x) +end -# mul!(x, y', A, α, β) -# the key here is that "adjoint" and "transpose" are lazy "views" -function LinearAlgebra.mul!(x::AdjointAbsVec, y::AdjointAbsVec, +Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::AdjointAbsVec, A::LinearMap, α::Number, β::Number) - check_dim_mul(x, y, A) - mul!(adjoint(x), A', y', α, β) + @boundscheck check_dim_mul(x, y, A) + @inbounds mul!(adjoint(x), A', y', α, β) return adjoint(x) end -# mul!(x, transpose(y), A, α, β) -function LinearAlgebra.mul!(x::TransposeAbsVec, y::TransposeAbsVec, +Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::TransposeAbsVec, A::LinearMap) + @boundscheck check_dim_mul(x, y, A) + @inbounds mul!(transpose(x), transpose(A), transpose(y)) + return transpose(x) +end + +Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::TransposeAbsVec, A::LinearMap, α::Number, β::Number) - check_dim_mul(x, y, A) - mul!(transpose(x), transpose(A), transpose(y), α, β) + @boundscheck check_dim_mul(x, y, A) + @inbounds mul!(transpose(x), transpose(A), transpose(y), α, β) return transpose(x) end diff --git a/test/left.jl b/test/left.jl index f38c3191..232d627d 100644 --- a/test/left.jl +++ b/test/left.jl @@ -5,7 +5,7 @@ using LinearAlgebra: mul! function left_tester(L::LinearMap{T}) where {T} - M,N = size(L) + M, N = size(L) A = Matrix(L) x = rand(T, N) @@ -26,12 +26,10 @@ function left_tester(L::LinearMap{T}) where {T} b1 = transpose(y) * A b2 = similar(b1) -@static if VERSION ≥ v"1.4.0" - mul!(b2, transpose(y), L) # 3-arg # todo: fails in Julia 1.0 + mul!(b2, transpose(y), L) # 3-arg @test b1 ≈ b2 mul!(b2, transpose(y), L, true, false) # 5-arg @test b1 ≈ b2 -end # version true end From 1fbcae13fd8735433ee85430cfaf5580dae55693 Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Fri, 31 Jul 2020 12:01:55 -0400 Subject: [PATCH 13/18] For codecov --- src/LinearMaps.jl | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index 2b94cad7..4281ce2c 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -52,12 +52,10 @@ function check_dim_mul(C, A, B) # @info "checked vector dimensions" # uncomment for testing mA, nA = size(A) # A always has two dimensions mB, nB = size(B, 1), size(B, 2) - if mB != nA - throw(DimensionMismatch("left factor has dimensions ($mA,$nA), right factor has dimensions ($mB,$nB)")) - end - if size(C, 1) != mA || size(C, 2) != nB + (mB == nA) || + throw(DimensionMismatch("left factor has dimensions ($mA,$nA), right factor has dimensions ($mB,$nB)") + (size(C, 1) != mA || size(C, 2) != nB) && throw(DimensionMismatch("result has dimensions $(size(C)), needs ($mA,$nB)")) - end return nothing end From 46489e1ce7192f41857a42f4008ddf3b9d7c242c Mon Sep 17 00:00:00 2001 From: Jeff Fessler Date: Fri, 31 Jul 2020 12:04:10 -0400 Subject: [PATCH 14/18] Oops, lost a paren --- src/LinearMaps.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index 4281ce2c..6e764f4a 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -53,7 +53,7 @@ function check_dim_mul(C, A, B) mA, nA = size(A) # A always has two dimensions mB, nB = size(B, 1), size(B, 2) (mB == nA) || - throw(DimensionMismatch("left factor has dimensions ($mA,$nA), right factor has dimensions ($mB,$nB)") + throw(DimensionMismatch("left factor has dimensions ($mA,$nA), right factor has dimensions ($mB,$nB)")) (size(C, 1) != mA || size(C, 2) != nB) && throw(DimensionMismatch("result has dimensions $(size(C)), needs ($mA,$nB)")) return nothing From 22b6fc950217ba891a6997e7b2ff65abecc0b53a Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sat, 1 Aug 2020 18:26:17 +0200 Subject: [PATCH 15/18] add dimension checking tests --- test/linearmaps.jl | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/test/linearmaps.jl b/test/linearmaps.jl index ea708f8c..b1af5e2f 100644 --- a/test/linearmaps.jl +++ b/test/linearmaps.jl @@ -47,6 +47,15 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools @test @inferred mul!(copy(W), M, V) ≈ AV @test typeof(M * V) <: LinearMap end + + @testset "dimension checking" begin + @test_throws DimensionMismatch M * similar(v, length(v) + 1) + @test_throws DimensionMismatch mul!(similar(w, length(w) + 1), M, v) + @test_throws DimensionMismatch similar(w, length(w) + 1)' * M + @test_throws DimensionMismatch mul!(copy(v)', similar(w, length(w) + 1)', M) + @test_throws DimensionMismatch mul!(similar(W, size(W).+(0,1)), M, V) + @test_throws DimensionMismatch mul!(copy(W), M, similar(V, size(V).+(0,1))) + end end # new type From ad2d33e0e27591d93ddf615979d625eb2033c3d3 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 7 Aug 2020 15:49:36 +0200 Subject: [PATCH 16/18] include left matrix multiply and some fixes --- src/left.jl | 49 ++++++++++++++++++++++++++----------------------- test/left.jl | 43 ++++++++++++++++++++++++++++++------------- 2 files changed, 56 insertions(+), 36 deletions(-) diff --git a/src/left.jl b/src/left.jl index 84f2cbcf..32e8d993 100644 --- a/src/left.jl +++ b/src/left.jl @@ -11,32 +11,35 @@ import LinearAlgebra: AdjointAbsVec, TransposeAbsVec # x = y'*A ⇐⇒ x' = (A'*y) -Base.:(*)(y::AdjointAbsVec, A::LinearMap) = adjoint(*(A', y')) +Base.:(*)(y::AdjointAbsVec, A::LinearMap) = adjoint(A' * y') Base.:(*)(y::TransposeAbsVec, A::LinearMap) = transpose(transpose(A) * transpose(y)) -# mul!(x, y', A) -Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::AdjointAbsVec, A::LinearMap) - @boundscheck check_dim_mul(x, y, A) - @inbounds mul!(adjoint(x), A', y') - return adjoint(x) -end - -Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::AdjointAbsVec, - A::LinearMap, α::Number, β::Number) - @boundscheck check_dim_mul(x, y, A) - @inbounds mul!(adjoint(x), A', y', α, β) - return adjoint(x) -end +# multiplication with vector/matrix +for Atype in (AdjointAbsVec, Adjoint{<:Any,<:AbstractMatrix}) + @eval Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::$Atype, A::LinearMap) + @boundscheck check_dim_mul(x, y, A) + @inbounds mul!(adjoint(x), A', y') + return x + end -Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::TransposeAbsVec, A::LinearMap) - @boundscheck check_dim_mul(x, y, A) - @inbounds mul!(transpose(x), transpose(A), transpose(y)) - return transpose(x) + @eval Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::$Atype, A::LinearMap, + α::Number, β::Number) + @boundscheck check_dim_mul(x, y, A) + @inbounds mul!(adjoint(x), A', y', conj(α), conj(β)) + return x + end end +for Atype in (TransposeAbsVec, Transpose{<:Any,<:AbstractMatrix}) + @eval Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::$Atype, A::LinearMap) + @boundscheck check_dim_mul(x, y, A) + @inbounds mul!(transpose(x), transpose(A), transpose(y)) + return x + end -Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::TransposeAbsVec, - A::LinearMap, α::Number, β::Number) - @boundscheck check_dim_mul(x, y, A) - @inbounds mul!(transpose(x), transpose(A), transpose(y), α, β) - return transpose(x) + @eval Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::$Atype, A::LinearMap, + α::Number, β::Number) + @boundscheck check_dim_mul(x, y, A) + @inbounds mul!(transpose(x), transpose(A), transpose(y), α, β) + return x + end end diff --git a/test/left.jl b/test/left.jl index 232d627d..b76cfbf4 100644 --- a/test/left.jl +++ b/test/left.jl @@ -17,34 +17,51 @@ function left_tester(L::LinearMap{T}) where {T} @test transpose(y) * A ≈ transpose(y) * L # mul! + α = rand(T); β = rand(T) b1 = y' * A b2 = similar(b1) - mul!(b2, y', L) # 3-arg - @test b1 ≈ b2 - mul!(b2, y', L, true, false) # 5-arg - @test b1 ≈ b2 + bt = copy(b1')' + # bm = Matrix(bt) # TODO: this requires a generalization of the output to AbstractVecOrMat + @test mul!(b2, y', L) ≈ mul!(bt, y', L)# ≈ mul!(bm, y', L)# 3-arg + @test mul!(b2, y', L) === b2 + @test b1 ≈ b2 ≈ bt + b3 = copy(b2) + mul!(b3, y', L, α, β) + @test b3 ≈ b2*β + b1*α b1 = transpose(y) * A b2 = similar(b1) - mul!(b2, transpose(y), L) # 3-arg - @test b1 ≈ b2 - mul!(b2, transpose(y), L, true, false) # 5-arg - @test b1 ≈ b2 + bt = transpose(copy(transpose(b1))) + @test mul!(b2, transpose(y), L) ≈ mul!(bt, transpose(y), L) + @test b1 ≈ b2 ≈ bt + b3 = copy(b2) + mul!(b3, transpose(y), L, α, β) + @test b3 ≈ b2*β + b1*α + + Y = rand(T, M, 3) + X = similar(Y, 3, N) + Xt = copy(X')' + @test Y'L isa LinearMap + @test Matrix(Y'L) ≈ Y'A + @test mul!(X, Y', L) ≈ mul!(Xt, Y', L) ≈ Y'A + @test mul!(Xt, Y', L) === Xt + @test mul!(copy(X), Y', L, α, β) ≈ X*β + Y'A*α + @test mul!(X, Y', L) === X + @test mul!(X, Y', L, α, β) === X true end - -@testset "left mul vec" begin +@testset "left multiplication" begin T = ComplexF32 - M,N = 5,5 + N = 5 L = LinearMap{T}(cumsum, reverse ∘ cumsum ∘ reverse, N) @test left_tester(L) # FunctionMap @test left_tester(L'*L) # CompositeMap @test left_tester(2L) # ScaledMap - @test left_tester(kron(L,L')) # KroneckerMap - @test left_tester(2L+3L') # LinearCombination + @test left_tester(kron(L, L')) # KroneckerMap + @test left_tester(2L + 3L') # LinearCombination @test left_tester([L L]) # BlockMap W = LinearMap(randn(T,5,4)) From b760a469271ecc54d8f149479bcaaea9b0ff204d Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 24 Aug 2020 11:41:52 +0200 Subject: [PATCH 17/18] final tweaks, don't bump version number (yet) --- Project.toml | 2 +- src/left.jl | 44 ++++++++++++++++++++++++-------------------- test/left.jl | 6 +----- 3 files changed, 26 insertions(+), 26 deletions(-) diff --git a/Project.toml b/Project.toml index f1af504e..7ed11620 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "LinearMaps" uuid = "7a12625a-238d-50fd-b39a-03d52299707e" -version = "2.8.0" +version = "2.7.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/left.jl b/src/left.jl index 32e8d993..7561bd1e 100644 --- a/src/left.jl +++ b/src/left.jl @@ -16,30 +16,34 @@ Base.:(*)(y::TransposeAbsVec, A::LinearMap) = transpose(transpose(A) * transpose # multiplication with vector/matrix for Atype in (AdjointAbsVec, Adjoint{<:Any,<:AbstractMatrix}) - @eval Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::$Atype, A::LinearMap) - @boundscheck check_dim_mul(x, y, A) - @inbounds mul!(adjoint(x), A', y') - return x - end + @eval begin + Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::$Atype, A::LinearMap) + @boundscheck check_dim_mul(x, y, A) + @inbounds mul!(x', A', y') + return x + end - @eval Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::$Atype, A::LinearMap, - α::Number, β::Number) - @boundscheck check_dim_mul(x, y, A) - @inbounds mul!(adjoint(x), A', y', conj(α), conj(β)) - return x + Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::$Atype, A::LinearMap, + α::Number, β::Number) + @boundscheck check_dim_mul(x, y, A) + @inbounds mul!(x', A', y', conj(α), conj(β)) + return x + end end end for Atype in (TransposeAbsVec, Transpose{<:Any,<:AbstractMatrix}) - @eval Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::$Atype, A::LinearMap) - @boundscheck check_dim_mul(x, y, A) - @inbounds mul!(transpose(x), transpose(A), transpose(y)) - return x - end + @eval begin + Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::$Atype, A::LinearMap) + @boundscheck check_dim_mul(x, y, A) + @inbounds mul!(transpose(x), transpose(A), transpose(y)) + return x + end - @eval Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::$Atype, A::LinearMap, - α::Number, β::Number) - @boundscheck check_dim_mul(x, y, A) - @inbounds mul!(transpose(x), transpose(A), transpose(y), α, β) - return x + Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::$Atype, A::LinearMap, + α::Number, β::Number) + @boundscheck check_dim_mul(x, y, A) + @inbounds mul!(transpose(x), transpose(A), transpose(y), α, β) + return x + end end end diff --git a/test/left.jl b/test/left.jl index b76cfbf4..8f3777ba 100644 --- a/test/left.jl +++ b/test/left.jl @@ -1,8 +1,4 @@ -using Test, LinearMaps -using LinearAlgebra: mul! - -# y'*L is an exception to the left multiplication rule that makes a WrappedMap - +using Test, LinearMaps, LinearAlgebra function left_tester(L::LinearMap{T}) where {T} M, N = size(L) From 7fff78da2d17c0f2306cd66c0e6250ac2180964e Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 24 Aug 2020 11:45:09 +0200 Subject: [PATCH 18/18] don't announce in README (yet) --- README.md | 4 ---- 1 file changed, 4 deletions(-) diff --git a/README.md b/README.md index 24f0cc5e..155b33bc 100644 --- a/README.md +++ b/README.md @@ -8,10 +8,6 @@ A Julia package for defining and working with linear maps, also known as linear transformations or linear operators acting on vectors. The only requirement for a LinearMap is that it can act on a vector (by multiplication) efficiently. -## What's new in v2.8 -* Left multiplying by a transpose or adjoint vector (e.g., `y'*A`) - produces a transpose or adjoint vector output, rather than a `LinearMap`. - ## What's new in v2.7 * Potential reduction of memory allocations in multiplication of `LinearCombination`s, `BlockMap`s, and real- or complex-scaled `LinearMap`s.