Support banded matrix with block structure (#109)

* Support banded matrix with block structure

* DiagonalBlock subarrays are BandedBlockBandedMatrix

* Update Project.toml

* v0.10.4

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "BlockBandedMatrices"
 uuid = "ffab5731-97b5-5995-9138-79e8c1846df0"
-version = "0.10.3"
+version = "0.10.4"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -13,9 +13,9 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-ArrayLayouts = "0.5, 0.6"
-BandedMatrices = "0.16"
-BlockArrays = "0.14, 0.15"
+ArrayLayouts = "0.6.3"
+BandedMatrices = "0.16.8"
+BlockArrays = "0.15"
 FillArrays = "0.11"
 MatrixFactorizations = "0.7.1, 0.8"
 julia = "1.5"

src/BandedBlockBandedMatrix.jl

@@ -294,8 +294,14 @@ bandwidths given by `subblockbandwidths(A)`.
 """
 subblockbandwidths(A::BandedBlockBandedMatrix) = (A.λ, A.μ)
 
-# default is to use bandwidth
-subblockbandwidths(A::AbstractMatrix) = bandwidths(A)
+# default is to use whole block
+_subblockbandwidths(A::AbstractMatrix, ::NTuple{2,OneTo{Int}}) = bandwidths(A)
+function _subblockbandwidths(A::AbstractMatrix, _)
+    M,N = map(maximum, blocksizes(A))
+    M-1,N-1
+end
+
+subblockbandwidths(A::AbstractMatrix) = _subblockbandwidths(A, axes(A))
 
 """
     subblockbandwidth(A, i)
@@ -438,6 +444,9 @@ sub_materialize(::AbstractBandedBlockBandedLayout, V, ::Tuple{<:AbstractUnitRang
 sub_materialize(::AbstractBandedBlockBandedLayout, V, ::Tuple{<:BlockedUnitRange,<:AbstractUnitRange}) = PseudoBlockArray(V)
 
 
+sub_materialize(::AbstractBandedLayout, V, ::Tuple{<:BlockedUnitRange,<:BlockedUnitRange}) = BandedMatrix(V)
+
+
 isbanded(A::SubArray{<:Any,2,<:BandedBlockBandedMatrix}) = MemoryLayout(A) isa AbstractBandedLayout
 isbandedblockbanded(A::SubArray{<:Any,2,<:BandedBlockBandedMatrix}) = MemoryLayout(A) isa AbstractBandedBlockBandedLayout
 

src/interfaceimpl.jl

@@ -12,6 +12,7 @@ function blockbandwidths(P::PseudoBlockMatrix{<:Any,<:Diagonal})
     (0,0)
 end
 
+blockbandwidths(::Diagonal) = (0,0)
 bandeddata(P::PseudoBlockMatrix) = bandeddata(P.blocks)
 bandwidths(P::PseudoBlockMatrix) = bandwidths(P.blocks)
 
@@ -123,4 +124,13 @@ end
 ###
 
 blockbandwidths(::Zeros) = (-1,-1)
-subblockbandwidths(::Zeros) = (-1,-1)
+subblockbandwidths(::Zeros) = (-1,-1)
+
+
+###
+# DiagonalBlockMatrix
+###
+
+sublayout(::DiagonalLayout{L}, inds::Type{<:NTuple{2,BS}}) where {L,BS<:BlockSlice{<:BlockRange1}} = bandedblockbandedcolumns(sublayout(L(),Tuple{BS}))
+subblockbandwidths(::Diagonal) = (0,0)
+bandedblockbandeddata(D::Diagonal) = permutedims(D.diag)

test/test_misc.jl

@@ -105,4 +105,23 @@ end
         @test A-I == mortar(Tridiagonal(fill([1 2; 1 2],3), fill([2 4; 3 3],4), fill([4 5; 4 5],3))) == Matrix(A) - I
         @test I-A == mortar(Tridiagonal(fill(-[1 2; 1 2],3), fill([-2 -4; -3 -3],4), fill(-[4 5; 4 5],3))) == I - Matrix(A)
     end
+
+    @testset "DiagonalBlock" begin
+        D = Diagonal(PseudoBlockVector(randn(6), 1:3))
+        D̃ = BandedBlockBandedMatrix(D)
+        @test D̃ == D
+        @test blockbandwidths(D) == blockbandwidths(D̃) == subblockbandwidths(D) == subblockbandwidths(D̃) == (0,0)
+
+        @test @inferred(D[Block(2,2)]) isa BandedMatrix
+        @test D[Block(2,3)] isa BandedMatrix
+        @test bandwidths(D[Block(2,2)]) == (0,0)
+        @test bandwidths(D[Block(2,3)]) == (-720,-720)
+        @test D[Block(2,2)] == D[2:3,2:3]
+        @test D[Block(2,3)] == zeros(2,3)
+
+        @test D[Block.(1:2),Block.(2:3)] isa BandedBlockBandedMatrix
+        @test D[Block.(1:2),Block.(2:3)] == D[1:3,2:6]
+        @test blockbandwidths(D[Block.(1:2),Block.(2:3)]) == (1,-1)
+        @test subblockbandwidths(D[Block.(1:2),Block.(2:3)]) == (0,0)
+    end
 end