Skip to content

Commit d02efe6

Browse files
dlfivefiftydependabot[bot]
dlfivefifty
and
dependabot[bot]
authored
Banded Matrix extension (#388)
* Move BandedMatrices+BlockArrays code in BlockBandedMatrices to extension * 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](julia-actions/setup-julia@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] <[email protected]> Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> * Move over blockbanded code * Add tests * Update Project.toml * add tests * Update Project.toml --------- Signed-off-by: dependabot[bot] <[email protected]> Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com>
1 parent b361279 commit d02efe6

File tree

10 files changed

+211
-5
lines changed

10 files changed

+211
-5
lines changed

‎.github/workflows/CompatHelper.yml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@ jobs:
1515
run: which julia
1616
continue-on-error: true
1717
- name: Install Julia, but only if it is not already available in the PATH
18-
uses: julia-actions/setup-julia@v1
18+
uses: julia-actions/setup-julia@v2
1919
with:
2020
version: '1'
2121
arch: ${{ runner.arch }}

‎.github/workflows/ci.yml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -46,7 +46,7 @@ jobs:
4646
- x64
4747
steps:
4848
- uses: actions/checkout@v4
49-
- uses: julia-actions/setup-julia@v1
49+
- uses: julia-actions/setup-julia@v2
5050
with:
5151
version: ${{ matrix.version }}
5252
arch: ${{ matrix.arch }}

‎.github/workflows/docs.yml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,7 @@ jobs:
2626
runs-on: ubuntu-latest
2727
steps:
2828
- uses: actions/checkout@v4
29-
- uses: julia-actions/setup-julia@v1
29+
- uses: julia-actions/setup-julia@v2
3030
with:
3131
version: '1.10'
3232
- uses: julia-actions/cache@v1

‎.github/workflows/downstream.yml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -41,7 +41,7 @@ jobs:
4141

4242
steps:
4343
- uses: actions/checkout@v4
44-
- uses: julia-actions/setup-julia@v1
44+
- uses: julia-actions/setup-julia@v2
4545
with:
4646
version: ${{ matrix.julia-version }}
4747
arch: x64

‎Project.toml‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,19 +4,23 @@ version = "1.0.0-dev"
44

55
[deps]
66
ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
7+
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
78
FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
89
LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
910
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
1011

1112
[weakdeps]
13+
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
1214
LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
1315

1416
[extensions]
17+
BlockArraysBandedMatricesExt = "BandedMatrices"
1518
BlockArraysLazyArraysExt = "LazyArrays"
1619

1720
[compat]
1821
Aqua = "0.8"
1922
ArrayLayouts = "1.0.8"
23+
BandedMatrices = "1.0"
2024
Documenter = "1"
2125
FillArrays = "1"
2226
InfiniteArrays = "0.13"
@@ -31,6 +35,7 @@ julia = "1.6"
3135

3236
[extras]
3337
Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
38+
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
3439
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
3540
InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
3641
OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
@@ -40,4 +45,4 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
4045
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
4146

4247
[targets]
43-
test = ["Aqua", "Documenter", "InfiniteArrays", "OffsetArrays", "SparseArrays", "StaticArrays", "Test", "Random"]
48+
test = ["Aqua", "BandedMatrices", "Documenter", "InfiniteArrays", "OffsetArrays", "SparseArrays", "StaticArrays", "Test", "Random"]

‎ext/BlockArraysBandedMatricesExt.jl‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
module BlockArraysBandedMatricesExt
2+
3+
using BandedMatrices, BlockArrays
4+
using BlockArrays.ArrayLayouts
5+
import BandedMatrices: isbanded, AbstractBandedLayout, bandeddata, bandwidths
6+
import BlockArrays: blockcolsupport, blockrowsupport, AbstractBlockedUnitRange
7+
import ArrayLayouts: sub_materialize
8+
9+
10+
bandeddata(P::PseudoBlockMatrix) = bandeddata(P.blocks)
11+
bandwidths(P::PseudoBlockMatrix) = bandwidths(P.blocks)
12+
13+
function blockcolsupport(::AbstractBandedLayout, B, j)
14+
m,n = axes(B)
15+
cs = colsupport(B,n[j])
16+
findblock(m,first(cs)):findblock(m,last(cs))
17+
end
18+
19+
function blockrowsupport(::AbstractBandedLayout, B, k)
20+
m,n = axes(B)
21+
rs = rowsupport(B,m[k])
22+
findblock(n,first(rs)):findblock(n,last(rs))
23+
end
24+
25+
# ambiguity
26+
sub_materialize(::AbstractBandedLayout, V, ::Tuple{AbstractBlockedUnitRange,Base.OneTo{Int}}) = BandedMatrix(V)
27+
sub_materialize(::AbstractBandedLayout, V, ::Tuple{Base.OneTo{Int},AbstractBlockedUnitRange}) = BandedMatrix(V)
28+
29+
30+
end

‎src/BlockArrays.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -69,13 +69,15 @@ include("show.jl")
6969
include("blockreduce.jl")
7070
include("blockdeque.jl")
7171
include("blockarrayinterface.jl")
72+
include("blockbanded.jl")
7273

7374
@deprecate getblock(A::AbstractBlockArray{T,N}, I::Vararg{Integer, N}) where {T,N} view(A, Block(I))
7475
@deprecate getblock!(X, A::AbstractBlockArray{T,N}, I::Vararg{Integer, N}) where {T,N} copyto!(X, view(A, Block(I)))
7576
@deprecate setblock!(A::AbstractBlockArray{T,N}, v, I::Vararg{Integer, N}) where {T,N} (A[Block(I...)] = v)
7677

7778
if !isdefined(Base, :get_extension)
7879
include("../ext/BlockArraysLazyArraysExt.jl")
80+
include("../ext/BlockArraysBandedMatricesExt.jl")
7981
end
8082

8183
end # module

‎src/blockbanded.jl‎

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
const BlockDiagonal{T,VT<:Matrix{T}} = BlockMatrix{T,<:Diagonal{VT}}
2+
3+
BlockDiagonal(A) = mortar(Diagonal(A))
4+
5+
6+
# Block Bi/Tridiagonal
7+
const BlockTridiagonal{T,VT<:Matrix{T}} = BlockMatrix{T,<:Tridiagonal{VT}}
8+
const BlockBidiagonal{T,VT<:Matrix{T}} = BlockMatrix{T,<:Bidiagonal{VT}}
9+
10+
BlockTridiagonal(A,B,C) = mortar(Tridiagonal(A,B,C))
11+
BlockBidiagonal(A, B, uplo) = mortar(Bidiagonal(A,B,uplo))
12+
13+
sizes_from_blocks(A::Diagonal, _) = (size.(A.diag, 1), size.(A.diag,2))
14+
15+
function sizes_from_blocks(A::Tridiagonal, _)
16+
for k = 1:length(A.du)
17+
size(A.du[k],1) == size(A.d[k],1) || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
18+
size(A.du[k],2) == size(A.d[k+1],2) || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
19+
size(A.dl[k],1) == size(A.d[k+1],1) || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
20+
size(A.dl[k],2) == size(A.d[k],2) || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
21+
end
22+
(size.(A.d, 1), size.(A.d,2))
23+
end
24+
25+
function sizes_from_blocks(A::Bidiagonal, _)
26+
if A.uplo == 'U'
27+
for k = 1:length(A.ev)
28+
size(A.ev[k],1) == size(A.dv[k],1) || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
29+
size(A.ev[k],2) == size(A.dv[k+1],2) || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
30+
end
31+
else
32+
for k = 1:length(A.ev)
33+
size(A.ev[k],1) == size(A.dv[k+1],1) || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
34+
size(A.ev[k],2) == size(A.dv[k],2) || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
35+
end
36+
end
37+
(size.(A.dv, 1), size.(A.dv,2))
38+
end
39+
40+
41+
# viewblock needs to handle zero blocks
42+
@inline function viewblock(block_arr::BlockBidiagonal{T,VT}, KJ::Block{2}) where {T,VT<:AbstractMatrix}
43+
K,J = KJ.n
44+
@boundscheck blockcheckbounds(block_arr, K, J)
45+
l,u = block_arr.blocks.uplo == 'U' ? (0,1) : (1,0)
46+
-l ≤ (J-K) ≤ u || return convert(VT, Zeros{T}(length.(getindex.(axes(block_arr),(Block(K),Block(J))))...))
47+
block_arr.blocks[K,J]
48+
end
49+
50+
@inline function viewblock(block_arr::BlockTridiagonal{T,VT}, KJ::Block{2}) where {T,VT<:AbstractMatrix}
51+
K,J = KJ.n
52+
@boundscheck blockcheckbounds(block_arr, K, J)
53+
abs(J-K) ≥ 2 && return convert(VT, Zeros{T}(length.(getindex.(axes(block_arr),(Block(K),Block(J))))...))
54+
block_arr.blocks[K,J]
55+
end
56+
57+
checksquareblocks(A) = blockisequal(axes(A)...) || throw(DimensionMismatch("blocks are not square: block dimensions are $(axes(A))"))
58+
59+
# special case UniformScaling arithmetic
60+
for op in (:-, :+)
61+
@eval begin
62+
function $op(A::BlockDiagonal, λ::UniformScaling)
63+
checksquareblocks(A)
64+
mortar(Diagonal(broadcast($op, A.blocks.diag, Ref(λ))))
65+
end
66+
function $op(λ::UniformScaling, A::BlockDiagonal)
67+
checksquareblocks(A)
68+
mortar(Diagonal(broadcast($op, Ref(λ), A.blocks.diag)))
69+
end
70+
71+
function $op(A::BlockTridiagonal, λ::UniformScaling)
72+
checksquareblocks(A)
73+
mortar(Tridiagonal(broadcast($op, A.blocks.dl, Ref(0λ)),
74+
broadcast($op, A.blocks.d, Ref(λ)),
75+
broadcast($op, A.blocks.du, Ref(0λ))))
76+
end
77+
function $op(λ::UniformScaling, A::BlockTridiagonal)
78+
checksquareblocks(A)
79+
mortar(Tridiagonal(broadcast($op, Ref(0λ), A.blocks.dl),
80+
broadcast($op, Ref(λ), A.blocks.d),
81+
broadcast($op, Ref(0λ), A.blocks.du)))
82+
end
83+
function $op(A::BlockBidiagonal, λ::UniformScaling)
84+
checksquareblocks(A)
85+
mortar(Bidiagonal(broadcast($op, A.blocks.dv, Ref(λ)), A.blocks.ev, A.blocks.uplo))
86+
end
87+
function $op(λ::UniformScaling, A::BlockBidiagonal)
88+
checksquareblocks(A)
89+
mortar(Bidiagonal(broadcast($op, Ref(λ), A.blocks.dv), broadcast($op,A.blocks.ev), A.blocks.uplo))
90+
end
91+
end
92+
end

‎test/runtests.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26,3 +26,4 @@ include("test_blockproduct.jl")
2626
include("test_blockreduce.jl")
2727
include("test_blockdeque.jl")
2828
include("test_blockcholesky.jl")
29+
include("test_blockbanded.jl")

‎test/test_blockbanded.jl‎

Lines changed: 76 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,76 @@
1+
module TestBlockArraysBandedMatrices
2+
3+
using BlockArrays, LinearAlgebra, BandedMatrices, Test
4+
using BlockArrays: BlockDiagonal, BlockBidiagonal, BlockTridiagonal, blockcolsupport, blockrowsupport
5+
using BandedMatrices: _BandedMatrix
6+
7+
8+
@testset "Block-Banded" begin
9+
@testset "Block Diagonal" begin
10+
A = BlockDiagonal(fill([1 2],3))
11+
@test A[Block(1,1)] == [1 2]
12+
@test @inferred(A[Block(1,2)]) == [0 0]
13+
@test_throws DimensionMismatch A+I
14+
A = BlockDiagonal(fill([1 2; 1 2],3))
15+
@test A+I == I+A == mortar(Diagonal(fill([2 2; 1 3],3))) == Matrix(A) + I
16+
end
17+
18+
@testset "Block Bidiagonal" begin
19+
Bu = BlockBidiagonal(fill([1 2],4), fill([3 4],3), :U)
20+
Bl = BlockBidiagonal(fill([1 2],4), fill([3 4],3), :L)
21+
@test Bu[Block(1,1)] == Bl[Block(1,1)] == [1 2]
22+
@test @inferred(Bu[Block(1,2)]) == @inferred(Bl[Block(2,1)]) == [3 4]
23+
@test @inferred(view(Bu,Block(1,3))) == @inferred(Bu[Block(1,3)]) == [0 0]
24+
@test_throws DimensionMismatch Bu+I
25+
Bu = BlockBidiagonal(fill([1 2; 1 2],4), fill([3 4; 3 4],3), :U)
26+
Bl = BlockBidiagonal(fill([1 2; 1 2],4), fill([3 4; 3 4],3), :L)
27+
@test Bu+I == I+Bu == mortar(Bidiagonal(fill([2 2; 1 3],4), fill([3 4; 3 4],3), :U)) == Matrix(Bu) + I
28+
@test Bl+I == I+Bl == mortar(Bidiagonal(fill([2 2; 1 3],4), fill([3 4; 3 4],3), :L)) == Matrix(Bl) + I
29+
@test Bu-I == mortar(Bidiagonal(fill([0 2; 1 1],4), fill([3 4; 3 4],3), :U)) == Matrix(Bu) - I
30+
@test I-Bu == mortar(Bidiagonal(fill([0 -2; -1 -1],4), fill(-[3 4; 3 4],3), :U)) == I - Matrix(Bu)
31+
end
32+
33+
@testset "Block Tridiagonal" begin
34+
A = BlockTridiagonal(fill([1 2],3), fill([3 4],4), fill([4 5],3))
35+
@test A[Block(1,1)] == [3 4]
36+
@test @inferred(A[Block(1,2)]) == [4 5]
37+
@test @inferred(view(A,Block(1,3))) == @inferred(A[Block(1,3)]) == [0 0]
38+
@test_throws DimensionMismatch A+I
39+
A = BlockTridiagonal(fill([1 2; 1 2],3), fill([3 4; 3 4],4), fill([4 5; 4 5],3))
40+
@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
41+
@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
42+
@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)
43+
end
44+
45+
46+
@testset "Block-BandedMatrix" begin
47+
a = blockedrange(1:5)
48+
B = _BandedMatrix(PseudoBlockArray(randn(5,length(a)),(Base.OneTo(5),a)), a, 3, 1)
49+
@test blockcolsupport(B,Block(1)) == Block.(1:3)
50+
@test blockcolsupport(B,Block(3)) == Block.(2:4)
51+
@test blockrowsupport(B,Block(1)) == Block.(1:2)
52+
@test blockrowsupport(B,Block(4)) == Block.(3:5)
53+
54+
Q = Eye((a,))[:,Block(2)]
55+
@test Q isa BandedMatrix
56+
@test blockcolsupport(Q,Block(1)) == Block.(2:2)
57+
58+
Q = Eye((a,))[Block(2),:]
59+
@test Q isa BandedMatrix
60+
@test blockrowsupport(Q,Block(1)) == Block.(2:2)
61+
62+
@testset "constant blocks" begin
63+
a = blockedrange(Fill(2,5))
64+
Q = Eye((a,))[:,Block(2)]
65+
@test Q isa BandedMatrix
66+
end
67+
end
68+
69+
@testset "Banded PseudoMatrix" begin
70+
A = PseudoBlockArray(brand(5,4,1,2), [3,2], [2,2])
71+
@test bandwidths(A) == (1,2)
72+
@test BandedMatrix(A) == A
73+
end
74+
end
75+
76+
end # module

0 commit comments

Comments
 (0)