1+ using LazyBandedMatrices, InfiniteArrays, ArrayLayouts, LazyArrays, BlockArrays, BandedMatrices, BlockBandedMatrices, LinearAlgebra, Test
2+ using InfiniteArrays: TridiagonalToeplitzLayout, BidiagonalToeplitzLayout, TridiagonalToeplitzLayout
3+ using Base: oneto
4+ using BlockArrays: blockcolsupport
5+ using LazyArrays: arguments
6+ using LazyBandedMatrices: BroadcastBandedBlockBandedLayout
7+
8+ const InfiniteArraysBlockArraysExt = Base. get_extension (InfiniteArrays, :InfiniteArraysBlockArraysExt )
9+ const LazyBandedMatricesInfiniteArraysExt = Base. get_extension (LazyBandedMatrices, :LazyBandedMatricesInfiniteArraysExt )
10+
11+ const OneToInfBlocks = InfiniteArraysBlockArraysExt. OneToInfBlocks
12+ const InfKronTravBandedBlockBandedLayout = LazyBandedMatricesInfiniteArraysExt. InfKronTravBandedBlockBandedLayout
13+
14+ @testset " ∞ LazyBandedMatrices" begin
15+ @test MemoryLayout (LazyBandedMatrices. Tridiagonal (Fill (1 ,∞), Zeros (∞), Fill (3 ,∞))) isa TridiagonalToeplitzLayout
16+ @test MemoryLayout (LazyBandedMatrices. Bidiagonal (Fill (1 ,∞), Zeros (∞), :U )) isa BidiagonalToeplitzLayout
17+ @test MemoryLayout (LazyBandedMatrices. SymTridiagonal (Fill (1 ,∞), Zeros (∞))) isa TridiagonalToeplitzLayout
18+
19+ T = LazyBandedMatrices. Tridiagonal (Fill (1 ,∞), Zeros (∞), Fill (3 ,∞))
20+ @test T[2 : ∞,3 : ∞] isa SubArray
21+ @test exp .(T) isa BroadcastMatrix
22+ @test exp .(T)[2 : ∞,3 : ∞][1 : 10 ,1 : 10 ] == exp .(T[2 : ∞,3 : ∞])[1 : 10 ,1 : 10 ] == exp .(T[2 : 11 ,3 : 12 ])
23+ @test exp .(T)[2 : ∞,3 : ∞] isa BroadcastMatrix
24+ @test exp .(T[2 : ∞,3 : ∞]) isa BroadcastMatrix
25+
26+ B = LazyBandedMatrices. Bidiagonal (Fill (1 ,∞), Zeros (∞), :U )
27+ @test B[2 : ∞,3 : ∞] isa SubArray
28+ @test exp .(B) isa BroadcastMatrix
29+ @test exp .(B)[2 : ∞,3 : ∞][1 : 10 ,1 : 10 ] == exp .(B[2 : ∞,3 : ∞])[1 : 10 ,1 : 10 ] == exp .(B[2 : 11 ,3 : 12 ])
30+ @test exp .(B)[2 : ∞,3 : ∞] isa BroadcastMatrix
31+
32+ @testset " Diagonal{Fill} * Bidiagonal" begin
33+ A, B = Diagonal (Fill (2 ,∞)) , LazyBandedMatrices. Bidiagonal (exp .(1 : ∞), exp .(1 : ∞), :L )
34+ @test (A* B)[1 : 10 ,1 : 10 ] ≈ (B* A)[1 : 10 ,1 : 10 ] ≈ 2 B[1 : 10 ,1 : 10 ]
35+ end
36+
37+ @testset " ∞-unit blocks" begin
38+ @test unitblocks (oneto (∞)) ≡ blockedrange (Ones {Int} (∞))
39+ @test unitblocks (2 : ∞) == 2 : ∞
40+
41+ @test unitblocks (oneto (∞))[Block .(2 : ∞)] == 2 : ∞
42+ end
43+
44+ @testset " concat" begin
45+ a = unitblocks (1 : ∞)
46+ b = exp .(a)
47+ c = BlockBroadcastArray (vcat, a, b)
48+ @test length (c) == ∞
49+ @test blocksize (c) == (∞,)
50+ @test c[Block (5 )] == [a[5 ], b[5 ]]
51+
52+ A = unitblocks (BandedMatrix (0 => 1 : ∞, 1 => Fill (2.0 , ∞), - 1 => Fill (3.0 , ∞)))
53+ B = BlockBroadcastArray (hvcat, 2 , A, Zeros (axes (A)), Zeros (axes (A)), A)
54+ @test B[Block (3 , 3 )] == [A[3 , 3 ] 0 ; 0 A[3 , 3 ]]
55+ @test B[Block (3 , 4 )] == [A[3 , 4 ] 0 ; 0 A[3 , 4 ]]
56+ @test B[Block (3 , 5 )] == [A[3 , 5 ] 0 ; 0 A[3 , 5 ]]
57+ @test blockbandwidths (B) == (1 , 1 )
58+ @test subblockbandwidths (B) == (0 , 0 )
59+ @test B[Block .(1 : 10 ), Block .(1 : 10 )] isa BlockSkylineMatrix
60+
61+ C = BlockBroadcastArray (hvcat, 2 , A, A, A, A)
62+ @test C[Block (3 , 3 )] == fill (A[3 , 3 ], 2 , 2 )
63+ @test C[Block (3 , 4 )] == fill (A[3 , 4 ], 2 , 2 )
64+ @test C[Block (3 , 5 )] == fill (A[3 , 5 ], 2 , 2 )
65+ @test blockbandwidths (C) == (1 , 1 )
66+ @test subblockbandwidths (C) == (1 , 1 )
67+ @test B[Block .(1 : 10 ), Block .(1 : 10 )] isa BlockSkylineMatrix
68+ end
69+
70+ @testset " DiagTrav" begin
71+ C = zeros (∞,∞);
72+ C[1 : 3 ,1 : 4 ] .= [1 2 3 4 ; 5 6 7 8 ; 9 10 11 12 ]
73+ A = DiagTrav (C)
74+ @test blockcolsupport (A) == Block .(1 : 6 )
75+ @test A[Block .(1 : 7 )] == [1 ; 5 ; 2 ; 9 ; 6 ; 3 ; 0 ; 10 ; 7 ; 4 ; 0 ; 0 ; 11 ; 8 ; 0 ; 0 ; 0 ; 0 ; 12 ; zeros (9 )]
76+
77+ C = zeros (∞,4 );
78+ C[1 : 3 ,1 : 4 ] .= [1 2 3 4 ; 5 6 7 8 ; 9 10 11 12 ]
79+ A = DiagTrav (C)
80+ @test A[Block .(1 : 7 )] == [1 ; 5 ; 2 ; 9 ; 6 ; 3 ; 0 ; 10 ; 7 ; 4 ; 0 ; 0 ; 11 ; 8 ; 0 ; 0 ; 0 ; 12 ; zeros (4 )]
81+
82+ C = zeros (3 ,∞);
83+ C[1 : 3 ,1 : 4 ] .= [1 2 3 4 ; 5 6 7 8 ; 9 10 11 12 ]
84+ A = DiagTrav (C)
85+ @test A[Block .(1 : 7 )] == [1 ; 5 ; 2 ; 9 ; 6 ; 3 ; 10 ; 7 ; 4 ; 11 ; 8 ; 0 ; 12 ; zeros (5 )]
86+ end
87+
88+ @testset " KronTrav" begin
89+ Δ = BandedMatrix (1 => Ones (∞), - 1 => Ones (∞)) / 2
90+ A = KronTrav (Δ - 2 I, Eye (∞))
91+ @test axes (A, 1 ) isa OneToInfBlocks
92+ V = view (A, Block .(Base. OneTo (3 )), Block .(Base. OneTo (3 )))
93+
94+ @test MemoryLayout (A) isa InfKronTravBandedBlockBandedLayout
95+ @test MemoryLayout (V) isa LazyBandedMatrices. KronTravBandedBlockBandedLayout
96+
97+ @test A[Block .(Base. OneTo (3 )), Block .(Base. OneTo (3 ))] isa KronTrav
98+
99+ u = A * [1 ; zeros (∞)]
100+ @test u[1 : 3 ] == A[1 : 3 , 1 ]
101+ @test bandwidths (view (A, Block (1 , 1 ))) == (1 , 1 )
102+
103+ @test A* A isa KronTrav
104+ @test (A* A)[Block .(Base. OneTo (3 )), Block .(Base. OneTo (3 ))] ≈ A[Block .(1 : 3 ), Block .(1 : 4 )]A[Block .(1 : 4 ), Block .(1 : 3 )]
105+ end
106+
107+ @testset " BlockHcat copyto!" begin
108+ n = mortar (Fill .(oneto (∞), oneto (∞)))
109+ k = mortar (Base. OneTo .(oneto (∞)))
110+
111+ a = b = c = 0.0
112+ dat = BlockHcat (
113+ BroadcastArray ((n, k, b, bc1) -> (k + b - 1 ) * (n + k + bc1) / (2 k + bc1), n, k, b, b + c - 1 ),
114+ BroadcastArray ((n, k, abc, bc, bc1) -> (n + k + abc) * (k + bc) / (2 k + bc1), n, k, a + b + c, b + c, b + c - 1 )
115+ )
116+ N = 1000
117+ KR = Block .(Base. OneTo (N))
118+ V = view (dat, Block .(Base. OneTo (N)), :)
119+ @test MemoryLayout (V) isa LazyArrays. ApplyLayout{typeof (hcat)}
120+ @test BlockedArray (V)[Block .(1 : 5 ), :] == dat[Block .(1 : 5 ), :]
121+ V = view (dat' , :, Block .(Base. OneTo (N)))
122+ @test MemoryLayout (V) isa LazyArrays. ApplyLayout{typeof (vcat)}
123+ a = dat. arrays[1 ]'
124+ N = 100
125+ KR = Block .(Base. OneTo (N))
126+ v = view (a, :, KR)
127+ @time r = BlockedArray (v)
128+ @test v == r
129+ end
130+
131+ @testset " Symmetric" begin
132+ k = mortar (Base. OneTo .(oneto (∞)))
133+ n = mortar (Fill .(oneto (∞), oneto (∞)))
134+
135+ dat = BlockHcat (
136+ BlockBroadcastArray (hcat, float .(k), Zeros ((axes (n, 1 ),)), float .(n)),
137+ Zeros ((axes (n, 1 ), Base. OneTo (3 ))),
138+ Zeros ((axes (n, 1 ), Base. OneTo (3 ))))
139+ M = BlockBandedMatrices. _BandedBlockBandedMatrix (dat' , axes (k, 1 ), (1 , 1 ), (1 , 1 ))
140+ Ms = Symmetric (M)
141+ @test blockbandwidths (M) == (1 , 1 )
142+ @test blockbandwidths (Ms) == (1 , 1 )
143+ @test Ms[Block .(1 : 5 ), Block .(1 : 5 )] == Symmetric (M[Block .(1 : 5 ), Block .(1 : 5 )])
144+ @test Ms[Block .(1 : 5 ), Block .(1 : 5 )] isa BandedBlockBandedMatrix
145+
146+ b = [ones (10 ); zeros (∞)]
147+ @test (Ms * b)[Block .(1 : 6 )] == Ms[Block .(1 : 6 ), Block .(1 : 4 )]* ones (10 )
148+ @test ((Ms * Ms) * b)[Block .(1 : 6 )] == (Ms * (Ms * b))[Block .(1 : 6 )]
149+ @test ((Ms + Ms) * b)[Block .(1 : 6 )] == (2 * (Ms * b))[Block .(1 : 6 )]
150+
151+ dat = BlockBroadcastArray (hcat, float .(k), Zeros ((axes (n, 1 ),)), float .(n))
152+ M = BlockBandedMatrices. _BandedBlockBandedMatrix (dat' , axes (k, 1 ), (- 1 , 1 ), (1 , 1 ))
153+ Ms = Symmetric (M)
154+ @test Symmetric ((M+ M)[Block .(1 : 10 ), Block .(1 : 10 )]) == (Ms+ Ms)[Block .(1 : 10 ), Block .(1 : 10 )]
155+ end
156+
157+ @testset " BlockBidiagonal" begin
158+ B = mortar (LazyBandedMatrices. Bidiagonal (Fill ([1 2 ; 3 4 ], ∞), Fill ([1 2 ; 3 4 ], ∞), :U ));
159+ # TODO : copy BlockBidiagonal code from BlockBandedMatrices to LazyBandedMatrices
160+ @test B[Block (2 , 3 )] == [1 2 ; 3 4 ]
161+ @test_broken B[Block (1 , 3 )] == Zeros (2 , 2 )
162+ end
163+
164+ @testset " KronTrav" begin
165+ Δ = BandedMatrix (1 => Ones (∞) / 2 , - 1 => Ones (∞))
166+ A = KronTrav (Δ, Eye (∞))
167+ @test A[Block (100 , 101 )] isa BandedMatrix
168+ @test A[Block (100 , 100 )] isa BandedMatrix
169+ @test A[Block .(1 : 5 ), Block .(1 : 5 )] isa BandedBlockBandedMatrix
170+ B = KronTrav (Eye (∞), Δ)
171+ @test B[Block (100 , 101 )] isa BandedMatrix
172+ @test B[Block (100 , 100 )] isa BandedMatrix
173+ V = view (A + B, Block .(1 : 5 ), Block .(1 : 5 ))
174+ @test MemoryLayout (typeof (V)) isa BroadcastBandedBlockBandedLayout{typeof (+ )}
175+ @test arguments (V) == (A[Block .(1 : 5 ), Block .(1 : 5 )], B[Block .(1 : 5 ), Block .(1 : 5 )])
176+ @test (A+ B)[Block .(1 : 5 ), Block .(1 : 5 )] == A[Block .(1 : 5 ), Block .(1 : 5 )] + B[Block .(1 : 5 ), Block .(1 : 5 )]
177+
178+ @test blockbandwidths (A + B) == (1 , 1 )
179+ @test blockbandwidths (2 A) == (1 , 1 )
180+ @test blockbandwidths (2 * (A + B)) == (1 , 1 )
181+
182+ @test subblockbandwidths (A + B) == (1 , 1 )
183+ @test subblockbandwidths (2 A) == (1 , 1 )
184+ @test subblockbandwidths (2 * (A + B)) == (1 , 1 )
185+ end
186+
187+ @testset " BlockTridiagonal" begin
188+ T = mortar (LazyBandedMatrices. Tridiagonal (Fill ([1 2 ; 3 4 ], ∞), Fill ([1 2 ; 3 4 ], ∞), Fill ([1 2 ; 3 4 ], ∞)));
189+ # TODO : copy BlockBidiagonal code from BlockBandedMatrices to LazyBandedMatrices
190+ @test T[Block (2 , 2 )] == [1 2 ; 3 4 ]
191+ @test_broken T[Block (1 , 3 )] == Zeros (2 , 2 )
192+ end
193+ end
0 commit comments