fixup tests, allow X'W-W*X by better colsupport

add some moment-based tests

Author: MikaelSlevinsky

src/FastTransforms.jl

@@ -19,9 +19,9 @@ import AbstractFFTs: Plan, ScaledPlan,
                      fftshift, ifftshift, rfft_output_size, brfft_output_size,
                      normalization
 
-import ArrayLayouts: colsupport, LayoutMatrix, MemoryLayout, AbstractBandedLayout
+import ArrayLayouts: rowsupport, colsupport, LayoutMatrix, MemoryLayout, AbstractBandedLayout
 
-import BandedMatrices: bandwidths
+import BandedMatrices: bandwidths, BandedLayout
 
 import FFTW: dct, dct!, idct, idct!, plan_dct!, plan_idct!,
              plan_dct, plan_idct, fftwNumber

src/GramMatrix.jl

@@ -2,7 +2,6 @@ abstract type AbstractGramMatrix{T} <: LayoutMatrix{T} end
 
 @inline issymmetric(G::AbstractGramMatrix) = true
 @inline isposdef(G::AbstractGramMatrix) = true
-@inline colsupport(G::AbstractGramMatrix, j) = colrange(G, j)
 
 struct GramMatrix{T, WT <: AbstractMatrix{T}, XT <: AbstractMatrix{T}} <: AbstractGramMatrix{T}
     W::WT
@@ -53,6 +52,8 @@ GramMatrix(W::WT, X::XT) where {T, WT <: AbstractMatrix{T}, XT <: AbstractMatrix
 @inline getindex(G::GramMatrix, i::Integer, j::Integer) = G.W[i, j]
 @inline bandwidths(G::GramMatrix) = bandwidths(G.W)
 @inline MemoryLayout(G::GramMatrix) = MemoryLayout(G.W)
+@inline rowsupport(G::GramMatrix, j) = rowsupport(MemoryLayout(G), G.W, j)
+@inline colsupport(G::GramMatrix, j) = colsupport(MemoryLayout(G), G.W, j)
 
 """
     GramMatrix(μ::AbstractVector, X::AbstractMatrix)
@@ -282,6 +283,7 @@ end
 @inline size(G::ChebyshevGramMatrix) = (G.n, G.n)
 @inline getindex(G::ChebyshevGramMatrix, i::Integer, j::Integer) = (G.μ[abs(i-j)+1] + G.μ[i+j-1])/2
 @inline bandwidths(G::ChebyshevGramMatrix{T, <: PaddedVector{T}}) where T = (length(G.μ.args[2])-1, length(G.μ.args[2])-1)
+@inline MemoryLayout(G::ChebyshevGramMatrix{T, <: PaddedVector{T}}) where T = BandedLayout()
 
 #
 # 2X'W-W*2X = G*J*G'

test/grammatrixtests.jl

@@ -5,38 +5,43 @@ using FastTransforms, BandedMatrices, LazyArrays, LinearAlgebra, Test
     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
+        W = GramMatrix(Symmetric(R'R), X)
+        F = cholesky(W)
+        @test F.L*F.L' ≈ Symmetric(R'R)
         @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
+        W = GramMatrix(Symmetric(R'R), X)
+        F = cholesky(W)
+        @test F.L*F.L' ≈ Symmetric(R'R)
         @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
+        M = Symmetric((I+X^2+X^4)[1:n, 1:n])
+        W = GramMatrix(M, X[1:n, 1:n])
+        @test bandwidths(W) == (b, b)
+        F = cholesky(W)
+        @test F.L*F.L' ≈ M
 
         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
+        M = Symmetric((I+X^2+X^4)[1:n, 1:n])
+        W = GramMatrix(M, X[1:n, 1:n])
+        @test bandwidths(W) == (b, b)
+        F = cholesky(W)
+        @test F.L*F.L' ≈ M
+
+        for μ in (PaddedVector([T(4)/3;0;-T(4)/15], 2n-1), # w(x) = 1-x^2
+                  PaddedVector([T(26)/15;0;-T(4)/105;0;T(16)/315], 2n-1), # w(x) = 1-x^2+x^4
+                  T(1) ./ (1:2n-1)) # Related to a log weight
+            X = Tridiagonal([T(n)/(2n-1) for n in 1:2n-2], zeros(T, 2n-1), [T(n)/(2n+1) for n in 1:2n-2]) # Legendre X
+            W = GramMatrix(μ, X)
+            X = Tridiagonal(X[1:n, 1:n])
+            G = FastTransforms.compute_skew_generators(W)
+            J = T[0 1; -1 0]
+            @test X'W-W*X ≈ G*J*G'
+        end
     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)
@@ -50,38 +55,43 @@ end
     n = 128
     for T in (Float32, Float64, BigFloat)
         μ = FastTransforms.chebyshevmoments1(T, 2n-1)
-        G = ChebyshevGramMatrix(μ)
-        F = cholesky(G)
-        @test F.L*F.L' ≈ G
+        W = ChebyshevGramMatrix(μ)
+        F = cholesky(W)
+        @test F.L*F.L' ≈ W
         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
+        W = ChebyshevGramMatrix(μ)
+        F = cholesky(W)
+        @test F.L*F.L' ≈ W
         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
+        W = ChebyshevGramMatrix(μ)
+        F = cholesky(W)
+        @test F.L*F.L' ≈ W
 
         μ = FastTransforms.chebyshevabsmoments1(T, 2n-1)
-        G = ChebyshevGramMatrix(μ)
-        F = cholesky(G)
-        @test F.L*F.L' ≈ G
+        W = ChebyshevGramMatrix(μ)
+        F = cholesky(W)
+        @test F.L*F.L' ≈ W
 
         μ = 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
+        W = ChebyshevGramMatrix(μ)
+        @test bandwidths(W) == (4, 4)
+        F = cholesky(W)
+        @test F.L*F.L' ≈ W
         μd = Vector{T}(μ)
-        Gd = ChebyshevGramMatrix(μd)
-        Fd = cholesky(Gd)
+        Wd = ChebyshevGramMatrix(μd)
+        Fd = cholesky(Wd)
         @test F.L ≈ Fd.L
+
+        X = Tridiagonal([T(1); ones(T, n-2)/2], zeros(T, n), ones(T, n-1)/2)
+        G = FastTransforms.compute_skew_generators(W)
+        J = T[0 1; -1 0]
+        @test 2*(X'W-W*X) ≈ G*J*G'
     end
 end