Skip to content

Commit ca34659

Browse files
dlfivefiftydlfivefifty
authored
Fix A*B where A is lower-constant diagonal (#87)
1 parent d53da6f commit ca34659

File tree

4 files changed

+130
-122
lines changed

4 files changed

+130
-122
lines changed

‎Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
name = "InfiniteLinearAlgebra"
22
uuid = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
3-
version = "0.5.12"
3+
version = "0.5.13"
44

55
[deps]
66
ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

‎src/banded/infbanded.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -416,7 +416,7 @@ function _bandedfill_mul(M::MulAdd, ::Tuple{InfAxes,InfAxes}, ::Tuple{InfAxes,In
416416
l,u = Al+Bl,Au+Bu
417417
m = min(Au+Al,Bl+Bu)+1
418418
λ = getindex_value(bandeddata(A))*getindex_value(bandeddata(B))
419-
ret = _BandedMatrix(Hcat(Array{typeof(λ)}(undef, l+u+1,u), [1:m-1; Fill(m,l+u-2m+3); m-1:-1:1]*Fill(λ,1,∞)), ℵ₀, l, u)
419+
ret = _BandedMatrix(Hcat(Array{typeof(λ)}(undef, l+u+1,max(0,u)), [1:m-1; Fill(m,l+u-2m+3); m-1:-1:1]*Fill(λ,1,∞)), ℵ₀, l, u)
420420
mul!(view(ret, 1:l+u,1:u), view(A,1:l+u,1:u+Bl), view(B,1:u+Bl,1:u))
421421
ret
422422
end

‎test/runtests.jl

Lines changed: 1 addition & 120 deletions
Original file line numberDiff line numberDiff line change
@@ -14,126 +14,7 @@ import InfiniteArrays: OneToInf, oneto, RealInfinity
1414
import LazyBandedMatrices: BroadcastBandedBlockBandedLayout, BroadcastBandedLayout, LazyBandedLayout
1515

1616
include("test_infconv.jl")
17-
18-
19-
@testset "∞-banded" begin
20-
@testset "Diagonal and BandedMatrix" begin
21-
D = Diagonal(Fill(2,∞))
22-
23-
B = D[1:∞,2:∞]
24-
@test B isa BandedMatrix
25-
@test B[1:10,1:10] == diagm(-1 => Fill(2,9))
26-
@test B[1:∞,2:∞] isa BandedMatrix
27-
28-
A = BandedMatrix(0 => 1:∞, 1=> Fill(2.0,∞), -1 => Fill(3.0,∞))
29-
x = [1; 2; zeros(∞)]
30-
@test A*x isa Vcat
31-
@test (A*x)[1:10] == A[1:10,1:10]*x[1:10]
32-
33-
@test InfBandCartesianIndices(0)[1:5] == CartesianIndex.(1:5,1:5)
34-
@test InfBandCartesianIndices(1)[1:5] == CartesianIndex.(1:5,2:6)
35-
@test InfBandCartesianIndices(-1)[1:5] == CartesianIndex.(2:6,1:5)
36-
37-
@test A[band(0)][2:10] == 2:10
38-
@test D[band(0)] ≡ Fill(2,∞)
39-
@test D[band(1)] ≡ Fill(0,∞)
40-
41-
@test B*A*x isa Vcat
42-
@test (B*A*x)[1:10] == [0; 10; 14; 12; zeros(6)]
43-
44-
@test _BandedMatrix((1:∞)', ∞, -1, 1) isa BandedMatrix
45-
end
46-
47-
@testset "∞-Toeplitz" begin
48-
A = BandedMatrix(1 => Fill(2im,∞), 2 => Fill(-1,∞), 3 => Fill(2,∞), -2 => Fill(-4,∞), -3 => Fill(-2im,∞))
49-
@test A isa InfToeplitz
50-
@test MemoryLayout(typeof(A.data)) == ConstRows()
51-
@test MemoryLayout(typeof(A)) == BandedToeplitzLayout()
52-
V = view(A,:,3:∞)
53-
@test MemoryLayout(typeof(bandeddata(V))) == ConstRows()
54-
@test MemoryLayout(typeof(V)) == BandedToeplitzLayout()
55-
@test BandedMatrix(V) isa InfToeplitz
56-
@test A[:,3:end] isa InfToeplitz
57-
58-
@test (A + 2I)[1:10,1:10] == (2I + A)[1:10,1:10] == A[1:10,1:10] + 2I
59-
60-
@test MemoryLayout(Tridiagonal(Fill(1,∞), Fill(2,∞), Fill(3,∞))) isa TridiagonalToeplitzLayout
61-
@test MemoryLayout(Bidiagonal(Fill(1,∞), Fill(2,∞), :U)) isa BidiagonalToeplitzLayout
62-
@test MemoryLayout(SymTridiagonal(Fill(1,∞), Fill(2,∞))) isa TridiagonalToeplitzLayout
63-
@test MemoryLayout(LazyBandedMatrices.Tridiagonal(Fill(1,∞), Zeros(∞), Fill(3,∞))) isa TridiagonalToeplitzLayout
64-
@test MemoryLayout(LazyBandedMatrices.Bidiagonal(Fill(1,∞), Zeros(∞), :U)) isa BidiagonalToeplitzLayout
65-
@test MemoryLayout(LazyBandedMatrices.SymTridiagonal(Fill(1,∞), Zeros(∞))) isa TridiagonalToeplitzLayout
66-
67-
T = LazyBandedMatrices.Tridiagonal(Fill(1,∞), Zeros(∞), Fill(3,∞))
68-
@test T[2:∞,3:∞] isa SubArray
69-
@test exp.(T) isa BroadcastMatrix
70-
@test exp.(T)[2:∞,3:∞] isa SubArray
71-
72-
B = LazyBandedMatrices.Bidiagonal(Fill(1,∞), Zeros(∞), :U)
73-
@test B[2:∞,3:∞] isa SubArray
74-
@test exp.(B) isa BroadcastMatrix
75-
@test exp.(B)[2:∞,3:∞] isa SubArray
76-
77-
@testset "algebra" begin
78-
T = Tridiagonal(Fill(1,∞), Fill(2,∞), Fill(3,∞))
79-
@test T isa InfiniteLinearAlgebra.TriToeplitz
80-
@test (T + 2I)[1:10,1:10] == (2I + T)[1:10,1:10] == T[1:10,1:10] + 2I
81-
end
82-
end
83-
84-
@testset "Pert-Toeplitz" begin
85-
@testset "Inf Pert" begin
86-
A = BandedMatrix(-2 => Vcat(Float64[], Fill(1/4,∞)), 0 => Vcat([1.0+im,2,3],Fill(0,∞)), 1 => Vcat(Float64[], Fill(1,∞)))
87-
@test A isa PertToeplitz
88-
@test MemoryLayout(typeof(A)) == PertToeplitzLayout()
89-
V = view(A,2:∞,2:∞)
90-
@test MemoryLayout(typeof(V)) == PertToeplitzLayout()
91-
@test BandedMatrix(V) isa PertToeplitz
92-
@test A[2:∞,2:∞] isa PertToeplitz
93-
94-
@test (A + 2I)[1:10,1:10] == (2I + A)[1:10,1:10] == A[1:10,1:10] + 2I
95-
end
96-
97-
@testset "TriPert" begin
98-
A = SymTridiagonal(Vcat([1,2.], Fill(2.,∞)), Vcat([3.,4.], Fill.(0.5,∞)))
99-
@test A isa InfiniteLinearAlgebra.SymTriPertToeplitz
100-
@test (A + 2I)[1:10,1:10] == (2I + A)[1:10,1:10] == A[1:10,1:10] + 2I
101-
102-
A = Tridiagonal(Vcat([3.,4.], Fill.(0.5,∞)), Vcat([1,2.], Fill(2.,∞)), Vcat([3.,4.], Fill.(0.5,∞)))
103-
@test A isa InfiniteLinearAlgebra.TriPertToeplitz
104-
@test Adjoint(A) isa InfiniteLinearAlgebra.AdjTriPertToeplitz
105-
@test (A + 2I)[1:10,1:10] == (2I + A)[1:10,1:10] == A[1:10,1:10] + 2I
106-
@test (Adjoint(A) + 2I)[1:10,1:10] == (2I + Adjoint(A))[1:10,1:10] == Adjoint(A)[1:10,1:10] + 2I
107-
end
108-
109-
110-
@testset "InfBanded" begin
111-
A = _BandedMatrix(Fill(2,4,∞),ℵ₀,2,1)
112-
B = _BandedMatrix(Fill(3,2,∞),ℵ₀,-1,2)
113-
@test mul(A,A) isa PertToeplitz
114-
@test A*A isa PertToeplitz
115-
@test (A*A)[1:20,1:20] == A[1:20,1:23]*A[1:23,1:20]
116-
@test (A*B)[1:20,1:20] == A[1:20,1:23]*B[1:23,1:20]
117-
end
118-
end
119-
120-
@testset "adjortrans" begin
121-
A = BandedMatrix(0 => 1:∞, 1=> Fill(2.0+im,∞), -1 => Fill(3.0,∞))
122-
@test copy(A')[1:10,1:10] == (A')[1:10,1:10]
123-
@test copy(transpose(A))[1:10,1:10] == transpose(A)[1:10,1:10]
124-
end
125-
126-
@testset "Eye subindex" begin
127-
@test Eye(∞)[:,1:3][1:5,:] == Eye(∞)[Base.Slice(oneto(∞)),1:3][1:5,:] == Eye(5,3)
128-
@test Eye(∞)[1:3,:][:,1:5] == Eye(∞)[1:3,Base.Slice(oneto(∞))][:,1:5] == Eye(3,5)
129-
@test Eye(∞)[:,:][1:5,1:3] == Eye(∞)[Base.Slice(oneto(∞)),Base.Slice(oneto(∞))][1:5,1:3] == Eye(5,3)
130-
end
131-
132-
@testset "band(0) indexing" begin
133-
D = ApplyArray(*, Diagonal(1:∞), Diagonal(1:∞))
134-
@test D[band(0)][1:10] == (1:10).^2
135-
end
136-
end
17+
include("test_infbanded.jl")
13718

13819
@testset "∞-block arrays" begin
13920
@testset "fixed block size" begin

‎test/test_infbanded.jl

Lines changed: 127 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,127 @@
1+
using InfiniteLinearAlgebra, InfiniteArrays, BandedMatrices, FillArrays, Test
2+
import BandedMatrices: _BandedMatrix
3+
4+
@testset "∞-banded" begin
5+
@testset "Diagonal and BandedMatrix" begin
6+
D = Diagonal(Fill(2,∞))
7+
8+
B = D[1:∞,2:∞]
9+
@test B isa BandedMatrix
10+
@test B[1:10,1:10] == diagm(-1 => Fill(2,9))
11+
@test B[1:∞,2:∞] isa BandedMatrix
12+
13+
A = BandedMatrix(0 => 1:∞, 1=> Fill(2.0,∞), -1 => Fill(3.0,∞))
14+
x = [1; 2; zeros(∞)]
15+
@test A*x isa Vcat
16+
@test (A*x)[1:10] == A[1:10,1:10]*x[1:10]
17+
18+
@test InfBandCartesianIndices(0)[1:5] == CartesianIndex.(1:5,1:5)
19+
@test InfBandCartesianIndices(1)[1:5] == CartesianIndex.(1:5,2:6)
20+
@test InfBandCartesianIndices(-1)[1:5] == CartesianIndex.(2:6,1:5)
21+
22+
@test A[band(0)][2:10] == 2:10
23+
@test D[band(0)] ≡ Fill(2,∞)
24+
@test D[band(1)] ≡ Fill(0,∞)
25+
26+
@test B*A*x isa Vcat
27+
@test (B*A*x)[1:10] == [0; 10; 14; 12; zeros(6)]
28+
29+
@test _BandedMatrix((1:∞)', ∞, -1, 1) isa BandedMatrix
30+
end
31+
32+
@testset "∞-Toeplitz" begin
33+
A = BandedMatrix(1 => Fill(2im,∞), 2 => Fill(-1,∞), 3 => Fill(2,∞), -2 => Fill(-4,∞), -3 => Fill(-2im,∞))
34+
@test A isa InfToeplitz
35+
@test MemoryLayout(typeof(A.data)) == ConstRows()
36+
@test MemoryLayout(typeof(A)) == BandedToeplitzLayout()
37+
V = view(A,:,3:∞)
38+
@test MemoryLayout(typeof(bandeddata(V))) == ConstRows()
39+
@test MemoryLayout(typeof(V)) == BandedToeplitzLayout()
40+
@test BandedMatrix(V) isa InfToeplitz
41+
@test A[:,3:end] isa InfToeplitz
42+
43+
@test (A + 2I)[1:10,1:10] == (2I + A)[1:10,1:10] == A[1:10,1:10] + 2I
44+
45+
@test MemoryLayout(Tridiagonal(Fill(1,∞), Fill(2,∞), Fill(3,∞))) isa TridiagonalToeplitzLayout
46+
@test MemoryLayout(Bidiagonal(Fill(1,∞), Fill(2,∞), :U)) isa BidiagonalToeplitzLayout
47+
@test MemoryLayout(SymTridiagonal(Fill(1,∞), Fill(2,∞))) isa TridiagonalToeplitzLayout
48+
@test MemoryLayout(LazyBandedMatrices.Tridiagonal(Fill(1,∞), Zeros(∞), Fill(3,∞))) isa TridiagonalToeplitzLayout
49+
@test MemoryLayout(LazyBandedMatrices.Bidiagonal(Fill(1,∞), Zeros(∞), :U)) isa BidiagonalToeplitzLayout
50+
@test MemoryLayout(LazyBandedMatrices.SymTridiagonal(Fill(1,∞), Zeros(∞))) isa TridiagonalToeplitzLayout
51+
52+
T = LazyBandedMatrices.Tridiagonal(Fill(1,∞), Zeros(∞), Fill(3,∞))
53+
@test T[2:∞,3:∞] isa SubArray
54+
@test exp.(T) isa BroadcastMatrix
55+
@test exp.(T)[2:∞,3:∞] isa SubArray
56+
57+
B = LazyBandedMatrices.Bidiagonal(Fill(1,∞), Zeros(∞), :U)
58+
@test B[2:∞,3:∞] isa SubArray
59+
@test exp.(B) isa BroadcastMatrix
60+
@test exp.(B)[2:∞,3:∞] isa SubArray
61+
62+
@testset "algebra" begin
63+
T = Tridiagonal(Fill(1,∞), Fill(2,∞), Fill(3,∞))
64+
@test T isa InfiniteLinearAlgebra.TriToeplitz
65+
@test (T + 2I)[1:10,1:10] == (2I + T)[1:10,1:10] == T[1:10,1:10] + 2I
66+
end
67+
end
68+
69+
@testset "Pert-Toeplitz" begin
70+
@testset "Inf Pert" begin
71+
A = BandedMatrix(-2 => Vcat(Float64[], Fill(1/4,∞)), 0 => Vcat([1.0+im,2,3],Fill(0,∞)), 1 => Vcat(Float64[], Fill(1,∞)))
72+
@test A isa PertToeplitz
73+
@test MemoryLayout(typeof(A)) == PertToeplitzLayout()
74+
V = view(A,2:∞,2:∞)
75+
@test MemoryLayout(typeof(V)) == PertToeplitzLayout()
76+
@test BandedMatrix(V) isa PertToeplitz
77+
@test A[2:∞,2:∞] isa PertToeplitz
78+
79+
@test (A + 2I)[1:10,1:10] == (2I + A)[1:10,1:10] == A[1:10,1:10] + 2I
80+
end
81+
82+
@testset "TriPert" begin
83+
A = SymTridiagonal(Vcat([1,2.], Fill(2.,∞)), Vcat([3.,4.], Fill.(0.5,∞)))
84+
@test A isa InfiniteLinearAlgebra.SymTriPertToeplitz
85+
@test (A + 2I)[1:10,1:10] == (2I + A)[1:10,1:10] == A[1:10,1:10] + 2I
86+
87+
A = Tridiagonal(Vcat([3.,4.], Fill.(0.5,∞)), Vcat([1,2.], Fill(2.,∞)), Vcat([3.,4.], Fill.(0.5,∞)))
88+
@test A isa InfiniteLinearAlgebra.TriPertToeplitz
89+
@test Adjoint(A) isa InfiniteLinearAlgebra.AdjTriPertToeplitz
90+
@test (A + 2I)[1:10,1:10] == (2I + A)[1:10,1:10] == A[1:10,1:10] + 2I
91+
@test (Adjoint(A) + 2I)[1:10,1:10] == (2I + Adjoint(A))[1:10,1:10] == Adjoint(A)[1:10,1:10] + 2I
92+
end
93+
94+
95+
@testset "InfBanded" begin
96+
A = _BandedMatrix(Fill(2,4,∞),ℵ₀,2,1)
97+
B = _BandedMatrix(Fill(3,2,∞),ℵ₀,-1,2)
98+
@test mul(A,A) isa PertToeplitz
99+
@test A*A isa PertToeplitz
100+
@test (A*A)[1:20,1:20] == A[1:20,1:23]*A[1:23,1:20]
101+
@test (A*B)[1:20,1:20] == A[1:20,1:23]*B[1:23,1:20]
102+
end
103+
end
104+
105+
@testset "adjortrans" begin
106+
A = BandedMatrix(0 => 1:∞, 1=> Fill(2.0+im,∞), -1 => Fill(3.0,∞))
107+
@test copy(A')[1:10,1:10] == (A')[1:10,1:10]
108+
@test copy(transpose(A))[1:10,1:10] == transpose(A)[1:10,1:10]
109+
end
110+
111+
@testset "Eye subindex" begin
112+
@test Eye(∞)[:,1:3][1:5,:] == Eye(∞)[Base.Slice(oneto(∞)),1:3][1:5,:] == Eye(5,3)
113+
@test Eye(∞)[1:3,:][:,1:5] == Eye(∞)[1:3,Base.Slice(oneto(∞))][:,1:5] == Eye(3,5)
114+
@test Eye(∞)[:,:][1:5,1:3] == Eye(∞)[Base.Slice(oneto(∞)),Base.Slice(oneto(∞))][1:5,1:3] == Eye(5,3)
115+
end
116+
117+
@testset "band(0) indexing" begin
118+
D = ApplyArray(*, Diagonal(1:∞), Diagonal(1:∞))
119+
@test D[band(0)][1:10] == (1:10).^2
120+
end
121+
122+
@testset "Fill * Banded" begin
123+
A = _BandedMatrix(Ones(1,∞), ∞, 1,-1)
124+
B = _BandedMatrix(Fill(1.0π,1,∞), ∞, 0,0)
125+
@test (A*B)[1:10,1:10] ≈ BandedMatrix(-1 => Fill(1.0π,9))
126+
end
127+
end

0 commit comments

Comments
 (0)