Add in-place left-multiplication for absmat (#131)

Author: dkarrasch

docs/src/types.md

@@ -102,6 +102,7 @@ Base.:*(::LinearMap,::AbstractMatrix)
 Base.:*(::AbstractMatrix,::LinearMap)
 LinearAlgebra.mul!(::AbstractVecOrMat,::LinearMap,::AbstractVector)
 LinearAlgebra.mul!(::AbstractVecOrMat,::LinearMap,::AbstractVector,::Number,::Number)
+LinearAlgebra.mul!(::AbstractMatrix,::AbstractMatrix,::LinearMap)
 *(::LinearAlgebra.AdjointAbsVec,::LinearMap)
 *(::LinearAlgebra.TransposeAbsVec,::LinearMap)
 ```

src/kronecker.jl

@@ -288,7 +288,7 @@ function _unsafe_mul!(y::AbstractVecOrMat, L::KroneckerSumMap, x::AbstractVector
     mb, nb = size(B)
     X = reshape(x, (nb, na))
     Y = reshape(y, (nb, na))
-    _unsafe_mul!(Y, X, convert(AbstractMatrix, transpose(A)))
+    _unsafe_mul!(Y, X, transpose(A))
     _unsafe_mul!(Y, B, X, true, true)
     return y
 end

src/left.jl

@@ -38,29 +38,61 @@ julia> A=LinearMap([1.0 2.0; 3.0 4.0]); x=[1.0, 1.0]; transpose(x)*A
 Base.:(*)(y::LinearAlgebra.TransposeAbsVec, A::LinearMap) = transpose(transpose(A) * transpose(y))
 
 # multiplication with vector/matrix
-const AdjointAbsVecOrMat{T} = Adjoint{T,<:AbstractVecOrMat}
 const TransposeAbsVecOrMat{T} = Transpose{T,<:AbstractVecOrMat}
 
-function mul!(x::AbstractMatrix, y::AdjointAbsVecOrMat, A::LinearMap)
-    check_dim_mul(x, y, A)
-    _unsafe_mul!(x', A', y')
-    return x
+# handles both y::AbstractMatrix and y::AdjointAbsVecOrMat
+"""
+    mul!(C::AbstractMatrix, A::AbstractMatrix, B::LinearMap) -> C
+
+Calculates the matrix representation of `A*B` and stores the result in `C`,
+overwriting the existing value of `C`. Note that `C` must not be aliased with
+either `A` or `B`. The computation `C = A*B` is performed via `C' = B'A'`.
+
+## Examples
+```jldoctest; setup=(using LinearAlgebra, LinearMaps)
+julia> A=[1.0 1.0; 1.0 1.0]; B=LinearMap([1.0 2.0; 3.0 4.0]); C = similar(A); mul!(C, A, B);
+
+julia> C
+2×2 Array{Float64,2}:
+ 4.0  6.0
+ 4.0  6.0
+```
+"""
+function mul!(X::AbstractMatrix, Y::AbstractMatrix, A::LinearMap)
+    check_dim_mul(X, Y, A)
+    _unsafe_mul!(X', A', Y')
+    return X
+end
+
+function mul!(X::AbstractMatrix, Y::TransposeAbsVecOrMat, A::LinearMap)
+    check_dim_mul(X, Y, A)
+    _unsafe_mul!(transpose(X), transpose(A), transpose(Y))
+    return X
 end
 
-function mul!(x::AbstractMatrix, y::AdjointAbsVecOrMat, A::LinearMap, α::Number, β::Number)
-    check_dim_mul(x, y, A)
-    _unsafe_mul!(x', conj(α)*A', y', true, conj(β))
-    return x
+# commutative case, handles both the abstract and adjoint case
+function mul!(X::AbstractMatrix{<:RealOrComplex}, Y::AbstractMatrix{<:RealOrComplex}, A::LinearMap{<:RealOrComplex},
+                α::RealOrComplex, β::RealOrComplex)
+    check_dim_mul(X, Y, A)
+    _unsafe_mul!(X', A', Y', conj(α), conj(β))
+    return X
 end
 
-function mul!(x::AbstractMatrix, y::TransposeAbsVecOrMat, A::LinearMap)
-    check_dim_mul(x, y, A)
-    _unsafe_mul!(transpose(x), transpose(A), transpose(y))
-    return x
+function mul!(X::AbstractMatrix{<:RealOrComplex}, Y::TransposeAbsVecOrMat{<:RealOrComplex}, A::LinearMap{<:RealOrComplex},
+                α::RealOrComplex, β::RealOrComplex)
+    check_dim_mul(X, Y, A)
+    _unsafe_mul!(transpose(X), transpose(A), transpose(Y), α, β)
+    return X
 end
 
-function mul!(x::AbstractMatrix, y::TransposeAbsVecOrMat, A::LinearMap, α::Number, β::Number)
-    check_dim_mul(x, y, A)
-    _unsafe_mul!(transpose(x), α*transpose(A), transpose(y), true, β)
-    return x
+# non-commutative case
+function mul!(X::AbstractMatrix, Y::AbstractMatrix, A::LinearMap, α::Number, β::Number)
+    check_dim_mul(X, Y, A)
+    if iszero(β)
+        _unsafe_mul!(X', conj(α)*A', Y')
+    else
+        !isone(β) && rmul!(X, β)
+        _unsafe_mul!(X', conj(α)*A', Y', true, true)
+    end
+    return X
 end

src/wrappedmap.jl

@@ -61,6 +61,10 @@ for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractM
     end
 end
 
+mul!(Y::AbstractMatrix, X::AbstractMatrix, A::MatrixMap) = mul!(Y, X, A.lmap)
+# the following method is needed for disambiguation with left-multiplication
+mul!(Y::AbstractMatrix, X::TransposeAbsVecOrMat, A::MatrixMap) = mul!(Y, X, A.lmap)
+
 if VERSION ≥ v"1.3.0-alpha.115"
     for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
         @eval begin
@@ -82,6 +86,18 @@ if VERSION ≥ v"1.3.0-alpha.115"
             end
         end
     end
+
+    mul!(X::AbstractMatrix, Y::AbstractMatrix, A::MatrixMap, α::Number, β::Number) =
+        mul!(X, Y, A.lmap, α, β)
+    # the following method is needed for disambiguation with left-multiplication
+    function mul!(Y::AbstractMatrix{<:RealOrComplex}, X::AbstractMatrix{<:RealOrComplex}, A::MatrixMap{<:RealOrComplex},
+                    α::RealOrComplex, β::RealOrComplex)
+        return mul!(Y, X, A.lmap, α, β)
+    end
+    function mul!(Y::AbstractMatrix{<:RealOrComplex}, X::TransposeAbsVecOrMat{<:RealOrComplex}, A::MatrixMap{<:RealOrComplex},
+                    α::RealOrComplex, β::RealOrComplex)
+        return mul!(Y, X, A.lmap, α, β)
+    end
 end # VERSION
 
 # combine LinearMap and Matrix objects: linear combinations and map composition
@@ -112,7 +128,8 @@ Base.:(*)(A₁::LinearMap, A₂::AbstractMatrix) = *(A₁, WrappedMap(A₂))
     *(X::AbstractMatrix, A::LinearMap)::CompositeMap
 
 Return the `CompositeMap` `LinearMap(X)*A`, interpreting the matrix `X` as a linear
-operator. To compute the right-action of `A` on each row of `X`, call `Matrix(X*A)`.
+operator. To compute the right-action of `A` on each row of `X`, call `Matrix(X*A)`
+or `mul!(Y, X, A)` for the in-place version.
 
 ## Examples
 ```jldoctest; setup=(using LinearMaps)

test/numbertypes.jl

@@ -1,28 +1,41 @@
 using Test, LinearMaps, LinearAlgebra, Quaternions
 
 # type piracy because Quaternions.jl doesn't have it right
-Base.:(*)(z::Complex{T}, q::Quaternion{T}) where {T<:Real} = quat(z) * q
-Base.:(*)(q::Quaternion{T}, z::Complex{T}) where {T<:Real} = q * quat(z)
+Base.:(*)(z::Complex, q::Quaternion) = quat(z) * q
+Base.:(*)(q::Quaternion, z::Complex) = q * quat(z)
+Base.:(+)(q::Quaternion, z::Complex) = q + quat(z)
 
 @testset "noncommutative number type" begin
     x = Quaternion.(rand(10), rand(10), rand(10), rand(10))
     v = rand(10)
     A = Quaternion.(rand(10,10), rand(10,10), rand(10,10), rand(10,10))
     B = rand(ComplexF64, 10, 10)
+    C = similar(A)
     γ = Quaternion.(rand(4)...) # "Number"
     α = UniformScaling(γ)
     β = UniformScaling(Quaternion.(rand(4)...))
     λ = rand(ComplexF64)
     L = LinearMap(A)
-    @test Array(L) == A
-    @test Array(L') == A'
-    @test Array(transpose(L)) == transpose(A)
-    @test Array(α * L) == α * A
-    @test Array(L * α) == A * α
-    @test Array(α * L) == α * A
-    @test Array(L * α ) == A * α
-    @test Array(α * L') == α * A'
-    @test Array((α * L')') ≈ (α * A')' ≈ A * conj(α)
+    F = LinearMap{eltype(A)}(x -> A*x, y -> A'y, 10)
+    @test Array(F) == A
+    @test Array(F') == A'
+    @test Array(transpose(F)) == transpose(A)
+    @test Array(α * F) == α * A
+    @test Array(F * α) == A * α
+    @test Array(α * F) == α * A
+    @test Array(F * α ) == A * α
+    @test Array(α * F') == α * A'
+    for M in (L, F)
+        @test mul!(C, transpose(A), M) ≈ transpose(A)*A
+        @test mul!(C, A', M) ≈ A'A
+        @test mul!(C, A, M) ≈ A*A
+        @test mul!(copy(C), M, A, γ, λ) ≈ A*A*γ + C*λ
+        @test mul!(copy(C), A, M, γ, λ) ≈ A*A*γ + C*λ
+        @test mul!(copy(C), A, M, γ, 0) ≈ A*A*γ
+        @test mul!(copy(C), transpose(A), M, γ, λ) ≈ transpose(A)*A*γ + C*λ
+        @test mul!(copy(C), adjoint(A), M, γ, λ) ≈ A'*A*γ + C*λ
+    end
+    @test Array((α * F')') ≈ (α * A')' ≈ A * conj(α)
     @test L * x ≈ A * x
     @test L' * x ≈ A' * x
     @test α * (L * x) ≈ α * (A * x)

test/wrappedmap.jl

@@ -72,4 +72,15 @@ using Test, LinearMaps, LinearAlgebra
     B = @inferred LinearMap(Hermitian(rand(ComplexF64, 10, 10)))
     @test adjoint(B) == B
     @test B == B'
+
+    N = 10
+    Id = LinearMap(identity, N; issymmetric=true)
+    A = rand(N, N)
+    B = similar(A)
+    @test mul!(copy(B), Id, A) == A
+    @test mul!(B, A, Id) == B == A
+    @test mul!(copy(B), A, LinearMap(Matrix(Id))) == A
+    @test mul!(B, Id, A, true, true) ≈ 2A
+    @test mul!(B, A, Id, true, true) == B == 3A
+    @test mul!(copy(B), A, LinearMap(Matrix(Id)), true, true) == 4A
 end