Allow BlockMap construction with matrices (#71)

Author: dkarrasch

Project.toml

@@ -11,10 +11,11 @@ julia = "1"
 
 [extras]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["LinearAlgebra", "SparseArrays", "Test", "BenchmarkTools", "Quaternions"]
+test = ["LinearAlgebra", "SparseArrays", "Test", "BenchmarkTools", "InteractiveUtils", "Quaternions"]

src/blockmap.jl

@@ -16,11 +16,6 @@ BlockMap{T}(maps::As, rows::S) where {T,As<:Tuple{Vararg{LinearMap}},S} = BlockM
 
 MulStyle(A::BlockMap) = MulStyle(A.maps...)
 
-function check_dim(A::LinearMap, dim, n)
-    n == size(A, dim) || throw(DimensionMismatch("Expected $n, got $(size(A, dim))"))
-    return nothing
-end
-
 """
     rowcolranges(maps, rows)
 
@@ -51,31 +46,14 @@ end
 
 Base.size(A::BlockMap) = (last(last(A.rowranges)), last(last(A.colranges)))
 
-############
-# concatenation
-############
-
-for k in 1:8 # is 8 sufficient?
-    Is = ntuple(n->:($(Symbol(:A,n))::UniformScaling), Val(k-1))
-    L = :($(Symbol(:A,k))::LinearMap)
-    args = ntuple(n->Symbol(:A,n), Val(k))
-
-    @eval begin
-        Base.hcat($(Is...), $L, As::Union{LinearMap,UniformScaling}...) = _hcat($(args...), As...)
-        Base.vcat($(Is...), $L, As::Union{LinearMap,UniformScaling}...) = _vcat($(args...), As...)
-        Base.hvcat(rows::Tuple{Vararg{Int}}, $(Is...), $L, As::Union{LinearMap,UniformScaling}...) = _hvcat(rows, $(args...), As...)
-    end
-end
-
 ############
 # hcat
 ############
 """
-    hcat(As::Union{LinearMap,UniformScaling}...)::BlockMap
+    hcat(As::Union{LinearMap,UniformScaling,AbstractVecOrMat}...)::BlockMap
 
 Construct a (lazy) representation of the horizontal concatenation of the arguments.
-`UniformScaling` objects are promoted to `LinearMap` automatically. To avoid fallback
-to the generic `Base.hcat`, there must be a `LinearMap` object among the first 8 arguments.
+All arguments are promoted to `LinearMap`s automatically.
 
 # Examples
 ```jldoctest; setup=(using LinearMaps)
@@ -90,9 +68,7 @@ julia> L * ones(Int, 6)
  6
 ```
 """
-Base.hcat
-
-function _hcat(As::Union{LinearMap,UniformScaling}...)
+function Base.hcat(As::Union{LinearMap,UniformScaling,AbstractVecOrMat}...)
     T = promote_type(map(eltype, As)...)
     nbc = length(As)
 
@@ -112,11 +88,10 @@ end
 # vcat
 ############
 """
-    vcat(As::Union{LinearMap,UniformScaling}...)::BlockMap
+    vcat(As::Union{LinearMap,UniformScaling,AbstractVecOrMat}...)::BlockMap
 
 Construct a (lazy) representation of the vertical concatenation of the arguments.
-`UniformScaling` objects are promoted to `LinearMap` automatically. To avoid fallback
-to the generic `Base.vcat`, there must be a `LinearMap` object among the first 8 arguments.
+All arguments are promoted to `LinearMap`s automatically.
 
 # Examples
 ```jldoctest; setup=(using LinearMaps)
@@ -134,9 +109,7 @@ julia> L * ones(Int, 3)
  3
 ```
 """
-Base.vcat
-
-function _vcat(As::Union{LinearMap,UniformScaling}...)
+function Base.vcat(As::Union{LinearMap,UniformScaling,AbstractVecOrMat}...)
     T = promote_type(map(eltype, As)...)
     nbr = length(As)
 
@@ -157,12 +130,11 @@ end
 # hvcat
 ############
 """
-    hvcat(rows::Tuple{Vararg{Int}}, As::Union{LinearMap,UniformScaling}...)::BlockMap
+    hvcat(rows::Tuple{Vararg{Int}}, As::Union{LinearMap,UniformScaling,AbstractVecOrMat}...)::BlockMap
 
 Construct a (lazy) representation of the horizontal-vertical concatenation of the arguments.
 The first argument specifies the number of arguments to concatenate in each block row.
-`UniformScaling` objects are promoted to `LinearMap` automatically. To avoid fallback
-to the generic `Base.hvcat`, there must be a `LinearMap` object among the first 8 arguments.
+All arguments are promoted to `LinearMap`s automatically.
 
 # Examples
 ```jldoctest; setup=(using LinearMaps)
@@ -185,7 +157,7 @@ julia> L * ones(Int, 6)
 """
 Base.hvcat
 
-function _hvcat(rows::Tuple{Vararg{Int}}, As::Union{LinearMap,UniformScaling}...)
+function Base.hvcat(rows::Tuple{Vararg{Int}}, As::Union{LinearMap,UniformScaling,AbstractVecOrMat}...)
     nr = length(rows)
     T = promote_type(map(eltype, As)...)
     sum(rows) == length(As) || throw(ArgumentError("mismatch between row sizes and number of arguments"))
@@ -237,6 +209,13 @@ function _hvcat(rows::Tuple{Vararg{Int}}, As::Union{LinearMap,UniformScaling}...
     return BlockMap{T}(promote_to_lmaps(n, 1, 1, As...), rows)
 end
 
+function check_dim(A, dim, n)
+    n == size(A, dim) || throw(DimensionMismatch("Expected $n, got $(size(A, dim))"))
+    return nothing
+end
+
+promote_to_lmaps_(n::Int, dim, A::AbstractMatrix) = (check_dim(A, dim, n); LinearMap(A))
+promote_to_lmaps_(n::Int, dim, A::AbstractVector) = (check_dim(A, dim, n); LinearMap(reshape(A, length(A), 1)))
 promote_to_lmaps_(n::Int, dim, J::UniformScaling) = UniformScalingMap(J.λ, n)
 promote_to_lmaps_(n::Int, dim, A::LinearMap) = (check_dim(A, dim, n); A)
 promote_to_lmaps(n, k, dim) = ()
@@ -292,11 +271,6 @@ end
 
 Base.:(==)(A::BlockMap, B::BlockMap) = (eltype(A) == eltype(B) && A.maps == B.maps && A.rows == B.rows)
 
-# special transposition behavior
-
-LinearAlgebra.transpose(A::BlockMap) = TransposeMap(A)
-LinearAlgebra.adjoint(A::BlockMap)  = AdjointMap(A)
-
 ############
 # multiplication helper functions
 ############
@@ -310,6 +284,7 @@ function _blockmul!(y, A::BlockMap, x, α, β)
     return __blockmul!(MulStyle(A), y, A, x, α, β)
 end
 
+# provide one global intermediate storage vector if necessary
 __blockmul!(::FiveArg, y, A, x, α, β)  = ___blockmul!(y, A, x, α, β, nothing)
 __blockmul!(::ThreeArg, y, A, x, α, β) = ___blockmul!(y, A, x, α, β, similar(y))
 
@@ -401,7 +376,6 @@ for (intype, outtype) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, A
             function _unsafe_mul!(y::$outtype, wrapA::$maptype, x::$intype,
                             α::Number, β::Number)
                 require_one_based_indexing(y, x)
-
                 return _transblockmul!(y, wrapA.lmap, x, α, β, $transform)
             end
         end
@@ -437,22 +411,24 @@ BlockDiagonalMap{T}(maps::As) where {T,As<:Tuple{Vararg{LinearMap}}} =
 BlockDiagonalMap(maps::LinearMap...) =
     BlockDiagonalMap{promote_type(map(eltype, maps)...)}(maps)
 
+# since the below methods are more specific than the Base method,
+# they would redefine Base/SparseArrays behavior
 for k in 1:8 # is 8 sufficient?
-    Is = ntuple(n->:($(Symbol(:A,n))::AbstractMatrix), Val(k-1))
+    Is = ntuple(n->:($(Symbol(:A,n))::AbstractVecOrMat), Val(k-1))
     # yields (:A1, :A2, :A3, ..., :A(k-1))
     L = :($(Symbol(:A,k))::LinearMap)
     # yields :Ak
-    mapargs = ntuple(n -> :(LinearMap($(Symbol(:A,n)))), Val(k-1))
+    mapargs = ntuple(n -> :($(Symbol(:A,n))), Val(k-1))
     # yields (:LinearMap(A1), :LinearMap(A2), ..., :LinearMap(A(k-1)))
 
     @eval begin
-        function SparseArrays.blockdiag($(Is...), $L, As::Union{LinearMap,AbstractMatrix}...)
-            return BlockDiagonalMap($(mapargs...), $(Symbol(:A,k)), convert_to_lmaps(As...)...)
+        function SparseArrays.blockdiag($(Is...), $L, As::Union{LinearMap,AbstractVecOrMat}...)
+            return BlockDiagonalMap(convert_to_lmaps($(mapargs...))..., $(Symbol(:A,k)), convert_to_lmaps(As...)...)
         end
 
-        function Base.cat($(Is...), $L, As::Union{LinearMap,AbstractMatrix}...; dims::Dims{2})
+        function Base.cat($(Is...), $L, As::Union{LinearMap,AbstractVecOrMat}...; dims::Dims{2})
             if dims == (1,2)
-                return BlockDiagonalMap($(mapargs...), $(Symbol(:A,k)), convert_to_lmaps(As...)...)
+                return BlockDiagonalMap(convert_to_lmaps($(mapargs...))..., $(Symbol(:A,k)), convert_to_lmaps(As...)...)
             else
                 throw(ArgumentError("dims keyword in cat of LinearMaps must be (1,2)"))
             end
@@ -461,7 +437,7 @@ for k in 1:8 # is 8 sufficient?
 end
 
 """
-    blockdiag(As::Union{LinearMap,AbstractMatrix}...)::BlockDiagonalMap
+    blockdiag(As::Union{LinearMap,AbstractVecOrMat}...)::BlockDiagonalMap
 
 Construct a (lazy) representation of the diagonal concatenation of the arguments.
 To avoid fallback to the generic `SparseArrays.blockdiag`, there must be a `LinearMap`
@@ -470,7 +446,7 @@ object among the first 8 arguments.
 SparseArrays.blockdiag
 
 """
-    cat(As::Union{LinearMap,AbstractMatrix}...; dims=(1,2))::BlockDiagonalMap
+    cat(As::Union{LinearMap,AbstractVecOrMat}...; dims=(1,2))::BlockDiagonalMap
 
 Construct a (lazy) representation of the diagonal concatenation of the arguments.
 To avoid fallback to the generic `Base.cat`, there must be a `LinearMap`

src/conversion.jl

@@ -132,7 +132,7 @@ function SparseArrays.sparse(A::BlockMap)
         convert.(AbstractMatrix, Base.tail(A.maps))...
     )
 end
-Base.Matrix{T}(A::BlockDiagonalMap) where {T} = cat(convert.(Matrix{T}, A.maps)...; dims=(1,2))
+Base.Matrix{T}(A::BlockDiagonalMap) where {T} = Base._cat((1,2), convert.(Matrix{T}, A.maps)...)
 Base.convert(::Type{AbstractMatrix}, A::BlockDiagonalMap) = sparse(A)
 function SparseArrays.sparse(A::BlockDiagonalMap)
     return blockdiag(convert.(SparseMatrixCSC, A.maps)...)

test/blockmap.jl

@@ -1,10 +1,11 @@
-using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools
+using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools, InteractiveUtils
 
 @testset "block maps" begin
     @testset "hcat" begin
         for elty in (Float32, ComplexF64), n2 = (0, 20)
             A11 = rand(elty, 10, 10)
             A12 = rand(elty, 10, n2)
+            v = rand(elty, 10)
             L = @inferred hcat(LinearMap(A11), LinearMap(A12))
             @test @inferred(LinearMaps.MulStyle(L)) === matrixstyle
             @test L isa LinearMaps.BlockMap{elty}
@@ -16,10 +17,14 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools
             L = @inferred hcat(LinearMap(A11), LinearMap(A12), LinearMap(A11))
             A = [A11 A12 A11]
             @test Matrix(L) ≈ A
-            A = [I I I A11 A11 A11]
-            L = @inferred hcat(I, I, I, LinearMap(A11), LinearMap(A11), LinearMap(A11))
-            @test L == [I I I LinearMap(A11) LinearMap(A11) LinearMap(A11)]
-            x = rand(elty, 60)
+            A = [I I I A11 A11 A11 v]
+            @test (@which [A11 A11 A11]).module != LinearMaps
+            @test (@which [I I I A11 A11 A11]).module != LinearMaps
+            @test (@which hcat(I, I, I)).module != LinearMaps
+            @test (@which hcat(I, I, I, LinearMap(A11), A11, A11)).module == LinearMaps
+            L = @inferred hcat(I, I, I, LinearMap(A11), A11, A11, v)
+            @test L == [I I I LinearMap(A11) LinearMap(A11) LinearMap(A11) LinearMap(reshape(v, :, 1))]
+            x = rand(elty, 61)
             @test L isa LinearMaps.BlockMap{elty}
             @test L * x ≈ A * x
             A11 = rand(elty, 11, 10)
@@ -31,21 +36,24 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools
     @testset "vcat" begin
         for elty in (Float32, ComplexF64)
             A11 = rand(elty, 10, 10)
+            v = rand(elty, 10)
             L = @inferred vcat(LinearMap(A11))
             @test L == [LinearMap(A11);]
             @test Matrix(L) ≈ A11
             A21 = rand(elty, 20, 10)
             L = @inferred vcat(LinearMap(A11), LinearMap(A21))
             @test L isa LinearMaps.BlockMap{elty}
             @test @inferred(LinearMaps.MulStyle(L)) === matrixstyle
+            @test (@which [A11; A21]).module != LinearMaps
             A = [A11; A21]
             x = rand(10)
             @test size(L) == size(A)
             @test Matrix(L) ≈ A
             @test L * x ≈ A * x
-            A = [I; I; I; A11; A11; A11]
-            L = @inferred vcat(I, I, I, LinearMap(A11), LinearMap(A11), LinearMap(A11))
-            @test L == [I; I; I; LinearMap(A11); LinearMap(A11); LinearMap(A11)]
+            A = [I; I; I; A11; A11; A11; v v v v v v v v v v]
+            @test (@which [I; I; I; A11; A11; A11; v v v v v v v v v v]).module != LinearMaps
+            L = @inferred vcat(I, I, I, LinearMap(A11), LinearMap(A11), LinearMap(A11), reduce(hcat, fill(v, 10)))
+            @test L == [I; I; I; LinearMap(A11); LinearMap(A11); LinearMap(A11); reduce(hcat, fill(v, 10))]
             x = rand(elty, 10)
             @test L isa LinearMaps.BlockMap{elty}
             @test L * x ≈ A * x
@@ -62,6 +70,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools
             A21 = rand(elty, 20, 10)
             A22 = rand(elty, 20, 20)
             A = [A11 A12; A21 A22]
+            @test (@which [A11 A12; A21 A22]).module != LinearMaps
             @inferred hvcat((2,2), LinearMap(A11), LinearMap(A12), LinearMap(A21), LinearMap(A22))
             L = [LinearMap(A11) LinearMap(A12); LinearMap(A21) LinearMap(A22)]
             @test @inferred(LinearMaps.MulStyle(L)) === matrixstyle
@@ -74,6 +83,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools
             @test Matrix(L) == A
             @test convert(AbstractMatrix, L) == A
             A = [I A12; A21 I]
+            @test (@which [I A12; A21 I]).module != LinearMaps
             @inferred hvcat((2,2), I, LinearMap(A12), LinearMap(A21), I)
             L = @inferred hvcat((2,2), I, LinearMap(A12), LinearMap(A21), I)
             @test L isa LinearMaps.BlockMap{elty}
@@ -173,6 +183,8 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools
 
             # Md = diag(M1, M2, M3, M2, M1) # unsupported so use sparse:
             Md = Matrix(blockdiag(sparse.((M1, M2, M3, M2, M1))...))
+            @test (@which blockdiag(sparse.((M1, M2, M3, M2, M1))...)).module != LinearMaps
+            @test (@which cat(M1, M2, M3, M2, M1; dims=(1,2))).module != LinearMaps
             x = randn(elty, size(Md, 2))
             Bd = @inferred blockdiag(L1, L2, L3, L2, L1)
             @test Matrix(Bd) == Md