Overhaul of conversion methods (#88)

Author: dkarrasch

src/conversion.jl

@@ -4,7 +4,7 @@ function Base.Matrix(A::LinearMap)
     T = eltype(A)
     mat = Matrix{T}(undef, (M, N))
     v = fill(zero(T), N)
-    @inbounds for i = 1:N
+    @inbounds for i in 1:N
         v[i] = one(T)
         mul!(view(mat, :, i), A, v)
         v[i] = zero(T)
@@ -17,35 +17,6 @@ Base.convert(::Type{Array}, A::LinearMap) = convert(Matrix, A)
 Base.convert(::Type{AbstractMatrix}, A::LinearMap) = convert(Matrix, A)
 Base.convert(::Type{AbstractArray}, A::LinearMap) = convert(AbstractMatrix, A)
 
-# special cases
-Base.convert(::Type{AbstractMatrix}, A::WrappedMap) = convert(AbstractMatrix, A.lmap)
-Base.convert(::Type{Matrix}, A::WrappedMap) = convert(Matrix, A.lmap)
-function Base.convert(::Type{Matrix}, ΣA::LinearCombination{<:Any,<:Tuple{Vararg{MatrixMap}}})
-    if length(ΣA.maps) <= 10
-        return (+).(map(A->getfield(A, :lmap), ΣA.maps)...)
-    else
-        S = zero(first(ΣA.maps).lmap)
-        for A in ΣA.maps
-            S .+= A.lmap
-        end
-        return S
-    end
-end
-function Base.convert(::Type{Matrix}, AB::CompositeMap{<:Any,<:Tuple{MatrixMap,MatrixMap}})
-    B, A = AB.maps
-    return A.lmap*B.lmap
-end
-function Base.convert(::Type{Matrix}, λA::CompositeMap{<:Any,<:Tuple{MatrixMap,UniformScalingMap}})
-    A, J = λA.maps
-    return J.λ*A.lmap
-end
-function Base.convert(::Type{Matrix}, Aλ::CompositeMap{<:Any,<:Tuple{UniformScalingMap,MatrixMap}})
-    J, A = Aλ.maps
-    return A.lmap*J.λ
-end
-
-Base.Matrix(A::BlockMap) = hvcat(A.rows, convert.(AbstractMatrix, A.maps)...)
-
 # sparse: create sparse matrix representation of LinearMap
 function SparseArrays.sparse(A::LinearMap{T}) where {T}
     M, N = size(A)
@@ -55,7 +26,7 @@ function SparseArrays.sparse(A::LinearMap{T}) where {T}
     v = fill(zero(T), N)
     Av = Vector{T}(undef, M)
 
-    for i = 1:N
+    @inbounds for i in 1:N
         v[i] = one(T)
         mul!(Av, A, v)
         js = findall(!iszero, Av)
@@ -70,7 +41,92 @@ function SparseArrays.sparse(A::LinearMap{T}) where {T}
 
     return SparseMatrixCSC(M, N, colptr, rowind, nzval)
 end
-
 Base.convert(::Type{SparseMatrixCSC}, A::LinearMap) = sparse(A)
+
+# special cases
+
+# UniformScalingMap
+Base.Matrix(J::UniformScalingMap) = Matrix(J.λ*I, size(J))
+Base.convert(::Type{AbstractMatrix}, J::UniformScalingMap) = Diagonal(fill(J.λ, size(J, 1)))
+
+# WrappedMap
+Base.Matrix(A::WrappedMap) = Matrix(A.lmap)
+Base.convert(::Type{AbstractMatrix}, A::WrappedMap) = convert(AbstractMatrix, A.lmap)
+Base.convert(::Type{Matrix}, A::WrappedMap) = convert(Matrix, A.lmap)
+SparseArrays.sparse(A::WrappedMap) = sparse(A.lmap)
 Base.convert(::Type{SparseMatrixCSC}, A::WrappedMap) = convert(SparseMatrixCSC, A.lmap)
-Base.convert(::Type{SparseMatrixCSC}, A::BlockMap) = hvcat(A.rows, convert.(SparseMatrixCSC, A.maps)...)
+
+# TransposeMap & AdjointMap
+for (TT, T) in ((AdjointMap, adjoint), (TransposeMap, transpose))
+    @eval Base.convert(::Type{AbstractMatrix}, A::$TT) = $T(convert(AbstractMatrix, A.lmap))
+    @eval SparseArrays.sparse(A::$TT) = $T(convert(SparseMatrixCSC, A.lmap))
+end
+
+# LinearCombination
+for (TT, T) in ((Type{Matrix}, Matrix), (Type{SparseMatrixCSC}, SparseMatrixCSC))
+    @eval function Base.convert(::$TT, ΣA::LinearCombination{<:Any,<:Tuple{Vararg{MatrixMap}}})
+        maps = ΣA.maps
+        mats = map(A->getfield(A, :lmap), maps)
+        return convert($T, sum(mats))
+    end
+end
+
+# CompositeMap
+for (TT, T) in ((Type{Matrix}, Matrix), (Type{SparseMatrixCSC}, SparseMatrixCSC))
+    @eval function Base.convert(::$TT, AB::CompositeMap{<:Any,<:Tuple{MatrixMap,MatrixMap}})
+        B, A = AB.maps
+        return convert($T, A.lmap*B.lmap)
+    end
+    @eval function Base.convert(::$TT, λA::CompositeMap{<:Any,<:Tuple{MatrixMap,UniformScalingMap}})
+        A, J = λA.maps
+        return convert($T, J.λ*A.lmap)
+    end
+    @eval function Base.convert(::$TT, Aλ::CompositeMap{<:Any,<:Tuple{UniformScalingMap,MatrixMap}})
+        J, A = Aλ.maps
+        return convert($T, A.lmap*J.λ)
+    end
+end
+
+# BlockMap & BlockDiagonalMap
+Base.Matrix(A::BlockMap) = hvcat(A.rows, convert.(Matrix, A.maps)...)
+Base.convert(::Type{AbstractMatrix}, A::BlockMap) = hvcat(A.rows, convert.(AbstractMatrix, A.maps)...)
+function Base.convert(::Type{SparseMatrixCSC}, A::BlockMap)
+    return hvcat(
+        A.rows,
+        convert(SparseMatrixCSC, first(A.maps)),
+        convert.(AbstractMatrix, Base.tail(A.maps))...
+    )
+end
+Base.Matrix(A::BlockDiagonalMap) = cat(convert.(Matrix, A.maps)...; dims=(1,2))
+Base.convert(::Type{AbstractMatrix}, A::BlockDiagonalMap) = sparse(A)
+function SparseArrays.sparse(A::BlockDiagonalMap)
+    return blockdiag(convert.(SparseMatrixCSC, A.maps)...)
+end
+
+# KroneckerMap & KroneckerSumMap
+Base.Matrix(A::KroneckerMap) = kron(convert.(Matrix, A.maps)...)
+Base.convert(::Type{AbstractMatrix}, A::KroneckerMap) = kron(convert.(AbstractMatrix, A.maps)...)
+function SparseArrays.sparse(A::KroneckerMap)
+    return kron(
+        convert(SparseMatrixCSC, first(A.maps)),
+        convert.(AbstractMatrix, Base.tail(A.maps))...
+    )
+end
+function Base.Matrix(L::KroneckerSumMap)
+    A, B = L.maps
+    IA = Diagonal(ones(Bool, size(A, 1)))
+    IB = Diagonal(ones(Bool, size(B, 1)))
+    return kron(Matrix(A), IB) + kron(IA, Matrix(B))
+end
+function Base.convert(::Type{AbstractMatrix}, L::KroneckerSumMap)
+    A, B = L.maps
+    IA = Diagonal(ones(Bool, size(A, 1)))
+    IB = Diagonal(ones(Bool, size(B, 1)))
+    return kron(convert(AbstractMatrix, A), IB) + kron(IA, convert(AbstractMatrix, B))
+end
+function SparseArrays.sparse(L::KroneckerSumMap)
+    A, B = L.maps
+    IA = sparse(Diagonal(ones(Bool, size(A, 1))))
+    IB = sparse(Diagonal(ones(Bool, size(B, 1))))
+    return kron(convert(AbstractMatrix, A), IB) + kron(IA, convert(AbstractMatrix, B))
+end

src/uniformscalingmap.jl

@@ -33,85 +33,57 @@ Base.:(*)(J::UniformScalingMap, α::Number) = UniformScalingMap(J.λ * α, size(
 
 # multiplication with vector
 Base.:(*)(A::UniformScalingMap, x::AbstractVector) =
-    length(x) == A.M ? A.λ * x : throw(DimensionMismatch("A_mul_B!"))
+    length(x) == A.M ? A.λ * x : throw(DimensionMismatch("*"))
 
-if VERSION < v"1.3.0-alpha.115"
-function A_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector)
-    (length(x) == length(y) == A.M || throw(DimensionMismatch("A_mul_B!")))
-    if iszero(A.λ)
-        return fill!(y, zero(eltype(y)))
-    elseif isone(A.λ)
-        return copyto!(y, x)
-    else
-        # call of LinearAlgebra.generic_mul! since order of arguments in mul! in
-        # stdlib/LinearAlgebra/src/generic.jl reversed
-        return LinearAlgebra.generic_mul!(y, A.λ, x)
+# multiplication with vector/matrix
+for Atype in (AbstractVector, AbstractMatrix)
+    @eval Base.@propagate_inbounds function LinearAlgebra.mul!(y::$Atype, J::UniformScalingMap, x::$Atype,
+                α::Number=true, β::Number=false)
+        @boundscheck check_dim_mul(y, J, x)
+        _scaling!(y, J.λ, x, α, β)
+        return y
     end
 end
-else # 5-arg mul! exists and order of arguments is corrected
-
-A_mul_B!(y::AbstractVector, J::UniformScalingMap, x::AbstractVector) = mul!(y, J, x)
-
-end # VERSION
-
-@inline function LinearAlgebra.mul!(y::AbstractVector, J::UniformScalingMap, x::AbstractVector, α::Number=true, β::Number=false)
-    @boundscheck (length(x) == length(y) == J.M || throw(DimensionMismatch("mul!")))
-    _scaling!(y, J, x, α, β)
-    return y
-end
 
-@inline function LinearAlgebra.mul!(Y::AbstractMatrix, J::UniformScalingMap, X::AbstractMatrix, α::Number=true, β::Number=false)
-    @boundscheck size(X) == size(Y) || throw(DimensionMismatch("mul!"))
-    @boundscheck size(X,1) == J.M || throw(DimensionMismatch("mul!"))
-    _scaling!(Y, J, X, α, β)
-    return Y
-end
-
-function _scaling!(y, J::UniformScalingMap, x, α::Number=true, β::Number=false)
-    λ = J.λ
-    if isone(α)
+function _scaling!(y, λ::Number, x, α::Number, β::Number)
+    if (iszero(α) || iszero(λ))
+        iszero(β) && (fill!(y, zero(eltype(y))); return y)
+        isone(β) && return y
+        # β != 0, 1
+        rmul!(y, β)
+        return y
+    elseif isone(α)
         if iszero(β)
-            iszero(λ) && return fill!(y, zero(eltype(y)))
             isone(λ) && return copyto!(y, x)
             y .= λ .* x
             return y
         elseif isone(β)
-            iszero(λ) && return y
             isone(λ) && return y .+= x
             y .+= λ .* x
             return y
         else # β != 0, 1
-            iszero(λ) && (rmul!(y, β); return y)
-            isone(λ) && (y .= y .* β .+ x; return y)
+            isone(λ) && (axpby!(one(eltype(x)), x, β, y); return y)
             y .= y .* β .+ λ .* x
             return y
         end
-    elseif iszero(α)
-        iszero(β) && (fill!(y, zero(eltype(y))); return y)
-        isone(β) && return y
-        # β != 0, 1
-        rmul!(y, β)
-        return y
     else # α != 0, 1
         if iszero(β)
-            iszero(λ) && return fill!(y, zero(eltype(y)))
-            isone(λ) && return y .= x .* α
+            isone(λ) && (y .= x .* α; return y)
             y .= λ .* x .* α
             return y
         elseif isone(β)
-            iszero(λ) && return y
-            isone(λ) && return y .+= x .* α
+            isone(λ) && (axpby!(α, x, β, y); return y)
             y .+= λ .* x .* α
             return y
         else # β != 0, 1
-            iszero(λ) && (rmul!(y, β); return y)
             isone(λ) && (y .= y .* β .+ x .* α; return y)
             y .= y .* β .+ λ .* x .* α
             return y
-        end
-    end
-end
+        end # β-cases
+    end # α-cases
+end # function _scaling!
 
+A_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector) = mul!(y, A, x)
 At_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x)
 Ac_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x)
 

test/blockmap.jl

@@ -71,7 +71,8 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays
             @test L isa LinearMaps.BlockMap{elty}
             @test size(L) == size(A)
             @test L * x ≈ A * x
-            @test Matrix(L) ≈ A
+            @test Matrix(L) == A
+            @test convert(AbstractMatrix, L) == A
             A = [I A12; A21 I]
             @inferred hvcat((2,2), I, LinearMap(A12), LinearMap(A21), I)
             L = @inferred hvcat((2,2), I, LinearMap(A12), LinearMap(A21), I)
@@ -136,13 +137,17 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays
             @test L isa LinearMaps.LinearMap{elty}
             @test size(L) == size(A)
             @test L * x ≈ A * x
-            @test Matrix(L) ≈ A
+            @test Matrix(L) == A
+            @test convert(AbstractMatrix, L) == A
+            @test sparse(L) == sparse(A)
             Lt = @inferred transform(L)
             @test Lt isa LinearMaps.LinearMap{elty}
             @test Lt * x ≈ transform(A) * x
+            @test convert(AbstractMatrix, Lt) == transform(A)
+            @test sparse(transform(L)) == transform(A)
             Lt = @inferred transform(LinearMap(L))
             @test Lt * x ≈ transform(A) * x
-            @test Matrix(Lt) ≈ Matrix(transform(A))
+            @test Matrix(Lt) == Matrix(transform(A))
             A21 = rand(elty, 10, 10)
             A = [I A12; A21 I]
             L = [I LinearMap(A12); LinearMap(A21) I]
@@ -170,6 +175,9 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays
             Md = Matrix(blockdiag(sparse.((M1, M2, M3, M2, M1))...))
             x = randn(elty, size(Md, 2))
             Bd = @inferred blockdiag(L1, L2, L3, L2, L1)
+            @test Matrix(Bd) == Md
+            @test convert(AbstractMatrix, Bd) isa SparseMatrixCSC
+            @test sparse(Bd) == Md
             @test Matrix(@inferred blockdiag(L1)) == M1
             @test Matrix(@inferred blockdiag(L1, L2)) == blockdiag(sparse.((M1, M2))...)
             Bd2 = @inferred cat(L1, L2, L3, L2, L1; dims=(1,2))

test/composition.jl

@@ -1,18 +1,5 @@
 using Test, LinearMaps, LinearAlgebra
 
-# new type
-struct SimpleFunctionMap <: LinearMap{Float64}
-    f::Function
-    N::Int
-end
-struct SimpleComplexFunctionMap <: LinearMap{Complex{Float64}}
-    f::Function
-    N::Int
-end
-Base.size(A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}) = (A.N, A.N)
-Base.:(*)(A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}, v::Vector) = A.f(v)
-LinearAlgebra.mul!(y::Vector, A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}, x::Vector) = copyto!(y, *(A, x))
-
 @testset "composition" begin
     F = @inferred LinearMap(cumsum, y -> reverse(cumsum(reverse(x))), 10; ismutating=false)
     FC = @inferred LinearMap{ComplexF64}(cumsum, y -> reverse(cumsum(reverse(x))), 10; ismutating=false)
@@ -66,24 +53,6 @@ LinearAlgebra.mul!(y::Vector, A::Union{SimpleFunctionMap,SimpleComplexFunctionMa
     mul!(w, L, v)
     @test w ≈ LF * v
 
-    # test new type
-    F = SimpleFunctionMap(cumsum, 10)
-    FC = SimpleComplexFunctionMap(cumsum, 10)
-    @test @inferred ndims(F) == 2
-    @test @inferred size(F, 1) == 10
-    @test @inferred length(F) == 100
-    @test @inferred !issymmetric(F)
-    @test @inferred !ishermitian(F)
-    @test @inferred !ishermitian(FC)
-    @test @inferred !isposdef(F)
-    w = similar(v)
-    mul!(w, F, v)
-    @test w == F * v
-    @test_throws MethodError F' * v
-    @test_throws MethodError transpose(F) * v
-    @test_throws MethodError mul!(w, adjoint(F), v)
-    @test_throws MethodError mul!(w, transpose(F), v)
-
     # test composition of several maps with shared data #31
     global sizes = ( (5, 2), (3, 3), (3, 2), (2, 2), (9, 2), (7, 1) )
     N = length(sizes) - 1

test/conversion.jl

@@ -24,6 +24,12 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays
     @test Matrix(transpose(M)) == copy(transpose(A))
     @test convert(AbstractMatrix, M') isa Adjoint
     @test convert(Matrix, M*3I) == A*3
+    @test convert(Matrix, M+M) == A + A
+
+    # UniformScalingMap
+    J = LinearMap(α*I, 10)
+    JM = convert(AbstractMatrix, J)
+    @test JM == Diagonal(fill(α, 10))
 
     # sparse matrix generation/conversion
     @test sparse(M) == sparse(Array(M))