BlockArrays.jl
BlockArrays.jl copied to clipboard
Indexing with `Vector{<:BlockIndexRange{1}}`
Related to #184, I think this would be very helpful to define:
julia> a = BlockArray(randn(6, 6), [3, 3], [3, 3])
2×2-blocked 6×6 BlockMatrix{Float64}:
1.107 0.791783 -0.445287 │ 0.04164 -1.53479 0.0830437
-0.192029 -0.56695 -1.10488 │ -0.727144 0.0331893 0.478663
0.5843 -1.18295 -3.09405 │ -0.0200601 0.663132 0.130971
──────────────────────────────────┼────────────────────────────────────
-2.38948 -1.62213 -0.0218833 │ -0.45099 0.596954 0.163427
0.598626 -0.672888 1.98241 │ -0.317914 1.12006 -0.0961988
1.07971 -0.127525 -0.13917 │ 1.83828 0.258948 1.1028
julia> a[[Block(1)[1:2], Block(2)[1:2]], [Block(1)[1:2], Block(2)[1:2]]]
2×2-blocked 4×4 BlockMatrix{Float64}:
1.107 0.791783 │ 0.04164 -1.53479
-0.192029 -0.56695 │ -0.727144 0.0331893
──────────────────────┼───────────────────────
-2.38948 -1.62213 │ -0.45099 0.596954
0.598626 -0.672888 │ -0.317914 1.12006
which currently errors. It would be equivalent to:
julia> mortar(reshape([a[Block(1, 1)[1:2, 1:2]], a[Block(2, 1)[1:2, 1:2]], a[Block(1, 2)[1:2, 1:2]], a[Block(2, 2)[1:2, 1:2]]], 2, 2))
2×2-blocked 4×4 BlockMatrix{Float64}:
1.107 0.791783 │ 0.04164 -1.53479
-0.192029 -0.56695 │ -0.727144 0.0331893
──────────────────────┼───────────────────────
-2.38948 -1.62213 │ -0.45099 0.596954
0.598626 -0.672888 │ -0.317914 1.12006
and related to:
julia> a[[1, 2, 4, 5], [1, 2, 4, 5]]
4×4 Matrix{Float64}:
1.107 0.791783 0.04164 -1.53479
-0.192029 -0.56695 -0.727144 0.0331893
-2.38948 -1.62213 -0.45099 0.596954
0.598626 -0.672888 -0.317914 1.12006
but it preserves the blocking structure.
I'm not sure if a[[Block(1)[1:2], Block(2)[1:2]], [Block(1)[1:2], Block(2)[1:2]]] is the correct "spelling" for this kind of operation, or if it warrants a new type, but the analogy with a[[Block(1), Block(2)], [Block(1), Block(2)]] makes it seem reasonable to me.
We would have some nice use cases for this, for example if you perform a low rank approximation of a block diagonal matrix which preserves the block structure, you want to perform a low rank approximation of the blocks separately (where potentially each block could have an associated rank). That corresponds to this kind of slicing operation, where you slice each block to the corresponding rank, or drop a block if the block rank goes to zero.
Closed by #1489.