Skip to content

Commit 50d30ad

Browse files
dlfivefiftydlfivefifty
committed
Allow general bandwiddths in R
1 parent 0325f28 commit 50d30ad

File tree

3 files changed

+176
-122
lines changed

3 files changed

+176
-122
lines changed

‎src/InfiniteLinearAlgebra.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -125,5 +125,6 @@ include("infqr.jl")
125125
include("inful.jl")
126126
include("infcholesky.jl")
127127
include("banded/bidiagonalconjugation.jl")
128+
include("banded/tridiagonalconjugation.jl")
128129

129130
end # module

‎src/banded/tridiagonalconjugation.jl‎

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
"""
2+
upper_mul_tri_triview(U, X) == Tridiagonal(U*X) where U is Upper triangular BandedMatrix and X is Tridiagonal
3+
"""
4+
function upper_mul_tri_triview(U::BandedMatrix, X::Tridiagonal)
5+
T = promote_type(eltype(U), eltype(X))
6+
n = size(U,1)
7+
upper_mul_tri_triview!(Tridiagonal(Vector{T}(undef, n-1), Vector{T}(undef, n), Vector{T}(undef, n-1)), U, X)
8+
end
9+
10+
function upper_mul_tri_triview!(UX::Tridiagonal, U::BandedMatrix, X::Tridiagonal)
11+
n = size(UX,1)
12+
13+
14+
Xdl, Xd, Xdu = X.dl, X.d, X.du
15+
UXdl, UXd, UXdu = UX.dl, UX.d, UX.du
16+
Udat = U.data
17+
18+
l,u = bandwidths(U)
19+
20+
@assert size(U) == (n,n)
21+
@assert l == 0 && u ≥ 2
22+
# Tridiagonal bands can be resized
23+
@assert length(Xdl)+1 == length(Xd) == length(Xdu)+1 == length(UXdl)+1 == length(UXd) == length(UXdu)+1 == n
24+
25+
j = 1
26+
aⱼ, cⱼ = Xd[1], Xdl[1]
27+
Uⱼⱼ, Uⱼⱼ₊₁, Uⱼⱼ₊₂ = Udat[u+1,1], Udat[u,2], Udat[u-1,3] # U[j,j], U[j,j+1], U[j,j+2]
28+
UXd[1] = Uⱼⱼ*aⱼ + Uⱼⱼ₊₁*cⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j]
29+
bⱼ, aⱼ, cⱼ, cⱼ₋₁ = Xdu[1], Xd[2], Xdl[2], cⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j]
30+
UXdu[1] = Uⱼⱼ*bⱼ + Uⱼⱼ₊₁*aⱼ + Uⱼⱼ₊₂*cⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+1]*X[j+1,j]
31+
32+
@inbounds for j = 2:n-2
33+
Uⱼⱼ, Uⱼⱼ₊₁, Uⱼⱼ₊₂ = Udat[u+1,j], Udat[u,j+1], Udat[u-1,j+2] # U[j,j], U[j,j+1], U[j,j+2]
34+
UXdl[j-1] = Uⱼⱼ*cⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1]
35+
UXd[j] = Uⱼⱼ*aⱼ + Uⱼⱼ₊₁*cⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j]
36+
bⱼ, aⱼ, cⱼ, cⱼ₋₁ = Xdu[j], Xd[j+1], Xdl[j+1], cⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j]
37+
UXdu[j] = Uⱼⱼ*bⱼ + Uⱼⱼ₊₁*aⱼ + Uⱼⱼ₊₂*cⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+2]*X[j+2,j+1]
38+
end
39+
40+
j = n-1
41+
Uⱼⱼ, Uⱼⱼ₊₁ = Udat[u+1,j], Udat[u,j+1] # U[j,j], U[j,j+1]
42+
UXdl[j-1] = Uⱼⱼ*cⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1]
43+
UXd[j] = Uⱼⱼ*aⱼ + Uⱼⱼ₊₁*cⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j]
44+
bⱼ, aⱼ, cⱼ₋₁ = Xdu[j], Xd[j+1], cⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j]
45+
UXdu[j] = Uⱼⱼ*bⱼ + Uⱼⱼ₊₁*aⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+2]*X[j+2,j+1]
46+
47+
j = n
48+
Uⱼⱼ = Udat[u+1,j] # U[j,j]
49+
UXdl[j-1] = Uⱼⱼ*cⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1]
50+
UXd[j] = Uⱼⱼ*aⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j]
51+
52+
UX
53+
end
54+
55+
56+
# X*R^{-1} = X*[1/R₁₁ -R₁₂/(R₁₁R₂₂) -R₁₃/R₂₂ …
57+
# 0 1/R₂₂ -R₂₃/R₃₃
58+
# 1/R₃₃
59+
60+
tri_mul_invupper_triview(X::Tridiagonal, R::BandedMatrix) = tri_mul_invupper_triview!(similar(X, promote_type(eltype(X), eltype(R))), X, R)
61+
62+
function tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::BandedMatrix)
63+
n = size(X,1)
64+
Xdl, Xd, Xdu = X.dl, X.d, X.du
65+
Ydl, Yd, Ydu = Y.dl, Y.d, Y.du
66+
Rdat = R.data
67+
68+
l,u = bandwidths(R)
69+
70+
@assert size(R) == (n,n)
71+
@assert l == 0 && u ≥ 2
72+
# Tridiagonal bands can be resized
73+
@assert length(Xdl)+1 == length(Xd) == length(Xdu)+1 == length(Ydl)+1 == length(Yd) == length(Ydu)+1 == n
74+
75+
76+
k = 1
77+
aₖ,bₖ = Xd[k], Xdu[k]
78+
Rₖₖ,Rₖₖ₊₁ = Rdat[u+1,k], Rdat[u,k+1] # R[1,1], R[1,2]
79+
Yd[k] = aₖ/Rₖₖ
80+
Ydu[k] = bₖ - aₖ * Rₖₖ₊₁/Rₖₖ
81+
82+
@inbounds for k = 2:n-1
83+
cₖ₋₁,aₖ,bₖ = Xdl[k-1], Xd[k], Xdu[k]
84+
Ydl[k-1] = cₖ₋₁/Rₖₖ
85+
Yd[k] = aₖ-cₖ₋₁*Rₖₖ₊₁/Rₖₖ
86+
Ydu[k] = cₖ₋₁/Rₖₖ
87+
Rₖₖ,Rₖₖ₊₁,Rₖ₋₁ₖ₊₁,Rₖ₋₁ₖ = Rdat[u+1,k], Rdat[u,k+1],Rdat[u-1,k+1],Rₖₖ₊₁ # R[2,2], R[2,3], R[1,3]
88+
Yd[k] /= Rₖₖ
89+
Ydu[k-1] /= Rₖₖ
90+
Ydu[k] *= Rₖ₋₁ₖ*Rₖₖ₊₁/Rₖₖ - Rₖ₋₁ₖ₊₁
91+
Ydu[k] += bₖ - aₖ * Rₖₖ₊₁ / Rₖₖ
92+
end
93+
94+
k = n
95+
cₖ₋₁,aₖ = Xdl[k-1], Xd[k]
96+
Ydl[k-1] = cₖ₋₁/Rₖₖ
97+
Yd[k] = aₖ-cₖ₋₁*Rₖₖ₊₁/Rₖₖ
98+
Rₖₖ = Rdat[u+1,k] # R[2,2], R[2,3], R[1,3]
99+
Yd[k] /= Rₖₖ
100+
Ydu[k-1] /= Rₖₖ
101+
102+
Y
103+
end
104+
"""
105+
TridiagonalConjugationData(U, X, V, Y)
106+
107+
caches the infinite dimensional Tridiagonal(U*X/V)
108+
in the tridiagonal matrix `Y`
109+
"""
110+
111+
mutable struct TridiagonalConjugationData{T}
112+
const U::AbstractMatrix{T}
113+
const X::AbstractMatrix{T}
114+
const V::AbstractMatrix{T}
115+
116+
const UX::Tridiagonal{T,Vector{T}} # cache Tridiagonal(U*X)
117+
const Y::Tridiagonal{T,Vector{T}} # cache Tridiagonal(U*X/V)
118+
119+
datasize::Int
120+
end
121+
122+
function TridiagonalConjugationData(U, X, V, uplo::Char)
123+
T = promote_type(typeof(inv(V[1, 1])), eltype(U), eltype(C)) # include inv so that we can't get Ints
124+
return BidiagonalConjugationData(U, X, V, Tridiagonal(T[], T[], T[]), Tridiagonal(T[], T[], T[]), 0)
125+
end
126+
127+
copy(data::TridiagonalConjugationData) = TridiagonalConjugationData(copy(data.U), copy(data.X), copy(data.V), copy(data.UX), copy(data.Y), data.datasize)
128+
129+
130+
function resizedata!(data::TridiagonalConjugationData, n)
131+
n ≤ 0 && return data
132+
n = max(v, n)
133+
dv, ev = data.dv, data.ev
134+
if n > length(ev) # Avoid O(n²) growing. Note min(length(dv), length(ev)) == length(ev)
135+
resize!(dv, 2n + 1)
136+
resize!(ev, 2n)
137+
end
138+
139+
140+
end
141+

‎test/test_bidiagonalconjugation.jl‎

Lines changed: 34 additions & 122 deletions
Original file line numberDiff line numberDiff line change
@@ -101,124 +101,6 @@ using LazyArrays: LazyLayout
101101
end
102102

103103

104-
"""
105-
upper_mul_tri_triview(U, X) == Tridiagonal(U*X) where U is Upper triangular BandedMatrix and X is Tridiagonal
106-
"""
107-
function upper_mul_tri_triview(U::BandedMatrix, X::Tridiagonal)
108-
T = promote_type(eltype(U), eltype(X))
109-
n = size(U,1)
110-
upper_mul_tri_triview!(Tridiagonal(Vector{T}(undef, n-1), Vector{T}(undef, n), Vector{T}(undef, n-1)), U, X)
111-
end
112-
113-
function upper_mul_tri_triview!(UX::Tridiagonal, U::BandedMatrix, X::Tridiagonal)
114-
n = size(U,1)
115-
116-
j = 1
117-
Xⱼⱼ, Xⱼ₊₁ⱼ = X.d[1], X.dl[1]
118-
Uⱼⱼ, Uⱼⱼ₊₁, Uⱼⱼ₊₂ = U.data[3,1], U.data[2,2], U.data[1,3] # U[j,j], U[j,j+1], U[j,j+2]
119-
UX.d[1] = Uⱼⱼ*Xⱼⱼ + Uⱼⱼ₊₁*Xⱼ₊₁ⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j]
120-
Xⱼⱼ₊₁, Xⱼⱼ, Xⱼ₊₁ⱼ, Xⱼⱼ₋₁ = X.du[1], X.d[2], X.dl[2], Xⱼ₊₁ⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j]
121-
UX.du[1] = Uⱼⱼ*Xⱼⱼ₊₁ + Uⱼⱼ₊₁*Xⱼⱼ + Uⱼⱼ₊₂*Xⱼ₊₁ⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+1]*X[j+1,j]
122-
123-
@inbounds for j = 2:n-2
124-
Uⱼⱼ, Uⱼⱼ₊₁, Uⱼⱼ₊₂ = U.data[3,j], U.data[2,j+1], U.data[1,j+2] # U[j,j], U[j,j+1], U[j,j+2]
125-
UX.dl[j-1] = Uⱼⱼ*Xⱼⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1]
126-
UX.d[j] = Uⱼⱼ*Xⱼⱼ + Uⱼⱼ₊₁*Xⱼ₊₁ⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j]
127-
Xⱼⱼ₊₁, Xⱼⱼ, Xⱼ₊₁ⱼ, Xⱼⱼ₋₁ = X.du[j], X.d[j+1], X.dl[j+1], Xⱼ₊₁ⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j]
128-
UX.du[j] = Uⱼⱼ*Xⱼⱼ₊₁ + Uⱼⱼ₊₁*Xⱼⱼ + Uⱼⱼ₊₂*Xⱼ₊₁ⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+2]*X[j+2,j+1]
129-
end
130-
131-
j = n-1
132-
Uⱼⱼ, Uⱼⱼ₊₁ = U.data[3,j], U.data[2,j+1] # U[j,j], U[j,j+1]
133-
UX.dl[j-1] = Uⱼⱼ*Xⱼⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1]
134-
UX.d[j] = Uⱼⱼ*Xⱼⱼ + Uⱼⱼ₊₁*Xⱼ₊₁ⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j]
135-
Xⱼⱼ₊₁, Xⱼⱼ, Xⱼⱼ₋₁ = X.du[j], X.d[j+1], Xⱼ₊₁ⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j]
136-
UX.du[j] = Uⱼⱼ*Xⱼⱼ₊₁ + Uⱼⱼ₊₁*Xⱼⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+2]*X[j+2,j+1]
137-
138-
j = n
139-
Uⱼⱼ = U.data[3,j] # U[j,j]
140-
UX.dl[j-1] = Uⱼⱼ*Xⱼⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1]
141-
UX.d[j] = Uⱼⱼ*Xⱼⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j]
142-
143-
UX
144-
end
145-
146-
147-
# X*R^{-1} = X*[1/R₁₁ -R₁₂/(R₁₁R₂₂) -R₁₃/R₂₂ …
148-
# 0 1/R₂₂ -R₂₃/R₃₃
149-
# 1/R₃₃
150-
151-
tri_mul_invupper_triview(X::Tridiagonal, R::BandedMatrix) = tri_mul_invupper_triview!(similar(X, promote_type(eltype(X), eltype(R))), X, R)
152-
153-
function tri_mul_invupper_triview!(Y, X, R)
154-
n = size(X,1)
155-
k = 1
156-
Xₖₖ,Xₖₖ₊₁ = X.d[k], X.du[k]
157-
Rₖₖ,Rₖₖ₊₁ = R.data[3,k], R.data[2,k+1] # R[1,1], R[1,2]
158-
Y.d[k] = Xₖₖ/Rₖₖ
159-
Y.du[k] = Xₖₖ₊₁ - Xₖₖ * Rₖₖ₊₁/Rₖₖ
160-
161-
@inbounds for k = 2:n-1
162-
Xₖₖ₋₁,Xₖₖ,Xₖₖ₊₁ = X.dl[k-1], X.d[k], X.du[k]
163-
Y.dl[k-1] = Xₖₖ₋₁/Rₖₖ
164-
Y.d[k] = Xₖₖ-Xₖₖ₋₁*Rₖₖ₊₁/Rₖₖ
165-
Y.du[k] = Xₖₖ₋₁/Rₖₖ
166-
Rₖₖ,Rₖₖ₊₁,Rₖ₋₁ₖ₊₁,Rₖ₋₁ₖ = R.data[3,k], R.data[2,k+1],R.data[1,k+1],Rₖₖ₊₁ # R[2,2], R[2,3], R[1,3]
167-
Y.d[k] /= Rₖₖ
168-
Y.du[k-1] /= Rₖₖ
169-
Y.du[k] *= Rₖ₋₁ₖ*Rₖₖ₊₁/Rₖₖ - Rₖ₋₁ₖ₊₁
170-
Y.du[k] += Xₖₖ₊₁ - Xₖₖ * Rₖₖ₊₁ / Rₖₖ
171-
end
172-
173-
k = n
174-
Xₖₖ₋₁,Xₖₖ = X.dl[k-1], X.d[k]
175-
Y.dl[k-1] = Xₖₖ₋₁/Rₖₖ
176-
Y.d[k] = Xₖₖ-Xₖₖ₋₁*Rₖₖ₊₁/Rₖₖ
177-
Rₖₖ = R.data[3,k] # R[2,2], R[2,3], R[1,3]
178-
Y.d[k] /= Rₖₖ
179-
Y.du[k-1] /= Rₖₖ
180-
181-
Y
182-
end
183-
"""
184-
TridiagonalConjugationData(U, X, V, Y)
185-
186-
caches the infinite dimensional Tridiagonal(U*X/V)
187-
in the tridiagonal matrix `Y`
188-
"""
189-
190-
mutable struct TridiagonalConjugationData{T}
191-
const U::AbstractMatrix{T}
192-
const X::AbstractMatrix{T}
193-
const V::AbstractMatrix{T}
194-
195-
const UX::Tridiagonal{T,Vector{T}} # cache Tridiagonal(U*X)
196-
const Y::Tridiagonal{T,Vector{T}} # cache Tridiagonal(U*X/V)
197-
198-
datasize::Int
199-
end
200-
201-
function TridiagonalConjugationData(U, X, V, uplo::Char)
202-
T = promote_type(typeof(inv(V[1, 1])), eltype(U), eltype(C)) # include inv so that we can't get Ints
203-
return BidiagonalConjugationData(U, X, V, Tridiagonal(T[], T[], T[]), Tridiagonal(T[], T[], T[]), 0)
204-
end
205-
206-
copy(data::TridiagonalConjugationData) = TridiagonalConjugationData(copy(data.U), copy(data.X), copy(data.V), copy(data.UX), copy(data.Y), data.datasize)
207-
208-
209-
function resizedata!(data::TridiagonalConjugationData, n)
210-
n ≤ 0 && return data
211-
n = max(v, n)
212-
dv, ev = data.dv, data.ev
213-
if n > length(ev) # Avoid O(n²) growing. Note min(length(dv), length(ev)) == length(ev)
214-
resize!(dv, 2n + 1)
215-
resize!(ev, 2n)
216-
end
217-
218-
219-
end
220-
221-
222104
@testset "TridiagonalConjugation" begin
223105
@testset "T -> U" begin
224106
R = BandedMatrices._BandedMatrix(Vcat(-Ones(1,∞)/2,
@@ -228,10 +110,10 @@ end
228110
n = 1000
229111
@time U = V = R[1:n,1:n];
230112
@time X = Tridiagonal(Vector(X_T.dl[1:n-1]), Vector(X_T.d[1:n]), Vector(X_T.du[1:n-1]));
231-
@time UX = upper_mul_tri_triview(U, X)
113+
@time UX = InfiniteLinearAlgebra.upper_mul_tri_triview(U, X)
232114
@test Tridiagonal(U*X) ≈ UX
233115
# U*X*inv(U) only depends on Tridiagonal(U*X)
234-
@time Y = tri_mul_invupper_triview(UX, U)
116+
@time Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U)
235117
@test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y
236118
end
237119
@testset "P -> Ultraspherical(3/2)" begin
@@ -242,10 +124,40 @@ end
242124
n = 1000
243125
@time U = V = R[1:n,1:n]
244126
@time X = Tridiagonal(Vector(X_P.dl[1:n-1]), Vector(X_P.d[1:n]), Vector(X_P.du[1:n-1]))
245-
@time UX = upper_mul_tri_triview(U, X)
127+
@time UX = InfiniteLinearAlgebra.upper_mul_tri_triview(U, X)
128+
@test Tridiagonal(U*X) ≈ UX
129+
# U*X*inv(U) only depends on Tridiagonal(U*X)
130+
@time Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U)
131+
@test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y
132+
end
133+
134+
@testset "Jacobi(1,0) -> Jacobi(2,0)" begin
135+
R = BandedMatrices._BandedMatrix(Vcat(Zeros(1,∞), # extra band since code assumes two bands
136+
(-(0:∞) ./ (2:2:∞))',
137+
((2:∞) ./ (2:2:∞))'), ℵ₀, 0,2)
138+
X_P = LazyBandedMatrices.Tridiagonal((2:∞) ./ (3:2:∞), -1 ./ ((1:2:∞) .* (3:2:∞)), (1:∞) ./ (3:2:∞))
139+
n = 1000
140+
@time U = V = R[1:n,1:n]
141+
@time X = Tridiagonal(Vector(X_P.dl[1:n-1]), Vector(X_P.d[1:n]), Vector(X_P.du[1:n-1]))
142+
@time UX = InfiniteLinearAlgebra.upper_mul_tri_triview(U, X)
143+
@test Tridiagonal(U*X) ≈ UX
144+
# U*X*inv(U) only depends on Tridiagonal(U*X)
145+
@time Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U)
146+
@test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y
147+
end
148+
149+
@testset "Legendre() -> Jacobi(5/2)" begin
150+
R = BandedMatrices._BandedMatrix(Vcat((-3 ./ (3:2:∞))', Zeros(1,∞), (3 ./ (3:2:∞))'), ℵ₀, 0,2) *
151+
BandedMatrices._BandedMatrix(Vcat((-1 ./ (1:2:∞))', Zeros(1,∞), (1 ./ (1:2:∞))'), ℵ₀, 0,2)
152+
X_P = LazyBandedMatrices.Tridiagonal((1:∞) ./ (1:2:∞), Zeros(∞), (1:∞) ./ (3:2:∞))
153+
154+
n = 1000
155+
@time U = V = R[1:n,1:n]
156+
@time X = Tridiagonal(Vector(X_P.dl[1:n-1]), Vector(X_P.d[1:n]), Vector(X_P.du[1:n-1]))
157+
@time UX = InfiniteLinearAlgebra.upper_mul_tri_triview(U, X)
246158
@test Tridiagonal(U*X) ≈ UX
247159
# U*X*inv(U) only depends on Tridiagonal(U*X)
248-
@time Y = tri_mul_invupper_triview(UX, U)
160+
@time Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U)
249161
@test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y
250162
end
251163
end

0 commit comments

Comments
 (0)