diff --git a/Project.toml b/Project.toml index fbc68ec..e5851ae 100644 --- a/Project.toml +++ b/Project.toml @@ -19,7 +19,7 @@ BlockBandedMatricesSparseArraysExt = "SparseArrays" [compat] Aqua = "0.8" ArrayLayouts = "1" -BandedMatrices = "1" +BandedMatrices = "1.9.4" BlockArrays = "1" Documenter = "1" FillArrays = "1" diff --git a/src/BandedBlockBandedMatrix.jl b/src/BandedBlockBandedMatrix.jl index 770efd1..ad2f403 100644 --- a/src/BandedBlockBandedMatrix.jl +++ b/src/BandedBlockBandedMatrix.jl @@ -39,7 +39,7 @@ struct BandedBlockBandedMatrix{T, BLOCKS, RAXIS<:AbstractUnitRange{Int}} <: Abst end end -const DefaultBandedBlockBandedMatrix{T} = BandedBlockBandedMatrix{T, BlockedMatrix{T, Matrix{T}, NTuple{2,DefaultBlockAxis}}, DefaultBlockAxis} +const DefaultBandedBlockBandedMatrix{T} = BandedBlockBandedMatrix{T, <:BlockedMatrix{T, Matrix{T}}} @inline _BandedBlockBandedMatrix(data::AbstractMatrix, axes::NTuple{2,AbstractUnitRange{Int}}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) = _BandedBlockBandedMatrix(BlockedArray(data,(blockedrange(Fill(sum(λμ)+1,sum(lu)+1)),axes[2])), axes[1], lu, λμ) @@ -167,7 +167,7 @@ BandedBlockBandedMatrix{T,B}(m::Union{AbstractMatrix, UniformScaling}, BandedBlockBandedMatrix{T}(m::Union{AbstractMatrix, UniformScaling}, axes::NTuple{2,AbstractUnitRange{Int}}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where T = - DefaultBandedBlockBandedMatrix{T}(m, axes, lu, λμ) + BandedBlockBandedMatrix{T,BlockedMatrix{T,Matrix{T},typeof(_bbb_data_axes(axes[2],lu,λμ))},typeof(axes[1])}(m, axes, lu, λμ) BandedBlockBandedMatrix{T,B,R}(m::Union{AbstractMatrix, UniformScaling}, rdims::AbstractVector{Int}, cdims::AbstractVector{Int}, @@ -334,6 +334,39 @@ zeroblock(A::BandedBlockBandedMatrix, K::Int, J::Int) = Base.size(arr::BandedBlockBandedMatrix) = @inbounds return map(length,axes(arr)) +@inline function inbands_viewblock(A::BandedBlockBandedMatrix, KJ::Block{2}) + l,u = blockbandwidths(A) + K,J = KJ.n + _BandedMatrix(view(A.data, Block(K-J+u+1, J)), length(axes(A,1)[Block(K)]), subblockbandwidths(A)...) +end + +@inline inbands_viewblock(A::BandedBlockBandedMatrix, K::Block{1}, J::Block{1}) = inbands_viewblock(A, Block(Int(K), Int(J))) + +emptyview(A::BlockedMatrix, J) = view(A.blocks, 1:0, axes(A,2)[J]) +function emptyview(A, J) + V = view(A, Block(1), J) + view(V, 1:0, 1:size(V,2)) +end + +@inline function outbands_viewblock(A::BandedBlockBandedMatrix, KJ::Block{2}) + l,u = blockbandwidths(A) + K,J = KJ.n + # Need to make an empty banded matrix in a type-stable way + _BandedMatrix(emptyview(A.data, Block(J)), length(axes(A,1)[Block(K)]), -720,-720) +end + + + +@inline function viewblock(A::BandedBlockBandedMatrix, KJ::Block{2}) + @boundscheck blockcheckbounds(A, KJ) + K,J = KJ.n + if -A.u ≤ K-J ≤ A.l + inbands_viewblock(A, KJ) + else + outbands_viewblock(A, KJ) + end +end + @inline function getindex(A::BandedBlockBandedMatrix{T}, i::Int, j::Int) where T @boundscheck checkbounds(A, i, j) @@ -346,7 +379,7 @@ end @boundscheck checkbounds(A, i, j) BI,BJ = findblockindex.(axes(A), (i,j)) if -A.l ≤ Int(block(BJ)-block(BI)) ≤ A.u - V = view(A, block(BI),block(BJ)) + V = inbands_viewblock(A, block(BI),block(BJ)) @inbounds V[blockindex(BI),blockindex(BJ)] = convert(T, v)::T elseif !iszero(v) throw(BandError(A)) diff --git a/src/BlockBandedMatrices.jl b/src/BlockBandedMatrices.jl index c9f950e..a05eaf3 100644 --- a/src/BlockBandedMatrices.jl +++ b/src/BlockBandedMatrices.jl @@ -24,7 +24,7 @@ import BlockArrays: AbstractBlockLayout, AbstractBlockedUnitRange, Block, BlockI _blockkron, _blocklengths2blocklasts, block, blockcheckbounds, blockcolstart, blockcolstop, blockcolsupport, blockindex, blockisequal, blockrowstart, blockrowstop, blockrowsupport, blocks, blocksize, checksquareblocks, - hasmatchingblocks, sizes_from_blocks + hasmatchingblocks, sizes_from_blocks, viewblock import FillArrays: Fill, Ones, Zeros