BlockArrays.jl icon indicating copy to clipboard operation
BlockArrays.jl copied to clipboard

Indexing with `Vector{<:BlockIndexRange{1}}`

Open mtfishman opened this issue 1 year ago • 1 comments

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.

mtfishman avatar Mar 28 '24 15:03 mtfishman

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.

mtfishman avatar Mar 28 '24 15:03 mtfishman

Closed by #1489.

mtfishman avatar Jun 18 '24 22:06 mtfishman