New interface blocks(a) for array-of-arrays view

[tkf]

This is an RFC to add:

blocks(a::AbstractArray{T,N}) :: AbstractArray{<:AbstractArray{T,N},N}

Return the array-of-arrays view to a such that

blocks(a)[i₁, i₂, ..., iₙ] == a[Block(i₁), Block(i₂), ..., Block(iₙ)]

This function does not copy the blocks and give a mutable viwe to the original
array. This is an "inverse" of mortar.

Examples

julia> bs1 = permutedims(reshape([
               1ones(1, 3), 2ones(1, 2),
               3ones(2, 3), 4ones(2, 2),
           ], (2, 2)))
2×2 Array{Array{Float64,2},2}:
 [1.0 1.0 1.0]               [2.0 2.0]
 [3.0 3.0 3.0; 3.0 3.0 3.0]  [4.0 4.0; 4.0 4.0]

julia> a = mortar(bs1)
2×2-blocked 3×5 BlockArray{Float64,2}:
 1.0  1.0  1.0  │  2.0  2.0
 ───────────────┼──────────
 3.0  3.0  3.0  │  4.0  4.0
 3.0  3.0  3.0  │  4.0  4.0

julia> bs2 = blocks(a)
2×2 Array{Array{Float64,2},2}:
 [1.0 1.0 1.0]               [2.0 2.0]
 [3.0 3.0 3.0; 3.0 3.0 3.0]  [4.0 4.0; 4.0 4.0]

julia> bs1 == bs2
true

julia> bs2[1, 1] .*= 100;

julia> a  # in-place mutation is reflected to the block array
2×2-blocked 3×5 BlockArray{Float64,2}:
 100.0  100.0  100.0  │  2.0  2.0
 ─────────────────────┼──────────
   3.0    3.0    3.0  │  4.0  4.0
   3.0    3.0    3.0  │  4.0  4.0

It is useful for passing block arrays to some generic function that does not know the block array interfaces. For example:

function sum_array_of_arrays(arrays)
    acc = zero(eltype(eltype(arrays)))
    for a in arrays
        for x in a
            acc += x
        end
    end
    return acc
end

After this PR, we can just do sum_array_of_arrays(blocks(block_array)).

[codecov]

Codecov Report

Merging #117 into master will increase coverage by 0.81%.
The diff coverage is 89.65%.

@@            Coverage Diff             @@
##           master     #117      +/-   ##
==========================================
+ Coverage   84.42%   85.24%   +0.81%     
==========================================
  Files          12       15       +3     
  Lines         822      935     +113     
==========================================
+ Hits          694      797     +103     
- Misses        128      138      +10     

Impacted Files	Coverage Δ
src/BlockArrays.jl	100.00% <ø> (ø)
src/blockarray.jl	97.18% <ø> (ø)
src/blocks.jl	86.95% <86.95%> (ø)
src/pseudo_blockarray.jl	89.28% <100.00%> (+0.60%)	⬆️
src/blockdeque.jl	91.13% <0.00%> (ø)
src/blockreduce.jl	100.00% <0.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 2424ea1...494a95c. Read the comment docs.

[dlfivefifty]

Haha, you literally read my mind: I coincidentally added blocks to #119

Let's get this PR merged first since it's more general. Though we need to add blocks(::Adjoint), blocks(::Transpose), blocks(::SubArray) as in PR #119

[tkf]

Do you want special handling for SubArrays? This currently already works because of the generic fallback.

[dlfivefifty]

Yes because the generic fallback is not efficient. My primary reason for adding this (which got muddled by a tangential desire for nice printing) is for fast paths to convert subviews of a BlockMatrix{T,<:Tridiagonal} to a BlockBandedMatrix see JuliaLinearAlgebra/BlockBandedMatrices.jl#73

[tkf]

OK. So I copied your code into this PR and add a few tests.

[dlfivefifty]

Looks good to me. Can you get the coverage up and we can merge?

[tkf]

OK now Codecov seems to be happy with it.

[tkf]

I can't hit this method with the test. But manually trying this branch with BlockBandedMatrix seems to hit this method correctly.

[dlfivefifty]

Happy for me to merge this?

[tkf]

I added everything I wanted to add. Please go ahead and merge this.

[dlfivefifty]

Wait, why was the Compat dependency added?

[tkf]

That's

BlockArrays.jl/src/BlockArrays.jl

Lines 42 to 44 in 42ae241

	if !@isdefined(only)
	using Compat: only
	end

Initially, I inlined the simple definition but I thought it'd be better to use Compat. Otherwise, I needed to test the error case to make codecov happy. I can do this but I was lazy :)

[dlfivefifty]

Ahh, sorry didn’t see the using statement. That’s fine