Add `blockvcat`, `blockhcat`, and `blockhvcat`

[dlfivefifty]

We should add blockvcat(A,B,...) that produces a blocked array. For example blockvcat(1:5,10:12)[Block(1)] == 1:5. But also supporting the case whereA and B are blocked, where the block structure is concatenated. An initial prototype is here

https://github.com/JuliaMatrices/LazyBandedMatrices.jl/blob/master/src/blockconcat.jl

Following blockkron we first produce lazy versions. We may also want to add a macro, something like:

@blocked [A; B] # lowers to blockvcat(A,B)
@blocked [A B] # lowers to blockhcat(A,B)
@blocked [A B; C D] # lowers toblockhvcat(2, A, B, C, D)

[putianyi889]

We could override Base.vcat , Base.hcat and Base.hvcat for concrete block array types. When block dimensions or types don't match, we may issue a warning and invoke the original Base function. For banded arrays, we should issue an error to avoid memory issues.

Each concrete type needs a set of these functions.

[dlfivefifty]

The problem with that is that we are treating Array as a block-array with one block. So whatever behaviour vcat and hcat have should be consistent.

[putianyi889]

Currently applying cats will convert the block arrays to concrete arrays.

julia> A=BlockArray{Float64}(zeros(7,9), [1,2,4], [3,2,4]);

julia> blocksizes(A)
([1, 2, 4], [3, 2, 4])

julia> blocksizes(vcat(A,A))
([14], [9])

If we implement cats for block arrays, then for blockcats we just need to do a conversion first.

blockvcat(1:5,10:12)=vcat(BlockArray(1:5),BlockArray(10:12))

[dlfivefifty]

But I don't think we want this, because we want consistency, that is, we want:

a = 1:5
b = 10:12
blocksizes(vcat(a, b)) == blocksizes(vcat(BlockArray(a), BlockArray(b)))

The reason this is important is that many block arrays are not <: AbstractBlockArray, e.g., subviews of block arrays.

[putianyi889]

True. So we may need both blockcat and cat. blockcat indicates that the arguments are in distinct blocks while cat performs the natural concatenation.

Examples that blockcat is needed:

• An array has block structures but is not typical AbstractBlockArray. E.g. views of block arrays (SubArray, Adjoint, etc.).
• The concatenation is on an extra dimension, e.g. in cat(BlockMatrix,BlockMatrix,dims=3), the user has to decide whether the two matrices belong to different blocks in the 3rd dimension.
• Concatenation of two common arrays, e.g. in cat(Matrix,Matrix), the user has to decide whether the two matrices are considered as two blocks. The same for block arrays which has only one block in the concatenation dimension. (The intention of this issue)

[putianyi889]

blockcat is more "blocked" than cat. blockcat always split the result when possible. cat also generates block arrays when it makes strong sense.

[putianyi889]

Should Base.OneTo and 1-blocked BlockedUnitRange, as axes of arrays, be treated equally?

For example

julia> A=BlockArray{Float64}(zeros(2,2),[1,1],[2])
2×1-blocked 2×2 BlockArray{Float64,2}:
 0.0  0.0
 ────────
 0.0  0.0

julia> axes(A,1)
2-blocked 2-element BlockedUnitRange{Array{Int64,1}}:
 1
 ─
 2

julia> axes(A,2)
1-blocked 2-element BlockedUnitRange{Array{Int64,1}}:
 1
 2

julia> axes(A,3)
Base.OneTo(1)

If I concatenate two BlockedUnitRange, I'd want two blocks. For Base.OneTo, I'd want no block.

[dlfivefifty]

yes. There's no such thing as "no block"

[putianyi889]

I mean unblocked.

[dlfivefifty]

unblocked == 1-blocked that's my point

[putianyi889]

My point was unblocked != 1-blocked. Anyway, when it comes to 1-blocked, we use cat to treat it as unblocked and blockcat to treat it as blocked. This should be the point we agree.

[putianyi889]

What result do you expect from vcat(blockrange([2,2]), blockrange([2,2])), blockrange([2,4,2]) or blockrange([2,2,2,2])? I suppose the former, and the latter will be from blockvcat.