conj and trans layouts are bandedblockbanded (#118)

* conj and trans layouts are bandedblockbanded

* Start blockbandwidths of banded matrix

* shifts properly

* Support Slice in sublayout

* tests pass

* Drop Julia v1.5

* add tests

* increase coverage

Author: dlfivefifty

.github/workflows/ci.yml

@@ -10,8 +10,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.5'
-          - '1'
+          - '1.6'
           - '^1.7.0-0'
         os:
           - ubuntu-latest

Project.toml

@@ -1,6 +1,6 @@
 name = "BlockBandedMatrices"
 uuid = "ffab5731-97b5-5995-9138-79e8c1846df0"
-version = "0.10.9"
+version = "0.11"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -13,12 +13,12 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-ArrayLayouts = "0.6.3, 0.7"
+ArrayLayouts = "0.7"
 BandedMatrices = "0.16.8"
-BlockArrays = "0.15.2, 0.16"
+BlockArrays = "0.16.6"
 FillArrays = "0.11, 0.12"
 MatrixFactorizations = "0.7.1, 0.8"
-julia = "1.5"
+julia = "1.6"
 
 [extras]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

src/AbstractBlockBandedMatrix.jl

@@ -25,29 +25,35 @@ struct BlockBandedColumns{LAY} <: AbstractBlockBandedLayout end
 struct BlockBandedRows{LAY} <: AbstractBlockBandedLayout end
 
 const BandedBlockBandedColumnMajor = BandedBlockBandedColumns{ColumnMajor}
-const BandedBlockBandedRowMajor = BandedBlockBandedColumns{RowMajor}
+const BandedBlockBandedRowMajor = BandedBlockBandedRows{ColumnMajor}
 const BlockBandedColumnMajor = BlockBandedColumns{ColumnMajor}
-const BlockBandedRowMajor = BlockBandedColumns{RowMajor}
+const BlockBandedRowMajor = BlockBandedRows{ColumnMajor}
 
-transposelayout(::BandedBlockBandedColumnMajor) = BandedBlockBandedRowMajor()
-transposelayout(::BandedBlockBandedRowMajor) = BandedBlockBandedColumnMajor()
-transposelayout(::BlockBandedColumnMajor) = BlockBandedRowMajor()
-transposelayout(::BlockBandedRowMajor) = BlockBandedColumnMajor()
+transposelayout(::BandedBlockBandedColumns{Lay}) where Lay = BandedBlockBandedRows{Lay}()
+transposelayout(::BandedBlockBandedRows{Lay}) where Lay = BandedBlockBandedColumns{Lay}()
+transposelayout(::BlockBandedColumns{Lay}) where Lay = BlockBandedRows{Lay}()
+transposelayout(::BlockBandedRows{Lay}) where Lay = BlockBandedColumns{Lay}()
 
+conjlayout(::Type{T}, ::BandedBlockBandedColumns{Lay}) where {T<:Complex,Lay} = BandedBlockBandedColumns{typeof(conjlayout(T,Lay))}()
+conjlayout(::Type{T}, ::BandedBlockBandedRows{Lay}) where {T<:Complex,Lay} = BandedBlockBandedRows{typeof(conjlayout(T,Lay))}()
+conjlayout(::Type{T}, ::BlockBandedColumns{Lay}) where {T<:Complex,Lay} = BlockBandedColumns{typeof(conjlayout(T,Lay))}()
+conjlayout(::Type{T}, ::BlockBandedRows{Lay}) where {T<:Complex,Lay} = BlockBandedRows{typeof(conjlayout(T,Lay))}()
 
 
 # AbstractBandedMatrix must implement
 
 # A BlockBandedMatrix is a BlockMatrix, but is not a BandedMatrix
 abstract type AbstractBlockBandedMatrix{T} <: AbstractBlockMatrix{T} end
+MemoryLayout(::Type{<:AbstractBlockBandedMatrix}) = BlockBandedLayout()
 
 
 """
     blockbandwidths(A)
 
 Returns a tuple containing the upper and lower blockbandwidth of `A`.
 """
-blockbandwidths(A::AbstractVecOrMat) = (blocksize(A,1)-1 , blocksize(A,2)-1)
+blockbandwidths(A::AbstractVecOrMat) = blockbandwidths(MemoryLayout(A), axes(A), A)
+blockbandwidths(_, _, A) = (blocksize(A,1)-1 , blocksize(A,2)-1)
 
 """
     blockbandwidth(A,i)

src/BlockBandedMatrices.jl

@@ -26,13 +26,14 @@ import ArrayLayouts: BlasMatLmulVec, MatLmulVec, MatLmulMat,
                     triangulardata, sublayout, sub_materialize,
                     AbstractColumnMajor, DenseColumnMajor, ColumnMajor,
                     DiagonalLayout, MulAdd, mul, colsupport, rowsupport,
-                    _qr, _factorize, _copyto!, zero!, layout_replace_in_print_matrix
+                    _qr, _factorize, _copyto!, zero!, layout_replace_in_print_matrix,
+                    transposelayout, conjlayout
 
 import BlockArrays: blocksize, blockcheckbounds, BlockedUnitRange, blockisequal, DefaultBlockAxis,
                         Block, BlockSlice, unblock, block, blockindex,
                         _blocklengths2blocklasts, BlockIndexRange, sizes_from_blocks, BlockSlice1,
                         blockcolsupport, blockrowsupport, blockcolstart, blockcolstop, blockrowstart, blockrowstop,
-                        AbstractBlockLayout, BlockLayout, blocks, hasmatchingblocks, BlockStyle
+                        AbstractBlockLayout, BlockLayout, blocks, hasmatchingblocks, BlockStyle, BlockSlices
 
 import BandedMatrices: isbanded, bandwidths, bandwidth, banded_getindex, colrange,
                         inbands_setindex!, inbands_getindex, banded_setindex!,

src/adjtransblockbanded.jl

@@ -1,3 +1,7 @@
 
 blockbandwidths(A::AdjOrTrans) = reverse(blockbandwidths(parent(A)))
 subblockbandwidths(A::AdjOrTrans) = reverse(subblockbandwidths(parent(A)))
+transposelayout(::AbstractBlockBandedLayout) = BlockBandedLayout()
+transposelayout(::AbstractBandedBlockBandedLayout) = BandedBlockBandedLayout()
+conjlayout(::Type{<:Complex}, ::M) where M<:AbstractBandedBlockBandedLayout = ConjLayout{M}()
+conjlayout(::Type{<:Complex}, ::M) where M<:AbstractBlockBandedLayout = ConjLayout{M}()

src/interfaceimpl.jl

@@ -151,6 +151,17 @@ function blockrowsupport(::AbstractBandedLayout, B, k)
     findblock(n,first(rs)):findblock(n,last(rs))
 end
 
+# fixed block sizes, we can figure out how far we encroach other blocks by looking at last column
+function blockbandwidths(::AbstractBandedLayout, (a,b)::Tuple{BlockedUnitRange{<:AbstractRange}, OneTo{Int}}, A)
+    l,u = bandwidths(A)
+    m = min(length(a), l + length(b))
+    if u ≥ 0
+        Int(findblock(a,m))-1,0
+    else
+        Int(findblock(a,m))-1,1-Int(findblock(a,min(length(a),1-u)))
+    end
+end
+
 # ambiguity
 sub_materialize(::AbstractBandedLayout, V, ::Tuple{BlockedUnitRange,Base.OneTo{Int}}) = BandedMatrix(V)
 sub_materialize(::AbstractBandedLayout, V, ::Tuple{Base.OneTo{Int},BlockedUnitRange}) = BandedMatrix(V)

src/linalg.jl

@@ -90,15 +90,14 @@ end
 ####
 # BlockIndexRange subblocks
 ####
+sublayout(::AbstractBlockBandedLayout, ::Type{<:Tuple{<:BlockSlices, <:BlockSlices}}) = BlockBandedLayout()
 
-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{<:BlockSlices, <:BlockSlices}}) = BlockBandedColumnMajor()
+sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlices, <:BlockSlice{Block1}}}) = ColumnMajor()
+sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlices, <: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{<:BlockSlices, <:BlockSlices}}) =
     isblockbanded(parent(V))
 
 sub_materialize(::AbstractBlockBandedLayout, V, _) = BlockBandedMatrix(V)

test/test_adjtransblockbanded.jl

@@ -1,16 +1,41 @@
-using BlockBandedMatrices
+using BlockBandedMatrices, ArrayLayouts, Test
+import BlockBandedMatrices: BandedBlockBandedRowMajor, BandedBlockBandedRows, BandedBlockBandedColumns,BlockBandedRows, BlockBandedColumns
 
 @testset "Adj/Trans" begin
-    A = BandedBlockBandedMatrix(randn(ComplexF64,10,14), 1:4,2:5, (1,2), (2,1))
+    @testset "BandedBlockBanded" begin
+        A = BandedBlockBandedMatrix(randn(ComplexF64,10,14), 1:4,2:5, (1,2), (2,1))
 
-    @test A'[Block(1,1)] == A[Block(1,1)]'
-    @test A'[Block(2,3)] == A[Block(3,2)]'
-    @test transpose(A)[Block(1,1)] == transpose(A[Block(1,1)])
-    @test transpose(A)[Block(2,3)] == transpose(A[Block(3,2)])
+        @test MemoryLayout(transpose(A)) isa BandedBlockBandedRowMajor
+        @test MemoryLayout(A') isa BandedBlockBandedRows
+        @test MemoryLayout(transpose(A)') isa BandedBlockBandedColumns
 
-    @test blockbandwidths(A') == blockbandwidths(transpose(A)) == (2,1)
-    @test subblockbandwidths(A') == subblockbandwidths(transpose(A)) == (1,2)
+        @test A'[Block(1,1)] == A[Block(1,1)]'
+        @test A'[Block(2,3)] == A[Block(3,2)]'
+        @test transpose(A)[Block(1,1)] == transpose(A[Block(1,1)])
+        @test transpose(A)[Block(2,3)] == transpose(A[Block(3,2)])
 
-    @test BandedBlockBandedMatrix(A') == A'
-    @test BandedBlockBandedMatrix(transpose(A)) == transpose(A)
+        @test blockbandwidths(A') == blockbandwidths(transpose(A)) == (2,1)
+        @test subblockbandwidths(A') == subblockbandwidths(transpose(A)) == (1,2)
+
+        @test BandedBlockBandedMatrix(A') == A'
+        @test BandedBlockBandedMatrix(transpose(A)) == transpose(A)
+    end
+
+    @testset "BlockBanded" begin
+        A = BlockBandedMatrix(randn(ComplexF64,10,14), 1:4,2:5, (1,2))
+
+        @test MemoryLayout(transpose(A)) isa BlockBandedRows
+        @test MemoryLayout(A') isa BlockBandedRows
+        @test MemoryLayout(transpose(A)') isa BlockBandedColumns
+
+        @test A'[Block(1,1)] == A[Block(1,1)]'
+        @test A'[Block(2,3)] == A[Block(3,2)]'
+        @test transpose(A)[Block(1,1)] == transpose(A[Block(1,1)])
+        @test transpose(A)[Block(2,3)] == transpose(A[Block(3,2)])
+
+        @test blockbandwidths(A') == blockbandwidths(transpose(A)) == (2,1)
+
+        @test BlockBandedMatrix(A') == A'
+        @test BlockBandedMatrix(transpose(A)) == transpose(A)
+    end
 end

test/test_misc.jl

@@ -1,5 +1,5 @@
 using ArrayLayouts, BlockBandedMatrices, BandedMatrices, BlockArrays, LinearAlgebra, Test
-import BlockBandedMatrices: AbstractBandedBlockBandedMatrix, block, blockindex, blockcolsupport, blockrowsupport
+import BlockBandedMatrices: AbstractBandedBlockBandedMatrix, AbstractBlockBandedMatrix, block, blockindex, blockcolsupport, blockrowsupport
 import BandedMatrices: bandwidths, AbstractBandedMatrix, BandedStyle, bandeddata, BandedColumns, _BandedMatrix
 
 struct MyBandedBlockBandedMatrix <: AbstractBandedBlockBandedMatrix{Float64}
@@ -17,6 +17,19 @@ function Base.getindex(A::MyBandedBlockBandedMatrix, k::Int, j::Int)
 end
 
 
+struct MyBlockBandedMatrix <: AbstractBlockBandedMatrix{Float64}
+    A::BlockMatrix{Float64}
+end
+
+BlockBandedMatrices.blockbandwidths(::MyBlockBandedMatrix) = (1,1)
+Base.axes(A::MyBlockBandedMatrix) = axes(A.A)
+function Base.getindex(A::MyBlockBandedMatrix, k::Int, j::Int)
+    Kk, Jj = findblockindex(axes(A,1), k), findblockindex(axes(A,2), j)
+    -1 ≤ Int(block(Kk)-block(Jj)) ≤ 1 || return 0.0
+    A.A[k,j]
+end
+
+
 @testset "Interfaces" begin
     @testset "MyBandedBlockBandedMatrix" begin
         A = MyBandedBlockBandedMatrix(BlockMatrix(randn(6,6), 1:3, 1:3))
@@ -35,6 +48,25 @@ end
         @test A[Block(2)[1:2],Block.(2:3)] isa PseudoBlockArray
     end
 
+    @testset "MyBlockBandedMatrix" begin
+        A = MyBlockBandedMatrix(BlockMatrix(randn(6,6), 1:3, 1:3))
+        @test MemoryLayout(A) isa BlockBandedMatrices.BlockBandedLayout
+        @test MemoryLayout(A') isa BlockBandedMatrices.BlockBandedLayout
+        @test MemoryLayout(transpose(A)) isa BlockBandedMatrices.BlockBandedLayout
+        @test A[Block(3,3)] == A.A[Block(3,3)]
+        @test A[Block.(1:3),Block.(1:3)] == A
+
+        @test A[Block(3,3)] isa Matrix
+        @test A[Block(3)[2:3],Block(3)[1:2]] isa Matrix
+        @test A[Block(3)[2:3],Block(3)] isa Matrix
+        @test A[Block(3),Block(3)[1:2]] isa Matrix
+        @test A[Block.(1:3),Block.(1:3)] isa BlockBandedMatrix
+        @test A[Block(1),Block.(2:3)] isa PseudoBlockArray
+        @test A[Block.(2:3),Block(1)] isa PseudoBlockArray
+        @test A[Block.(2:3),Block(2)[1:2]] isa PseudoBlockArray
+        @test A[Block(2)[1:2],Block.(2:3)] isa PseudoBlockArray
+    end
+
     @testset "Zeros" begin
         Z = Zeros(blockedrange(1:3), blockedrange(1:3))
         B = BandedBlockBandedMatrix(Z)
@@ -140,5 +172,14 @@ end
         Q = Eye((a,))[Block(2),:]
         @test Q isa BandedMatrix
         @test blockrowsupport(Q,1) == Block.(2:2)
+
+        @testset "constant blocks" begin
+            a = blockedrange(Fill(2,5))
+            Q = Eye((a,))[:,Block(2)]
+            @test blockbandwidths(Q) == (1,-1)
+
+            B = _BandedMatrix(randn(5,length(a)), a, 3, 1)
+            @test blockbandwidths(B) == (4,0)
+        end
     end
 end