diff --git a/Project.toml b/Project.toml index eeb9112..429ac13 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ClassicalOrthogonalPolynomials" uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd" authors = ["Sheehan Olver "] -version = "0.14.2" +version = "0.14.3" [deps] ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" diff --git a/src/choleskyQR.jl b/src/choleskyQR.jl index 362e464..8f6b527 100644 --- a/src/choleskyQR.jl +++ b/src/choleskyQR.jl @@ -1,5 +1,5 @@ """ -Represent an Orthogonal polynomial which has a conversion operator from P, that is, Q = P * inv(U). +Represents orthonormal polynomials defined via a conversion operator from P, that is, Q = P * inv(U). """ struct ConvertedOrthogonalPolynomial{T, WW<:AbstractQuasiVector{T}, XX, UU, PP} <: OrthonormalPolynomial{T} weight::WW @@ -8,15 +8,22 @@ struct ConvertedOrthogonalPolynomial{T, WW<:AbstractQuasiVector{T}, XX, UU, PP} P::PP end -_p0(Q::ConvertedOrthogonalPolynomial) = _p0(Q.P) +_p0(Q::ConvertedOrthogonalPolynomial) = _p0(Q.P)/Q.U[1,1] axes(Q::ConvertedOrthogonalPolynomial) = axes(Q.P) + + +struct ConvertedOPLayout <: AbstractNormalizedOPLayout end MemoryLayout(::Type{<:ConvertedOrthogonalPolynomial}) = ConvertedOPLayout() + + + + jacobimatrix(Q::ConvertedOrthogonalPolynomial) = Q.X orthogonalityweight(Q::ConvertedOrthogonalPolynomial) = Q.weight -# transform to P * U if needed for differentiation, etc. +# transform to P * inv(U) if needed for differentiation, etc. arguments(::ApplyLayout{typeof(*)}, Q::ConvertedOrthogonalPolynomial) = Q.P, ApplyArray(inv, Q.U) OrthogonalPolynomial(w::AbstractQuasiVector) = OrthogonalPolynomial(w, orthogonalpolynomial(singularities(w))) diff --git a/src/normalized.jl b/src/normalized.jl index b0cf762..efc89e8 100644 --- a/src/normalized.jl +++ b/src/normalized.jl @@ -35,7 +35,6 @@ isnormalized(_) = false represents OPs that are of the form P * R where P is another family of OPs and R is upper-triangular. """ abstract type AbstractNormalizedOPLayout <: AbstractOPLayout end -struct ConvertedOPLayout <: AbstractNormalizedOPLayout end struct NormalizedOPLayout{LAY<:AbstractBasisLayout} <: AbstractNormalizedOPLayout end MemoryLayout(::Type{<:Normalized{<:Any, OPs}}) where OPs = NormalizedOPLayout{typeof(MemoryLayout(OPs))}() diff --git a/test/test_choleskyQR.jl b/test/test_choleskyQR.jl index 839a0a7..905c82a 100644 --- a/test/test_choleskyQR.jl +++ b/test/test_choleskyQR.jl @@ -1,5 +1,5 @@ using Test, ClassicalOrthogonalPolynomials, BandedMatrices, LinearAlgebra, LazyArrays, ContinuumArrays, LazyBandedMatrices, InfiniteLinearAlgebra -import ClassicalOrthogonalPolynomials: cholesky_jacobimatrix, qr_jacobimatrix, orthogonalpolynomial +import ClassicalOrthogonalPolynomials: cholesky_jacobimatrix, qr_jacobimatrix, orthogonalpolynomial, _p0 import LazyArrays: AbstractCachedMatrix, resizedata! @testset "CholeskyQR" begin @@ -228,9 +228,8 @@ import LazyArrays: AbstractCachedMatrix, resizedata! @test Q == Q̃ @test Q̃ == Q - @test Q[0.1,1] ≈ 1/sqrt(2) + @test Q[0.1,1] ≈ _p0(Q) ≈ 1/sqrt(2) @test Q[0.1,1:10] ≈ Q̃[0.1,1:10] - # AWESOME, thanks TSGUT!! @test Q[0.1,10_000] ≈ Q̃[0.1,10_000] R = P \ Q @@ -247,6 +246,9 @@ import LazyArrays: AbstractCachedMatrix, resizedata! @testset "Chebyshev" begin U = ChebyshevU() Q = orthogonalpolynomial(x -> (1+x^2)*sqrt(1-x^2), U) + x = axes(U,1) + + @test Q[0.1,1] ≈ _p0(Q) ≈ 1/sqrt(sum(expand(Weighted(U),x -> (1+x^2)*sqrt(1-x^2)))) @test bandwidths(Q\U) == (0,2) Q̃ = OrthogonalPolynomial(x -> (1+x^2)*sqrt(1-x^2), U) diff --git a/test/test_gram.jl b/test/test_gram.jl new file mode 100644 index 0000000..82039da --- /dev/null +++ b/test/test_gram.jl @@ -0,0 +1,10 @@ +using ClassicalOrthogonalPolynomials, FastTransforms + +P = Legendre() +x = axes(P,1) +w = @.(1-x^2) + +μ = P'w +X = jacobimatrix(P) +n = 20 +@test GramMatrix(μ[1:2n], X[1:2n,1:2n]) ≈ (P' * (w .* P))[1:n,1:n] \ No newline at end of file