include left matrix multiply and some fixes

Author: dkarrasch

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

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))