Support non-UnitRange based ranges (#94)

* Support non-UnitRange based ranges

* BlockArrays new view

* blockbandwidths and subblockbandwidths for Zeros

* Add ==, sublayout

* tests pass

* v0.10

* Update Project.toml

* simplify sublayout

* add simple implementation for tests

Author: dlfivefifty

Project.toml

@@ -1,13 +1,14 @@
 name = "BlockBandedMatrices"
 uuid = "ffab5731-97b5-5995-9138-79e8c1846df0"
-version = "0.9.6"
+version = "0.10.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MatrixFactorizations = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
@@ -16,11 +17,11 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-ArrayLayouts = "0.4.1"
-BandedMatrices = "0.15.16"
-BlockArrays = "0.13"
-FillArrays = "0.9.2, 0.10"
-MatrixFactorizations = "0.5.2, 0.6, 0.7, 0.8"
+ArrayLayouts = "0.5"
+BandedMatrices = "0.16"
+BlockArrays = "0.14"
+FillArrays = "0.11"
+MatrixFactorizations = "0.7.1, 0.8"
 julia = "1.5"
 
 [extras]

src/AbstractBlockBandedMatrix.jl

@@ -172,3 +172,15 @@ function fill!(A::AbstractBlockBandedMatrix, val::Any)
   fill!(A.data, val)
   A
 end
+
+
+function ==(A::AbstractBlockBandedMatrix, B::AbstractBlockBandedMatrix)
+    axes(A) == axes(B) || return false
+    blockisequal(axes(A), axes(B)) || return Base.invoke(==, Tuple{AbstractArray,AbstractArray}, A, B)
+    l,u = max.(blockbandwidths(A), blockbandwidths(B))
+    N = blocksize(A,1)
+    for J = blockaxes(A,2), K = max(Block(1),J-u):min(J+l,Block(N))
+        view(A, K, J) == view(B, K, J) || return false
+    end
+    return true 
+end

src/BandedBlockBandedMatrix.jl

@@ -29,7 +29,7 @@ import ArrayLayouts: BlasMatLmulVec, MatLmulVec, MatLmulMat,
                     _qr, _factorize, _copyto!, zero!, layout_replace_in_print_matrix
 
 import BlockArrays: blocksize, blockcheckbounds, BlockedUnitRange, blockisequal, DefaultBlockAxis,
-                        Block, BlockSlice, getblock, unblock, setblock!, block, blockindex,
+                        Block, BlockSlice, unblock, block, blockindex,
                         _blocklengths2blocklasts, BlockIndexRange, sizes_from_blocks, BlockSlice1,
                         blockcolsupport, blockrowsupport, blockcolstart, blockcolstop, blockrowstart, blockrowstop,
                         AbstractBlockLayout, BlockLayout, blocks, hasmatchingblocks, BlockStyle
@@ -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,
-                        BandedLayout, BandedColumnMajor, BandedColumns,
+                        BandedLayout, BandedColumnMajor, BandedColumns, bandedcolumns,
                         BandedSubBandedMatrix, bandeddata,
                         _BandedMatrix, colstart, colstop, rowstart, rowstop,
                         BandedStyle, _fill_lmul!, bandshift,
@@ -49,8 +49,8 @@ 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 BlockRange1{R<:AbstractUnitRange{Int}} = BlockRange{1,Tuple{R}}
+const BlockIndexRange1{R<:AbstractUnitRange{Int}} = BlockIndexRange{1,Tuple{R}}
 
 blockcolrange(A...) = blockcolsupport(A...)
 blockrowrange(A...) = blockrowsupport(A...)

src/BlockBandedMatrices.jl

@@ -300,23 +300,6 @@ BroadcastStyle(::Type{<:BlockBandedMatrix}) = BlockBandedStyle()
 zeroblock(A::BlockSkylineMatrix, K::Int, J::Int) =
     Matrix(Zeros{eltype(A)}(length.(getindex.(axes(A),(Block(K),Block(J))))))
 
-@inline function getblock(A::BlockSkylineMatrix, K::Int, J::Int)
-    @boundscheck blockcheckbounds(A, K, J)
-    if -A.block_sizes.l[J] ≤ J - K ≤ A.block_sizes.u[J]
-        convert(Matrix, view(A, Block(K, J)))
-    else
-        zeroblock(A, K, J)
-    end
-end
-
-# @inline function Base.getindex(block_arr::BlockArray{T,N}, blockindex::BlockIndex{N}) where {T,N}
-#     @boundscheck checkbounds(block_arr.blocks, blockindex.I...)
-#     @inbounds block = block_arr.blocks[blockindex.I...]
-#     @boundscheck checkbounds(block, blockindex.α...)
-#     @inbounds v = block[blockindex.α...]
-#     return v
-# end
-
 
 ###########################
 # AbstractArray Interface #
@@ -366,48 +349,7 @@ end
     return A
 end
 
-# @propagate_inbounds function Base.setindex!(block_array::BlockArray{T, N}, v, block_index::BlockIndex{N}) where {T,N}
-#     getblock(block_array, block_index.I...)[block_index.α...] = v
-# end
 
-########
-# Misc #
-########
-
-# @generated function Base.Array(block_array::BlockArray{T, N, R}) where {T,N,R}
-#     # TODO: This will fail for empty block array
-#     return quote
-#         block_sizes = block_array.block_sizes
-#         arr = similar(block_array.blocks[1], size(block_array)...)
-#         @nloops $N i i->(1:nblocks(block_sizes, i)) begin
-#             block_index = @ntuple $N i
-#             indices = globalrange(block_sizes, block_index)
-#             arr[indices...] = getblock(block_array, block_index...)
-#         end
-#
-#         return arr
-#     end
-# end
-#
-# @generated function Base.copyto!(block_array::BlockArray{T, N, R}, arr::R) where {T,N,R <: AbstractArray}
-#     return quote
-#         block_sizes = block_array.block_sizes
-#
-#         @nloops $N i i->(1:nblocks(block_sizes, i)) begin
-#             block_index = @ntuple $N i
-#             indices = globalrange(block_sizes, block_index)
-#             copyto!(getblock(block_array, block_index...), arr[indices...])
-#         end
-#
-#         return block_array
-#     end
-# end
-#
-# function Base.fill!(block_array::BlockArray, v)
-#     for block in block_array.blocks
-#         fill!(block, v)
-#     end
-# end
 
 
 ##################
@@ -417,7 +359,7 @@ end
 #   with StridedMatrix.
 ##################
 
-const BlockBandedBlock{T} = SubArray{T,2,<:BlockSkylineMatrix,<:Tuple{<:BlockSlice1,<:BlockSlice1},false}
+const BlockBandedBlock{T} = SubArray{T,2,<:BlockSkylineMatrix,<:Tuple{BlockSlice1,BlockSlice1},false}
 
 
 

src/BlockSkylineMatrix.jl

@@ -214,6 +214,15 @@ function _copyto!(_, ::BlockLayout{<:AbstractBandedLayout}, dest::AbstractMatrix
     dest
 end
 
+function _copyto!(::BandedBlockBandedColumns, ::BandedBlockBandedColumns, dest::AbstractMatrix, src::AbstractMatrix)
+    if blockbandwidths(dest) == blockbandwidths(src) && subblockbandwidths(dest) == subblockbandwidths(src)
+        copyto!(bandedblockbandeddata(dest), bandedblockbandeddata(src))
+        dest
+    else # TODO: if subblockbandwidths match we can still do something better
+        _copyto!(BandedBlockBandedLayout(), BandedBlockBandedLayout(), dest, src)
+    end
+end
+
 
 function copyto!(dest::AbstractArray, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(identity)})
     (A,) = bc.args

src/broadcast.jl

@@ -66,14 +66,16 @@ end
 blockbandwidths(A::BlockArray) = bandwidths(A.blocks)
 isblockbanded(A::BlockArray) = isbanded(A.blocks)
 
-@inline function getblock(block_arr::BlockBidiagonal{T,VT}, K::Int, J::Int) where {T,VT<:AbstractMatrix}
+@inline function BlockArrays.viewblock(block_arr::BlockBidiagonal{T,VT}, KJ::Block{2}) where {T,VT<:AbstractMatrix}
+    K,J = KJ.n
     @boundscheck blockcheckbounds(block_arr, K, J)
     l,u = blockbandwidths(block_arr)
     -l ≤ (J-K) ≤ u || return convert(VT, Zeros{T}(length.(getindex.(axes(block_arr),(Block(K),Block(J))))...))
     block_arr.blocks[K,J]
 end
 
-@inline function getblock(block_arr::BlockTridiagonal{T,VT}, K::Int, J::Int) where {T,VT<:AbstractMatrix}
+@inline function BlockArrays.viewblock(block_arr::BlockTridiagonal{T,VT}, KJ::Block{2}) where {T,VT<:AbstractMatrix}
+    K,J = KJ.n
     @boundscheck blockcheckbounds(block_arr, K, J)
     abs(J-K) ≥ 2 && return convert(VT, Zeros{T}(length.(getindex.(axes(block_arr),(Block(K),Block(J))))...))
     block_arr.blocks[K,J]
@@ -114,3 +116,11 @@ for op in (:-, :+)
         end
     end
 end
+
+
+###
+# Zeros
+###
+
+blockbandwidths(::Zeros) = (-1,-1)
+subblockbandwidths(::Zeros) = (-1,-1)

src/interfaceimpl.jl

@@ -1,6 +1,6 @@
 # BlockBandedMatrix with block range indexes is also block-banded
 const SubBlockSkylineMatrix{T,LL,UU,R1,R2} =
-    SubArray{T,2,BlockSkylineMatrix{T,LL,UU},<:Tuple{<:BlockSlice{R1},<:BlockSlice{R2}}}
+    SubArray{T,2,BlockSkylineMatrix{T,LL,UU},<:Tuple{BlockSlice{<:R1},BlockSlice{<:R2}}}
 
 
 
@@ -71,7 +71,7 @@ similar(M::MulAdd{<:AbstractBlockBandedLayout,<:DiagonalLayout}, ::Type{T}) wher
 
 
 
-function blockbandwidths(V::SubBlockSkylineMatrix{<:Any,LL,UU,BlockRange1,BlockRange1}) where {LL,UU}
+function blockbandwidths(V::SubBlockSkylineMatrix{<:Any,LL,UU,<:BlockRange1,<:BlockRange1}) where {LL,UU}
     A = parent(V)
 
     KR = parentindices(V)[1].block.indices[1]
@@ -87,13 +87,14 @@ end
 # BlockIndexRange subblocks
 ####
 
+sublayout(::AbstractBlockBandedLayout, ::Type{<:Tuple{<:BlockSlice{<:BlockRange1}, <:BlockSlice{<:BlockRange1}}}) = BlockBandedLayout()
 
-sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:BlockSlice{BlockRange1}}}) = BlockBandedColumnMajor()
-sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:BlockSlice{Block1}}}) = ColumnMajor()
-sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:BlockSlice{BlockIndexRange1}}}) = ColumnMajor()
-sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockIndexRange1}, <:BlockSlice{BlockIndexRange1}}}) = ColumnMajor()
+sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{<:BlockRange1}, <:BlockSlice{<:BlockRange1}}}) = BlockBandedColumnMajor()
+sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{<:BlockRange1}, <:BlockSlice{Block1}}}) = ColumnMajor()
+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(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{<:BlockRange1}, <:BlockSlice{<:BlockRange1}}}) =
     isblockbanded(parent(V))
 
 sub_materialize(::AbstractBlockBandedLayout, V, _) = BlockBandedMatrix(V)
@@ -113,10 +114,10 @@ function unsafe_convert(::Type{Ptr{T}}, V::SubBlockSkylineMatrix{T,LL,UU,<:Union
     p = unsafe_convert(Ptr{T}, view(A, first(KR), J))
 end
 
-strides(V::SubBlockSkylineMatrix{<:Any,LL,UU,BlockRange1,BlockIndexRange1}) where {LL,UU} =
+strides(V::SubBlockSkylineMatrix{<:Any,LL,UU,<:BlockRange1,<:BlockIndexRange1}) where {LL,UU} =
     (1,parent(V).block_sizes.block_strides[Int(Block(parentindices(V)[2]))])
 
-function unsafe_convert(::Type{Ptr{T}}, V::SubBlockSkylineMatrix{T,LL,UU,BlockRange1,BlockIndexRange1}) where {T,LL,UU}
+function unsafe_convert(::Type{Ptr{T}}, V::SubBlockSkylineMatrix{T,LL,UU,<:BlockRange1,<:BlockIndexRange1}) where {T,LL,UU}
     A = parent(V)
     JR = parentindices(V)[2]
     K = first(parentindices(V)[1].block)

src/linalg.jl

@@ -1,6 +1,7 @@
 using BlockArrays, BandedMatrices, BlockBandedMatrices, FillArrays, SparseArrays, Test, ArrayLayouts , LinearAlgebra
 import BlockBandedMatrices: _BandedBlockBandedMatrix, blockcolsupport, blockrowsupport, colsupport, rowsupport, 
                             isbandedblockbanded, bandeddata, BandedBlockBandedColumns
+import ArrayLayouts: RangeCumsum
 
 @testset "BandedBlockBandedMatrix" begin
     @testset "constructors" begin
@@ -181,7 +182,7 @@ import BlockBandedMatrices: _BandedBlockBandedMatrix, blockcolsupport, blockrows
         A = _BandedBlockBandedMatrix(data, rows,cols, (l,u), (λ,μ))
 
         @test A[Block(1), Block(1)] isa BandedMatrix
-        @test A[Block(1), Block(1)] == A[Block(1,1)] == BlockArrays.getblock(A, 1, 1) == BandedMatrix(view(A, Block(1,1)))
+        @test A[Block(1), Block(1)] == A[Block(1,1)] == view(A, Block(1, 1)) == BandedMatrix(view(A, Block(1,1)))
         @test A[1,1] == view(A,Block(1),Block(1))[1,1] == view(A,Block(1,1))[1,1] == A[Block(1,1)][1,1]  == A[Block(1),Block(1)][1,1] == 5
         @test A[2,1] == view(A,Block(2),Block(1))[1,1] == view(A,Block(2,1))[1,1] == 8
         @test A[3,1] == view(A,Block(2),Block(1))[2,1] == 9
@@ -454,20 +455,20 @@ import BlockBandedMatrices: _BandedBlockBandedMatrix, blockcolsupport, blockrows
         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 MemoryLayout(typeof(V)) == BandedMatrices.BandedColumns{ColumnMajor}()
         @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 MemoryLayout(typeof(V2)) == BandedMatrices.BandedColumns{ColumnMajor}()
         @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 MemoryLayout(typeof(V2)) == BandedMatrices.BandedColumns{ColumnMajor}()
         @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 MemoryLayout(typeof(V2)) == BandedMatrices.BandedColumns{ColumnMajor}()
         @test bandwidths(V2) == (-1,1)
         @test BandedMatrix(V2) == A[3:3,5:6]
     end
@@ -487,6 +488,22 @@ import BlockBandedMatrices: _BandedBlockBandedMatrix, blockcolsupport, blockrows
         A = _BandedBlockBandedMatrix(PseudoBlockVector([1,2,3],[1,2])', blockedrange([1,2]), (-1,1), (-1,1))
         @test MemoryLayout(A) isa BandedBlockBandedColumns{RowMajor}
     end
+
+    @testset "1:N blocks" begin
+        N = 10
+        A = BandedBlockBandedMatrix{Float64}(undef, 1:N,1:N, (1,1), (1,1))
+        @test axes(A) isa NTuple{2,BlockedUnitRange{<:RangeCumsum}}
+    end
+
+    @testset "change bandwidths" begin
+        l , u = 1,1
+        λ , μ = 1,1
+        N = M = 4
+        cols = rows = 1:N
+        data = reshape(collect(1:(λ+μ+1)*(l+u+1)*sum(cols)), ((λ+μ+1)*(l+u+1), sum(cols)))
+        A = _BandedBlockBandedMatrix(data, rows,cols, (l,u), (λ,μ))
+        @test A == BandedBlockBandedMatrix(A, (2,1), (2,1)) == BandedBlockBandedMatrix{Float64}(A, (2,1), (2,1))
+    end
 end
 
 if false # turned off since tests have check-bounds=yes

test/test_bandedblockbanded.jl

@@ -47,7 +47,7 @@ import Base.Broadcast: materialize!
 
 
         @test A[Block(1), Block(1)] isa Matrix
-        @test A[Block(1), Block(1)] == A[Block(1,1)] == BlockArrays.getblock(A, 1, 1) == Matrix(view(A, Block(1,1)))
+        @test A[Block(1), Block(1)] == A[Block(1,1)] == view(A, Block(1, 1)) == Matrix(view(A, Block(1,1)))
         @test A[1,1] == view(A,Block(1),Block(1))[1,1] == view(A,Block(1,1))[1,1] == A[Block(1,1)][1,1]  == A[Block(1),Block(1)][1,1] == 1
         @test A[2,1] == view(A,Block(2),Block(1))[1,1] == view(A,Block(2,1))[1,1] == 2
         @test A[3,1] == view(A,Block(2),Block(1))[2,1] == 3

test/test_blockbanded.jl

@@ -150,14 +150,17 @@ using BandedMatrices, BlockBandedMatrices, BlockArrays, LinearAlgebra, ArrayLayo
         @test B == C
 
         N = 10
-        A = BandedBlockBandedMatrix{Float64}(undef, 1:N,1:N, (1,1), (1,1))
+        A = BandedBlockBandedMatrix{Float64}(undef, Base.OneTo(N),Base.OneTo(N), (1,1), (1,1))
             A.data .= randn.()
-        B = BandedBlockBandedMatrix{Float64}(undef, 1:N,1:N, (2,2), (2,2))
+        B = BandedBlockBandedMatrix{Float64}(undef, Base.OneTo(N),Base.OneTo(N), (2,2), (2,2))
             B.data .= randn.()
-        C = BandedBlockBandedMatrix{Float64}(undef, 1:N,1:N, (3,3), (3,3))
+        C = BandedBlockBandedMatrix{Float64}(undef, Base.OneTo(N),Base.OneTo(N), (3,3), (3,3))
         @time C .= A .+ B
         @test C == A + B == A .+ B
 
+        bc = Base.broadcasted(+, A, B)
+        @test @inferred(axes(bc)) === axes(A)
+
         @test A + B isa typeof(A)
         @test A .+ B isa typeof(A)
         @test blockbandwidths(A+B) == blockbandwidths(A.+B) == (2,2)