Banded Matrix extension

.github/workflows/CompatHelper.yml

@@ -15,7 +15,7 @@ jobs:
         run: which julia
         continue-on-error: true
       - name: Install Julia, but only if it is not already available in the PATH
-        uses: julia-actions/setup-julia@v1
+        uses: julia-actions/setup-julia@v2
         with:
           version: '1'
           arch: ${{ runner.arch }}

.github/workflows/ci.yml

@@ -46,7 +46,7 @@ jobs:
           - x64
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}

.github/workflows/docs.yml

@@ -26,7 +26,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: '1.10'
       - uses: julia-actions/cache@v1

.github/workflows/downstream.yml

@@ -41,7 +41,7 @@ jobs:
 
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.julia-version }}
           arch: x64

Project.toml

@@ -4,19 +4,23 @@ version = "1.0.0-dev"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
+BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [weakdeps]
+BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 
 [extensions]
+BlockArraysBandedMatricesExt = "BandedMatrices"
 BlockArraysLazyArraysExt = "LazyArrays"
 
 [compat]
 Aqua = "0.8"
 ArrayLayouts = "1.0.8"
+BandedMatrices = "1.0"
 Documenter = "1"
 FillArrays = "1"
 InfiniteArrays = "0.13"
@@ -31,6 +35,7 @@ julia = "1.6"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
@@ -40,4 +45,4 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Documenter", "InfiniteArrays", "OffsetArrays", "SparseArrays", "StaticArrays", "Test", "Random"]
+test = ["Aqua", "BandedMatrices", "Documenter", "InfiniteArrays", "OffsetArrays", "SparseArrays", "StaticArrays", "Test", "Random"]

ext/BlockArraysBandedMatricesExt.jl

@@ -0,0 +1,30 @@
+module BlockArraysBandedMatricesExt
+
+using BandedMatrices, BlockArrays
+using BlockArrays.ArrayLayouts
+import BandedMatrices: isbanded, AbstractBandedLayout, bandeddata, bandwidths
+import BlockArrays: blockcolsupport, blockrowsupport, AbstractBlockedUnitRange
+import ArrayLayouts: sub_materialize
+
+
+bandeddata(P::PseudoBlockMatrix) = bandeddata(P.blocks)
+bandwidths(P::PseudoBlockMatrix) = bandwidths(P.blocks)
+
+function blockcolsupport(::AbstractBandedLayout, B, j)
+    m,n = axes(B)
+    cs = colsupport(B,n[j])
+    findblock(m,first(cs)):findblock(m,last(cs))
+end
+
+function blockrowsupport(::AbstractBandedLayout, B, k)
+    m,n = axes(B)
+    rs = rowsupport(B,m[k])
+    findblock(n,first(rs)):findblock(n,last(rs))
+end
+
+# ambiguity
+sub_materialize(::AbstractBandedLayout, V, ::Tuple{AbstractBlockedUnitRange,Base.OneTo{Int}}) = BandedMatrix(V)
+sub_materialize(::AbstractBandedLayout, V, ::Tuple{Base.OneTo{Int},AbstractBlockedUnitRange}) = BandedMatrix(V)
+
+
+end

src/BlockArrays.jl

@@ -69,13 +69,15 @@ include("show.jl")
 include("blockreduce.jl")
 include("blockdeque.jl")
 include("blockarrayinterface.jl")
+include("blockbanded.jl")
 
 @deprecate getblock(A::AbstractBlockArray{T,N}, I::Vararg{Integer, N}) where {T,N} view(A, Block(I))
 @deprecate getblock!(X, A::AbstractBlockArray{T,N}, I::Vararg{Integer, N}) where {T,N} copyto!(X, view(A, Block(I)))
 @deprecate setblock!(A::AbstractBlockArray{T,N}, v, I::Vararg{Integer, N}) where {T,N} (A[Block(I...)] = v)
 
 if !isdefined(Base, :get_extension)
     include("../ext/BlockArraysLazyArraysExt.jl")
+    include("../ext/BlockArraysBandedMatricesExt.jl")
 end
 
 end # module

src/blockbanded.jl

@@ -0,0 +1,92 @@
+const BlockDiagonal{T,VT<:Matrix{T}} = BlockMatrix{T,<:Diagonal{VT}}
+
+BlockDiagonal(A) = mortar(Diagonal(A))
+
+
+# Block Bi/Tridiagonal
+const BlockTridiagonal{T,VT<:Matrix{T}} = BlockMatrix{T,<:Tridiagonal{VT}}
+const BlockBidiagonal{T,VT<:Matrix{T}} = BlockMatrix{T,<:Bidiagonal{VT}}
+
+BlockTridiagonal(A,B,C) = mortar(Tridiagonal(A,B,C))
+BlockBidiagonal(A, B, uplo) = mortar(Bidiagonal(A,B,uplo))
+
+sizes_from_blocks(A::Diagonal, _) = (size.(A.diag, 1), size.(A.diag,2))
+
+function sizes_from_blocks(A::Tridiagonal, _)
+    for k = 1:length(A.du)
+        size(A.du[k],1) == size(A.d[k],1) || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
+        size(A.du[k],2) == size(A.d[k+1],2) || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
+        size(A.dl[k],1) == size(A.d[k+1],1) || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
+        size(A.dl[k],2) == size(A.d[k],2) || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
+    end
+    (size.(A.d, 1), size.(A.d,2))
+end
+
+function sizes_from_blocks(A::Bidiagonal, _)
+    if A.uplo == 'U'
+        for k = 1:length(A.ev)
+            size(A.ev[k],1) == size(A.dv[k],1) || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
+            size(A.ev[k],2) == size(A.dv[k+1],2) || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
+        end
+    else
+        for k = 1:length(A.ev)
+            size(A.ev[k],1) == size(A.dv[k+1],1) || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
+            size(A.ev[k],2) == size(A.dv[k],2) || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
+        end
+    end
+    (size.(A.dv, 1), size.(A.dv,2))
+end
+
+
+# viewblock needs to handle zero blocks
+@inline function 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 = block_arr.blocks.uplo == 'U' ? (0,1) : (1,0)
+    -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 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]
+end
+
+checksquareblocks(A) = blockisequal(axes(A)...) || throw(DimensionMismatch("blocks are not square: block dimensions are $(axes(A))"))
+
+# special case UniformScaling arithmetic
+for op in (:-, :+)
+    @eval begin
+        function $op(A::BlockDiagonal, λ::UniformScaling)
+            checksquareblocks(A)
+            mortar(Diagonal(broadcast($op, A.blocks.diag, Ref(λ))))
+        end
+        function $op(λ::UniformScaling, A::BlockDiagonal)
+            checksquareblocks(A)
+            mortar(Diagonal(broadcast($op, Ref(λ), A.blocks.diag)))
+        end
+
+        function $op(A::BlockTridiagonal, λ::UniformScaling)
+            checksquareblocks(A)
+            mortar(Tridiagonal(broadcast($op, A.blocks.dl, Ref(0λ)),
+                               broadcast($op, A.blocks.d, Ref(λ)),
+                               broadcast($op, A.blocks.du, Ref(0λ))))
+        end
+        function $op(λ::UniformScaling, A::BlockTridiagonal)
+            checksquareblocks(A)
+            mortar(Tridiagonal(broadcast($op, Ref(0λ), A.blocks.dl),
+                               broadcast($op, Ref(λ), A.blocks.d),
+                               broadcast($op, Ref(0λ), A.blocks.du)))
+        end
+        function $op(A::BlockBidiagonal, λ::UniformScaling)
+            checksquareblocks(A)
+            mortar(Bidiagonal(broadcast($op, A.blocks.dv, Ref(λ)), A.blocks.ev, A.blocks.uplo))
+        end
+        function $op(λ::UniformScaling, A::BlockBidiagonal)
+            checksquareblocks(A)
+            mortar(Bidiagonal(broadcast($op, Ref(λ), A.blocks.dv), broadcast($op,A.blocks.ev), A.blocks.uplo))
+        end
+    end
+end

test/runtests.jl

@@ -26,3 +26,4 @@ include("test_blockproduct.jl")
 include("test_blockreduce.jl")
 include("test_blockdeque.jl")
 include("test_blockcholesky.jl")
+include("test_blockbanded.jl")

test/test_blockbanded.jl

@@ -0,0 +1,76 @@
+module TestBlockArraysBandedMatrices
+
+using BlockArrays, LinearAlgebra, BandedMatrices, Test
+using BlockArrays: BlockDiagonal, BlockBidiagonal, BlockTridiagonal, blockcolsupport, blockrowsupport
+using BandedMatrices: _BandedMatrix
+
+
+@testset "Block-Banded" begin
+    @testset "Block Diagonal" begin
+        A = BlockDiagonal(fill([1 2],3))
+        @test A[Block(1,1)] == [1 2]
+        @test @inferred(A[Block(1,2)]) == [0 0]
+        @test_throws DimensionMismatch A+I
+        A = BlockDiagonal(fill([1 2; 1 2],3))
+        @test A+I == I+A == mortar(Diagonal(fill([2 2; 1 3],3))) == Matrix(A) + I
+    end
+
+    @testset "Block Bidiagonal" begin
+        Bu = BlockBidiagonal(fill([1 2],4), fill([3 4],3), :U)
+        Bl = BlockBidiagonal(fill([1 2],4), fill([3 4],3), :L)
+        @test Bu[Block(1,1)] == Bl[Block(1,1)] == [1 2]
+        @test @inferred(Bu[Block(1,2)]) == @inferred(Bl[Block(2,1)]) == [3 4]
+        @test @inferred(view(Bu,Block(1,3))) == @inferred(Bu[Block(1,3)]) == [0 0]
+        @test_throws DimensionMismatch Bu+I
+        Bu = BlockBidiagonal(fill([1 2; 1 2],4), fill([3 4; 3 4],3), :U)
+        Bl = BlockBidiagonal(fill([1 2; 1 2],4), fill([3 4; 3 4],3), :L)
+        @test Bu+I == I+Bu == mortar(Bidiagonal(fill([2 2; 1 3],4), fill([3 4; 3 4],3), :U)) == Matrix(Bu) + I
+        @test Bl+I == I+Bl == mortar(Bidiagonal(fill([2 2; 1 3],4), fill([3 4; 3 4],3), :L)) == Matrix(Bl) + I
+        @test Bu-I == mortar(Bidiagonal(fill([0 2; 1 1],4), fill([3 4; 3 4],3), :U)) == Matrix(Bu) - I
+        @test I-Bu == mortar(Bidiagonal(fill([0 -2; -1 -1],4), fill(-[3 4; 3 4],3), :U)) == I - Matrix(Bu)
+    end
+
+    @testset "Block Tridiagonal" begin
+        A = BlockTridiagonal(fill([1 2],3), fill([3 4],4), fill([4 5],3))
+        @test A[Block(1,1)] == [3 4]
+        @test @inferred(A[Block(1,2)]) == [4 5]
+        @test @inferred(view(A,Block(1,3))) == @inferred(A[Block(1,3)]) == [0 0]
+        @test_throws DimensionMismatch A+I
+        A = BlockTridiagonal(fill([1 2; 1 2],3), fill([3 4; 3 4],4), fill([4 5; 4 5],3))
+        @test A+I == I+A == mortar(Tridiagonal(fill([1 2; 1 2],3), fill([4 4; 3 5],4), fill([4 5; 4 5],3))) == Matrix(A) + I
+        @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 "Block-BandedMatrix" begin
+        a = blockedrange(1:5)
+        B = _BandedMatrix(PseudoBlockArray(randn(5,length(a)),(Base.OneTo(5),a)), a, 3, 1)
+        @test blockcolsupport(B,Block(1)) == Block.(1:3)
+        @test blockcolsupport(B,Block(3)) == Block.(2:4)
+        @test blockrowsupport(B,Block(1)) == Block.(1:2)
+        @test blockrowsupport(B,Block(4)) == Block.(3:5)
+
+        Q = Eye((a,))[:,Block(2)]
+        @test Q isa BandedMatrix
+        @test blockcolsupport(Q,Block(1)) == Block.(2:2)
+
+        Q = Eye((a,))[Block(2),:]
+        @test Q isa BandedMatrix
+        @test blockrowsupport(Q,Block(1)) == Block.(2:2)
+
+        @testset "constant blocks" begin
+            a = blockedrange(Fill(2,5))
+            Q = Eye((a,))[:,Block(2)]
+            @test Q isa BandedMatrix
+        end
+    end
+
+    @testset "Banded PseudoMatrix" begin
+        A = PseudoBlockArray(brand(5,4,1,2), [3,2], [2,2])
+        @test bandwidths(A) == (1,2)
+        @test BandedMatrix(A) == A
+    end
+end
+
+end # module