Return a BandedMatrix from a view of a BandedBlockBandedMatrix

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"

src/BandedBlockBandedMatrix.jl

@@ -39,7 +39,7 @@
     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}(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 @@
 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 @@
     @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))

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