diff --git a/Project.toml b/Project.toml index 1e225fc..4af185f 100644 --- a/Project.toml +++ b/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,14 +29,15 @@ 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] @@ -43,10 +45,11 @@ 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"] diff --git a/ext/InfiniteArraysBlockBandedMatricesExt.jl b/ext/InfiniteArraysBlockBandedMatricesExt.jl new file mode 100644 index 0000000..79fd77d --- /dev/null +++ b/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 \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 9eed2bc..6a95be4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1264,4 +1264,5 @@ end include("test_infconv.jl") include("test_infblock.jl") -include("test_infbanded.jl") \ No newline at end of file +include("test_infbanded.jl") +include("test_infblockbanded.jl") \ No newline at end of file diff --git a/test/test_infblockbanded.jl b/test/test_infblockbanded.jl index 01e6331..523a179 100644 --- a/test/test_infblockbanded.jl +++ b/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 \ No newline at end of file + + @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 \ No newline at end of file