JuliaLinearAlgebra / BlockBandedMatrices.jl

A Julia package for representing block-banded matrices and banded-block-banded matrices
https://julialinearalgebra.github.io/BlockBandedMatrices.jl/
MIT License
56 stars 13 forks source link

StackOverFlow Error #106

Open putianyi889 opened 3 years ago

putianyi889 commented 3 years ago
julia> A=BlockBandedMatrix(Zeros(1,1),[1],[1],(0,0))
1×1-blocked 1×1 BlockSkylineMatrix{Float64,Array{Float64,1},BlockBandedMatrices.BlockSkylineSizes{Tuple{BlockedUnitRange{Array{Int64,1}},BlockedUnitRange{Array{Int64,1}}},Fill{Int64,1,Tuple{Base.OneTo{Int64}}},Fill{Int64,1,Tuple{Base.OneTo{Int64}}},BandedMatrices.BandedMatrix{Int64,Array{Int64,2},Base.OneTo{Int64}},Array{Int64,1}}}:
 0.0

julia> A[Block(1)]
ERROR: StackOverflowError:
Stacktrace:
 [1] view(::BlockSkylineMatrix{Float64,Array{Float64,1},BlockBandedMatrices.BlockSkylineSizes{Tuple{BlockedUnitRange{Array{Int64,1}},BlockedUnitRange{Array{Int64,1}}},Fill{Int64,1,Tuple{Base.OneTo{Int64}}},Fill{Int64,1,Tuple{Base.OneTo{Int64}}},BandedMatrices.BandedMatrix{Int64,Array{Int64,2},Base.OneTo{Int64}},Array{Int64,1}}}, ::Block{1,Int64}) at C:\Users\pty\.julia\packages\BlockArrays\g01fJ\src\abstractblockarray.jl:129 (repeats 79984 times)

An error is expected, but not stack overflow maybe?

dlfivefifty commented 3 years ago

Same as https://github.com/JuliaArrays/BlockArrays.jl/issues/120

would love a PR fixing this

putianyi889 commented 3 years ago

What does this code serve for? In which case does the right hand side not refer back again?

@inline @propagate_inbounds Base.view(block_arr::AbstractBlockArray, block::Block{1}...) = view(block_arr, Block(block))

https://github.com/JuliaArrays/BlockArrays.jl/blob/master/src/abstractblockarray.jl#L129

putianyi889 commented 3 years ago

It would be super easy to achieve only A[Block(1)] but I'm not sure if it breaks flexibility of other codes.

dlfivefifty commented 3 years ago

What does this code serve for?

This converts, for example, view(A,Block(1), Block(2)) to view(A,Block(1,2)). Note if there is only one block and its a AbstractBlockVector then the line above is called.

This seems to work (leveraging the code for standard arrays using linear indexing):

julia> A = BlockArray(randn(3,3), [2,1], [2,1])
2×2-blocked 3×3 BlockArray{Float64,2}:
 -0.115769   0.690882  │  -1.04413 
  0.661802  -0.217254  │  -0.272209
 ──────────────────────┼───────────
  0.433854   0.501771  │  -0.884769

julia> Base.view(A::AbstractBlockArray, B::Block{1}) = view(A, Block(Base._ind2sub(map(Base.OneTo, blocksize(A)), Int(B))...))

julia> view(A, Block(3))
2×1 Matrix{Float64}:
 -1.0441288027748936
 -0.2722086284730829