Add support for multiple layers of block arrays

[dlfivefifty]

It's sometimes useful to have multiple layers of blocks, see conversation in

SciML/MultiScaleArrays.jl#25

I think one way to do this would be to support some sort of BlockView that converts an AbstractBlockVector{T} into an AbstractVector{Vector{T}}:

struct BlockView{AT, p, BA <: AbstractBlockArray} <: AbstractArray{AT, p}
   data::BA
end

_BlockView(data::BA) where BA<:AbstractBlockArray{T,p} where T = BlockView{Array{T,p},p,BA}(data)

getindex(A::BlockView, i...) = A[Block(i...)]
size(A::BlockView) = nblocks(A)

Then one create a hierarchical block array view (here we have sub-blocks of length 10 and super-blocks of length 5):

A_sub = PseudoBlockArray(randn(10*5*n), fill(10,5n))
A = PseudoBlockArray(BlockView(A_sub), fill(5,n))
A[Block(2)][3]  # returns the 3rd sub-block of the 2nd super-block

I'm still working out how to get it to work with A[Block(2)][Block(3)]...

[dlfivefifty]

I'm still working out how to get it to work with A[Block(2)][Block(3)]

This part is worked out already: see example of a triple kronecker product

BlockArrays.jl/test/test_blockproduct.jl

Line 116 in ac2fbec

@test k̄[Block(1)][Block(1)] == a[1]*b[1]*c

So in some sense we have blue-print for what to do but at the moment it's a bit of a hack. I think this example should be modified to incorporate the nested block structure in the axes.

One potentially elegant idea: have a BlockedUnitRange whose lasts are a blocked vector. That is, to represent a double-blocked vector, we would want something like this to work:

julia> a = blockedrange(Fill(3,(blockedrange(Fill(2,2)),)))
4-blocked 12-element BlockedUnitRange{StepRange{Int64,Int64}}:
  1
  2
  3
 ┄┄
  4
  5
  6
 ──
  7
  8
  9
 ┄┄
 10
 11
 12

julia> a[Block(1)]
julia> blockedrange(Fill(3,(blockedrange(Fill(2,1)),)))
2-blocked 6-element BlockedUnitRange{StepRange{Int64,Int64}}:
 1
 2
 3
 ─
 4
 5
 6

julia> a[Block(1)][Block(2)]
4:6

@ChrisRackauckas Are you still interested in something like BandedBlockBandedBlockBandedMatrix for 3D PDEs? I think getting it working is not too far off now since we already have:

julia> n = 5; Δ = SymTridiagonal(zeros(n), ones(n-1))/2
5×5 SymTridiagonal{Float64,Array{Float64,1}}:
 0.0  0.5   ⋅    ⋅    ⋅ 
 0.5  0.0  0.5   ⋅    ⋅ 
  ⋅   0.5  0.0  0.5   ⋅ 
  ⋅    ⋅   0.5  0.0  0.5
  ⋅    ⋅    ⋅   0.5  0.0

julia> Dx = blockkron(D²,Eye(n),Eye(n));

julia> Dy = blockkron(Eye(n),D²,Eye(n));

julia> Dz = blockkron(Eye(n),Eye(n),D²);

julia> Δ = Dx + Dy + Dz
3×3-blocked 27×27 BlockArray{Float64,2}:
 0.0  0.5  0.0  0.5  0.0  0.0  0.0  0.0  0.0  │  0.5  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  │  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.5  0.0  0.5  0.0  0.5  0.0  0.0  0.0  0.0  │  0.0  0.5  0.0  0.0  0.0  0.0  0.0  0.0  0.0  │  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.5  0.0  0.0  0.0  0.5  0.0  0.0  0.0  │  0.0  0.0  0.5  0.0  0.0  0.0  0.0  0.0  0.0  │  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.5  0.0  0.0  0.0  0.5  0.0  0.5  0.0  0.0  │  0.0  0.0  0.0  0.5  0.0  0.0  0.0  0.0  0.0  │  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.5  0.0  0.5  0.0  0.5  0.0  0.5  0.0  │  0.0  0.0  0.0  0.0  0.5  0.0  0.0  0.0  0.0  │  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.5  0.0  0.5  0.0  0.0  0.0  0.5  │  0.0  0.0  0.0  0.0  0.0  0.5  0.0  0.0  0.0  │  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.5  0.0  0.0  0.0  0.5  0.0  │  0.0  0.0  0.0  0.0  0.0  0.0  0.5  0.0  0.0  │  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.5  0.0  0.5  0.0  0.5  │  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.5  0.0  │  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.5  0.0  0.5  0.0  │  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.5  │  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 ─────────────────────────────────────────────┼───────────────────────────────────────────────┼─────────────────────────────────────────────
 0.5  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  │  0.0  0.5  0.0  0.5  0.0  0.0  0.0  0.0  0.0  │  0.5  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.5  0.0  0.0  0.0  0.0  0.0  0.0  0.0  │  0.5  0.0  0.5  0.0  0.5  0.0  0.0  0.0  0.0  │  0.0  0.5  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.5  0.0  0.0  0.0  0.0  0.0  0.0  │  0.0  0.5  0.0  0.0  0.0  0.5  0.0  0.0  0.0  │  0.0  0.0  0.5  0.0  0.0  0.0  0.0  0.0  0.0
 ⋮                        ⋮                        ⋮                        ⋮                        ⋮                        ⋮    
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.5  0.0  │  0.0  0.0  0.0  0.0  0.5  0.0  0.5  0.0  0.5  │  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.5  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.5  │  0.0  0.0  0.0  0.0  0.0  0.5  0.0  0.5  0.0  │  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.5
 ─────────────────────────────────────────────┼───────────────────────────────────────────────┼─────────────────────────────────────────────
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  │  0.5  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  │  0.0  0.5  0.0  0.5  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  │  0.0  0.5  0.0  0.0  0.0  0.0  0.0  0.0  0.0  │  0.5  0.0  0.5  0.0  0.5  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  │  0.0  0.0  0.5  0.0  0.0  0.0  0.0  0.0  0.0  │  0.0  0.5  0.0  0.0  0.0  0.5  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  │  0.0  0.0  0.0  0.5  0.0  0.0  0.0  0.0  0.0  │  0.5  0.0  0.0  0.0  0.5  0.0  0.5  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  │  0.0  0.0  0.0  0.0  0.5  0.0  0.0  0.0  0.0  │  0.0  0.5  0.0  0.5  0.0  0.5  0.0  0.5  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  │  0.0  0.0  0.0  0.0  0.0  0.5  0.0  0.0  0.0  │  0.0  0.0  0.5  0.0  0.5  0.0  0.0  0.0  0.5
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  │  0.0  0.0  0.0  0.0  0.0  0.0  0.5  0.0  0.0  │  0.0  0.0  0.0  0.5  0.0  0.0  0.0  0.5  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  │  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.5  0.0  │  0.0  0.0  0.0  0.0  0.5  0.0  0.5  0.0  0.5
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  │  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.5  │  0.0  0.0  0.0  0.0  0.0  0.5  0.0  0.5  0.0

So with a few tweeks we'll have the dense construction working with nested block structure, and it shouldn't be too hard to get the sparse banded-block-banded-block-banded matrix working. (I think the above proposal "stacks" so we could potentially get (banded-block)^3-banded-matrix for 4D PDEs)