From 73ebb102ec1d246f9437c85732954a99ec68e051 Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Thu, 12 Dec 2024 14:35:52 -0600 Subject: [PATCH 1/3] add O(bn) GramMatrix constructor from moments --- Project.toml | 2 +- src/GramMatrix.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index be1349a..a7e938a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "FastTransforms" uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c" -version = "0.16.6" +version = "0.16.7" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" diff --git a/src/GramMatrix.jl b/src/GramMatrix.jl index e7e8a70..8bf6dc4 100644 --- a/src/GramMatrix.jl +++ b/src/GramMatrix.jl @@ -64,7 +64,7 @@ function GramMatrix(μ::AbstractVector{T}, X::XT, p0::T) where {T, XT <: Abstrac n = (N+1)÷2 @assert N == size(X, 1) == size(X, 2) @assert bandwidths(X) == (1, 1) - W = Matrix{T}(undef, N, N) + W = LowerTriangular(Matrix{T}(undef, N, N)) if n > 0 @inbounds for m in 1:N W[m, 1] = p0*μ[m] From 0c5cc86c8feaef0f039387f7b63d13c06ef9a632 Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Fri, 13 Dec 2024 08:46:21 -0600 Subject: [PATCH 2/3] delete tests --- test/GramMatrixtests.jl | 87 ----------------------------------------- 1 file changed, 87 deletions(-) delete mode 100644 test/GramMatrixtests.jl diff --git a/test/GramMatrixtests.jl b/test/GramMatrixtests.jl deleted file mode 100644 index 3ac9d62..0000000 --- a/test/GramMatrixtests.jl +++ /dev/null @@ -1,87 +0,0 @@ -using FastTransforms, BandedMatrices, LazyArrays, LinearAlgebra, Test - -@testset "GramMatrix" begin - n = 128 - for T in (Float32, Float64, BigFloat) - R = plan_leg2cheb(T, n; normcheb=true)*I - X = Tridiagonal([T(n)/(2n-1) for n in 1:n-1], zeros(T, n), [T(n)/(2n+1) for n in 1:n-1]) # Legendre X - W = Symmetric(R'R) - G = GramMatrix(W, X) - F = cholesky(G) - @test F.L*F.L' ≈ W - @test F.U ≈ R - - R = plan_leg2cheb(T, n; normcheb=true, normleg=true)*I - X = SymTridiagonal(zeros(T, n), [sqrt(T(n)^2/(4*n^2-1)) for n in 1:n-1]) # normalized Legendre X - W = Symmetric(R'R) - G = GramMatrix(W, X) - F = cholesky(G) - @test F.L*F.L' ≈ W - @test F.U ≈ R - - b = 4 - X = BandedMatrix(SymTridiagonal(zeros(T, n+b), [sqrt(T(n)^2/(4*n^2-1)) for n in 1:n+b-1])) # normalized Legendre X - W = I+X^2+X^4 - W = Symmetric(W[1:n, 1:n]) - X = BandedMatrix(SymTridiagonal(zeros(T, n), [sqrt(T(n)^2/(4*n^2-1)) for n in 1:n-1])) # normalized Legendre X - G = GramMatrix(W, X) - @test bandwidths(G) == (b, b) - F = cholesky(G) - @test F.L*F.L' ≈ W - - X = BandedMatrix(SymTridiagonal(T[2n-1 for n in 1:n+b], T[-n for n in 1:n+b-1])) # Laguerre X, tests nonzero diagonal - W = I+X^2+X^4 - W = Symmetric(W[1:n, 1:n]) - X = BandedMatrix(SymTridiagonal(T[2n-1 for n in 1:n], T[-n for n in 1:n-1])) # Laguerre X - G = GramMatrix(W, X) - @test bandwidths(G) == (b, b) - F = cholesky(G) - @test F.L*F.L' ≈ W - end - W = reshape([i for i in 1.0:n^2], n, n) - X = reshape([i for i in 1.0:4n^2], 2n, 2n) - @test_throws "different sizes" GramMatrix(W, X) - X = X[1:n, 1:n] - @test_throws "nonsymmetric" GramMatrix(W, X) - @test_throws "nontridiagonal" GramMatrix(Symmetric(W), X) -end - -@testset "ChebyshevGramMatrix" begin - n = 128 - for T in (Float32, Float64, BigFloat) - μ = FastTransforms.chebyshevmoments1(T, 2n-1) - G = ChebyshevGramMatrix(μ) - F = cholesky(G) - @test F.L*F.L' ≈ G - R = plan_cheb2leg(T, n; normleg=true)*I - @test F.U ≈ R - - α, β = (T(0.123), T(0.456)) - μ = FastTransforms.chebyshevjacobimoments1(T, 2n-1, α, β) - G = ChebyshevGramMatrix(μ) - F = cholesky(G) - @test F.L*F.L' ≈ G - R = plan_cheb2jac(T, n, α, β; normjac=true)*I - @test F.U ≈ R - - μ = -FastTransforms.chebyshevlogmoments1(T, 2n-1) - G = ChebyshevGramMatrix(μ) - F = cholesky(G) - @test F.L*F.L' ≈ G - - μ = FastTransforms.chebyshevabsmoments1(T, 2n-1) - G = ChebyshevGramMatrix(μ) - F = cholesky(G) - @test F.L*F.L' ≈ G - - μ = PaddedVector(T(1) ./ [1,2,3,4,5], 2n-1) - G = ChebyshevGramMatrix(μ) - @test bandwidths(G) == (4, 4) - F = cholesky(G) - @test F.L*F.L' ≈ G - μd = Vector{T}(μ) - Gd = ChebyshevGramMatrix(μd) - Fd = cholesky(Gd) - @test F.L ≈ Fd.L - end -end From 5260921ef953b357e7eb9a0fd4b0567b4484ef61 Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Fri, 13 Dec 2024 08:47:06 -0600 Subject: [PATCH 3/3] reinstate tests --- Project.toml | 2 +- test/grammatrixtests.jl | 87 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 88 insertions(+), 1 deletion(-) create mode 100644 test/grammatrixtests.jl diff --git a/Project.toml b/Project.toml index a7e938a..11bc836 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "FastTransforms" uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c" -version = "0.16.7" +version = "0.16.8" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" diff --git a/test/grammatrixtests.jl b/test/grammatrixtests.jl new file mode 100644 index 0000000..3ac9d62 --- /dev/null +++ b/test/grammatrixtests.jl @@ -0,0 +1,87 @@ +using FastTransforms, BandedMatrices, LazyArrays, LinearAlgebra, Test + +@testset "GramMatrix" begin + n = 128 + for T in (Float32, Float64, BigFloat) + R = plan_leg2cheb(T, n; normcheb=true)*I + X = Tridiagonal([T(n)/(2n-1) for n in 1:n-1], zeros(T, n), [T(n)/(2n+1) for n in 1:n-1]) # Legendre X + W = Symmetric(R'R) + G = GramMatrix(W, X) + F = cholesky(G) + @test F.L*F.L' ≈ W + @test F.U ≈ R + + R = plan_leg2cheb(T, n; normcheb=true, normleg=true)*I + X = SymTridiagonal(zeros(T, n), [sqrt(T(n)^2/(4*n^2-1)) for n in 1:n-1]) # normalized Legendre X + W = Symmetric(R'R) + G = GramMatrix(W, X) + F = cholesky(G) + @test F.L*F.L' ≈ W + @test F.U ≈ R + + b = 4 + X = BandedMatrix(SymTridiagonal(zeros(T, n+b), [sqrt(T(n)^2/(4*n^2-1)) for n in 1:n+b-1])) # normalized Legendre X + W = I+X^2+X^4 + W = Symmetric(W[1:n, 1:n]) + X = BandedMatrix(SymTridiagonal(zeros(T, n), [sqrt(T(n)^2/(4*n^2-1)) for n in 1:n-1])) # normalized Legendre X + G = GramMatrix(W, X) + @test bandwidths(G) == (b, b) + F = cholesky(G) + @test F.L*F.L' ≈ W + + X = BandedMatrix(SymTridiagonal(T[2n-1 for n in 1:n+b], T[-n for n in 1:n+b-1])) # Laguerre X, tests nonzero diagonal + W = I+X^2+X^4 + W = Symmetric(W[1:n, 1:n]) + X = BandedMatrix(SymTridiagonal(T[2n-1 for n in 1:n], T[-n for n in 1:n-1])) # Laguerre X + G = GramMatrix(W, X) + @test bandwidths(G) == (b, b) + F = cholesky(G) + @test F.L*F.L' ≈ W + end + W = reshape([i for i in 1.0:n^2], n, n) + X = reshape([i for i in 1.0:4n^2], 2n, 2n) + @test_throws "different sizes" GramMatrix(W, X) + X = X[1:n, 1:n] + @test_throws "nonsymmetric" GramMatrix(W, X) + @test_throws "nontridiagonal" GramMatrix(Symmetric(W), X) +end + +@testset "ChebyshevGramMatrix" begin + n = 128 + for T in (Float32, Float64, BigFloat) + μ = FastTransforms.chebyshevmoments1(T, 2n-1) + G = ChebyshevGramMatrix(μ) + F = cholesky(G) + @test F.L*F.L' ≈ G + R = plan_cheb2leg(T, n; normleg=true)*I + @test F.U ≈ R + + α, β = (T(0.123), T(0.456)) + μ = FastTransforms.chebyshevjacobimoments1(T, 2n-1, α, β) + G = ChebyshevGramMatrix(μ) + F = cholesky(G) + @test F.L*F.L' ≈ G + R = plan_cheb2jac(T, n, α, β; normjac=true)*I + @test F.U ≈ R + + μ = -FastTransforms.chebyshevlogmoments1(T, 2n-1) + G = ChebyshevGramMatrix(μ) + F = cholesky(G) + @test F.L*F.L' ≈ G + + μ = FastTransforms.chebyshevabsmoments1(T, 2n-1) + G = ChebyshevGramMatrix(μ) + F = cholesky(G) + @test F.L*F.L' ≈ G + + μ = PaddedVector(T(1) ./ [1,2,3,4,5], 2n-1) + G = ChebyshevGramMatrix(μ) + @test bandwidths(G) == (4, 4) + F = cholesky(G) + @test F.L*F.L' ≈ G + μd = Vector{T}(μ) + Gd = ChebyshevGramMatrix(μd) + Fd = cholesky(Gd) + @test F.L ≈ Fd.L + end +end