BlockBandedMatrices extension

Project.toml

@@ -4,7 +4,6 @@ version = "0.15.0-dev"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
-BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 Infinities = "e1ba4f0e-776d-440f-acd9-e1d2e9742647"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
@@ -13,12 +12,14 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 [weakdeps]
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [extensions]
 InfiniteArraysBandedMatricesExt = "BandedMatrices"
 InfiniteArraysBlockArraysExt = "BlockArrays"
+InfiniteArraysBlockBandedMatricesExt = "BlockBandedMatrices"
 InfiniteArraysDSPExt = "DSP"
 InfiniteArraysStatisticsExt = "Statistics"
 
@@ -28,25 +29,27 @@ ArrayLayouts = "1.8"
 BandedMatrices = "1.7.6"
 Base64 = "1"
 BlockArrays = "1.0"
+BlockBandedMatrices = "0.13"
 DSP = "0.7"
 FillArrays = "1.0"
 Infinities = "0.1.1"
 LazyArrays = "2.2.3"
 LinearAlgebra = "1.6"
-SparseArrays = "1"
-Statistics = "1"
-Test = "1"
+SparseArrays = "1.0"
+Statistics = "1.0"
+Test = "1.0"
 julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Test", "BandedMatrices", "BlockArrays", "Statistics", "SparseArrays", "Base64", "DSP"]
+test = ["Aqua", "Test", "BandedMatrices", "BlockArrays", "BlockBandedMatrices", "Statistics", "SparseArrays", "Base64", "DSP"]

ext/InfiniteArraysBlockBandedMatricesExt.jl

@@ -0,0 +1,65 @@
+module InfiniteArraysBlockBandedMatricesExt
+
+using InfiniteArrays, BlockBandedMatrices
+
+using InfiniteArrays.LinearAlgebra, BlockBandedMatrices.BlockArrays, BlockBandedMatrices.BandedMatrices
+
+import BlockBandedMatrices: BlockSkylineSizes, BlockBandedMatrix, _BlockSkylineMatrix
+import InfiniteArrays: InfiniteCardinal, OneToInf, InfStepRange
+import BlockBandedMatrices.BandedMatrices: _BandedMatrix
+
+const BlockTriPertToeplitz{T} = BlockMatrix{T,Tridiagonal{Matrix{T},Vcat{Matrix{T},1,Tuple{Vector{Matrix{T}},Fill{Matrix{T},1,Tuple{OneToInf{Int}}}}}},
+                                        NTuple{2,BlockedOneTo{Int,Vcat{Int,1,Tuple{Vector{Int},InfStepRange{Int,Int}}}}}}
+
+
+BlockBandedMatrices.blockbanded_colstop(A, x::InfiniteCardinal{0}) = x
+BlockBandedMatrices.blockbanded_rowstop(A, x::InfiniteCardinal{0}) = x
+
+function BlockSkylineSizes(A::BlockTriPertToeplitz, (l,u)::NTuple{2,Int})
+    N = max(length(A.blocks.du.args[1])+1,length(A.blocks.d.args[1]),length(A.blocks.dl.args[1]))
+    block_sizes = Vector{Int}(undef, N) # assume square
+    block_starts = BandedMatrix{Int}(undef, (N+l,N),  (l,u))
+    block_strides = Vector{Int}(undef, N)
+    for J=1:N
+        block_starts[max(1,J-u),J] = J == 1 ? 1 :
+                            block_starts[max(1,J-1-u),J-1]+block_sizes[J-1]*block_strides[J-1]
+
+        for K=max(1,J-u)+1:J+l
+            block_starts[K,J] = block_starts[K-1,J]+size(A[Block(K-1,J)],1)
+        end
+        block_strides[J] = block_starts[J+l,J] + size(A[Block(J+l,J)],1) - block_starts[max(1,J-u),J]
+        block_sizes[J] = size(A[Block(J,J)],2)
+    end
+
+    block_stride∞ = 0
+    for K=max(1,N+1-u):N+1+l
+        block_stride∞ += size(A[Block(K,N+1)],1)
+    end
+    block_size∞ = size(A[Block(N+1,N+1)],2)
+
+    bs∞ = fill(block_starts[max(1,N-u),N]+block_strides[N]*size(A[Block(N,N)],2):block_stride∞*block_size∞:∞, l+u+1)
+    for k=2:l+u+1
+        bs∞[k] = bs∞[k-1] .+ size(A[Block(N+1-u+k-1,N+1)],1)
+    end
+
+    BlockSkylineSizes(axes(A),
+                        _BandedMatrix(Hcat(block_starts.data, Vcat(adjoint.(bs∞)...)), ℵ₀, l, u),
+                        Vcat(block_strides, Fill(block_stride∞,∞)),
+                        Fill(l,∞),Fill(u,∞))
+end
+
+function BlockBandedMatrix(A::BlockTriPertToeplitz{T}, (l,u)::NTuple{2,Int}) where T
+    data = T[]
+    append!(data,vec(A[Block.(1:1+l),Block(1)]))
+    N = max(length(A.blocks.du.args[1])+1,length(A.blocks.d.args[1]),length(A.blocks.dl.args[1]))
+    for J=2:N
+        append!(data, vec(A[Block.(max(1,J-u):J+l),Block(J)]))
+    end
+    tl = mortar(Fill(vec(A[Block.(max(1,N+1-u):N+1+l),Block(N+1)]),∞))
+
+    B = _BlockSkylineMatrix(Vcat(data,tl), BlockSkylineSizes(A, (l,u)))
+end
+
+BlockBandedMatrix(A::BlockTriPertToeplitz) = BlockBandedMatrix(A, blockbandwidths(A))
+
+end #module

test/runtests.jl

@@ -1264,4 +1264,5 @@ end
 
 include("test_infconv.jl")
 include("test_infblock.jl")
-include("test_infbanded.jl")
+include("test_infbanded.jl")
+include("test_infblockbanded.jl")

test/test_infblockbanded.jl

@@ -1,22 +1,26 @@
+using InfiniteArrays, BlockBandedMatrices, ArrayLayouts, FillArrays, BlockArrays, Test
+using Base: oneto
 
 @testset "BlockBanded" begin
-            a = b = c = 0.0
-            n = mortar(Fill.(oneto(∞), oneto(∞)))
-            k = mortar(Base.OneTo.(oneto(∞)))
-            Dy = BlockBandedMatrices._BandedBlockBandedMatrix((k .+ (b + c))', axes(k, 1), (-1, 1), (-1, 1))
-            N = 100
-            @test Dy[Block.(1:N), Block.(1:N)] == BlockBandedMatrices._BandedBlockBandedMatrix((k.+(b+c))[Block.(1:N)]', axes(k, 1)[Block.(1:N)], (-1, 1), (-1, 1))
-            @test colsupport(Dy, axes(Dy,2)) == 1:∞
-            @test rowsupport(Dy, axes(Dy,1)) == 2:∞
-        end
+    @testset "Triangle Recurrence" begin
+        a = b = c = 0.0
+        n = mortar(Fill.(oneto(∞), oneto(∞)))
+        k = mortar(Base.OneTo.(oneto(∞)))
+        Dy = BlockBandedMatrices._BandedBlockBandedMatrix((k .+ (b + c))', axes(k, 1), (-1, 1), (-1, 1))
+        N = 100
+        @test Dy[Block.(1:N), Block.(1:N)] == BlockBandedMatrices._BandedBlockBandedMatrix((k.+(b+c))[Block.(1:N)]', axes(k, 1)[Block.(1:N)], (-1, 1), (-1, 1))
+        @test colsupport(Dy, axes(Dy,2)) == 1:∞
+        @test rowsupport(Dy, axes(Dy,1)) == 2:∞
+    end
 
-
-
-        @testset "BlockTridiagonal" begin
-            A = BlockTridiagonal(Vcat([fill(1.0, 2, 1), Matrix(1.0I, 2, 2), Matrix(1.0I, 2, 2), Matrix(1.0I, 2, 2)], Fill(Matrix(1.0I, 2, 2), ∞)),
-                Vcat([zeros(1, 1)], Fill(zeros(2, 2), ∞)),
-                Vcat([fill(1.0, 1, 2), Matrix(1.0I, 2, 2)], Fill(Matrix(1.0I, 2, 2), ∞)))
 
-            @test isblockbanded(A)    
-            @test BlockBandedMatrix(A)[1:100, 1:100] == BlockBandedMatrix(A, (2, 1))[1:100, 1:100] == BlockBandedMatrix(A, (1, 1))[1:100, 1:100] == A[1:100, 1:100]
-        end
+
+    @testset "BlockTridiagonal" begin
+        A = BlockTridiagonal(Vcat([fill(1.0, 2, 1), Matrix(1.0I, 2, 2), Matrix(1.0I, 2, 2), Matrix(1.0I, 2, 2)], Fill(Matrix(1.0I, 2, 2), ∞)),
+            Vcat([zeros(1, 1)], Fill(zeros(2, 2), ∞)),
+            Vcat([fill(1.0, 1, 2), Matrix(1.0I, 2, 2)], Fill(Matrix(1.0I, 2, 2), ∞)))
+
+        @test isblockbanded(A)    
+        @test BlockBandedMatrix(A)[1:100, 1:100] == BlockBandedMatrix(A, (2, 1))[1:100, 1:100] == BlockBandedMatrix(A, (1, 1))[1:100, 1:100] == A[1:100, 1:100]
+    end
+end