More general printing

Project.toml

@@ -5,10 +5,11 @@ version = "0.12.7"
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [compat]
-ArrayLayouts = "0.3"
+ArrayLayouts = "0.3.2"
 Compat = "2.2, 3"
 julia = "1.1"
 

src/BlockArrays.jl

@@ -21,23 +21,24 @@ export blockappend!, blockpush!, blockpushfirst!, blockpop!, blockpopfirst!
 import Base: @propagate_inbounds, Array, to_indices, to_index,
             unsafe_indices, first, last, size, length, unsafe_length,
             unsafe_convert,
-            getindex, show,
+            getindex, ndims, show,
             step, 
             broadcast, eltype, convert, similar,
             @_inline_meta, _maybetail, tail, @_propagate_inbounds_meta, reindex,
             RangeIndex, Int, Integer, Number,
-            +, -, min, max, *, isless, in, copy, copyto!, axes, @deprecate,
+            +, -, *, /, \, min, max, isless, in, copy, copyto!, axes, @deprecate,
             BroadcastStyle, checkbounds, throw_boundserror
 using Base: ReshapedArray, dataids
 
 
-import Base: (:), IteratorSize, iterate, axes1, strides
-import Base.Broadcast: broadcasted, DefaultArrayStyle, AbstractArrayStyle, Broadcasted
+import Base: (:), IteratorSize, iterate, axes1, strides, isempty
+import Base.Broadcast: broadcasted, DefaultArrayStyle, AbstractArrayStyle, Broadcasted, broadcastable
 import LinearAlgebra: lmul!, rmul!, AbstractTriangular, HermOrSym, AdjOrTrans,
                         StructuredMatrixStyle
 import ArrayLayouts: _fill_lmul!, MatMulVecAdd, MatMulMatAdd, MatLmulVec, MatLdivVec,
                         materialize!, MemoryLayout, sublayout, transposelayout, conjlayout, 
-                        triangularlayout, triangulardata, _inv
+                        triangularlayout, triangulardata, _inv, _copyto!, axes_print_matrix_row,
+                        colsupport, rowsupport
 
 if !@isdefined(only)
     using Compat: only
@@ -49,12 +50,12 @@ include("abstractblockarray.jl")
 include("blockarray.jl")
 include("pseudo_blockarray.jl")
 include("views.jl")
-include("show.jl")
 include("blocks.jl")
 include("blockarrayinterface.jl")
 include("blockbroadcast.jl")
 include("blocklinalg.jl")
 include("blockproduct.jl")
+include("show.jl")
 include("blockreduce.jl")
 include("blockdeque.jl")
 

src/blockarray.jl

@@ -62,12 +62,12 @@ struct BlockArray{T, N, R <: AbstractArray{<:AbstractArray{T,N},N}, BS<:NTuple{N
     blocks::R
     axes::BS
 
-    global _BlockArray(blocks::R, block_sizes::BS) where {T, N, R<:AbstractArray{<:AbstractArray{T,N},N}, BS<:NTuple{N,AbstractUnitRange{Int}}} =
+    global @inline _BlockArray(blocks::R, block_sizes::BS) where {T, N, R<:AbstractArray{<:AbstractArray{T,N},N}, BS<:NTuple{N,AbstractUnitRange{Int}}} =
         new{T, N, R, BS}(blocks, block_sizes)
 end
 
 # Auxilary outer constructors
-_BlockArray(blocks::R, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N, R<:AbstractArray{<:AbstractArray{T,N},N}} =
+@inline _BlockArray(blocks::R, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N, R<:AbstractArray{<:AbstractArray{T,N},N}} =
     _BlockArray(blocks, map(blockedrange, block_sizes))
 
 # support non-concrete eltypes in blocks

src/blockarrayinterface.jl

@@ -3,12 +3,3 @@ getindex(a::Number, ::Block{0}) = a
 
 axes(A::AbstractTriangular{<:Any,<:AbstractBlockMatrix}) = axes(parent(A))
 axes(A::HermOrSym{<:Any,<:AbstractBlockMatrix}) = axes(parent(A))
-
-Base.print_matrix_row(io::IO,
-        X::Union{AbstractTriangular{<:Any,<:AbstractBlockMatrix},
-                 Symmetric{<:Any,<:AbstractBlockMatrix},
-                 Hermitian{<:Any,<:AbstractBlockMatrix},
-                 Adjoint{<:Any,<:AbstractBlockMatrix},
-                 Transpose{<:Any,<:AbstractBlockMatrix}}, A::Vector,
-        i::Integer, cols::AbstractVector, sep::AbstractString) =
-        _blockarray_print_matrix_row(io, X, A, i, cols, sep)

src/blockbroadcast.jl

@@ -147,10 +147,8 @@ end
     return copyto!(dest, Base.Broadcast.instantiate(Base.Broadcast.Broadcasted{BS}(bc.f, bc.args, combine_blockaxes.(axes(dest),axes(bc)))))
 end
 
-@generated function copyto!(
-        dest::AbstractArray,
-        bc::Broadcasted{<:AbstractBlockStyle{NDims}, <:Any, <:Any, Args},
-        ) where {NDims, Args <: Tuple}
+@generated function copyto!(dest::AbstractArray,
+                            bc::Broadcasted{<:AbstractBlockStyle{NDims}, <:Any, <:Any, Args}) where {NDims, Args <: Tuple}
 
     NArgs = length(Args.parameters)
 
@@ -190,29 +188,47 @@ end
     end
 end
 
-@inline function Broadcast.instantiate(
-        bc::Broadcasted{Style}) where {Style <:AbstractBlockStyle}
-    bcf = Broadcast.flatten(Broadcasted{Nothing}(bc.f, bc.args, bc.axes))
+@inline function Broadcast.instantiate(bc::Broadcasted{Style}) where {Style <:AbstractBlockStyle}
+    bcf = Broadcast.instantiate(Broadcast.flatten(Broadcasted{Nothing}(bc.f, bc.args, bc.axes)))
     return Broadcasted{Style}(bcf.f, bcf.args, bcf.axes)
 end
 
 
 for op in (:+, :-, :*)
-    @eval function copy(bc::Broadcasted{BlockStyle{N},<:Any,typeof($op),<:Tuple{<:BlockArray{<:Number,N}}}) where N
+    @eval function copy(bc::Broadcasted{BlockStyle{N},<:Any,typeof($op),<:Tuple{<:AbstractArray{<:Number,N}}}) where N
         (A,) = bc.args
-        _BlockArray(broadcast(a -> broadcast($op, a), A.blocks), axes(A))
+        _BlockArray(broadcast(a -> broadcast($op, a), blocks(A)), axes(A))
     end
 end
 
 for op in (:+, :-, :*, :/, :\)
     @eval begin
-        function copy(bc::Broadcasted{BlockStyle{N},<:Any,typeof($op),<:Tuple{<:Number,<:BlockArray{<:Number,N}}}) where N
+        function copy(bc::Broadcasted{BlockStyle{N},<:Any,typeof($op),<:Tuple{<:Number,<:AbstractArray{<:Number,N}}}) where N
             x,A = bc.args
-            _BlockArray(broadcast((x,a) -> broadcast($op, x, a), x, A.blocks), axes(A))
+            _BlockArray(broadcast((x,a) -> broadcast($op, x, a), x, blocks(A)), axes(A))
         end
-        function copy(bc::Broadcasted{BlockStyle{N},<:Any,typeof($op),<:Tuple{<:BlockArray{<:Number,N},<:Number}}) where N
+        function copy(bc::Broadcasted{BlockStyle{N},<:Any,typeof($op),<:Tuple{<:AbstractArray{<:Number,N},<:Number}}) where N
             A,x = bc.args
-            _BlockArray(broadcast((a,x) -> broadcast($op, a, x), A.blocks,x), axes(A))
+            _BlockArray(broadcast((a,x) -> broadcast($op, a, x), blocks(A),x), axes(A))
         end
     end
 end
+
+# exploit special cases for *, for example, *(::Number, ::Diagonal)
+for op in (:*, :/) 
+    @eval @inline $op(A::BlockArray, x::Number) = _BlockArray($op(blocks(A),x), axes(A))
+end
+for op in (:*, :\)
+    @eval @inline $op(x::Number, A::BlockArray) = _BlockArray($op(x,blocks(A)), axes(A))
+end
+
+###
+# SubViews
+###
+
+_blocktype(::Type{<:BlockArray{<:Any,N,<:AbstractArray{R,N}}}) where {N,R} = R
+
+BroadcastStyle(::Type{<:SubArray{T,N,Arr,<:NTuple{N,BlockSlice1},false}}) where {T,N,Arr<:BlockArray} = 
+    BroadcastStyle(_blocktype(Arr))
+BroadcastStyle(::Type{<:SubArray{T,N,Arr,<:NTuple{N,BlockSlice{BlockRange{1,Tuple{II}}}},false}}) where {T,N,Arr<:BlockArray,II} = 
+    BroadcastStyle(Arr)

src/blockindices.jl

@@ -35,6 +35,17 @@ Block(n::NTuple{N, T}) where {N,T} = Block{N, T}(n)
     Block{N, T}(ntuple(i -> blocks[i].n[1], Val(N)))
 end
 
+# iterate and broadcast like Number
+length(b::Block) = 1
+size(b::Block) = ()
+iterate(x::Block) = (x, nothing)
+iterate(x::Block, ::Any) = nothing
+isempty(x::Block) = false
+broadcastable(x::Block) = x
+ndims(::Type{<:Block}) = 0
+ndims(::Block) = 0
+eltype(::Type{B}) where B<:Block = B
+getindex(B::Block, ::CartesianIndex{0}) = B
 
 # The following code is taken from CartesianIndex
 @inline (+)(index::Block{N}) where {N} = Block{N}(map(+, index.n))
@@ -295,6 +306,7 @@ BlockRange(inds::Vararg{AbstractUnitRange{Int},N}) where {N} =
 Base.BroadcastStyle(::Type{<:BlockRange{1}}) = DefaultArrayStyle{1}()
 broadcasted(::DefaultArrayStyle{1}, ::Type{Block}, r::AbstractUnitRange) = Block(first(r)):Block(last(r))
 broadcasted(::DefaultArrayStyle{1}, ::Type{Int}, block_range::BlockRange{1}) = first(block_range.indices)
+broadcasted(::DefaultArrayStyle{0}, ::Type{Int}, block::Block{1}) = Int(block)
 
 
 # AbstractArray implementation

src/blocklinalg.jl

@@ -7,7 +7,7 @@ gives an iterator containing the possible non-zero blocks in the k-th block-row
 blockrowsupport(A, k) = blockrowsupport(MemoryLayout(A), A, k)
 blockrowsupport(A) = blockrowsupport(A, blockaxes(A,1))
 
-blockcolsupport(_, A, j) = blockaxes(A,1)
+blockcolsupport(_, A, j) = Block.(colsupport(blocks(A), Int.(j)))
 
 """"
     blockcolsupport(A, j)
@@ -28,24 +28,41 @@ end
 
 
 abstract type AbstractBlockLayout <: MemoryLayout end
-struct BlockLayout{LAY} <: AbstractBlockLayout end
+struct BlockLayout{ArrLay,BlockLay} <: AbstractBlockLayout end
 
+function colsupport(::BlockLayout, A, j)
+    JR = Int(findblock(axes(A,2), minimum(j))):Int(findblock(axes(A,2), maximum(j)))
+    KR = colsupport(blocks(A), JR)
+    axes(A,1)[Block.(KR)]
+end
+
+function rowsupport(::BlockLayout, A, k)
+    KR = Int(findblock(axes(A,1), minimum(k))):Int(findblock(axes(A,1), maximum(k)))
+    JR = rowsupport(blocks(A), KR)
+    axes(A,2)[Block.(JR)]
+end
 
-similar(M::MulAdd{<:AbstractBlockLayout,<:AbstractBlockLayout}, ::Type{T}, axes) where {T,N} = 
+similar(M::MulAdd{<:AbstractBlockLayout,<:AbstractBlockLayout}, ::Type{T}, axes) where {T,N} =
     similar(BlockArray{T}, axes)
 
 MemoryLayout(::Type{<:PseudoBlockArray{T,N,R}}) where {T,N,R} = MemoryLayout(R)
-MemoryLayout(::Type{<:BlockArray{T,N,R}}) where {T,N,R} = BlockLayout{typeof(MemoryLayout(R))}()
+MemoryLayout(::Type{<:BlockArray{T,N,R}}) where {T,N,D,R<:AbstractArray{D,N}} =
+    BlockLayout{typeof(MemoryLayout(R)),typeof(MemoryLayout(D))}()
 
-sublayout(::BlockLayout{LAY}, ::Type{NTuple{N,BlockSlice1}}) where {LAY,N} = LAY()
-sublayout(BL::BlockLayout, ::Type{<:NTuple{N,BlockSlice}}) where N = BL
+sublayout(::BlockLayout{MLAY,BLAY}, ::Type{NTuple{N,BlockSlice1}}) where {MLAY,BLAY,N} = BLAY()
+sublayout(BL::BlockLayout{MLAY,BLAY}, ::Type{<:NTuple{N,<:BlockSlice{BlockRange{1,Tuple{II}}}}}) where {N,MLAY,BLAY,II} =
+    BlockLayout{typeof(sublayout(MLAY(),NTuple{N,II})), BLAY}()
+sublayout(BL::BlockLayout{MLAY,BLAY}, ::Type{<:Tuple{BlockSlice1,<:BlockSlice{BlockRange{1,Tuple{II}}}}}) where {MLAY,BLAY,II} =
+    BlockLayout{typeof(sublayout(MLAY(),Tuple{Int,II})), BLAY}()
+sublayout(BL::BlockLayout{MLAY,BLAY}, ::Type{<:Tuple{<:BlockSlice{BlockRange{1,Tuple{II}}},BlockSlice1}}) where {MLAY,BLAY,II} =
+    BlockLayout{typeof(sublayout(MLAY(),Tuple{II,Int})), BLAY}()
 
-conjlayout(::Type{T}, ::BlockLayout{LAY}) where {T<:Complex,LAY} = BlockLayout{typeof(conjlayout(T,LAY()))}()
-conjlayout(::Type{T}, ::BlockLayout{LAY}) where {T<:Real,LAY} = BlockLayout{LAY}()
-
-transposelayout(::BlockLayout{LAY}) where LAY = BlockLayout{typeof(transposelayout(LAY()))}()
 
+conjlayout(::Type{T}, ::BlockLayout{MLAY,BLAY}) where {T<:Complex,MLAY,BLAY} = BlockLayout{MLAY,typeof(conjlayout(T,BLAY()))}()
+conjlayout(::Type{T}, ::BlockLayout{MLAY,BLAY}) where {T<:Real,MLAY,BLAY} = BlockLayout{MLAY,BLAY}()
 
+transposelayout(::BlockLayout{MLAY,BLAY}) where {MLAY,BLAY} =
+    BlockLayout{typeof(transposelayout(MLAY())),typeof(transposelayout(BLAY()))}()
 
 #############
 # BLAS overrides
@@ -78,10 +95,10 @@ function _block_muladd!(α, A, X, β, Y)
         muladd!(α, view(A,K,N), view(X,N,J), one(α), view(Y,K,J))
     end
     Y
-end 
+end
 
-mul_blockscompatible(A, B, C) = blockisequal(axes(A,2), axes(B,1)) && 
-    blockisequal(axes(A,1), axes(C,1)) && 
+mul_blockscompatible(A, B, C) = blockisequal(axes(A,2), axes(B,1)) &&
+    blockisequal(axes(A,1), axes(C,1)) &&
     blockisequal(axes(B,2), axes(C,2))
 
 function materialize!(M::MatMulMatAdd{<:AbstractBlockLayout,<:AbstractBlockLayout,<:AbstractBlockLayout})
@@ -193,7 +210,7 @@ for UNIT in ('U', 'N')
 
             A = triangulardata(U)
             if !hasmatchingblocks(A) # Use default for now
-                return materialize!(Ldiv{TriangularLayout{'U',$UNIT,UnknownLayout}, 
+                return materialize!(Ldiv{TriangularLayout{'U',$UNIT,UnknownLayout},
                                          typeof(MemoryLayout(dest))}(U, dest))
             end
 
@@ -226,7 +243,7 @@ for UNIT in ('U', 'N')
             T = eltype(dest)
             A = triangulardata(L)
             if !hasmatchingblocks(A) # Use default for now
-                return materialize!(Ldiv{TriangularLayout{'L',$UNIT,UnknownLayout}, 
+                return materialize!(Ldiv{TriangularLayout{'L',$UNIT,UnknownLayout},
                                          typeof(MemoryLayout(dest))}(L, dest))
             end
 

src/blocks.jl

@@ -52,8 +52,10 @@ blocks(A::Adjoint) = adjoint(blocks(parent(A)))
 blocks(A::Transpose) = transpose(blocks(parent(A)))
 
 # convert a tuple of BlockRange to a tuple of `AbstractUnitRange{Int}`
+_block2int(B::Block{1}) = Int(B):Int(B)
+_block2int(B::BlockRange{1}) = Int.(B)
 _blockrange2int() = ()
-_blockrange2int(A, B...) = tuple(Int.(A.block), _blockrange2int(B...)...)
+_blockrange2int(A, B...) = tuple(_block2int(A.block), _blockrange2int(B...)...)
 
 blocks(A::SubArray{<:Any,N,<:Any,<:NTuple{N,BlockSlice}}) where N =
     view(blocks(parent(A)), _blockrange2int(parentindices(A)...)...)

src/pseudo_blockarray.jl

@@ -321,3 +321,11 @@ end
 Base.strides(A::PseudoBlockArray) = strides(A.blocks)
 Base.stride(A::PseudoBlockArray, i::Integer) = stride(A.blocks, i)
 Base.unsafe_convert(::Type{Ptr{T}}, A::PseudoBlockArray) where T = Base.unsafe_convert(Ptr{T}, A.blocks)
+
+
+###
+# col/rowsupport
+###
+
+colsupport(A::PseudoBlockArray, j) = colsupport(A.blocks, j)
+rowsupport(A::PseudoBlockArray, j) = rowsupport(A.blocks, j)

src/show.jl

@@ -88,12 +88,6 @@ function _blockarray_print_matrix_row(io::IO,
     end
 end
 
-
-Base.print_matrix_row(io::IO,
-        X::AbstractBlockVecOrMat, A::Vector,
-        i::Integer, cols::AbstractVector, sep::AbstractString) =
-        _blockarray_print_matrix_row(io, X, A, i, cols, sep)
-
 function _show_typeof(io::IO, a::BlockArray{T,N,Array{Array{T,N},N},NTuple{N,DefaultBlockAxis}}) where {T,N}
     Base.show_type_name(io, typeof(a).name)
     print(io, '{')
@@ -103,6 +97,20 @@ function _show_typeof(io::IO, a::BlockArray{T,N,Array{Array{T,N},N},NTuple{N,Def
     print(io, '}')
 end
 
+# LayoutArray with blocked axes will dispatch to here
+axes_print_matrix_row(::Tuple{BlockedUnitRange}, io, X, A, i, cols, sep) =
+        _blockarray_print_matrix_row(io, X, A, i, cols, sep)
+axes_print_matrix_row(::NTuple{2,BlockedUnitRange}, io, X, A, i, cols, sep) =
+        _blockarray_print_matrix_row(io, X, A, i, cols, sep)
+axes_print_matrix_row(::Tuple{AbstractUnitRange,BlockedUnitRange}, io, X, A, i, cols, sep) =
+        _blockarray_print_matrix_row(io, X, A, i, cols, sep)
+axes_print_matrix_row(::Tuple{BlockedUnitRange,AbstractUnitRange}, io, X, A, i, cols, sep) =
+        _blockarray_print_matrix_row(io, X, A, i, cols, sep)
+
+# Need to handled BlockedUnitRange, which is not a LayoutVector
+Base.print_matrix_row(io::IO, X::BlockedUnitRange, A::Vector, i::Integer, cols::AbstractVector, sep::AbstractString) =
+        _blockarray_print_matrix_row(io, X, A, i, cols, sep)
+
 function _show_typeof(io::IO, a::PseudoBlockArray{T,N,Array{T,N},NTuple{N,DefaultBlockAxis}}) where {T,N}
     Base.show_type_name(io, typeof(a).name)
     print(io, '{')
@@ -112,13 +120,7 @@ function _show_typeof(io::IO, a::PseudoBlockArray{T,N,Array{T,N},NTuple{N,Defaul
     print(io, '}')
 end
 
-
 ## Cumsum
 
-Base.print_matrix_row(io::IO,
-        X::BlockedUnitRange, A::Vector,
-        i::Integer, cols::AbstractVector, sep::AbstractString) =
-        _blockarray_print_matrix_row(io, X, A, i, cols, sep)
-
 Base.show(io::IO, mimetype::MIME"text/plain", a::BlockedUnitRange) = 
     Base.invoke(show, Tuple{typeof(io),MIME"text/plain",AbstractArray},io, mimetype, a)

test/test_blockarrayinterface.jl

@@ -1,4 +1,4 @@
-using BlockArrays, LinearAlgebra, Base64
+using BlockArrays, LinearAlgebra, FillArrays, Base64
 
 struct PartiallyImplementedBlockVector <: AbstractBlockArray{Float64,1} end
 
@@ -79,3 +79,26 @@ end
     @test A isa BlockMatrix{Int,Matrix{Matrix{Int}},Tuple{BlockedUnitRange{StepRange{Int64,Int64}},BlockedUnitRange{Vector{Int}}}}
 end
 
+@testset "block Fill" begin
+    A = Fill(2,(blockedrange([1,2,2]),))
+    @test A[Block(1)] == [2]
+    @test A[Block.(1:2)] == [2,2,2]
+    @test_broken A[Block(1)] isa Fill
+    @test_broken A[Block.(1:2)] isa Fill
+    @test_broken 2A
+
+
+    B = Eye((blockedrange([1,2]),))
+    @test B[Block(2,2)] == Matrix(I,2,2)
+
+    C = Eye((blockedrange([1,2,2]),blockedrange([2,2])))
+    @test C[Block(2,2)] == [0 0; 1.0 0]
+
+    U = UpperTriangular(Ones((blockedrange([1,2]),blockedrange([2,1]))))
+
+    if VERSION ≥ v"1.2"
+        @test stringmime("text/plain", A) == "5-element Fill{Int64,1,Tuple{BlockedUnitRange{Array{Int64,1}}}} with indices 1:1:5:\n 2\n ─\n 2\n 2\n ─\n 2\n 2"
+        @test stringmime("text/plain", B) == "3×3 Diagonal{Float64,Ones{Float64,1,Tuple{BlockedUnitRange{Array{Int64,1}}}}} with indices 1:1:3×1:1:3:\n 1.0  │   ⋅    ⋅ \n ─────┼──────────\n  ⋅   │  1.0   ⋅ \n  ⋅   │   ⋅   1.0"
+        @test stringmime("text/plain", U) == "3×3 UpperTriangular{Float64,Ones{Float64,2,Tuple{BlockedUnitRange{Array{Int64,1}},BlockedUnitRange{Array{Int64,1}}}}} with indices 1:1:3×1:1:3:\n 1.0  1.0  │  1.0\n ──────────┼─────\n  ⋅   1.0  │  1.0\n  ⋅    ⋅   │  1.0"
+    end
+end

test/test_blockbroadcast.jl

@@ -119,4 +119,16 @@ using BlockArrays, FillArrays, Test
         C = BlockArray(randn(6), (BlockArrays._BlockedUnitRange(1,2:6),))
         @test axes(A+C,1) === BlockArrays._BlockedUnitRange(1,1:6)
     end
+
+    @testset "Views" begin
+        A = BlockArray(randn(6), 1:3)
+        @test Base.BroadcastStyle(typeof(view(A, Block(2)))) isa Base.Broadcast.DefaultArrayStyle{1}
+        V = view(A, Block.(2:3))
+        @test Base.BroadcastStyle(typeof(V)) isa BlockArrays.BlockStyle{1}
+        @test V .+ 1 isa BlockArray
+        @test 1 .+ V isa BlockArray
+        @test V .+ 1 == 1 .+ V == A[Block.(2:3)] .+ 1
+        @test -V isa BlockArray
+        @test -V == -A[Block.(2:3)]
+    end
 end