Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Linear indexing with Block broken? #120

Closed
tkf opened this issue May 21, 2020 · 2 comments
Closed

Linear indexing with Block broken? #120

tkf opened this issue May 21, 2020 · 2 comments

Comments

@tkf
Copy link
Member

tkf commented May 21, 2020

For a matrix M, we have M[1, 1] == M[1]. However, similar indexing with Block seems to be broken:

julia> A = PseudoBlockArray(rand(2,3), [1,1], [2,1])
2×2-blocked 2×3 PseudoBlockArray{Float64,2}:
 0.969606  0.5692780.767362
 ────────────────────┼──────────
 0.331377  0.404040.549861

julia> A[Block(1)]
1-element Array{Float64,1}:
 0.9696062405354238

julia> A[Block(1), Block(1)]
1×2 view(::Array{Float64,2}, 1:1, 1:2) with eltype Float64:
 0.969606  0.569278

I expect A[Block(1)] and A[Block(1), Block(1)] above to return equivalent matrices. (Or throwing with A[Block(1)] might be OK.)


This is also true with BlockArray:

julia> bs = 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(bs)
2×2-blocked 3×5 BlockArray{Float64,2}:
 1.0  1.0  1.02.0  2.0
 ───────────────┼──────────
 3.0  3.0  3.04.0  4.0
 3.0  3.0  3.04.0  4.0

julia> a[Block(1), Block(1)]
1×3 Array{Float64,2}:
 1.0  1.0  1.0

julia> a[Block(1)]
1-element Array{Float64,1}:
 1.0
@dlfivefifty
Copy link
Member

Yes this is not intentional, didn't think about this case when designing the code, and surprised it doesn't error out.

The way I would fix this is to first get it working on views and then use ArrayLayouts.jl approach of materializing a view (adding getindex(a::AbstractArray, K::Block{1}) = layout_getindex(a, K) will materialize a view).

What goes wrong with views is that it tries to reshape to a vector:

julia> view(a,Block(1))
1-element view(reshape(::PseudoBlockArray{Float64,2,Array{Float64,2},Tuple{BlockedUnitRange{Array{Int64,1}},BlockedUnitRange{Array{Int64,1}}}}, 15), BlockSlice(Block(1),1:1)) with eltype Float64:
 1.0

Probably an overload of _maybe_reshape_parent will help, though we'll need to add support for block ReshapedArrays.

@dlfivefifty
Copy link
Member

Fixed by #156

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants