Fixes for new BlockArrays

Project.toml

@@ -1,6 +1,6 @@
 name = "BlockBandedMatrices"
 uuid = "ffab5731-97b5-5995-9138-79e8c1846df0"
-version = "0.7"
+version = "0.7.1"
 
 [deps]
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"

src/AbstractBlockBandedMatrix.jl

@@ -117,7 +117,7 @@ const BlockSlice1 = BlockSlice{Block{1,Int}}
 ######################################
 
 @inline colstart(A::AbstractBlockBandedMatrix, i::Integer) =
-    first(axes(A,1)[blockcolstart(A,findblock(axes(A,2),i))])
+    isempty(axes(A,1)) ? 1 : first(axes(A,1)[blockcolstart(A,findblock(axes(A,2),i))])
 
 @inline function colstop(A::AbstractBlockBandedMatrix, i::Integer)
     CS = blockcolstop(A,findblock(axes(A,2),i))

src/BandedBlockBandedMatrix.jl

@@ -85,6 +85,9 @@ function convert(::Type{<:BandedBlockBandedMatrix}, B::BandedMatrix)
     end
 end
 
+convert(::Type{BandedBlockBandedMatrix{T,BLOCKS,RAXIS}}, A::BandedBlockBandedMatrix) where {T,BLOCKS,RAXIS} =
+    _BandedBlockBandedMatrix(convert(BLOCKS, A.data), convert(RAXIS, A.raxis), (A.l, A.u), (A.λ, A.μ))
+
 function BandedBlockBandedMatrix{T,B,R}(Z::Zeros, axes::NTuple{2,AbstractUnitRange{Int}},
                                        lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B,R<:AbstractUnitRange{Int}}
    if size(Z) ≠ map(length,axes)
@@ -191,10 +194,7 @@ similar(A::BandedBlockBandedMatrix, ::Type{T}, axes::NTuple{2,AbstractUnitRange{
 
 
 similar(A::BandedBlockBandedMatrix{T}, axes::NTuple{2,AbstractUnitRange{Int}}) where T =
-      BandedBlockBandedMatrix{T}(undef, axes, blockbandwidths(A), subblockbandwidths(A))   
-
-similar(A::BandedBlockBandedMatrix{T}, axes::NTuple{2,OneTo{Int}}) where T =
-      Matrix{T}(undef, map(length,axes))
+    similar(Matrix{T}, map(length,axes)...)
 
 axes(A::BandedBlockBandedMatrix) = (A.raxis, axes(A.data,2))
 
@@ -398,21 +398,91 @@ end
 #   with BLASBandedMatrix.
 ##################
 
+const SubBandedBlockBandedMatrix{T,R1,R2} =
+    SubArray{T,2,<:BandedBlockBandedMatrix{T},<:Tuple{<:BlockSlice{R1},<:BlockSlice{R2}}}
+
+
+sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{Block1}}}) = BandedColumnMajor()
+sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) = BandedColumnMajor()
+sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) = BandedColumnMajor()
+sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) = BandedColumnMajor()
+sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{BlockRange1}}}) = BandedBlockBandedColumnMajor()
+sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{BlockRange1}}}) = BandedBlockBandedColumnMajor()
+sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{Block1}}}) = BandedBlockBandedColumnMajor()
+sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{BlockIndexRange1}}}) = BandedBlockBandedColumnMajor()
+
+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))
+
+function blockbandwidths(V::SubBandedBlockBandedMatrix{<:Any,BlockRange1,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})
+    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
+end
+
+function blockbandwidths(V::SubBandedBlockBandedMatrix{<:Any,BlockRange1,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])
+
+    blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
+end
+
+
 const BandedBlockBandedBlock{T, BLOCKS, RAXIS} = SubArray{T,2,BandedBlockBandedMatrix{T, BLOCKS, RAXIS},<:Tuple{<:BlockSlice1,<:BlockSlice1},false}
 
 
-isbanded(::BandedBlockBandedBlock) = true
-MemoryLayout(::BandedBlockBandedBlock) = BandedColumnMajor()
 BroadcastStyle(::Type{<: BandedBlockBandedBlock}) = BandedStyle()
 
 
-function inblockbands(V::BandedBlockBandedBlock)
+function inblockbands(V::SubArray{<:Any,2,<:AbstractMatrix,<:Tuple{<:BlockSlice1,<:BlockSlice1},false})
     A = parent(V)
     K_sl, J_sl = parentindices(V)
     K, J = K_sl.block, J_sl.block
-    -A.l ≤ Int(J-K) ≤ A.u
+    l,u = blockbandwidths(A)
+    -l ≤ Int(J-K) ≤ u
 end
 
+function parentblock(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) where T
+    A = parent(V)
+    K_sl, J_sl = parentindices(V)
+    view(A, K_sl.block.block, J_sl.block.block)
+end
+
+function parentblock(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) where T
+    A = parent(V)
+    K_sl, J_sl = parentindices(V)
+    view(A, K_sl.block, J_sl.block.block)
+end
+
+function parentblock(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) where T
+    A = parent(V)
+    K_sl, J_sl = parentindices(V)
+    view(A, K_sl.block.block, J_sl.block)
+end
+
+# gives the columns of parent(V).data that encode the block
+parentblocks2Int(V::BandedBlockBandedBlock)::Tuple{Int,Int} = Int(first(parentindices(V)).block),
+                                                          Int(last(parentindices(V)).block)
+
 
 ######################################
 # BandedMatrix interface  for Blocks #
@@ -422,11 +492,23 @@ end
     (-720,-720)
 end
 
+function bandwidths(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) where T
+   B = parentblock(V)
+   K_sl, J_sl = parentindices(V)
+   bandwidths(B) .+ (-1,1) .* bandshift(K_sl.block.indices[1],J_sl.block.indices[1])
+end
 
+function bandwidths(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) where T
+    B = parentblock(V)
+    K_sl, J_sl = parentindices(V)
+    bandwidths(B) .+ (-1,1) .* bandshift(Base.OneTo(1),J_sl.block.indices[1])
+ end
 
-# gives the columns of parent(V).data that encode the block
-blocks(V::BandedBlockBandedBlock)::Tuple{Int,Int} = Int(first(parentindices(V)).block),
-                                                    Int(last(parentindices(V)).block)
+function bandwidths(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) where T
+    B = parentblock(V)
+    K_sl, J_sl = parentindices(V)
+    bandwidths(B) .+ (-1,1) .* bandshift(K_sl.block.indices[1],Base.OneTo(1))
+ end
 
 
 function bandeddata(V::BandedBlockBandedBlock{T}) where T
@@ -438,6 +520,27 @@ function bandeddata(V::BandedBlockBandedBlock{T}) where T
     view(A.data, u + K - J + 1, J)
 end
 
+function bandeddata(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) where T
+    A = parent(V)
+    K_sl, J_sl = parentindices(V)
+    view(bandeddata(parentblock(V)), :, J_sl.block.indices[1])
+end
+
+function bandeddata(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) where T
+    A = parent(V)
+    K_sl, J_sl = parentindices(V)
+    view(bandeddata(parentblock(V)), :, J_sl.block.indices[1])
+end
+
+function bandeddata(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) where T
+    A = parent(V)
+    K_sl, J_sl = parentindices(V)
+    bandeddata(parentblock(V))
+end
+
+
+
+
 
 @inline function inbands_getindex(V::BandedBlockBandedBlock, k::Int, j::Int)
     A = parent(V)
@@ -453,7 +556,7 @@ end
 @propagate_inbounds function getindex(V::BandedBlockBandedBlock, k::Int, j::Int)
     @boundscheck checkbounds(V, k, j)
     A = parent(V)
-    K,J = blocks(V)
+    K,J = parentblocks2Int(V)
     if -A.l ≤ J-K ≤ A.u && -A.λ ≤ j-k ≤ A.μ
         inbands_getindex(V, k, j)
     else
@@ -464,7 +567,7 @@ end
 @propagate_inbounds function setindex!(V::BandedBlockBandedBlock, v, k::Int, j::Int)
     @boundscheck checkbounds(V, k, j)
     A = parent(V)
-    K,J = blocks(V)
+    K,J = parentblocks2Int(V)
     if -A.l ≤ J-K ≤ A.u && -A.λ ≤ j-k ≤ A.μ
         inbands_setindex!(V, v, k, j)
     elseif iszero(v) # allow setindex for 0 datya

src/BlockBandedMatrices.jl

@@ -38,13 +38,17 @@ import BandedMatrices: isbanded, bandwidths, bandwidth, banded_getindex, colrang
                         BandedColumnMajor,
                         BandedSubBandedMatrix, bandeddata, tribandeddata,
                         _BandedMatrix, colstart, colstop, rowstart, rowstop,
-                        BandedStyle, _fill_lmul!,
+                        BandedStyle, _fill_lmul!, bandshift,
                         _banded_colval, _banded_rowval, _banded_nzval # for sparse
 
 export BandedBlockBandedMatrix, BlockBandedMatrix, BlockSkylineMatrix, blockbandwidth, blockbandwidths,
         subblockbandwidth, subblockbandwidths, Ones, Zeros, Fill, Block, BlockTridiagonal, isblockbanded
 
 
+const Block1 = Block{1,Int}
+const BlockRange1 = BlockRange{1,Tuple{UnitRange{Int}}}
+const BlockIndexRange1 = BlockIndexRange{1,Tuple{UnitRange{Int}}}        
+
 include("AbstractBlockBandedMatrix.jl")
 include("broadcast.jl")
 include("BlockSkylineMatrix.jl")

src/broadcast.jl

@@ -212,7 +212,7 @@ function blockbanded_axpy!(a, X::AbstractMatrix, Y::AbstractMatrix)
 end
 
 function _combine_blockaxes(A, B)
-    blockisequal(axes(A), axes(B)) || throw(DimensionMismatch())
+    blockisequal(axes(A), axes(B)) || throw(DimensionMismatch("Block sizes do not agree"))
     axes(A)
 end
 

src/linalg.jl

@@ -1,14 +1,7 @@
 # BlockBandedMatrix with block range indexes is also block-banded
-const Block1 = Block{1,Int}
-const BlockRange1 = BlockRange{1,Tuple{UnitRange{Int}}}
-const BlockIndexRange1 = BlockIndexRange{1,Tuple{UnitRange{Int}}}
 const SubBlockSkylineMatrix{T,LL,UU,R1,R2} =
     SubArray{T,2,BlockSkylineMatrix{T,LL,UU},<:Tuple{<:BlockSlice{R1},<:BlockSlice{R2}}}
 
-const SubBandedBlockBandedMatrix{T,R1,R2} =
-    SubArray{T,2,<:BandedBlockBandedMatrix{T},<:Tuple{<:BlockSlice{R1},<:BlockSlice{R2}}}
-
-
 
 getindex(A::BandedBlockBandedMatrix, KR::BlockRange1, JR::BlockRange1) = BandedBlockBandedMatrix(view(A, KR, JR))
 getindex(A::BandedBlockBandedMatrix, KR::BlockRange1, J::Block1) = BandedBlockBandedMatrix(view(A, KR, J))
@@ -180,48 +173,6 @@ sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:
 sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:BlockSlice{BlockIndexRange1}}}) = ColumnMajor()
 sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockIndexRange1}, <:BlockSlice{BlockIndexRange1}}}) = ColumnMajor()
 
-sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{Block1}, <:BlockSlice{Block1}}}) = BandedColumnMajor()
-sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:BlockSlice{BlockRange1}}}) = BandedBlockBandedColumnMajor()
-sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{Block1}, <:BlockSlice{BlockRange1}}}) = BandedBlockBandedColumnMajor()
-sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:BlockSlice{Block1}}}) = BandedBlockBandedColumnMajor()
-sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:BlockSlice{BlockIndexRange1}}}) = BandedBlockBandedColumnMajor()
-sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockIndexRange1}, <:BlockSlice{BlockIndexRange1}}}) = BandedBlockBandedColumnMajor()
-
-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))
-
-function blockbandwidths(V::SubBandedBlockBandedMatrix{<:Any,BlockRange1,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})
-    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
-end
-
-function blockbandwidths(V::SubBandedBlockBandedMatrix{<:Any,BlockRange1,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])
-
-    blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
-end
-
 
 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)])

test/test_bandedblockbanded.jl

@@ -444,6 +444,30 @@ import BlockBandedMatrices: _BandedBlockBandedMatrix, blockcolrange, blockrowran
         A = BandedBlockBandedMatrix{Float64}(undef, 1:5,1:5, (-1,1), (-1,1))
         @test size(sparse(A)) == size(A) == (15,15)
     end
+
+    @testset "non-standard blocks" begin
+        A = BandedBlockBandedMatrix{Float64}(undef, Int[], 1:5,(-1,1), (-1,1))
+        @test BlockBandedMatrices.colstart(A,1) == 1
+        A = BandedBlockBandedMatrix{Float64}(undef, 1:2, 1:5,(-1,1), (-1,1))
+        A.data .= randn.()
+        V = view(A, Block(2,3))
+        @test MemoryLayout(typeof(V)) == BandedMatrices.BandedColumnMajor()
+        @test isbanded(V)
+        @test bandwidths(V) == (-1,1)
+        @test BandedMatrix(V) == A[2:3,4:6]
+        V2 = view(V, :, 2:3)
+        @test MemoryLayout(typeof(V2)) == BandedMatrices.BandedColumnMajor()
+        @test bandwidths(V2) == (0,0)
+        @test BandedMatrix(V2) == A[2:3,5:6]
+        V2 = view(V, 2:2, :)
+        @test MemoryLayout(typeof(V2)) == BandedMatrices.BandedColumnMajor()
+        @test bandwidths(V2) == (-2,2)
+        @test BandedMatrix(V2) == A[3:3,4:6]
+        V2 = view(V, 2:2, 2:3)
+        @test MemoryLayout(typeof(V2)) == BandedMatrices.BandedColumnMajor()
+        @test bandwidths(V2) == (-1,1)
+        @test BandedMatrix(V2) == A[3:3,5:6]
+    end
 end
 
 if false # turned off since tests have check-bounds=yes

test/test_linalg.jl

@@ -136,6 +136,19 @@ end
 
     @time BLAS.axpy!(1.0, A, B)
     @test B ≈ AB
+
+    @testset "degenerate" begin
+        A = BandedBlockBandedMatrix{Float64}(undef, 1:5, 1:5,(-1,1), (-1,1))
+        B = BandedBlockBandedMatrix{Float64}(undef, 1:5, 1:5,(1,-1), (1,-1))
+        A.data .= randn.()
+        B.data .= randn.()
+        @test blockbandwidths(A^2) == subblockbandwidths(A^2) == (-2,2)
+        @test A^2 == Matrix(A)^2
+        @test blockbandwidths(B^2) == subblockbandwidths(B^2) == (2,-2)
+        @test B^2 == Matrix(B)^2        
+        @test blockbandwidths(A*B) == subblockbandwidths(A*B) == (0,0)
+        @test A*B == Matrix(A)*Matrix(B)
+    end
 end
 
 @testset "Rectangular block *" begin