1
1
using Test, ClassicalOrthogonalPolynomials, BandedMatrices, LinearAlgebra, LazyArrays, ContinuumArrays, LazyBandedMatrices
2
- import ClassicalOrthogonalPolynomials: symmjacobim
2
+ import ClassicalOrthogonalPolynomials: symmjacobim, CholeskyJacobiBands
3
+ import LazyArrays: AbstractCachedMatrix
3
4
4
- @testset begin " Test the Q&D conversion to BandedMatrix format"
5
- # Legendre
6
- P = Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])
7
- x = axes (P,1 )
8
- J = jacobimatrix (P)
9
- Jx = symmjacobim (J)
10
- @test J[1 : 100 ,1 : 100 ] == Jx[1 : 100 ,1 : 100 ]
11
- # Jacobi
12
- P = Normalized (Jacobi (1 ,2 )[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])
13
- x = axes (P,1 )
14
- J = jacobimatrix (P)
15
- Jx = symmjacobim (J)
16
- @test J[1 : 100 ,1 : 100 ] == Jx[1 : 100 ,1 : 100 ]
17
- end
5
+ @testset " Basic properties" begin
6
+ @testset " Test the Q&D conversion to BandedMatrix format" begin
7
+ # Legendre
8
+ P = Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])
9
+ x = axes (P,1 )
10
+ J = jacobimatrix (P)
11
+ Jx = symmjacobim (J)
12
+ @test J[1 : 100 ,1 : 100 ] == Jx[1 : 100 ,1 : 100 ]
13
+ # Jacobi
14
+ P = Normalized (Jacobi (1 ,2 )[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])
15
+ x = axes (P,1 )
16
+ J = jacobimatrix (P)
17
+ Jx = symmjacobim (J)
18
+ @test J[1 : 100 ,1 : 100 ] == Jx[1 : 100 ,1 : 100 ]
19
+ end
18
20
19
- @testset begin " Basic properties of Cholesky based Jacobi operators"
20
- P = Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])
21
- x = axes (P,1 )
22
- J = jacobimatrix (P)
23
- Jx = symmjacobim (J)
24
- w = (I - Jx^ 2 )
25
- # banded cholesky for symmetric-tagged W
26
- @test cholesky (w). U isa UpperTriangular
27
- # compute Jacobi matrix via cholesky
28
- Jchol = cholesky_jacobimatrix (w)
29
- @test Jchol isa LazyBandedMatrices. SymTridiagonal
21
+ @testset " Basic types and getindex variants" begin
22
+ # basis
23
+ P = Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])
24
+ x = axes (P,1 )
25
+ J = jacobimatrix (P)
26
+ Jx = symmjacobim (J)
27
+ # example weight
28
+ w = (I - Jx^ 2 )
29
+ # banded cholesky for symmetric-tagged W
30
+ @test cholesky (w). U isa UpperTriangular
31
+ # compute Jacobi matrix via cholesky
32
+ Jchol = cholesky_jacobimatrix (w)
33
+ @test Jchol isa LazyBandedMatrices. SymTridiagonal
34
+ # CholeskyJacobiBands object
35
+ Cbands = CholeskyJacobiBands (w)
36
+ @test Cbands isa CholeskyJacobiBands
37
+ @test Cbands isa AbstractCachedMatrix
38
+ @test Cbands[1 ,100 ] == (Cbands[1 ,1 : 100 ])[100 ]
39
+ @test Cbands[:,50 ] == [Cbands[1 ,50 ],Cbands[2 ,50 ]]
40
+ end
30
41
end
31
42
32
- @testset begin " Comparison with Lanczos and Classical for w(x) = x^2*(1-x)"
33
- P = Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])
34
- x = axes (P,1 )
35
- J = jacobimatrix (P)
36
- Jx = symmjacobim (J)
37
- w = (Jx^ 2 - Jx^ 3 )
38
- # compute Jacobi matrix via cholesky
39
- Jchol = cholesky_jacobimatrix (w)
40
- # compute Jacobi matrix via classical recurrence
41
- Q = Normalized (Jacobi (1 ,2 )[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])
42
- Jclass = jacobimatrix (Q)
43
- # compute Jacobi matrix via Lanczos
44
- wf = x.^ 2 .* (1 .- x)
45
- Jlanc = jacobimatrix (LanczosPolynomial (@. (wf),Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])))
46
- # Comparison with Lanczos
47
- @test Jchol[1 : 500 ,1 : 500 ] ≈ Jlanc[1 : 500 ,1 : 500 ]
48
- # Comparison with Classical
49
- @test Jchol[1 : 500 ,1 : 500 ] ≈ Jclass[1 : 500 ,1 : 500 ]
50
- end
43
+ @testset " Comparison with Lanczos and Classical" begin
44
+ @testset " w(x) = x^2*(1-x)" begin
45
+ P = Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])
46
+ x = axes (P,1 )
47
+ J = jacobimatrix (P)
48
+ Jx = symmjacobim (J)
49
+ w = (Jx^ 2 - Jx^ 3 )
50
+ # compute Jacobi matrix via cholesky
51
+ Jchol = cholesky_jacobimatrix (w)
52
+ # compute Jacobi matrix via classical recurrence
53
+ Q = Normalized (Jacobi (1 ,2 )[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])
54
+ Jclass = jacobimatrix (Q)
55
+ # compute Jacobi matrix via Lanczos
56
+ wf = x.^ 2 .* (1 .- x)
57
+ Jlanc = jacobimatrix (LanczosPolynomial (@. (wf),Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])))
58
+ # Comparison with Lanczos
59
+ @test Jchol[1 : 500 ,1 : 500 ] ≈ Jlanc[1 : 500 ,1 : 500 ]
60
+ # Comparison with Classical
61
+ @test Jchol[1 : 500 ,1 : 500 ] ≈ Jclass[1 : 500 ,1 : 500 ]
62
+ end
51
63
52
- @testset begin " Comparison with Lanczos for w(x) = (1-x^2)"
53
- P = Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])
54
- x = axes (P,1 )
55
- J = jacobimatrix (P)
56
- Jx = symmjacobim (J)
57
- w = (I - Jx^ 2 )
58
- # compute Jacobi matrix via cholesky
59
- Jchol = cholesky_jacobimatrix (w)
60
- # compute Jacobi matrix via Lanczos
61
- wf = (1 .- x.^ 2 )
62
- Jlanc = jacobimatrix (LanczosPolynomial (@. (wf),Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])))
63
- # Comparison
64
- @test Jchol[1 : 500 ,1 : 500 ] ≈ Jlanc[1 : 500 ,1 : 500 ]
65
- end
64
+ @testset " w(x) = (1-x^2)" begin
65
+ P = Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])
66
+ x = axes (P,1 )
67
+ J = jacobimatrix (P)
68
+ Jx = symmjacobim (J)
69
+ w = (I - Jx^ 2 )
70
+ # compute Jacobi matrix via cholesky
71
+ Jchol = cholesky_jacobimatrix (w)
72
+ # compute Jacobi matrix via Lanczos
73
+ wf = (1 .- x.^ 2 )
74
+ Jlanc = jacobimatrix (LanczosPolynomial (@. (wf),Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])))
75
+ # Comparison
76
+ @test Jchol[1 : 500 ,1 : 500 ] ≈ Jlanc[1 : 500 ,1 : 500 ]
77
+ end
66
78
67
- @testset begin " Comparison with Lanczos for w(x) = (1-x^4)"
68
- P = Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])
69
- x = axes (P,1 )
70
- J = jacobimatrix (P)
71
- Jx = symmjacobim (J)
72
- w = (I - Jx^ 4 )
73
- # compute Jacobi matrix via cholesky
74
- Jchol = cholesky_jacobimatrix (w)
75
- # compute Jacobi matrix via Lanczos
76
- wf = (1 .- x.^ 4 )
77
- Jlanc = jacobimatrix (LanczosPolynomial (@. (wf),Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])))
78
- # Comparison
79
- @test Jchol[1 : 500 ,1 : 500 ] ≈ Jlanc[1 : 500 ,1 : 500 ]
80
- end
79
+ @testset " w(x) = (1-x^4)" begin
80
+ P = Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])
81
+ x = axes (P,1 )
82
+ J = jacobimatrix (P)
83
+ Jx = symmjacobim (J)
84
+ w = (I - Jx^ 4 )
85
+ # compute Jacobi matrix via cholesky
86
+ Jchol = cholesky_jacobimatrix (w)
87
+ # compute Jacobi matrix via Lanczos
88
+ wf = (1 .- x.^ 4 )
89
+ Jlanc = jacobimatrix (LanczosPolynomial (@. (wf),Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])))
90
+ # Comparison
91
+ @test Jchol[1 : 500 ,1 : 500 ] ≈ Jlanc[1 : 500 ,1 : 500 ]
92
+ end
81
93
82
- @testset begin " Comparison with Lanczos for w(x) = (1.014-x^3)"
83
- P = Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])
84
- x = axes (P,1 )
85
- J = jacobimatrix (P)
86
- Jx = symmjacobim (J)
87
- t = 1.014
88
- w = (t* I - Jx^ 3 )
89
- # compute Jacobi matrix via cholesky
90
- Jchol = cholesky_jacobimatrix (w)
91
- # compute Jacobi matrix via Lanczos
92
- wf = (t .- x.^ 3 )
93
- Jlanc = jacobimatrix (LanczosPolynomial (@. (wf),Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])))
94
- # Comparison
95
- @test Jchol[1 : 500 ,1 : 500 ] ≈ Jlanc[1 : 500 ,1 : 500 ]
94
+ @testset " w(x) = (1.014-x^3)" begin
95
+ P = Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])
96
+ x = axes (P,1 )
97
+ J = jacobimatrix (P)
98
+ Jx = symmjacobim (J)
99
+ t = 1.014
100
+ w = (t* I - Jx^ 3 )
101
+ # compute Jacobi matrix via cholesky
102
+ Jchol = cholesky_jacobimatrix (w)
103
+ # compute Jacobi matrix via Lanczos
104
+ wf = (t .- x.^ 3 )
105
+ Jlanc = jacobimatrix (LanczosPolynomial (@. (wf),Normalized (Legendre ()[affine (0 .. 1 ,Inclusion (- 1 .. 1 )),:])))
106
+ # Comparison
107
+ @test Jchol[1 : 500 ,1 : 500 ] ≈ Jlanc[1 : 500 ,1 : 500 ]
108
+ end
96
109
end
0 commit comments