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

Diagonal(BlockVector)[Block(i, i)] is not Diagonal #427

Open
mtfishman opened this issue Oct 31, 2024 · 2 comments
Open

Diagonal(BlockVector)[Block(i, i)] is not Diagonal #427

mtfishman opened this issue Oct 31, 2024 · 2 comments

Comments

@mtfishman
Copy link
Collaborator

mtfishman commented Oct 31, 2024

For example:

julia> using BlockArrays

julia> d = Diagonal(BlockedVector([1, 2, 3, 4], [2, 2]))
4×4 Diagonal{Int64, BlockedVector{Int64, Vector{Int64}, Tuple{BlockedOneTo{Int64, Vector{Int64}}}}} with indices BlockedOneTo([2, 4])×BlockedOneTo([2, 4]):
 1    
   2  
 ──────┼──────
   3  
     4

julia> d[Block(1, 1)]
2×2 Matrix{Int64}:
 1  0
 0  2

Shows up in #426.

@jishnub
Copy link
Member

jishnub commented Nov 1, 2024

In general, the off-diagonal blocks won't be square, so we won't have a type-stable result.

julia> d = Diagonal(BlockedVector(1:6, [2, 4]))
6×6 Diagonal{Int64, BlockedVector{Int64, UnitRange{Int64}, Tuple{BlockedOneTo{Int64, Vector{Int64}}}}} with indices BlockedOneTo([2, 6])×BlockedOneTo([2, 6]):
 1        
   2      
 ──────┼────────────
   3      
     4    
       5  
         6

@mtfishman
Copy link
Collaborator Author

That's the same as when you make a BlockDiagonal with Diagonal blocks:

julia> D = BlockDiagonal([Diagonal([1, 2]), Diagonal([3, 4])])
2×2-blocked 4×4 BlockMatrix{Int64, Diagonal{Diagonal{Int64, Vector{Int64}}, Vector{Diagonal{Int64, Vector{Int64}}}}, Tuple{BlockedOneTo{Int64, Vector{Int64}}, BlockedOneTo{Int64, Vector{Int64}}}}:
 1    
   2  
 ──────┼──────
   3  
     4

julia> D[Block(1, 1)]
2×2 Diagonal{Int64, Vector{Int64}}:
 1  
   2

julia> D[Block(1, 2)]
2×2 Matrix{Int64}:
 0  0
 0  0

julia> @code_warntype D[Block(1, 1)]
MethodInstance for getindex(::BlockMatrix{Int64, Diagonal{Diagonal{Int64, Vector{Int64}}, Vector{Diagonal{Int64, Vector{Int64}}}}, Tuple{BlockedOneTo{Int64, Vector{Int64}}, BlockedOneTo{Int64, Vector{Int64}}}}, ::Block{2, Int64})
  from getindex(A::AbstractArray{T, N}, block::Block{N}) where {T, N} @ BlockArrays ~/.julia/packages/BlockArrays/VccC2/src/abstractblockarray.jl:229
Static Parameters
  T = Int64
  N = 2
Arguments
  #self#::Core.Const(getindex)
  A::BlockMatrix{Int64, Diagonal{Diagonal{Int64, Vector{Int64}}, Vector{Diagonal{Int64, Vector{Int64}}}}, Tuple{BlockedOneTo{Int64, Vector{Int64}}, BlockedOneTo{Int64, Vector{Int64}}}}
  block::Block{2, Int64}
Body::Union{Diagonal{Int64, Vector{Int64}}, Matrix{Int64}}
1nothing%2 = BlockArrays.ArrayLayouts::Core.Const(ArrayLayouts)
│   %3 = Base.getproperty(%2, :layout_getindex)::Core.Const(ArrayLayouts.layout_getindex)
│   %4 = (%3)(A, block)::Union{Diagonal{Int64, Vector{Int64}}, Matrix{Int64}}
└──      return %4

It seems like BlockMatrix{T,<:Diagonal{<:Diagonal}} should act like Diagonal{T,<:BlockVector}. Also that seems like a pretty mild form of type instability.

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