Skip to content

Commit 52bb87b

Browse files
TSGutTSGut
authored
Merge branch 'main' into main
2 parents ad552d7 + d7b2997 commit 52bb87b

File tree

4 files changed

+208
-0
lines changed

4 files changed

+208
-0
lines changed

‎src/ClassicalOrthogonalPolynomials.jl

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -324,5 +324,9 @@ include("classical/laguerre.jl")
324324
include("classical/fourier.jl")
325325
include("stieltjes.jl")
326326
include("roots.jl")
327+
include("decompOPs.jl")
328+
329+
export cholesky_jacobimatrix
330+
327331

328332
end # module

‎src/decompOPs.jl

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
# This currently takes the weight multiplication operator as input.
2+
# I will probably change this to take the weight function instead.
3+
function cholesky_jacobimatrix(w, P::OrthogonalPolynomial)
4+
W = Symmetric(P \ (w.(axes(P,1)) .* P)) # Compute weight multiplication via Clenshaw
5+
bands = CholeskyJacobiBands(W, P) # the cached array only needs to store two bands bc of symmetry
6+
return SymTridiagonal(bands[1,:],bands[2,:])
7+
end
8+
9+
# The generated Jacobi operators are symmetric tridiagonal, so we store their data as two bands
10+
mutable struct CholeskyJacobiBands{T} <: AbstractCachedMatrix{T}
11+
data::Matrix{T}
12+
U::UpperTriangular
13+
X::SymTridiagonal{T}
14+
datasize::Int
15+
array
16+
end
17+
18+
# Computes the initial data for the Jacobi operator bands
19+
function CholeskyJacobiBands(W::Symmetric{T}, P::OrthogonalPolynomial) where T
20+
@assert P isa Normalized
21+
U = cholesky(W).U
22+
X = jacobimatrix(P)
23+
dat = zeros(T,2,10)
24+
for k in 1:10
25+
dat[1,k] = (U * (X * (U \ [zeros(k-1); 1; zeros(∞)])))[k]
26+
dat[2,k] = (U * (X * (U \ [zeros(k); 1; zeros(∞)])))[k]
27+
end
28+
return CholeskyJacobiBands{T}(dat, U, X, 10, dat)
29+
end
30+
31+
size(::CholeskyJacobiBands) = (2,ℵ₀) # Stored as two infinite bands
32+
33+
# Resize and filling functions for cached implementation
34+
function resizedata!(K::CholeskyJacobiBands, nm::Integer)
35+
νμ = K.datasize
36+
if nm > νμ
37+
olddata = copy(K.data)
38+
K.data = similar(K.data, 2, nm)
39+
K.data[axes(olddata)...] = olddata
40+
inds = νμ:nm
41+
cache_filldata!(K, inds)
42+
K.datasize = nm
43+
end
44+
K
45+
end
46+
function cache_filldata!(J::CholeskyJacobiBands, inds::UnitRange{Int})
47+
for k in inds
48+
J.data[1,k] = (J.U * (J.X * (J.U \ [zeros(k-1); 1; zeros(∞)])))[k]
49+
J.data[2,k] = (J.U * (J.X * (J.U \ [zeros(k); 1; zeros(∞)])))[k]
50+
end
51+
end
52+
53+
function getindex(K::CholeskyJacobiBands, k::Integer, j::Integer)
54+
resizedata!(K, max(k,j))
55+
K.data[k, j]
56+
end
57+
function getindex(K::CholeskyJacobiBands, kr::Integer, jr::UnitRange{Int})
58+
resizedata!(K, maximum(jr))
59+
K.data[kr, jr]
60+
end
61+
function getindex(K::CholeskyJacobiBands, I::Vararg{Int,2})
62+
resizedata!(K,maximum(I))
63+
getindex(K.data,I[1],I[2])
64+
end

‎test/runtests.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -43,6 +43,7 @@ include("test_lanczos.jl")
4343
include("test_interlace.jl")
4444
include("test_stieltjes.jl")
4545
include("test_roots.jl")
46+
include("test_decompOPs.jl")
4647

4748
@testset "Auto-diff" begin
4849
U = Ultraspherical(1)

‎test/test_decompOPs.jl

Lines changed: 139 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,139 @@
1+
using Test, ClassicalOrthogonalPolynomials, BandedMatrices, LinearAlgebra, LazyArrays, ContinuumArrays, LazyBandedMatrices
2+
import ClassicalOrthogonalPolynomials: CholeskyJacobiBands
3+
import LazyArrays: AbstractCachedMatrix
4+
import LazyBandedMatrices: SymTridiagonal
5+
6+
@testset "Basic properties" begin
7+
# basis
8+
P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
9+
x = axes(P,1)
10+
J = jacobimatrix(P)
11+
# example weight
12+
w(x) = (1 - x^2)
13+
W = Symmetric(P \ (w.(x) .* P))
14+
# banded cholesky for symmetric-tagged W
15+
@test cholesky(W).U isa UpperTriangular
16+
# compute Jacobi matrix via cholesky
17+
Jchol = cholesky_jacobimatrix(w,P)
18+
@test Jchol isa SymTridiagonal
19+
# CholeskyJacobiBands object
20+
Cbands = CholeskyJacobiBands(W,P)
21+
@test Cbands isa CholeskyJacobiBands
22+
@test Cbands isa AbstractCachedMatrix
23+
@test getindex(Cbands,1,100) == getindex(Cbands,1,1:100)[100]
24+
end
25+
26+
@testset "Comparison with Lanczos and Classical" begin
27+
@testset "Using Clenshaw for polynomial weights" begin
28+
@testset "w(x) = x^2*(1-x)" begin
29+
P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
30+
x = axes(P,1)
31+
J = jacobimatrix(P)
32+
wf(x) = x^2*(1-x)
33+
# compute Jacobi matrix via cholesky
34+
Jchol = cholesky_jacobimatrix(wf, P)
35+
# compute Jacobi matrix via classical recurrence
36+
Q = Normalized(Jacobi(1,2)[affine(0..1,Inclusion(-1..1)),:])
37+
Jclass = jacobimatrix(Q)
38+
# compute Jacobi matrix via Lanczos
39+
Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
40+
# Comparison with Lanczos
41+
@test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
42+
# Comparison with Classical
43+
@test Jchol[1:500,1:500] ≈ Jclass[1:500,1:500]
44+
end
45+
46+
@testset "w(x) = (1-x^2)" begin
47+
P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
48+
x = axes(P,1)
49+
J = jacobimatrix(P)
50+
wf(x) = (1-x^2)
51+
# compute Jacobi matrix via cholesky
52+
Jchol = cholesky_jacobimatrix(wf, P)
53+
# compute Jacobi matrix via Lanczos
54+
Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
55+
# Comparison with Lanczos
56+
@test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
57+
end
58+
59+
@testset "w(x) = (1-x^4)" begin
60+
P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
61+
x = axes(P,1)
62+
J = jacobimatrix(P)
63+
wf(x) = (1-x^4)
64+
# compute Jacobi matrix via cholesky
65+
Jchol = cholesky_jacobimatrix(wf, P)
66+
# compute Jacobi matrix via Lanczos
67+
Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
68+
# Comparison with Lanczos
69+
@test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
70+
end
71+
72+
@testset "w(x) = (1.014-x^3)" begin
73+
P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
74+
x = axes(P,1)
75+
J = jacobimatrix(P)
76+
wf(x) = 1.014-x^4
77+
# compute Jacobi matrix via cholesky
78+
Jchol = cholesky_jacobimatrix(wf, P)
79+
# compute Jacobi matrix via Lanczos
80+
Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
81+
# Comparison with Lanczos
82+
@test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
83+
end
84+
end
85+
86+
@testset "Using Clenshaw with exponential weights" begin
87+
@testset "w(x) = exp(x)" begin
88+
P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
89+
x = axes(P,1)
90+
J = jacobimatrix(P)
91+
wf(x) = exp(x)
92+
# compute Jacobi matrix via cholesky
93+
Jchol = cholesky_jacobimatrix(wf, P)
94+
# compute Jacobi matrix via Lanczos
95+
Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
96+
# Comparison with Lanczos
97+
@test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
98+
end
99+
100+
@testset "w(x) = (1-x)*exp(x)" begin
101+
P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
102+
x = axes(P,1)
103+
J = jacobimatrix(P)
104+
wf(x) = (1-x)*exp(x)
105+
# compute Jacobi matrix via cholesky
106+
Jchol = cholesky_jacobimatrix(wf, P)
107+
# compute Jacobi matrix via Lanczos
108+
Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
109+
# Comparison with Lanczos
110+
@test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
111+
end
112+
113+
@testset "w(x) = (1-x^2)*exp(x^2)" begin
114+
P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
115+
x = axes(P,1)
116+
J = jacobimatrix(P)
117+
wf(x) = (1-x^2)*exp(x^2)
118+
# compute Jacobi matrix via cholesky
119+
Jchol = cholesky_jacobimatrix(wf, P)
120+
# compute Jacobi matrix via Lanczos
121+
Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
122+
# Comparison with Lanczos
123+
@test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
124+
end
125+
126+
@testset "w(x) = x*(1-x^2)*exp(-x^2)" begin
127+
P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
128+
x = axes(P,1)
129+
J = jacobimatrix(P)
130+
wf(x) = x*(1-x^2)*exp(-x^2)
131+
# compute Jacobi matrix via cholesky
132+
Jchol = cholesky_jacobimatrix(wf, P)
133+
# compute Jacobi matrix via Lanczos
134+
Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
135+
# Comparison with Lanczos
136+
@test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
137+
end
138+
end
139+
end

0 commit comments

Comments
 (0)