diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index 30419ca5..6e764f4a 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -47,17 +47,15 @@ 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::AbstractVector, A::LinearMap, x::AbstractVector) +# 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 == 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!")) + 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)")) + (size(C, 1) != mA || size(C, 2) != nB) && + throw(DimensionMismatch("result has dimensions $(size(C)), needs ($mA,$nB)")) return nothing end @@ -118,6 +116,7 @@ A_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, A, 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 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/left.jl b/src/left.jl new file mode 100644 index 00000000..7561bd1e --- /dev/null +++ b/src/left.jl @@ -0,0 +1,49 @@ +# 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 + +# x = y'*A ⇐⇒ x' = (A'*y) +Base.:(*)(y::AdjointAbsVec, A::LinearMap) = adjoint(A' * y') +Base.:(*)(y::TransposeAbsVec, A::LinearMap) = transpose(transpose(A) * transpose(y)) + +# multiplication with vector/matrix +for Atype in (AdjointAbsVec, Adjoint{<:Any,<:AbstractMatrix}) + @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 + + 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 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 + + 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 new file mode 100644 index 00000000..8f3777ba --- /dev/null +++ b/test/left.jl @@ -0,0 +1,65 @@ +using Test, LinearMaps, LinearAlgebra + +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! + α = rand(T); β = rand(T) + b1 = y' * A + b2 = similar(b1) + 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) + 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 multiplication" begin + T = ComplexF32 + 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([L L]) # BlockMap + + W = LinearMap(randn(T,5,4)) + @test left_tester(W) # WrappedMap +end 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 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")