fix tests

Author: dlfivefifty

src/choleskyQR.jl

@@ -56,7 +56,7 @@ function cholesky_jacobimatrix(w::AbstractQuasiVector, P)
     W = Symmetric(Q \ ((w ./ w_P) .* Q)) # Compute weight multiplication via Clenshaw
     return cholesky_jacobimatrix(W, jacobimatrix(Q))
 end
-function cholesky_jacobimatrix(W::AbstractMatrix, X)
+function cholesky_jacobimatrix(W::AbstractMatrix, X::AbstractMatrix)
     U = cholesky(W).U
     return SymTridiagonalConjugation(U, X), U
 end
@@ -79,7 +79,7 @@ function qr_jacobimatrix(w::AbstractQuasiVector, P)
     sqrtW = Symmetric(Q \ (sqrt.(w ./ w_P) .* Q)) # Compute weight multiplication via Clenshaw
     return qr_jacobimatrix(sqrtW, jacobimatrix(Q))
 end
-function qr_jacobimatrix(sqrtW::AbstractMatrix{T}, X) where T
+function qr_jacobimatrix(sqrtW::AbstractMatrix{T}, X::AbstractMatrix) where T
     R = qr(sqrtW).R
     return SymTridiagonalConjugation(R, X), R
 end

test/test_choleskyQR.jl

@@ -138,20 +138,16 @@ import LazyArrays: AbstractCachedMatrix, resizedata!
             wf = x-> (1-x)^2
             sqrtwf = x-> (1-x)
             # compute Jacobi matrix via decomp
-            Jchol = cholesky_jacobimatrix(wf, P)
-            JqrQ = qr_jacobimatrix(wf, P)
-            JqrR = qr_jacobimatrix(wf, P, :R)
+            Jchol,_ = cholesky_jacobimatrix(wf, P)
+            JqrR,_ = qr_jacobimatrix(wf, P)
             # use alternative inputs
             sqrtW = (P \ (sqrtwf.(x) .* P))
-            JqrQalt = qr_jacobimatrix(sqrtW, P)
-            JqrRalt = qr_jacobimatrix(sqrtW, P, :R)
+            JqrRalt,_ = qr_jacobimatrix(sqrtW, J)
             # compute Jacobi matrix via Lanczos
             Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(legendre(0..1))))
             # Comparison with Lanczos
             @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
-            @test JqrQ[1:500,1:500] ≈ Jlanc[1:500,1:500]
             @test JqrR[1:500,1:500] ≈ Jlanc[1:500,1:500]
-            @test JqrQalt[1:500,1:500] ≈ Jlanc[1:500,1:500]
             @test JqrRalt[1:500,1:500] ≈ Jlanc[1:500,1:500]
         end
         @testset "QR case, w(x) = (1-x)^4" begin
@@ -161,19 +157,17 @@ import LazyArrays: AbstractCachedMatrix, resizedata!
             wf = x -> (1-x)^4
             sqrtwf = x -> (1-x)^2
             # compute Jacobi matrix via decomp
-            JqrQ = qr_jacobimatrix(wf, P)
-            JqrR = qr_jacobimatrix(wf, P, :R)
+            JqrR,R = qr_jacobimatrix(wf, P)
             # use alternative inputs
             sqrtW = (P \ (sqrtwf.(x) .* P))
-            JqrQalt = qr_jacobimatrix(sqrtW, P)
-            JqrRalt = qr_jacobimatrix(sqrtW, P, :R)
+            JqrRalt,Ralt = qr_jacobimatrix(sqrtW, J)
             # compute Jacobi matrix via Lanczos
             Jclass = jacobimatrix(Normalized(jacobi(4,0,0..1)))
             # Comparison with Lanczos
-            @test JqrQ[1:10,1:10] ≈ Jclass[1:10,1:10]
-            @test JqrR[1:10,1:10] ≈ Jclass[1:10,1:10]
-            @test JqrQalt[1:10,1:10] ≈ Jclass[1:10,1:10]
-            @test JqrRalt[1:10,1:10] ≈ Jclass[1:10,1:10]
+            S = Diagonal(sign.(diag(R[1:10,1:10])))
+            @test JqrR[1:10,1:10] ≈ S*Jclass[1:10,1:10]*S
+            S = Diagonal(sign.(diag(Ralt[1:10,1:10])))
+            @test JqrRalt[1:10,1:10] ≈ S*Jclass[1:10,1:10]*S
         end
         @testset "QR case, w(x) = (x)^2*(1-x)^4" begin
             P = Normalized(legendre(0..1))
@@ -182,35 +176,24 @@ import LazyArrays: AbstractCachedMatrix, resizedata!
             wf = x-> (x)^2*(1-x)^4
             sqrtwf = x-> (x)*(1-x)^2
             # compute Jacobi matrix via decomp
-            JqrQ = qr_jacobimatrix(wf, P)
-            JqrR = qr_jacobimatrix(wf, P, :R)
+            JqrR,R = qr_jacobimatrix(wf, P)
             # use alternative inputs
             sqrtW = (P \ (sqrtwf.(x) .* P))
-            JqrQalt = qr_jacobimatrix(sqrtW, P)
-            JqrRalt = qr_jacobimatrix(sqrtW, P, :R)
+            JqrRalt,_ = qr_jacobimatrix(sqrtW, J)
             # compute Jacobi matrix via Lanczos
             Jclass = jacobimatrix(Normalized(jacobi(4,2,0..1)))
             # Comparison with Lanczos
-            @test JqrQ[1:10,1:10] ≈ Jclass[1:10,1:10]
-            @test JqrR[1:10,1:10] ≈ Jclass[1:10,1:10]
-            @test JqrQalt[1:10,1:10] ≈ Jclass[1:10,1:10]
-            @test JqrRalt[1:10,1:10] ≈ Jclass[1:10,1:10]
-            # test consistency of resizing in succession
-            F = qr_jacobimatrix(wf, P);
-            resizedata!(JqrQ.dv,70)
-            resizedata!(JqrQ.ev,70)
-            @test JqrQ[1:5,1:5] ≈ F[1:5,1:5]
-            @test JqrQ[1:20,1:20] ≈ F[1:20,1:20]
-            @test JqrQ[50:70,50:70] ≈ F[50:70,50:70]
+
+            S = Diagonal(sign.(diag(R[1:10,1:10])))
+            @test JqrR[1:10,1:10] ≈ S*Jclass[1:10,1:10]*S
+            @test JqrRalt[1:10,1:10] ≈ S*Jclass[1:10,1:10]*S
         end
         @testset "BigFloat returns correct values" begin
             t = BigFloat("1.1")
             P = Normalized(legendre(big(0)..big(1)))
             X = jacobimatrix(P)
-            Xq = qr_jacobimatrix(t*I-X, P, :Q)
-            Xr = qr_jacobimatrix(t*I-X, P, :R)
-            @test Xq[1:20,1:20] ≈ Xr[1:20,1:20]
-            @test_broken Xq[1:20,1:20] ≈ cholesky_jacobimatrix(Symmetric((t*I-X)^2), P)[1:20,1:20]
+            Xr,R = qr_jacobimatrix(t*I-X, X)
+            @test_broken Xr[1:20,1:20] ≈ cholesky_jacobimatrix(Symmetric((t*I-X)^2), X)[1:20,1:20]
         end
     end