basic Clenshaw tests

Author: TSGut

test/test_decompOPs.jl

@@ -1,6 +1,7 @@
 using Test, ClassicalOrthogonalPolynomials, BandedMatrices, LinearAlgebra, LazyArrays, ContinuumArrays, LazyBandedMatrices
 import ClassicalOrthogonalPolynomials: symmjacobim, CholeskyJacobiBands
 import LazyArrays: AbstractCachedMatrix
+import LazyBandedMatrices: SymTridiagonal
 
 @testset "Basic properties" begin
     @testset "Test the Q&D conversion to BandedMatrix format" begin
@@ -35,12 +36,11 @@ import LazyArrays: AbstractCachedMatrix
         Cbands = CholeskyJacobiBands(w)
         @test Cbands isa CholeskyJacobiBands
         @test Cbands isa AbstractCachedMatrix
-        @test Cbands[1,100] == (Cbands[1,1:100])[100]
-        @test Cbands[:,50] == [Cbands[1,50],Cbands[2,50]]
+        @test getindex(Cbands,1,100) == getindex(Cbands,1,1:100)[100]
     end
 end
 
-@testset "Comparison with Lanczos and Classical" begin
+@testset "Comparison with Lanczos and Classical, non-Clenshaw" begin
     @testset "w(x) = x^2*(1-x)" begin
         P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
         x = axes(P,1)
@@ -106,4 +106,71 @@ end
         # Comparison
         @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
     end
+end
+
+@testset "Comparison with Lanczos and Classical, with Clenshaw, basics" begin
+    @testset "w(x) = x^2*(1-x)" begin
+        P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
+        x = axes(P,1)
+        J = jacobimatrix(P)
+        Jx = symmjacobim(J)
+        wf(x) = x^2*(1-x)
+        # compute Jacobi matrix via cholesky
+        W = P \ (wf.(x) .* P)
+        Jchol = cholesky_jacobimatrix(Symmetric(W))
+        # compute Jacobi matrix via classical recurrence
+        Q = Normalized(Jacobi(1,2)[affine(0..1,Inclusion(-1..1)),:])
+        Jclass = jacobimatrix(Q)
+        # compute Jacobi matrix via Lanczos
+        Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
+        # Comparison with Lanczos
+        @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
+        # Comparison with Classical
+        @test Jchol[1:500,1:500] ≈ Jclass[1:500,1:500]
+    end
+
+    @testset "w(x) = (1-x^2)" begin
+        P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
+        x = axes(P,1)
+        J = jacobimatrix(P)
+        Jx = symmjacobim(J)
+        wf(x) = (1-x^2)
+        # compute Jacobi matrix via cholesky
+        W = P \ (wf.(x) .* P)
+        Jchol = cholesky_jacobimatrix(Symmetric(W))
+        # compute Jacobi matrix via Lanczos
+        Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
+        # Comparison with Lanczos
+        @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
+    end
+
+    @testset "w(x) = (1-x^4)" begin
+        P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
+        x = axes(P,1)
+        J = jacobimatrix(P)
+        Jx = symmjacobim(J)
+        wf(x) = (1-x^4)
+        # compute Jacobi matrix via cholesky
+        W = P \ (wf.(x) .* P)
+        Jchol = cholesky_jacobimatrix(Symmetric(W))
+        # compute Jacobi matrix via Lanczos
+        Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
+        # Comparison with Lanczos
+        @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
+    end
+
+    @testset "w(x) = (1.014-x^3)" begin
+        P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
+        x = axes(P,1)
+        J = jacobimatrix(P)
+        Jx = symmjacobim(J)
+        wf(x) = 1.014-x^4
+        # compute Jacobi matrix via cholesky
+        W = P \ (wf.(x) .* P)
+        Jchol = cholesky_jacobimatrix(Symmetric(W))
+        # compute Jacobi matrix via Lanczos
+        Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
+        # Comparison with Lanczos
+        @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
+    end
 end