Merge branch 'master' into generalize_blockindex

Author: mtfishman

docs/src/lib/public.md

@@ -37,6 +37,7 @@ blockfirsts
 blocklasts
 blocklengths
 blocksizes
+eachblockaxes
 blocks
 eachblock
 blockcheckbounds

src/BlockArrays.jl

@@ -5,7 +5,7 @@ using LinearAlgebra, ArrayLayouts, FillArrays
 export AbstractBlockArray, AbstractBlockMatrix, AbstractBlockVector, AbstractBlockVecOrMat
 export Block, getblock, getblock!, setblock!, eachblock, blocks
 export blockaxes, blocksize, blocklength, blockcheckbounds, BlockBoundsError, BlockIndex, BlockIndexRange, BlockIndices
-export blocksizes, blocklengths, blocklasts, blockfirsts, blockisequal
+export blocksizes, blocklengths, eachblockaxes, blocklasts, blockfirsts, blockisequal
 export BlockRange, blockedrange, BlockedUnitRange, BlockedOneTo
 
 export BlockArray, BlockMatrix, BlockVector, BlockVecOrMat, mortar
@@ -20,16 +20,16 @@ export blockappend!, blockpush!, blockpushfirst!, blockpop!, blockpopfirst!
 import Base: @propagate_inbounds, Array, AbstractArray, to_indices, to_index,
             unsafe_indices, first, last, size, length, unsafe_length,
             unsafe_convert,
-            getindex, setindex!, ndims, show, view,
+            getindex, setindex!, ndims, show, print_array, view,
             step,
-            broadcast, eltype, convert, similar,
+            broadcast, eltype, convert, similar, collect,
             tail, reindex,
             RangeIndex, Int, Integer, Number, Tuple,
             +, -, *, /, \, min, max, isless, in, copy, copyto!, axes, @deprecate,
-            BroadcastStyle, checkbounds,
+            BroadcastStyle, checkbounds, checkindex, ensure_indexable,
             oneunit, ones, zeros, intersect, Slice, resize!
 
-using Base: ReshapedArray, dataids, oneto
+using Base: ReshapedArray, LogicalIndex, dataids, oneto
 
 import Base: (:), IteratorSize, iterate, axes1, strides, isempty
 import Base.Broadcast: broadcasted, DefaultArrayStyle, AbstractArrayStyle, Broadcasted, broadcastable

src/blockaxis.jl

@@ -105,6 +105,16 @@ function Base.AbstractUnitRange{T}(r::BlockedUnitRange) where {T}
     return _BlockedUnitRange(convert(T,first(r)), convert.(T,blocklasts(r)))
 end
 
+# See: https://github.com/JuliaLang/julia/blob/b06d26075bf7b3f4e7f1b64b120f5665d8ed76f9/base/range.jl#L991-L1004
+function Base.getindex(r::AbstractUnitRange, s::AbstractBlockedUnitRange{T}) where {T<:Integer}
+    @boundscheck checkbounds(r, s)
+
+    f = first(r)
+    start = oftype(f, f + first(s) - firstindex(r))
+    lens = map(Base.Fix1(oftype, f), blocklengths(s))
+    return blockedrange(start, lens)
+end
+
 """
     BlockedOneTo{T, <:Union{AbstractVector{T}, NTuple{<:Any,T}}} where {T}
 
@@ -164,6 +174,16 @@ function Base.AbstractUnitRange{T}(r::BlockedOneTo) where {T}
     return BlockedOneTo(convert.(T,blocklasts(r)))
 end
 
+# See: https://github.com/JuliaLang/julia/blob/b06d26075bf7b3f4e7f1b64b120f5665d8ed76f9/base/range.jl#L1006-L1010
+function getindex(r::Base.OneTo{T}, s::BlockedOneTo) where T
+    @inline
+    @boundscheck checkbounds(r, s)
+    return BlockedOneTo(convert(AbstractVector{T}, blocklasts(s)))
+end
+function getindex(r::BlockedOneTo{T}, s::BlockedOneTo) where T
+    return Base.oneto(r)[s]
+end
+
 """
     blockedrange(blocklengths::Union{Tuple, AbstractVector})
     blockedrange(first::Integer, blocklengths::Union{Tuple, AbstractVector})
@@ -406,9 +426,13 @@ end
 @propagate_inbounds getindex(b::AbstractBlockedUnitRange, KR::BlockSlice) = b[KR.block]
 
 getindex(b::AbstractBlockedUnitRange, KR::AbstractVector{<:Block{1}}) = mortar([b[K] for K in KR])
-getindex(b::AbstractBlockedUnitRange, KR::AbstractVector{<:BlockIndices{1}}) = mortar([b[K] for K in KR])
-getindex(b::AbstractBlockedUnitRange, KR::AbstractVector{<:BlockIndex{1}}) = [b[K] for K in KR]
+getindex(b::AbstractBlockedUnitRange, KR::AbstractVector{<:AbstractVector{<:Block{1}}}) = mortar([b[K] for K in KR])
 getindex(b::AbstractBlockedUnitRange, Kkr::BlockIndices{1}) =  b[block(Kkr)][Kkr.indices...]
+getindex(b::AbstractBlockedUnitRange, KR::AbstractVector{<:BlockIndex{1}}) = [b[K] for K in KR]
+getindex(b::AbstractBlockedUnitRange, KR::AbstractVector{<:BlockIndices{1}}) = mortar([b[K] for K in KR])
+getindex(b::AbstractBlockedUnitRange, KR::AbstractVector{<:AbstractVector{<:BlockIndex{1}}}) = mortar([b[K] for K in KR])
+getindex(b::AbstractBlockedUnitRange, KR::AbstractVector{<:AbstractVector{<:BlockIndices{1}}}) = mortar([b[K] for K in KR])
+getindex(b::AbstractBlockedUnitRange, KR::AbstractVector{<:AbstractVector{<:AbstractVector{<:BlockIndex{1}}}}) = mortar([b[K] for K in KR])
 
 _searchsortedfirst(a::AbstractVector, k) = searchsortedfirst(a, k)
 function _searchsortedfirst(a::Tuple, k)

src/blockedarray.jl

@@ -193,6 +193,10 @@ AbstractArray{T,N}(A::BlockedArray) where {T,N} = BlockedArray(AbstractArray{T,N
 
 copy(A::BlockedArray) = BlockedArray(copy(A.blocks), A.axes)
 
+# Blocked version of `collect(::AbstractArray)` that preserves the
+# block structure.
+blockcollect(a::AbstractArray) = BlockedArray(collect(a), axes(a))
+
 Base.dataids(A::BlockedArray) = Base.dataids(A.blocks)
 
 ###########################

src/blocks.jl

@@ -93,6 +93,21 @@ This is broken for now. See: https://github.com/JuliaArrays/BlockArrays.jl/issue
     a
 end
 
+# AbstractArray version of `Iterators.product`.
+# https://en.wikipedia.org/wiki/Cartesian_product
+# https://github.com/lazyLibraries/ProductArrays.jl
+# https://github.com/JuliaData/SplitApplyCombine.jl#productviewf-a-b
+# https://github.com/JuliaArrays/MappedArrays.jl/pull/42
+struct ProductArray{T,N,V<:Tuple{Vararg{AbstractVector,N}}} <: AbstractArray{T,N}
+    vectors::V
+end
+ProductArray(vectors::Vararg{AbstractVector,N}) where {N} =
+    ProductArray{Tuple{map(eltype, vectors)...},N,typeof(vectors)}(vectors)
+Base.size(p::ProductArray) = map(length, p.vectors)
+Base.axes(p::ProductArray) = map(Base.axes1, p.vectors)
+@propagate_inbounds getindex(p::ProductArray{T,N}, I::Vararg{Int,N}) where {T,N} =
+    map((v, i) -> v[i], p.vectors, I)
+
 """
     blocksizes(A::AbstractArray)
     blocksizes(A::AbstractArray, d::Integer)
@@ -110,7 +125,7 @@ julia> A = BlockArray(ones(3,3),[2,1],[1,1,1])
  1.0  │  1.0  │  1.0
 
 julia> blocksizes(A)
-2×3 BlockArrays.BlockSizes{Tuple{Int64, Int64}, 2, BlockMatrix{Float64, Matrix{Matrix{Float64}}, Tuple{BlockedOneTo{Int64, Vector{Int64}}, BlockedOneTo{Int64, Vector{Int64}}}}}:
+2×3 BlockArrays.ProductArray{Tuple{Int64, Int64}, 2, Tuple{Vector{Int64}, Vector{Int64}}}:
  (2, 1)  (2, 1)  (2, 1)
  (1, 1)  (1, 1)  (1, 1)
 
@@ -124,16 +139,80 @@ julia> blocksizes(A,2)
  1
 ```
 """
-blocksizes(A::AbstractArray) = BlockSizes(A)
+blocksizes(A::AbstractArray) = ProductArray(map(blocklengths, axes(A))...)
 @inline blocksizes(A::AbstractArray, d::Integer) = blocklengths(axes(A, d))
 
-struct BlockSizes{T,N,A<:AbstractArray{<:Any,N}} <: AbstractArray{T,N}
+"""
+    blocklengths(A::AbstractArray)
+
+Return an iterator over the lengths of each block.
+See also blocksizes.
+
+# Examples
+```jldoctest
+julia> A = BlockArray(ones(3,3),[2,1],[1,1,1])
+2×3-blocked 3×3 BlockMatrix{Float64}:
+ 1.0  │  1.0  │  1.0
+ 1.0  │  1.0  │  1.0
+ ─────┼───────┼─────
+ 1.0  │  1.0  │  1.0
+
+julia> blocklengths(A)
+2×3 BlockArrays.BlockLengths{Int64, 2, BlockMatrix{Float64, Matrix{Matrix{Float64}}, Tuple{BlockedOneTo{Int64, Vector{Int64}}, BlockedOneTo{Int64, Vector{Int64}}}}}:
+ 2  2  2
+ 1  1  1
+
+julia> blocklengths(A)[1,2]
+2
+```
+"""
+blocklengths(A::AbstractArray) = BlockLengths(A)
+blocklengths(A::AbstractVector) = map(length, blocks(A))
+
+struct BlockLengths{T,N,A<:AbstractArray{<:Any,N}} <: AbstractArray{T,N}
     array::A
 end
-BlockSizes(a::AbstractArray{<:Any,N}) where {N} =
-    BlockSizes{Tuple{eltype.(axes(a))...},N,typeof(a)}(a)
+BlockLengths(a::AbstractArray{<:Any,N}) where {N} =
+    BlockLengths{typeof(length(a)),N,typeof(a)}(a)
 
-size(bs::BlockSizes) = blocksize(bs.array)
-axes(bs::BlockSizes) = map(br -> Int.(br), blockaxes(bs.array))
-@propagate_inbounds getindex(a::BlockSizes{T,N}, i::Vararg{Int,N}) where {T,N} =
-    size(view(a.array, Block.(i)...))
+size(bs::BlockLengths) = blocksize(bs.array)
+axes(bs::BlockLengths) = map(br -> Int.(br), blockaxes(bs.array))
+@propagate_inbounds getindex(a::BlockLengths{T,N}, i::Vararg{Int,N}) where {T,N} =
+    length(view(a.array, Block.(i)...))
+
+"""
+    eachblockaxes(A::AbstractArray)
+    eachblockaxes(A::AbstractArray, d::Integer)
+
+Return an iterator over the axes of each block.
+See also blocksizes and blocklengths.
+
+# Examples
+```jldoctest
+julia> A = BlockArray(ones(3,3),[2,1],[1,1,1])
+2×3-blocked 3×3 BlockMatrix{Float64}:
+ 1.0  │  1.0  │  1.0
+ 1.0  │  1.0  │  1.0
+ ─────┼───────┼─────
+ 1.0  │  1.0  │  1.0
+
+julia> eachblockaxes(A)
+2×3 BlockArrays.ProductArray{Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, 2, Tuple{Vector{Base.OneTo{Int64}}, Vector{Base.OneTo{Int64}}}}:
+ (Base.OneTo(2), Base.OneTo(1))  …  (Base.OneTo(2), Base.OneTo(1))
+ (Base.OneTo(1), Base.OneTo(1))     (Base.OneTo(1), Base.OneTo(1))
+
+julia> eachblockaxes(A)[1,2]
+(Base.OneTo(2), Base.OneTo(1))
+
+julia> eachblockaxes(A,2)
+3-element Vector{Base.OneTo{Int64}}:
+ Base.OneTo(1)
+ Base.OneTo(1)
+ Base.OneTo(1)
+```
+"""
+eachblockaxes(A::AbstractArray) =
+    ProductArray(map(ax -> map(Base.axes1, blocks(ax)), axes(A))...)
+eachblockaxes(A::AbstractVector) = map(axes, blocks(Base.axes1(A)))
+eachblockaxes(A::AbstractArray, d::Integer) = map(Base.axes1, blocks(axes(A, d)))
+eachblockaxes1(A::AbstractArray) = map(Base.axes1, blocks(Base.axes1(A)))

src/views.jl

@@ -26,7 +26,13 @@ to_index(::BlockRange) = throw(ArgumentError("BlockRange must be converted by to
 
 @inline to_indices(A, inds, I::Tuple{Block{1}, Vararg{Any}}) =
     (unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)
-@inline to_indices(A, inds, I::Tuple{BlockRange{1,R}, Vararg{Any}}) where R =
+@inline to_indices(A, inds, I::Tuple{BlockRange{1}, Vararg{Any}}) =
+    (unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)
+@inline to_indices(A, inds, I::Tuple{AbstractVector{<:Block{1}}, Vararg{Any}}) =
+    (unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)
+@inline to_indices(A, inds, I::Tuple{AbstractVector{<:BlockRange{1}}, Vararg{Any}}) =
+    (unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)
+@inline to_indices(A, inds, I::Tuple{AbstractVector{<:AbstractVector{<:Block{1}}}, Vararg{Any}}) =
     (unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)
 @inline to_indices(A, inds, I::Tuple{BlockIndex{1}, Vararg{Any}}) =
     (inds[1][I[1]], to_indices(A, _maybetail(inds), tail(I))...)
@@ -38,12 +44,20 @@ to_index(::BlockRange) = throw(ArgumentError("BlockRange must be converted by to
     (unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)
 @inline to_indices(A, inds, I::Tuple{AbstractVector{<:BlockIndices{1}}, Vararg{Any}}) =
     (unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)
+@inline to_indices(A, inds, I::Tuple{AbstractVector{<:AbstractVector{<:BlockIndex{1}}}, Vararg{Any}}) =
+    (unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)
+@inline to_indices(A, inds, I::Tuple{AbstractVector{<:AbstractVector{<:BlockIndices{1}}}, Vararg{Any}}) =
+    (unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)
+@inline to_indices(A, inds, I::Tuple{AbstractVector{<:AbstractVector{<:AbstractVector{<:BlockIndex{1}}}}, Vararg{Any}}) =
+    (unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)
 
 
 # splat out higher dimensional blocks
 # this mimics view of a CartesianIndex
 @inline to_indices(A, inds, I::Tuple{Block, Vararg{Any}}) =
     to_indices(A, inds, (Block.(I[1].n)..., tail(I)...))
+@inline to_indices(A, inds, I::Tuple{BlockRange, Vararg{Any}}) =
+    to_indices(A, inds, (BlockRange.(tuple.(I[1].indices))..., tail(I)...))
 @inline to_indices(A, inds, I::Tuple{BlockIndex, Vararg{Any}}) =
     to_indices(A, inds, (BlockIndex.(I[1].I, I[1].α)..., tail(I)...))
 @inline to_indices(A, inds, I::Tuple{BlockIndices, Vararg{Any}}) =
@@ -57,8 +71,57 @@ to_index(::BlockRange) = throw(ArgumentError("BlockRange must be converted by to
 @inline to_indices(A, I::Tuple{Block, Vararg{Any}}) = to_indices(A, axes(A), I)
 @inline to_indices(A, I::Tuple{BlockRange, Vararg{Any}}) = to_indices(A, axes(A), I)
 @inline to_indices(A, I::Tuple{AbstractVector{<:Block{1}}, Vararg{Any}}) = to_indices(A, axes(A), I)
+@inline to_indices(A, I::Tuple{AbstractVector{<:BlockRange{1}}, Vararg{Any}}) = to_indices(A, axes(A), I)
+@inline to_indices(A, I::Tuple{AbstractVector{<:AbstractVector{<:Block{1}}}, Vararg{Any}}) = to_indices(A, axes(A), I)
 @inline to_indices(A, I::Tuple{AbstractVector{<:BlockIndex{1}}, Vararg{Any}}) = to_indices(A, axes(A), I)
 @inline to_indices(A, I::Tuple{AbstractVector{<:BlockIndices{1}}, Vararg{Any}}) = to_indices(A, axes(A), I)
+@inline to_indices(A, I::Tuple{AbstractVector{<:AbstractVector{<:BlockIndex{1}}}, Vararg{Any}}) = to_indices(A, axes(A), I)
+@inline to_indices(A, I::Tuple{AbstractVector{<:AbstractVector{<:BlockIndices{1}}}, Vararg{Any}}) = to_indices(A, axes(A), I)
+@inline to_indices(A, I::Tuple{AbstractVector{<:AbstractVector{<:AbstractVector{<:BlockIndex{1}}}}, Vararg{Any}}) = to_indices(A, axes(A), I)
+
+## BlockedLogicalIndex
+# Blocked version of `LogicalIndex`:
+# https://github.com/JuliaLang/julia/blob/3e2f90fbb8f6b0651f2601d7599c55d4e3efd496/base/multidimensional.jl#L819-L831
+const BlockedLogicalIndex{T,R<:LogicalIndex{T},BS<:Tuple{AbstractUnitRange{<:Integer}}} = BlockedVector{T,R,BS}
+function BlockedLogicalIndex(I::AbstractVector{Bool})
+    blocklengths = map(b -> count(view(I, b)), BlockRange(I))
+    return BlockedVector(LogicalIndex(I), blocklengths)
+end
+# https://github.com/JuliaLang/julia/blob/3e2f90fbb8f6b0651f2601d7599c55d4e3efd496/base/multidimensional.jl#L838-L839
+show(io::IO, r::BlockedLogicalIndex) = print(io, blockcollect(r))
+print_array(io::IO, X::BlockedLogicalIndex) = print_array(io, blockcollect(X))
+
+# Blocked version of `to_index(::AbstractArray{Bool})`:
+# https://github.com/JuliaLang/julia/blob/3e2f90fbb8f6b0651f2601d7599c55d4e3efd496/base/indices.jl#L309
+function to_index(I::AbstractBlockVector{Bool})
+    return BlockedLogicalIndex(I)
+end
+
+# Blocked version of `collect(::LogicalIndex)`:
+# https://github.com/JuliaLang/julia/blob/3e2f90fbb8f6b0651f2601d7599c55d4e3efd496/base/multidimensional.jl#L837
+# Without this definition, `collect` will try to call `getindex` on the `LogicalIndex`
+# which isn't defined.
+collect(I::BlockedLogicalIndex) = collect(I.blocks)
+
+# Iteration of BlockedLogicalIndex is just iteration over the underlying
+# LogicalIndex, which is implemented here:
+# https://github.com/JuliaLang/julia/blob/3e2f90fbb8f6b0651f2601d7599c55d4e3efd496/base/multidimensional.jl#L840-L890
+@inline iterate(I::BlockedLogicalIndex) = iterate(I.blocks)
+@inline iterate(I::BlockedLogicalIndex, s) = iterate(I.blocks, s)
+
+## Boundscheck for BlockLogicalindex
+# Like for LogicalIndex, map all calls to mask:
+# https://github.com/JuliaLang/julia/blob/3e2f90fbb8f6b0651f2601d7599c55d4e3efd496/base/multidimensional.jl#L892-L897
+checkbounds(::Type{Bool}, A::AbstractArray, i::BlockedLogicalIndex) = checkbounds(Bool, A, i.blocks.mask)
+# `checkbounds_indices` has been handled via `I::AbstractArray` fallback
+checkindex(::Type{Bool}, inds::AbstractUnitRange, i::BlockedLogicalIndex) = checkindex(Bool, inds, i.blocks.mask)
+checkindex(::Type{Bool}, inds::Tuple, i::BlockedLogicalIndex) = checkindex(Bool, inds, i.blocks.mask)
+
+# Instantiate the BlockedLogicalIndex when constructing a SubArray, similar to
+# `ensure_indexable(I::Tuple{LogicalIndex,Vararg{Any}})`:
+# https://github.com/JuliaLang/julia/blob/3e2f90fbb8f6b0651f2601d7599c55d4e3efd496/base/multidimensional.jl#L918
+@inline ensure_indexable(I::Tuple{BlockedLogicalIndex,Vararg{Any}}) =
+    (blockcollect(I[1]), ensure_indexable(tail(I))...)
 
 @propagate_inbounds reindex(idxs::Tuple{BlockSlice{<:BlockRange}, Vararg{Any}},
         subidxs::Tuple{BlockSlice{<:BlockIndices}, Vararg{Any}}) =

test/test_blockarrays.jl

@@ -1,7 +1,7 @@
 module TestBlockArrays
 
 using SparseArrays, BlockArrays, FillArrays, LinearAlgebra, Test, OffsetArrays, Images
-import BlockArrays: _BlockArray
+import BlockArrays: _BlockArray, blockcollect
 
 const Fill = FillArrays.Fill
 
@@ -255,6 +255,32 @@ end
             @test zero(b) isa typeof(b)
         end
 
+        @testset "blockcollect" begin
+            a = randn(6, 6)
+            @test blockcollect(a) == a
+            @test blockcollect(a) ≢ a
+            @test blockcollect(a).blocks ≢ a
+            # TODO: Maybe special case this to call `collect` and return a `Matrix`?
+            @test blockcollect(a) isa BlockedMatrix{Float64,Matrix{Float64}}
+            @test blockisequal(axes(blockcollect(a)), axes(a))
+            @test blocksize(blockcollect(a)) == (1, 1)
+
+            b = BlockedArray(randn(6, 6), [3, 3], [3, 3])
+            @test blockcollect(b) == b
+            @test blockcollect(b) ≢ b
+            @test blockcollect(b).blocks ≢ b
+            @test blockcollect(b) isa BlockedMatrix{Float64,Matrix{Float64}}
+            @test blockisequal(axes(blockcollect(b)), axes(b))
+            @test blocksize(blockcollect(b)) == (2, 2)
+
+            c = BlockArray(randn(6, 6), [3, 3], [3, 3])
+            @test blockcollect(c) == c
+            @test blockcollect(c) ≢ c
+            @test blockcollect(c) isa BlockedMatrix{Float64,Matrix{Float64}}
+            @test blockisequal(axes(blockcollect(c)), axes(c))
+            @test blocksize(blockcollect(c)) == (2, 2)
+        end
+
         @test_throws DimensionMismatch BlockArray([1,2,3],[1,1])
 
         @testset "mortar" begin
@@ -1085,6 +1111,21 @@ end
         @test a[[Block(1)[1:2], Block(2)[1:2]], [Block(1)[1:2], Block(2)[1:2]]] == [a[Block(1,1)[1:2,1:2]] a[Block(1,2)[1:2,1:2]]; a[Block(2,1)[1:2,1:2]] a[Block(2,2)[1:2,1:2]]]
         @test a[[Block(1)[1], Block(2)[2]], [Block(1)[1:2], Block(2)[1:2]]] == [a[Block(1)[1],Block(1)[1:2]]' a[Block(1)[1], Block(2)[1:2]]'; a[Block(2)[2],Block(1)[1:2]]' a[Block(2)[2], Block(2)[1:2]]']
     end
+    @testset "Blocked block-vector indexing (#359)" begin
+        for a in (BlockArray(randn(14, 14), 2:5, 2:5), BlockedArray(randn(14, 14), 2:5, 2:5))
+            for I in (
+                [Block.(1:2), Block.(3:4)],
+                [[Block(1), Block(3)], [Block(2), Block(4)]],
+                [[Block(1)[1:2], Block(3)[1:2]], [Block(2)[1:2], Block(4)[1:2]]],
+                [[[Block(1)[1], Block(1)[2]], [Block(3)[1], Block(3)[2]]], [[Block(2)[1], Block(2)[2]], [Block(4)[1], Block(4)[2]]]],
+            )
+                b = a[I, I]
+                for (i, j) in Iterators.product(1:length(I), 1:length(I))
+                    @test a[I[i], I[j]] == b[Block(i), Block(j)]
+                end
+            end
+        end
+    end
 end
 
 end # module