From 8accc9dd805e533ff89d405c273cbd73d6d1c50d Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Wed, 20 Dec 2023 11:19:39 +0000 Subject: [PATCH 1/7] Move BandedMatrices+BlockArrays code in BlockBandedMatrices to extension --- Project.toml | 8 ++++++++ ext/BlockArraysBandedMatricesExt.jl | 28 ++++++++++++++++++++++++++++ 2 files changed, 36 insertions(+) create mode 100644 ext/BlockArraysBandedMatricesExt.jl diff --git a/Project.toml b/Project.toml index 6443f52e..6fb7afa5 100644 --- a/Project.toml +++ b/Project.toml @@ -7,9 +7,17 @@ ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +[weakdeps] +BandedMatrices = "aae01518-5342-5314-be14-df237901396f" + +[extensions] +BlockArraysBandedMatricesExt = "SparseArrays" + + [compat] Aqua = "0.8" ArrayLayouts = "1.0.8" +BandedMatrices = "1.0" Documenter = "0.27" FillArrays = "1" LinearAlgebra = "1.6" diff --git a/ext/BlockArraysBandedMatricesExt.jl b/ext/BlockArraysBandedMatricesExt.jl new file mode 100644 index 00000000..21810d29 --- /dev/null +++ b/ext/BlockArraysBandedMatricesExt.jl @@ -0,0 +1,28 @@ +module BlockArraysBandedMatricesExt + +using BandedMatrices, BlockArrays +import BandedMatrices: isbanded, AbstractBandedLayout, bandeddata, bandwidths +import BlockArrays: blockcolsupport, blockrowsupport, 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{BlockedUnitRange,Base.OneTo{Int}}) = BandedMatrix(V) +sub_materialize(::AbstractBandedLayout, V, ::Tuple{Base.OneTo{Int},BlockedUnitRange}) = BandedMatrix(V) + + +end From 56ffbf1b470958d6f079a9307c63faca291da473 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 8 Apr 2024 16:48:03 +0530 Subject: [PATCH 2/7] Bump julia-actions/setup-julia from 1 to 2 (#387) Bumps [julia-actions/setup-julia](https://github.com/julia-actions/setup-julia) from 1 to 2. - [Release notes](https://github.com/julia-actions/setup-julia/releases) - [Commits](https://github.com/julia-actions/setup-julia/compare/v1...v2) --- updated-dependencies: - dependency-name: julia-actions/setup-julia dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/CompatHelper.yml | 2 +- .github/workflows/ci.yml | 2 +- .github/workflows/docs.yml | 2 +- .github/workflows/downstream.yml | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml index 09181610..8ad0284c 100644 --- a/.github/workflows/CompatHelper.yml +++ b/.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 }} diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index bb3e44c7..992c41d9 100644 --- a/.github/workflows/ci.yml +++ b/.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 }} diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 1570fc3f..5cf49ab6 100644 --- a/.github/workflows/docs.yml +++ b/.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 diff --git a/.github/workflows/downstream.yml b/.github/workflows/downstream.yml index 894646aa..b0a79864 100644 --- a/.github/workflows/downstream.yml +++ b/.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 From f8350d7f3282c8bdde17fb3bc116a0c59907f374 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 15 Apr 2024 13:43:42 +0100 Subject: [PATCH 3/7] Move over blockbanded code --- src/BlockArrays.jl | 1 + src/blockbanded.jl | 92 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 93 insertions(+) create mode 100644 src/blockbanded.jl diff --git a/src/BlockArrays.jl b/src/BlockArrays.jl index fcc98b38..e5dc86ec 100644 --- a/src/BlockArrays.jl +++ b/src/BlockArrays.jl @@ -69,6 +69,7 @@ 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))) diff --git a/src/blockbanded.jl b/src/blockbanded.jl new file mode 100644 index 00000000..5c7e23d5 --- /dev/null +++ b/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 \ No newline at end of file From fb99633a03c5efc5ec0220295450c40404973b72 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 15 Apr 2024 16:03:13 +0100 Subject: [PATCH 4/7] Add tests --- ext/BlockArraysBandedMatricesExt.jl | 8 ++-- test/runtests.jl | 1 + test/test_blockbanded.jl | 66 +++++++++++++++++++++++++++++ 3 files changed, 72 insertions(+), 3 deletions(-) create mode 100644 test/test_blockbanded.jl diff --git a/ext/BlockArraysBandedMatricesExt.jl b/ext/BlockArraysBandedMatricesExt.jl index 21810d29..a2d9d2ee 100644 --- a/ext/BlockArraysBandedMatricesExt.jl +++ b/ext/BlockArraysBandedMatricesExt.jl @@ -1,8 +1,10 @@ module BlockArraysBandedMatricesExt using BandedMatrices, BlockArrays +using BlockArrays.ArrayLayouts import BandedMatrices: isbanded, AbstractBandedLayout, bandeddata, bandwidths -import BlockArrays: blockcolsupport, blockrowsupport, sub_materialize +import BlockArrays: blockcolsupport, blockrowsupport, AbstractBlockedUnitRange +import ArrayLayouts: sub_materialize bandeddata(P::PseudoBlockMatrix) = bandeddata(P.blocks) @@ -21,8 +23,8 @@ function blockrowsupport(::AbstractBandedLayout, B, k) end # ambiguity -sub_materialize(::AbstractBandedLayout, V, ::Tuple{BlockedUnitRange,Base.OneTo{Int}}) = BandedMatrix(V) -sub_materialize(::AbstractBandedLayout, V, ::Tuple{Base.OneTo{Int},BlockedUnitRange}) = BandedMatrix(V) +sub_materialize(::AbstractBandedLayout, V, ::Tuple{AbstractBlockedUnitRange,Base.OneTo{Int}}) = BandedMatrix(V) +sub_materialize(::AbstractBandedLayout, V, ::Tuple{Base.OneTo{Int},AbstractBlockedUnitRange}) = BandedMatrix(V) end diff --git a/test/runtests.jl b/test/runtests.jl index a1873d8d..f6ce9823 100644 --- a/test/runtests.jl +++ b/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") diff --git a/test/test_blockbanded.jl b/test/test_blockbanded.jl new file mode 100644 index 00000000..e2c88744 --- /dev/null +++ b/test/test_blockbanded.jl @@ -0,0 +1,66 @@ +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 +end \ No newline at end of file From a1dcef706096dff47ce1fdeac7dbdb96b77b246d Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 15 Apr 2024 16:18:31 +0100 Subject: [PATCH 5/7] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 8b270b14..75206d0d 100644 --- a/Project.toml +++ b/Project.toml @@ -44,4 +44,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"] From ba9a70a5c7fc5ce70b75b5ebadcc9eacb776c757 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 15 Apr 2024 17:50:04 +0100 Subject: [PATCH 6/7] add tests --- src/BlockArrays.jl | 1 + test/test_blockbanded.jl | 12 +++++++++++- 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/src/BlockArrays.jl b/src/BlockArrays.jl index e5dc86ec..5f12a465 100644 --- a/src/BlockArrays.jl +++ b/src/BlockArrays.jl @@ -77,6 +77,7 @@ include("blockbanded.jl") if !isdefined(Base, :get_extension) include("../ext/BlockArraysLazyArraysExt.jl") + include("../ext/BlockArraysBandedMatricesExt.jl") end end # module diff --git a/test/test_blockbanded.jl b/test/test_blockbanded.jl index e2c88744..4e69e431 100644 --- a/test/test_blockbanded.jl +++ b/test/test_blockbanded.jl @@ -1,3 +1,5 @@ +module TestBlockArraysBandedMatrices + using BlockArrays, LinearAlgebra, BandedMatrices, Test using BlockArrays: BlockDiagonal, BlockBidiagonal, BlockTridiagonal, blockcolsupport, blockrowsupport using BandedMatrices: _BandedMatrix @@ -63,4 +65,12 @@ using BandedMatrices: _BandedMatrix @test Q isa BandedMatrix end end -end \ No newline at end of file + + @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 \ No newline at end of file From 4256d3565b20ec18378de2a50537a64c26c92431 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 15 Apr 2024 18:10:31 +0100 Subject: [PATCH 7/7] Update Project.toml --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 75206d0d..d256d3ab 100644 --- a/Project.toml +++ b/Project.toml @@ -35,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"