Support AbstractQLayout (#70)

* Revert "Revert "abstractq""

This reverts commit 31b1fce.

* Revert "Revert "Merge branch 'master' of https://github.com/JuliaMatrices/BlockBandedMatrices.jl""

This reverts commit 08d6dd4.

* Update BlockBandedMatrices.jl

* Use MatrixFactorizations.QR

* update qr

* Update Project.toml

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "BlockBandedMatrices"
 uuid = "ffab5731-97b5-5995-9138-79e8c1846df0"
-version = "0.8.1"
+version = "0.8.2"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -18,9 +18,9 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-ArrayLayouts = "0.2.3"
-BandedMatrices = "0.15.2"
-BlockArrays = "0.12"
-FillArrays = "0.8"
-MatrixFactorizations = "0.3.1"
+ArrayLayouts = "0.2.6"
+BandedMatrices = "0.15.6"
+BlockArrays = "0.12.4"
+FillArrays = "0.8.8"
+MatrixFactorizations = "0.4"
 julia = "1.2"

src/AbstractBlockBandedMatrix.jl

@@ -1,14 +1,24 @@
-####
-# Matrix memory layout traits
-#
-# if MemoryLayout(A) returns BandedColumnMajor, you must override
-# pointer and leadingdimension
-# in addition to the banded matrix interface
-####
 
+"""
+   AbstractBlockBandedLayout
+
+isa a `MemoryLayout` that indicates that the array implements the block-banded
+interface.
+"""
 abstract type AbstractBlockBandedLayout <: AbstractBlockLayout end
+
+"""
+   AbstractBandedBlockBandedLayout
+
+isa a `MemoryLayout` that indicates that the array implements the banded-block-banded
+interface.
+"""
 abstract type AbstractBandedBlockBandedLayout <: AbstractBlockBandedLayout end
 
+
+struct BandedBlockBandedLayout <: AbstractBandedBlockBandedLayout end
+struct BlockBandedLayout <: AbstractBlockBandedLayout end
+
 struct BandedBlockBandedColumns{LAY} <: AbstractBandedBlockBandedLayout end
 struct BandedBlockBandedRows{LAY} <: AbstractBandedBlockBandedLayout end
 struct BlockBandedColumns{LAY} <: AbstractBlockBandedLayout end
@@ -64,6 +74,12 @@ for Func in (:blockbanded_blockcolstart, :blockbanded_blockcolstop, :blockbanded
     @eval $Func(A, i::Block{1}) = $Func(A, Int(i))
 end
 
+blockbanded_blockcolstart(A, i::BlockRange) = blockbanded_blockcolstart(A, minimum(i))
+blockbanded_blocktowstart(A, i::BlockRange) = blockbanded_blockrowstart(A, minimum(i))
+blockbanded_blockcolstop(A, i::BlockRange) = blockbanded_blockcolstop(A, maximum(i))
+blockbanded_blocktowstop(A, i::BlockRange) = blockbanded_blockrowstop(A, maximum(i))
+
+
 @inline blockcolsupport(::AbstractBlockBandedLayout, A, i) = blockbanded_blockcolstart(A,i):blockbanded_blockcolstop(A,i)
 @inline blockrowsupport(::AbstractBlockBandedLayout, A, i) = blockbanded_blockrowstart(A,i):blockbanded_blockrowstop(A,i)
 

src/BandedBlockBandedMatrix.jl

@@ -175,8 +175,10 @@ function BandedBlockBandedMatrix{T,Blocks,RR}(A::AbstractMatrix, axes::NTuple{2,
 end
 
 
-BandedBlockBandedMatrix(A::AbstractMatrix) =
-    BandedBlockBandedMatrix(A, axes(A), blockbandwidths(A), subblockbandwidths(A))
+BandedBlockBandedMatrix{T}(A::AbstractMatrix) where T =
+    copyto!(BandedBlockBandedMatrix{T}(undef, axes(A), blockbandwidths(A), subblockbandwidths(A)), A)
+
+BandedBlockBandedMatrix(A::AbstractMatrix{T}) where T = BandedBlockBandedMatrix{T}(A)
 
 BandedBlockBandedMatrix(A::AbstractMatrix, rdims::AbstractVector{Int}, cdims::AbstractVector{Int}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) =
     BandedBlockBandedMatrix(A, (blockedrange(rdims), blockedrange(cdims)), lu, λμ)
@@ -316,13 +318,18 @@ end
 end
 
 ## structured matrix methods ##
-function Base.replace_in_print_matrix(A::BandedBlockBandedMatrix, i::Integer, j::Integer, s::AbstractString)
+function _bandedblockbanded_replace_in_print_matrix(A, i, j, s)
     bi = findblockindex.(axes(A), (i,j))
     I,J = block.(bi)
     i,j = blockindex.(bi)
-    -A.l ≤ Int(J-I) ≤ A.u && -A.λ ≤ j-i ≤ A.μ ? s : Base.replace_with_centered_mark(s)
+    l,u = blockbandwidths(A)
+    λ,μ = subblockbandwidths(A)
+    -l ≤ Int(J-I) ≤ u && -λ ≤ j-i ≤ μ ? s : Base.replace_with_centered_mark(s)
 end
 
+Base.replace_in_print_matrix(A::BandedBlockBandedMatrix, i::Integer, j::Integer, s::AbstractString) = 
+    _bandedblockbanded_replace_in_print_matrix(A, i, j, s)
+
 
 ############
 # Indexing #
@@ -411,37 +418,48 @@ sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{Block1},Bloc
 sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{Block1}}}) = BandedBlockBandedColumnMajor()
 sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{BlockIndexRange1}}}) = BandedBlockBandedColumnMajor()
 
+sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{Block1}}}) = BandedLayout()
+sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) = BandedLayout()
+sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) = BandedLayout()
+sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) = BandedLayout()
+sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{BlockRange1}}}) = BandedBlockBandedLayout()
+sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{BlockRange1}}}) = BandedBlockBandedLayout()
+sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{Block1}}}) = BandedBlockBandedLayout()
+sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{BlockIndexRange1}}}) = BandedBlockBandedLayout()
+
+
+sub_materialize(::AbstractBandedBlockBandedLayout, V, _) = BandedBlockBandedMatrix(V)
+
 isbanded(A::SubArray{<:Any,2,<:BandedBlockBandedMatrix}) = MemoryLayout(typeof(A)) == BandedColumnMajor()
 isbandedblockbanded(A::SubArray{<:Any,2,<:BandedBlockBandedMatrix}) = MemoryLayout(typeof(A)) == BandedBlockBandedColumnMajor()
 
 
-subblockbandwidths(V::SubBandedBlockBandedMatrix) = subblockbandwidths(parent(V))
+subblockbandwidths(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{BlockRange1},<:BlockSlice{BlockRange1}}}) = 
+    subblockbandwidths(parent(V))
 
-function blockbandwidths(V::SubBandedBlockBandedMatrix{<:Any,BlockRange1,Block1})
+function blockbandwidths(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{BlockRange1},<:BlockSlice{Block1}}})
     A = parent(V)
-
     KR = parentindices(V)[1].block.indices[1]
     J = parentindices(V)[2].block
     shift = Int(KR[1])-Int(J)
     blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
 end
 
-function blockbandwidths(V::SubBandedBlockBandedMatrix{<:Any,Block1,BlockRange1})
+function blockbandwidths(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{<:Block1},<:BlockSlice{<:BlockRange1}}})
     A = parent(V)
-
     K = parentindices(V)[1].block
     JR = parentindices(V)[2].block.indices[1]
     shift = Int(K)-Int(JR[1])
-
-    blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
+    l,u = blockbandwidths(A)
+    l - shift, u + shift
 end
 
-function blockbandwidths(V::SubBandedBlockBandedMatrix{<:Any,BlockRange1,BlockRange1})
+function blockbandwidths(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{<:BlockRange1},<:BlockSlice{<:BlockRange1}}})
     A = parent(V)
 
     KR = parentindices(V)[1].block.indices[1]
     JR = parentindices(V)[2].block.indices[1]
-    shift = Int(KR[1])-Int(JR[1])
+    shift = Int(first(KR))-Int(first(JR))
 
     blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
 end
@@ -487,7 +505,7 @@ parentblocks2Int(V::BandedBlockBandedBlock)::Tuple{Int,Int} = Int(first(parentin
 ######################################
 # BandedMatrix interface  for Blocks #
 ######################################
-@inline function bandwidths(V::BandedBlockBandedBlock) 
+@inline function bandwidths(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice1,BlockSlice1}}) where T
     inblockbands(V) && return subblockbandwidths(parent(V))
     (-720,-720)
 end

src/BlockBandedMatrices.jl

@@ -14,19 +14,19 @@ import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, DefaultArrayStyle, Br
 
 import LinearAlgebra: UniformScaling, isdiag, rmul!, lmul!, ldiv!, rdiv!,
                         AbstractTriangular, AdjOrTrans, HermOrSym, StructuredMatrixStyle,
-                        qr, qr!, QRPackedQ
+                        qr, qr!
 import LinearAlgebra.BLAS: BlasInt, BlasFloat, @blasfunc, libblas, BlasComplex, BlasReal
 import LinearAlgebra.LAPACK: chktrans, chkdiag, liblapack, chklapackerror, checksquare, chkstride1,
                     chkuplo
-import MatrixFactorizations: ql, ql!, _ql, QLPackedQ
+import MatrixFactorizations: ql, ql!, _ql, QLPackedQ, AdjQRPackedQLayout, AdjQLPackedQLayout, QR, QRPackedQ
 import SparseArrays: sparse
 
 import ArrayLayouts: BlasMatLmulVec,
                     triangularlayout, UpperTriangularLayout, TriangularLayout, MatLdivVec,
-                    triangulardata, sublayout, 
+                    triangulardata, sublayout, sub_materialize,
                     AbstractColumnMajor, DenseColumnMajor, ColumnMajor,
                     DiagonalLayout, MulAdd, mul, colsupport, rowsupport,
-                    _qr, _factorize, _copyto!
+                    _qr, _factorize, _copyto!, zero!
 
 import BlockArrays: blocksize, blockcheckbounds, BlockedUnitRange, blockisequal, DefaultBlockAxis,
                         Block, BlockSlice, getblock, unblock, setblock!, block, blockindex,
@@ -38,7 +38,7 @@ import BandedMatrices: isbanded, bandwidths, bandwidth, banded_getindex, colrang
                         inbands_setindex!, inbands_getindex, banded_setindex!,
                         banded_generic_axpy!,
                         BlasFloat, banded_dense_axpy!, MemoryLayout,
-                        BandedColumnMajor,
+                        BandedLayout, BandedColumnMajor,
                         BandedSubBandedMatrix, bandeddata, tribandeddata,
                         _BandedMatrix, colstart, colstop, rowstart, rowstop,
                         BandedStyle, _fill_lmul!, bandshift,
@@ -50,7 +50,7 @@ export BandedBlockBandedMatrix, BlockBandedMatrix, BlockSkylineMatrix, blockband
 
 const Block1 = Block{1,Int}
 const BlockRange1 = BlockRange{1,Tuple{UnitRange{Int}}}
-const BlockIndexRange1 = BlockIndexRange{1,Tuple{UnitRange{Int}}}    
+const BlockIndexRange1 = BlockIndexRange{1,Tuple{UnitRange{Int}}}
 
 blockcolrange(A...) = blockcolsupport(A...)
 blockrowrange(A...) = blockrowsupport(A...)

src/BlockSkylineMatrix.jl

@@ -65,11 +65,14 @@ BlockSkylineSizes(rows::AbstractVector{Int}, cols::AbstractVector{Int}, l::Abstr
     BlockSkylineSizes((blockedrange(rows),blockedrange(cols)), l, u)
 
 BlockBandedSizes(b_axes::NTuple{2,AbstractUnitRange{Int}}, l::Int, u::Int) =
-    BlockSkylineSizes(b_axes, Fill(l, blocklength(b_axes[2])), Fill(u, blocklength(b_axes[2])))
+    convert(BlockBandedSizes, BlockSkylineSizes(b_axes, Fill(l, blocklength(b_axes[2])), Fill(u, blocklength(b_axes[2]))))
 
 
 BlockBandedSizes(rows::AbstractVector{Int}, cols::AbstractVector{Int}, l::Int, u::Int) =
-    BlockSkylineSizes(rows, cols, Fill(l, length(cols)), Fill(u, length(cols)))
+    convert(BlockBandedSizes, BlockSkylineSizes(rows, cols, Fill(l, length(cols)), Fill(u, length(cols))))
+
+convert(::Type{BlockSkylineSizes{BS,LL,UU,BStarts,BStrides}}, b::BlockSkylineSizes) where {BS,LL,UU,BStarts,BStrides} =
+    BlockSkylineSizes(convert(BS, b.axes), convert(BStarts, b.block_starts), convert(BStrides,b.block_strides), convert(LL,b.l), convert(UU,b.u))
 
 colblockbandwidths(bs::BlockSkylineSizes) = (bs.l, bs.u)
 
@@ -227,7 +230,7 @@ BlockBandedMatrix(A::Union{AbstractMatrix,UniformScaling},
                         rdims::AbstractVector{Int}, cdims::AbstractVector{Int},
                         lu::NTuple{2,Int}) = BlockBandedMatrix{eltype(A)}(A, rdims, cdims, lu)
 
-BlockBandedMatrix(A::BlockBandedMatrix, lu::NTuple{2,Int}) = BlockBandedMatrix(A, BlockBandedSizes(axes(A), lu...))
+BlockBandedMatrix(A::AbstractMatrix, lu::NTuple{2,Int}) = BlockBandedMatrix(A, BlockBandedSizes(axes(A), lu...))
 
 function convert(::Type{BlockSkylineMatrix}, A::AbstractMatrix)
     @assert isblockbanded(A)

src/blockskylineqr.jl

@@ -1,36 +1,47 @@
-qrf!(A,τ) = _qrf!(MemoryLayout(typeof(A)),MemoryLayout(typeof(τ)),A,τ)
+qrf!(A,τ) = _qrf!(MemoryLayout(A),MemoryLayout(τ),A,τ)
 _qrf!(::AbstractColumnMajor,::AbstractStridedLayout,A::AbstractMatrix{T},τ::AbstractVector{T}) where T<:BlasFloat =
     LAPACK.geqrf!(A,τ)
 _apply_qr!(::AbstractColumnMajor, ::AbstractStridedLayout, ::AbstractStridedLayout, A::AbstractMatrix{T}, τ::AbstractVector{T}, B::AbstractVecOrMat{T}) where T<:BlasReal =
     LAPACK.ormqr!('L','T',A,τ,B)
 _apply_qr!(::AbstractColumnMajor, ::AbstractStridedLayout, ::AbstractStridedLayout, A::AbstractMatrix{T}, τ::AbstractVector{T}, B::AbstractVecOrMat{T}) where T<:BlasComplex =
     LAPACK.ormqr!('L','C',A,τ,B)
-apply_qr!(A, τ, B) = _apply_qr!(MemoryLayout(typeof(A)), MemoryLayout(typeof(τ)), MemoryLayout(typeof(B)), A, τ, B)
+apply_qr!(A, τ, B) = _apply_qr!(MemoryLayout(A), MemoryLayout(τ), MemoryLayout(B), A, τ, B)
 
-qlf!(A,τ) = _qlf!(MemoryLayout(typeof(A)),MemoryLayout(typeof(τ)),A,τ)
+qlf!(A,τ) = _qlf!(MemoryLayout(A),MemoryLayout(τ),A,τ)
 _qlf!(::AbstractColumnMajor,::AbstractStridedLayout,A::AbstractMatrix{T},τ::AbstractVector{T}) where T<:BlasFloat =
     LAPACK.geqlf!(A,τ)
 _apply_ql!(::AbstractColumnMajor, ::AbstractStridedLayout, ::AbstractStridedLayout, A::AbstractMatrix{T}, τ::AbstractVector{T}, B::AbstractVecOrMat{T}) where T<:BlasReal =
     LAPACK.ormql!('L','T',A,τ,B)
 _apply_ql!(::AbstractColumnMajor, ::AbstractStridedLayout, ::AbstractStridedLayout, A::AbstractMatrix{T}, τ::AbstractVector{T}, B::AbstractVecOrMat{T}) where T<:BlasComplex =
     LAPACK.ormql!('L','C',A,τ,B)
-apply_ql!(A, τ, B) = _apply_ql!(MemoryLayout(typeof(A)), MemoryLayout(typeof(τ)), MemoryLayout(typeof(B)), A, τ, B)
+apply_ql!(A, τ, B) = _apply_ql!(MemoryLayout(A), MemoryLayout(τ), MemoryLayout(B), A, τ, B)
 
-function qr!(A::BlockBandedMatrix{T}) where T
+function _blockbanded_qr!(A::AbstractMatrix, τ::AbstractVector)
+    M,N = Block.(blocksize(A))
+    (M < N ? axes(A,1) : axes(A,2)) == axes(τ,1) || throw(DimensionMismatch(""))
+    _blockbanded_qr!(A, τ, min(N,M))
+    QR(A,τ.blocks)
+end
+
+function _blockbanded_qr!(A::AbstractMatrix, τ::AbstractVector, NCOLS::Block{1})
     l,u = blockbandwidths(A)
-    M,N = blocksize(A)
-    ax1 = M < N ? axes(A,1) : axes(A,2)
-    τ = PseudoBlockVector{T}(undef, (ax1,))
-    for K = 1:min(N,M)
-        KR = Block.(K:min(K+l,M))
-        V = view(A,KR,Block(K))
-        t = view(τ,Block(K))
+    M,N = Block.(blocksize(A))
+    for K = Block(1):NCOLS
+        KR = K:min(K+l,M)
+        V = view(A,KR,K)
+        t = view(τ,K)
         qrf!(V,t)
         for J = K+1:min(K+u,N)
-            apply_qr!(V, t, view(A,KR,Block(J)))
+            apply_qr!(V, t, view(A,KR,J))
         end
     end
-    QR(A,τ.blocks)
+    A,τ
+end
+
+function qr!(A::BlockBandedMatrix{T}) where T 
+    M,N = blocksize(A)
+    ax1 = M < N ? axes(A,1) : axes(A,2)
+    _blockbanded_qr!(A, PseudoBlockVector{T}(undef, (ax1,)))
 end
 
 function ql!(A::BlockBandedMatrix{T}) where T
@@ -57,11 +68,13 @@ function ql!(A::BlockBandedMatrix{T}) where T
     QL(A,τ.blocks)
 end
 
-_qr(::AbstractBlockBandedLayout, _, A) = qr!(BlockBandedMatrix(A, (blockbandwidth(A,1), blockbandwidth(A,1)+blockbandwidth(A,2))))
+_qr(::AbstractBlockBandedLayout, ::Tuple{Integer,Integer}, A) = qr!(BlockBandedMatrix(A, (blockbandwidth(A,1), blockbandwidth(A,1)+blockbandwidth(A,2))))
+_qr(lay::AbstractBlockBandedLayout, ax::Tuple{AbstractUnitRange{Int},AbstractUnitRange{Int}}, A) = _qr(lay, map(length, ax), A)
 _ql(::AbstractBlockBandedLayout, _, A) = ql!(BlockBandedMatrix(A, (blockbandwidth(A,1)+blockbandwidth(A,2),blockbandwidth(A,2))))
 _factorize(::AbstractBlockBandedLayout, _, A) = qr(A)
 
-function lmul!(adjQ::Adjoint{<:Any,<:QRPackedQ{<:Any,<:BlockSkylineMatrix}}, Bin::AbstractVector)
+function materialize!(Mul::Lmul{<:AdjQRPackedQLayout{<:AbstractBlockBandedLayout}})
+    adjQ,Bin = Mul.A,Mul.B
     Q = parent(adjQ)
     A = Q.factors
     l,u = blockbandwidths(A)
@@ -78,23 +91,22 @@ function lmul!(adjQ::Adjoint{<:Any,<:QRPackedQ{<:Any,<:BlockSkylineMatrix}}, Bin
     end
     Bin
 end
-
-function lmul!(adjQ::Adjoint{<:Any,<:QLPackedQ{<:Any,<:BlockSkylineMatrix}}, Bin::AbstractVector)
-    Q = parent(adjQ)
+function materialize!(Mul::Lmul{<:AdjQLPackedQLayout{<:AbstractBlockBandedLayout}})
+    Q = parent(Mul.A)
     A = Q.factors
     l,u = blockbandwidths(A)
     N,M = blocksize(A)
     # impose block structure
     ax1 = (axes(A,1),)
     τ  = PseudoBlockArray(Q.τ, ax1)
-    B = PseudoBlockArray(Bin, ax1)
+    B = PseudoBlockArray(Mul.B, ax1)
     for K = N:-1:1
         KR = Block.(max(1,K-u):K)
         V = view(A,KR,Block(K))
         t = view(τ,Block(K))
         apply_ql!(V, t, view(B,KR))
     end
-    Bin
+    Mul.B
 end
 
 # avoid LinearALgebra Strided obsession 

src/broadcast.jl

@@ -46,7 +46,7 @@ similar(bc::Broadcasted{<:AbstractBlockBandedStyle}, ::Type{T}) where T =
 
 checkblocks(A, B) = blockisequal(axes(A), axes(B)) || throw(DimensionMismatch("*"))
 
-function _blockbanded_copyto!(dest::AbstractMatrix{T}, src::AbstractMatrix) where T
+function _blockbanded_copyto!(dest::AbstractMatrix, src::AbstractMatrix)
     @boundscheck checkblocks(dest, src)
 
     dl, du = colblockbandwidths(dest)
@@ -57,13 +57,13 @@ function _blockbanded_copyto!(dest::AbstractMatrix{T}, src::AbstractMatrix) wher
 
     for J = 1:N
         for K = max(1,J-du[J]):min(J-su[J]-1,M)
-            view(dest,Block(K),Block(J)) .= zero(T)
+            zero!(view(dest,Block(K),Block(J)))
         end
         for K = max(1,J-su[J]):min(J+sl[J],M)
-            view(dest,Block(K),Block(J)) .= view(src,Block(K),Block(J))
+            copyto!(view(dest,Block(K),Block(J)), view(src,Block(K),Block(J)))
         end
         for K = max(1,J+sl[J]+1):min(J+dl[J],M)
-            view(dest,Block(K),Block(J)) .= zero(T)
+            zero!(view(dest,Block(K),Block(J)))
         end
     end
     dest

src/linalg.jl

@@ -77,8 +77,9 @@ function blockbandwidths(V::SubBlockSkylineMatrix{<:Any,LL,UU,BlockRange1,BlockR
     KR = parentindices(V)[1].block.indices[1]
     JR = parentindices(V)[2].block.indices[1]
 
-    @assert KR[1] == JR[1] == 1
-    blockbandwidths(A)
+    shift =  first(KR) - first(JR)
+    l,u = blockbandwidths(A)
+    l-shift, u+shift
 end
 
 
@@ -92,6 +93,10 @@ sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:
 sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:BlockSlice{BlockIndexRange1}}}) = ColumnMajor()
 sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockIndexRange1}, <:BlockSlice{BlockIndexRange1}}}) = ColumnMajor()
 
+isblockbanded(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{BlockRange1}, <:BlockSlice{BlockRange1}}}) =
+    isblockbanded(parent(V))
+
+sub_materialize(::AbstractBlockBandedLayout, V, _) = BlockBandedMatrix(V)
 
 strides(V::SubBlockSkylineMatrix{<:Any,LL,UU,<:Union{BlockRange1,Block1},Block1}) where {LL,UU} =
     (1,parent(V).block_sizes.block_strides[Int(parentindices(V)[2].block)])
@@ -133,5 +138,3 @@ strides(V::SubBlockSkylineMatrix{T,LL,UU,BlockIndexRange1,BlockIndexRange1}) whe
     (1,parent(V).block_sizes.block_strides[Int(parentindices(V)[2].block.block)])
 
 MemoryLayout(V::SubBlockSkylineMatrix{T,LL,UU,BlockIndexRange1,BlockIndexRange1}) where {T,LL,UU} = ColumnMajor()
-
-