diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index c3c6f813..a76e2c2a 100644 --- a/src/LinearMaps.jl +++ b/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 diff --git a/src/linearcombination.jl b/src/linearcombination.jl index 24a7b3b8..55e2e6d3 100644 --- a/src/linearcombination.jl +++ b/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) diff --git a/src/wrappedmap.jl b/src/wrappedmap.jl index 7ba88424..5f919be5 100644 --- a/src/wrappedmap.jl +++ b/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 diff --git a/test/linearcombination.jl b/test/linearcombination.jl index e9767474..d3f0aaca 100644 --- a/test/linearcombination.jl +++ b/test/linearcombination.jl @@ -15,7 +15,7 @@ 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 @@ -23,6 +23,19 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools 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 diff --git a/test/transpose.jl b/test/transpose.jl index 711b64c1..784fefd7 100644 --- a/test/transpose.jl +++ b/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