Use 5-arg mul! for matrices (#68)

* intro MatrixMap, 5-arg mul! for WrappedMaps -> LinAlg

* make adjoint/transpose invariant for MatrixMaps

* inherit properties of MatrixMap, add comparison

* test TransposeMap/AdjointMap based on FunctionMaps

* add specialized matrix-muls for linearcombinations

* introduce FreeMap, test for zero allocations

Author: dkarrasch

src/LinearMaps.jl

@@ -105,11 +105,11 @@ function SparseArrays.sparse(A::LinearMap{T}) where {T}
     return SparseMatrixCSC(M, N, colptr, rowind, nzval)
 end
 
+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
 include("linearcombination.jl") # defining linear combinations of linear maps
 include("composition.jl") # composition of linear maps
-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("functionmap.jl") # using a function as linear map
 include("blockmap.jl") # block linear maps
 

src/linearcombination.jl

@@ -46,6 +46,8 @@ LinearAlgebra.transpose(A::LinearCombination) = LinearCombination{eltype(A)}(map
 LinearAlgebra.adjoint(A::LinearCombination)   = LinearCombination{eltype(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)
@@ -60,6 +62,67 @@ function A_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector)
     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
+        end
+    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(α)
+        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)
+    end
+    return y
+end
+
+end # VERSION
+
 At_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = A_mul_B!(y, transpose(A), x)
 
 Ac_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = A_mul_B!(y, adjoint(A), x)

src/wrappedmap.jl

@@ -17,6 +17,24 @@ function WrappedMap{T}(lmap::Union{AbstractMatrix, LinearMap};
     WrappedMap{T, typeof(lmap)}(lmap, issymmetric, ishermitian, isposdef)
 end
 
+const MatrixMap{T} = WrappedMap{T,<:AbstractMatrix}
+
+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} =
+    WrappedMap{T}(adjoint(A.lmap); issymmetric=A._issymmetric, ishermitian=A._ishermitian, isposdef=A._isposdef)
+
+Base.:(==)(A::MatrixMap, B::MatrixMap) =
+    (eltype(A)==eltype(B) && A.lmap==B.lmap && A._issymmetric==B._issymmetric &&
+     A._ishermitian==B._ishermitian && A._isposdef==B._isposdef)
+
+if VERSION ≥ v"1.3.0-alpha.115"
+
+LinearAlgebra.mul!(y::AbstractVector, A::WrappedMap, x::AbstractVector, α::Number=true, β::Number=false) =
+    mul!(y, A.lmap, x, α, β)
+
+end # VERSION
+
 # properties
 Base.size(A::WrappedMap) = size(A.lmap)
 LinearAlgebra.issymmetric(A::WrappedMap) = A._issymmetric

test/linearcombination.jl

@@ -15,14 +15,27 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools
     @test run(b, samples=5).allocs <= 1
 
     A = 2 * rand(ComplexF64, (10, 10)) .- 1
-    B = rand(size(A)...)
+    B = rand(ComplexF64, size(A)...)
     M = @inferred LinearMap(A)
     N = @inferred LinearMap(B)
     LC = @inferred M + N
     v = rand(ComplexF64, 10)
     w = similar(v)
     b = @benchmarkable mul!($w, $M, $v)
     @test run(b, samples=3).allocs == 0
+    if VERSION >= v"1.3.0-alpha.115"
+        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
+            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*β
+        end
+    end
     # @test_throws ErrorException LinearMaps.LinearCombination{ComplexF64}((M, N), (1, 2, 3))
     @test @inferred size(3M + 2.0N) == size(A)
     # addition
@@ -31,16 +44,16 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools
     @test @inferred convert(Matrix, M + LC) ≈ 2A + B
     @test @inferred convert(Matrix, M + M) ≈ 2A
     # subtraction
-    @test @inferred Matrix(LC - LC) == zeros(size(LC))
-    @test @inferred Matrix(LC - M) == B
-    @test @inferred Matrix(N - LC) == -A
-    @test @inferred Matrix(M - M) == zeros(size(M))
+    @test (@inferred Matrix(LC - LC)) ≈ zeros(eltype(LC), size(LC)) atol=10eps()
+    @test (@inferred Matrix(LC - M)) ≈ B
+    @test (@inferred Matrix(N - LC)) ≈ -A
+    @test (@inferred Matrix(M - M)) ≈ zeros(size(M)) atol=10eps()
     # scalar multiplication
     @test @inferred Matrix(-M) == -A
     @test @inferred Matrix(-LC) == -A - B
     @test @inferred Matrix(3 * M) == 3 * A
     @test @inferred Matrix(M * 3) == 3 * A
-    @test Matrix(3.0 * LC) ≈ Matrix(LC * 3) == 3A + 3B
+    @test Matrix(3.0 * LC) ≈ Matrix(LC * 3) ≈ 3A + 3B
     @test @inferred Matrix(3 \ M) ≈ A/3
     @test @inferred Matrix(M / 3) ≈ A/3
     @test @inferred Matrix(3 \ LC) ≈ (A + B) / 3

test/transpose.jl

@@ -49,4 +49,24 @@ using Test, LinearMaps, LinearAlgebra
     B = @inferred LinearMap(Hermitian(rand(ComplexF64, 10, 10)))
     @test adjoint(B) == B
     @test B == B'
+
+    CS = @inferred LinearMap{ComplexF64}(cumsum, x -> reverse(cumsum(reverse(x))), 10; ismutating=false)
+    for transform in (adjoint, transpose)
+        @test transform(transform(CS)) == CS
+    end
+    @test transpose(CS) != adjoint(CS)
+    @test adjoint(CS) != transpose(CS)
+    M = Matrix(CS)
+    x = rand(ComplexF64, 10)
+    for transform1 in (adjoint, transpose), transform2 in (adjoint, transpose)
+        @test transform2(LinearMap(transform1(CS))) * x ≈ transform2(transform1(M))*x
+    end
+
+    id = @inferred LinearMap(identity, identity, 10; issymmetric=true, ishermitian=true, isposdef=true)
+    for transform in (adjoint, transpose)
+        @test transform(id) == id
+        for prop in (issymmetric, ishermitian, isposdef)
+            @test prop(transform(id))
+        end
+    end
 end