add a test for nonzero diagonal jacobi matrices

JuliaApproximation · Nov 4, 2024 · 4947845 · 4947845
1 parent c0532d3
commit 4947845
Show file tree

Hide file tree

Showing 2 changed files with 11 additions and 2 deletions.
diff --git a/src/GramMatrix.jl b/src/GramMatrix.jl
@@ -56,10 +56,10 @@ GramMatrix(W::WT, X::XT) where {T, WT <: AbstractMatrix{T}, XT <: AbstractMatrix
 #
 function compute_skew_generators(W::GramMatrix{T}) where T
     X = W.X
-    m, n = size(W)
+    n = size(W, 1)
     G = zeros(T, n, 2)
     G[n, 1] = one(T)
-    G[:, 2] .= W[n-1, :]*X[n-1, n] - X'W[:, n]
+    G[:, 2] .= W[n-1, :]*X[n-1, n] + W[n, :]*X[n, n] - X'W[:, n]
     return G
 end
 

diff --git a/test/GramMatrixtests.jl b/test/GramMatrixtests.jl
@@ -28,6 +28,15 @@ using FastTransforms, BandedMatrices, LazyArrays, LinearAlgebra, Test
         @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)