add the code

Author: dkarrasch

src/LinearMaps.jl

@@ -16,6 +16,10 @@ const MapOrVecOrMat{T} = Union{LinearMap{T}, AbstractVecOrMat{T}}
 const MapOrMatrix{T} = Union{LinearMap{T}, AbstractMatrix{T}}
 const RealOrComplex = Union{Real, Complex}
 
+const LinearMapTuple = Tuple{Vararg{LinearMap}}
+const LinearMapVector = AbstractVector{<:LinearMap}
+const LinearMapTupleOrVector = Union{LinearMapTuple,LinearMapVector}
+
 Base.eltype(::LinearMap{T}) where {T} = T
 
 abstract type MulStyle end
@@ -69,6 +73,21 @@ convert_to_lmaps(A) = (convert_to_lmaps_(A),)
 @inline convert_to_lmaps(A, B, Cs...) =
     (convert_to_lmaps_(A), convert_to_lmaps_(B), convert_to_lmaps(Cs...)...)
 
+_front(As::Tuple) = Base.front(As)
+_front(As::AbstractVector) = @inbounds @views As[1:end-1]
+_tail(As::Tuple) = Base.tail(As)
+_tail(As::AbstractVector) = @inbounds @views As[2:end]
+
+_combine(A::LinearMap, B::LinearMap) = tuple(A, B)
+_combine(A::LinearMap, Bs::LinearMapTuple) = tuple(A, Bs...)
+_combine(As::LinearMapTuple, B::LinearMap) = tuple(As..., B)
+_combine(As::LinearMapTuple, Bs::LinearMapTuple) = tuple(As..., Bs...)
+_combine(A::LinearMap, Bs::LinearMapVector) = vcat(A, Bs...)
+_combine(As::LinearMapVector, B::LinearMap) = vcat(As..., B)
+_combine(As::LinearMapVector, Bs::LinearMapTuple) = vcat(As..., Bs...)
+_combine(As::LinearMapTuple, Bs::LinearMapVector) = vcat(As..., Bs...)
+_combine(As::LinearMapVector, Bs::LinearMapVector) = vcat(As..., Bs...)
+
 # The (internal) multiplication logic is as follows:
 #  - `*(A, x)` calls `mul!(y, A, x)` for appropriately-sized y
 #  - `mul!` checks consistency of the sizes, and calls `_unsafe_mul!`,
@@ -233,8 +252,6 @@ function _unsafe_mul!(y::AbstractMatrix, A::LinearMap, x::AbstractMatrix, α, β
     return _generic_mapmat_mul!(y, A, x, α, β)
 end
 
-const LinearMapTuple = Tuple{Vararg{LinearMap}}
-
 include("left.jl") # left multiplication by a transpose or adjoint vector
 include("transpose.jl") # transposing linear maps
 include("wrappedmap.jl") # wrap a matrix of linear map in a new type, thereby allowing to alter its properties

src/blockmap.jl

@@ -1,14 +1,14 @@
 struct BlockMap{T,
-                As<:LinearMapTuple,
+                As<:LinearMapTupleOrVector,
                 Rs<:Tuple{Vararg{Int}},
-                Rranges<:Tuple{Vararg{UnitRange{Int}}},
-                Cranges<:Tuple{Vararg{UnitRange{Int}}}} <: LinearMap{T}
+                Rranges<:Vector{UnitRange{Int}},
+                Cranges<:Vector{UnitRange{Int}}} <: LinearMap{T}
     maps::As
     rows::Rs
     rowranges::Rranges
     colranges::Cranges
     function BlockMap{T,As,Rs}(maps::As, rows::Rs) where
-                {T, As<:LinearMapTuple, Rs<:Tuple{Vararg{Int}}}
+                {T, As<:LinearMapTupleOrVector, Rs<:Tuple{Vararg{Int}}}
         for TA in Base.Generator(eltype, maps)
             promote_type(T, TA) == T ||
                 error("eltype $TA cannot be promoted to $T in BlockMap constructor")
@@ -19,15 +19,17 @@ struct BlockMap{T,
     end
 end
 
-BlockMap{T}(maps::As, rows::Rs) where {T, As<:LinearMapTuple, Rs} =
+BlockMap{T}(maps::As, rows::Rs) where {T, As<:LinearMapTupleOrVector, Rs} =
     BlockMap{T, As, Rs}(maps, rows)
+BlockMap(maps::As, rows::Rs) where {As<:LinearMapTupleOrVector, Rs} =
+    BlockMap{promote_type(map(eltype, maps)...), As, Rs}(maps, rows)
 
 MulStyle(A::BlockMap) = MulStyle(A.maps...)
 
-function _getranges(maps, dim, inds=ntuple(identity, Val(length(maps))))
-    sizes = map(i -> size(maps[i], dim)::Int, inds)
-    ends = cumsum(sizes)
-    starts = (1, (1 .+ Base.front(ends))...)
+function _getranges(maps, dim, inds=1:length(maps))
+    ends = map(i -> size(maps[i], dim)::Int, inds)
+    cumsum!(ends, ends)
+    starts = vcat(1, 1 .+ @views ends[1:end-1])
     return UnitRange.(starts, ends)
 end
 
@@ -40,14 +42,15 @@ block linear map obtained from `hvcat(rows, maps...)`.
 """
 function rowcolranges(maps, rows)
     # find indices of the row-wise first maps
-    firstmapinds = cumsum((1, Base.front(rows)...))
+    firstmapinds = vcat(1, Base.front(rows)...)
+    cumsum!(firstmapinds, firstmapinds)
     # compute rowranges from first dimension of the row-wise first maps
     rowranges = _getranges(maps, 1, firstmapinds)
 
     # compute ranges from second dimension as if all in one row
     temp = _getranges(maps, 2)
     # introduce "line breaks"
-    colranges = ntuple(Val(length(maps))) do i
+    colranges = map(1:length(maps)) do i
         # for each map find the index of the respective row-wise first map
         # something-trick just to assure the compiler that the index is an Int
         @inbounds firstmapind = firstmapinds[something(findlast(<=(i), firstmapinds), 1)]
@@ -277,7 +280,7 @@ end
 ############
 
 Base.:(==)(A::BlockMap, B::BlockMap) =
-    (eltype(A) == eltype(B) && A.maps == B.maps && A.rows == B.rows)
+    (eltype(A) == eltype(B) && all(A.maps .== B.maps) && all(A.rows .== B.rows))
 
 ############
 # multiplication helper functions
@@ -350,9 +353,9 @@ function _transblockmul!(y, A::BlockMap, x, α, β, transform)
         # subsequent block rows of A (block columns of A'),
         # add results to corresponding parts of y
         # TODO: think about multithreading
-        for (row, xi) in zip(Base.tail(rows), Base.tail(xinds))
-            xrow = selectdim(x, 1, xi)
-            for _ in 1:row
+        @inbounds for i in 2:length(rows)
+            xrow = selectdim(x, 1, xinds[i])
+            for _ in 1:rows[i]
                 mapind +=1
                 yrow = selectdim(y, 1, yinds[mapind])
                 _unsafe_mul!(yrow, transform(maps[mapind]), xrow, α, true)
@@ -397,12 +400,12 @@ end
 # BlockDiagonalMap
 ############
 struct BlockDiagonalMap{T,
-                        As<:LinearMapTuple,
-                        Ranges<:Tuple{Vararg{UnitRange{Int}}}} <: LinearMap{T}
+                        As<:LinearMapTupleOrVector,
+                        Ranges<:Vector{UnitRange{Int}}} <: LinearMap{T}
     maps::As
     rowranges::Ranges
     colranges::Ranges
-    function BlockDiagonalMap{T, As}(maps::As) where {T, As<:LinearMapTuple}
+    function BlockDiagonalMap{T, As}(maps::As) where {T, As<:LinearMapTupleOrVector}
         for TA in Base.Generator(eltype, maps)
             promote_type(T, TA) == T ||
                 error("eltype $TA cannot be promoted to $T in BlockDiagonalMap constructor")
@@ -413,7 +416,7 @@ struct BlockDiagonalMap{T,
     end
 end
 
-BlockDiagonalMap{T}(maps::As) where {T, As<:LinearMapTuple} =
+BlockDiagonalMap{T}(maps::As) where {T, As<:LinearMapTupleOrVector} =
     BlockDiagonalMap{T,As}(maps)
 BlockDiagonalMap(maps::LinearMap...) =
     BlockDiagonalMap{promote_type(map(eltype, maps)...)}(maps)
@@ -478,7 +481,7 @@ LinearAlgebra.transpose(A::BlockDiagonalMap{T}) where {T} =
     BlockDiagonalMap{T}(map(transpose, A.maps))
 
 Base.:(==)(A::BlockDiagonalMap, B::BlockDiagonalMap) =
-    (eltype(A) == eltype(B) && A.maps == B.maps)
+    (eltype(A) == eltype(B) && all(A.maps .== B.maps))
 
 for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
     @eval begin
@@ -496,7 +499,7 @@ end
 function _blockscaling!(y, A::BlockDiagonalMap, x, α, β)
     maps, yinds, xinds = A.maps, A.rowranges, A.colranges
     # TODO: think about multi-threading here
-    @inbounds for i in eachindex(yinds, maps, xinds)
+    @inbounds for i in 1:length(maps)
         _unsafe_mul!(selectdim(y, 1, yinds[i]), maps[i], selectdim(x, 1, xinds[i]), α, β)
     end
     return y

src/composition.jl

@@ -1,4 +1,4 @@
-struct CompositeMap{T, As<:LinearMapTuple} <: LinearMap{T}
+struct CompositeMap{T, As<:LinearMapTupleOrVector} <: LinearMap{T}
     maps::As # stored in order of application to vector
     function CompositeMap{T, As}(maps::As) where {T, As}
         N = length(maps)
@@ -12,7 +12,7 @@ struct CompositeMap{T, As<:LinearMapTuple} <: LinearMap{T}
         new{T, As}(maps)
     end
 end
-CompositeMap{T}(maps::As) where {T, As<:LinearMapTuple} = CompositeMap{T, As}(maps)
+CompositeMap{T}(maps::As) where {T, As<:LinearMapTupleOrVector} = CompositeMap{T, As}(maps)
 
 # basic methods
 Base.size(A::CompositeMap) = (size(A.maps[end], 1), size(A.maps[1], 2))
@@ -26,8 +26,18 @@ for (f, _f, g) in ((:issymmetric, :_issymmetric, :transpose),
         LinearAlgebra.$f(A::CompositeMap) = $_f(A.maps)
         $_f(maps::Tuple{}) = true
         $_f(maps::Tuple{<:LinearMap}) = $f(maps[1])
-        $_f(maps::Tuple{Vararg{LinearMap}}) =
+        $_f(maps::LinearMapTuple) =
             maps[end] == $g(maps[1]) && $_f(Base.front(Base.tail(maps)))
+        function $_f(maps::LinearMapVector)
+            n = length(maps)
+            if n == 0
+                return true
+            elseif n == 1
+                return $f(maps[1])
+            else
+                return maps[end] == $g(maps[1]) && $_f(maps[2:end-1])
+            end
+        end
         # since the introduction of ScaledMap, the following cases cannot occur
         # function $_f(maps::Tuple{Vararg{LinearMap}}) # length(maps) >= 2
             # if maps[1] isa UniformScalingMap{<:RealOrComplex}
@@ -45,23 +55,34 @@ end
 LinearAlgebra.isposdef(A::CompositeMap) = _isposdef(A.maps)
 _isposdef(maps::Tuple{}) = true # empty product is equivalent to "I" which is pos. def.
 _isposdef(maps::Tuple{<:LinearMap}) = isposdef(maps[1])
-function _isposdef(maps::Tuple{Vararg{LinearMap}})
+function _isposdef(maps::LinearMapTuple)
     (maps[end] == adjoint(maps[1]) || maps[end] == maps[1]) && 
-    isposdef(maps[1]) && _isposdef(Base.front(Base.tail(maps)))
+        isposdef(maps[1]) && _isposdef(Base.front(Base.tail(maps)))
+end
+function _isposdef(maps::LinearMapVector)
+    n = length(maps)
+    if n == 0
+        return true
+    elseif n == 1
+        return isposdef(maps[1])
+    else
+        return (maps[end] == adjoint(maps[1]) || maps[end] == maps[1]) && 
+            isposdef(maps[1]) && _isposdef(maps[2:end-1])
+    end
 end
 
 # scalar multiplication and division (non-commutative case)
 function Base.:(*)(α::Number, A::LinearMap)
     T = promote_type(typeof(α), eltype(A))
-    return CompositeMap{T}(tuple(A, UniformScalingMap(α, size(A, 1))))
+    return CompositeMap{T}(_combine(A, UniformScalingMap(α, size(A, 1))))
 end
 function Base.:(*)(α::Number, A::CompositeMap)
     T = promote_type(typeof(α), eltype(A))
     Alast = last(A.maps)
     if Alast isa UniformScalingMap
-        return CompositeMap{T}(tuple(Base.front(A.maps)..., α * Alast))
+        return CompositeMap{T}(_combine(_front(A.maps), α * Alast))
     else
-        return CompositeMap{T}(tuple(A.maps..., UniformScalingMap(α, size(A, 1))))
+        return CompositeMap{T}(_combine(A.maps, UniformScalingMap(α, size(A, 1))))
     end
 end
 # needed for disambiguation
@@ -71,15 +92,15 @@ function Base.:(*)(α::RealOrComplex, A::CompositeMap{<:RealOrComplex})
 end
 function Base.:(*)(A::LinearMap, α::Number)
     T = promote_type(typeof(α), eltype(A))
-    return CompositeMap{T}(tuple(UniformScalingMap(α, size(A, 2)), A))
+    return CompositeMap{T}(_combine(UniformScalingMap(α, size(A, 2)), A))
 end
 function Base.:(*)(A::CompositeMap, α::Number)
     T = promote_type(typeof(α), eltype(A))
     Afirst = first(A.maps)
     if Afirst isa UniformScalingMap
-        return CompositeMap{T}(tuple(Afirst * α, Base.tail(A.maps)...))
+        return CompositeMap{T}(_combine(Afirst * α, _tail(A.maps)))
     else
-        return CompositeMap{T}(tuple(UniformScalingMap(α, size(A, 2)), A.maps...))
+        return CompositeMap{T}(_combine(UniformScalingMap(α, size(A, 2)), A.maps))
     end
 end
 # needed for disambiguation
@@ -110,19 +131,19 @@ julia> LinearMap(ones(Int, 3, 3)) * CS * I * rand(3, 3);
 """
 function Base.:(*)(A₁::LinearMap, A₂::LinearMap)
     T = promote_type(eltype(A₁), eltype(A₂))
-    return CompositeMap{T}(tuple(A₂, A₁))
+    return CompositeMap{T}(_combine(A₂, A₁))
 end
 function Base.:(*)(A₁::LinearMap, A₂::CompositeMap)
     T = promote_type(eltype(A₁), eltype(A₂))
-    return CompositeMap{T}(tuple(A₂.maps..., A₁))
+    return CompositeMap{T}(_combine(A₂.maps, A₁))
 end
 function Base.:(*)(A₁::CompositeMap, A₂::LinearMap)
     T = promote_type(eltype(A₁), eltype(A₂))
-    return CompositeMap{T}(tuple(A₂, A₁.maps...))
+    return CompositeMap{T}(_combine(A₂, A₁.maps))
 end
 function Base.:(*)(A₁::CompositeMap, A₂::CompositeMap)
     T = promote_type(eltype(A₁), eltype(A₂))
-    return CompositeMap{T}(tuple(A₂.maps..., A₁.maps...))
+    return CompositeMap{T}(_combine(A₂.maps, A₁.maps))
 end
 # needed for disambiguation
 Base.:(*)(A₁::ScaledMap, A₂::CompositeMap) = A₁.λ * (A₁.lmap * A₂)
@@ -135,7 +156,8 @@ LinearAlgebra.adjoint(A::CompositeMap{T}) where {T} =
     CompositeMap{T}(map(adjoint, reverse(A.maps)))
 
 # comparison of CompositeMap objects
-Base.:(==)(A::CompositeMap, B::CompositeMap) = (eltype(A) == eltype(B) && A.maps == B.maps)
+Base.:(==)(A::CompositeMap, B::CompositeMap) =
+    (eltype(A) == eltype(B) && all(A.maps .== B.maps))
 
 # multiplication with vectors/matrices
 _unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap, x::AbstractVector) =

src/conversion.jl

@@ -136,7 +136,7 @@ function SparseArrays.sparse(A::BlockMap)
     return hvcat(
         A.rows,
         convert(SparseMatrixCSC, first(A.maps)),
-        convert.(AbstractMatrix, Base.tail(A.maps))...
+        convert.(AbstractArray, _tail(A.maps))...
     )
 end
 Base.Matrix{T}(A::BlockDiagonalMap) where {T} = Base._cat((1,2), convert.(Matrix{T}, A.maps)...)
@@ -152,7 +152,7 @@ Base.convert(::Type{AbstractMatrix}, A::KroneckerMap) =
 function SparseArrays.sparse(A::KroneckerMap)
     return kron(
         convert(SparseMatrixCSC, first(A.maps)),
-        convert.(AbstractMatrix, Base.tail(A.maps))...
+        convert.(AbstractMatrix, _tail(A.maps))...
     )
 end