BigFloat compatibility for qr_jacobimatrix (#146)

* fix QR-X BigFloat issue

* add test and broken test for Cholesky

* fix test

* fix isa vs ==

* use zeros for length 2

* up version a bit

* use isprimitivetype

Author: TSGut

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.11"
+version = "0.11.1"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/choleskyQR.jl

@@ -66,8 +66,9 @@ end
 
 # Computes the initial data for the Jacobi operator bands
 function CholeskyJacobiData(U::AbstractMatrix{T}, UX) where T
-    dv = Vector{T}(undef,2) # compute a length 2 vector on first go
-    ev = Vector{T}(undef,2)
+    # compute a length 2 vector on first go and circumvent BigFloat issue
+    dv = zeros(T,2) 
+    ev = zeros(T,2)
     dv[1] = UX[1,1]/U[1,1] # this is dot(view(UX,1,1), U[1,1] \ [one(T)])
     dv[2] = -U[1,2]*UX[2,1]/(U[1,1]*U[2,2])+UX[2,2]/U[2,2] # this is dot(view(UX,2,1:2), U[1:2,1:2] \ [zero(T); one(T)])
     ev[1] = -UX[1,1]*U[1,2]/(U[1,1]*U[2,2])+UX[1,2]/U[2,2] # this is dot(view(UX,1,1:2), U[1:2,1:2] \ [zero(T); one(T)])
@@ -148,10 +149,10 @@ end
 function QRJacobiData{:Q,T}(F, P) where T
     b = 3+bandwidths(F.R)[2]÷2
     X = jacobimatrix(P)
-        # we fill 1 entry on the first run
-    dv = zeros(T,2)
+    # we fill 1 entry on the first run and circumvent BigFloat issue
+    dv = zeros(T,2) 
     ev = zeros(T,1)
-        # fill first entry (special case)
+    # fill first entry (special case)
     M = Matrix(X[1:b,1:b])
     resizedata!(F.factors,b,b)
     # special case for first entry double Householder product
@@ -176,8 +177,9 @@ function QRJacobiData{:R,T}(F, P) where T
     U = ApplyArray(*,Diagonal(sign.(view(U,band(0)))),U)  # QR decomposition does not force positive diagonals on R by default
     X = jacobimatrix(P)
     UX = ApplyArray(*,U,X)
-    dv = Vector{T}(undef,2) # compute a length 2 vector on first go
-    ev = Vector{T}(undef,2)
+    # compute a length 2 vector on first go and circumvent BigFloat issue
+    dv = zeros(T,2) 
+    ev = zeros(T,2)
     dv[1] = UX[1,1]/U[1,1] # this is dot(view(UX,1,1), U[1,1] \ [one(T)])
     dv[2] = -U[1,2]*UX[2,1]/(U[1,1]*U[2,2])+UX[2,2]/U[2,2] # this is dot(view(UX,2,1:2), U[1:2,1:2] \ [zero(T); one(T)])
     ev[1] = -UX[1,1]*U[1,2]/(U[1,1]*U[2,2])+UX[1,2]/U[2,2] # this is dot(view(UX,1,1:2), U[1:2,1:2] \ [zero(T); one(T)])
@@ -216,7 +218,12 @@ function _fillqrbanddata!(J::QRJacobiData{:Q,T}, inds::UnitRange{Int}) where T
     resizedata!(J.U.τ,m)
     K, τ, F, dv, ev = J.UX, J.U.τ, J.U.factors, J.dv, J.ev
     D = sign.(view(J.U.R,band(0)).*view(J.U.R,band(0))[2:end])
-    M = Matrix{T}(undef,b+3,b+3)
+    M = zeros(T,b+3,b+3)
+    if isprimitivetype(T)
+        M = Matrix{T}(undef,b+3,b+3) 
+    else
+        M = zeros(T,b+3,b+3)
+    end
     @inbounds for n in jj
         dv[n] = K[1,1] # no sign correction needed on diagonal entry due to cancellation
         # doublehouseholderapply!(K,τ[n+1],view(F,n+2:n+b+2,n+1),w)
@@ -243,4 +250,4 @@ function _fillqrbanddata!(J::QRJacobiData{:R,T}, inds::UnitRange{Int}) where T
         dv[k] = -U[k-1,k]*UX[k,k-1]/(U[k-1,k-1]*U[k,k])+UX[k,k]./U[k,k] # this is dot(view(UX,k,k-1:k), U[k-1:k,k-1:k] \ ek)
         ev[k] = UX[k,k-1]/U[k+1,k+1]*(-U[k-1,k+1]/U[k-1,k-1]+U[k-1,k]*U[k,k+1]/(U[k-1,k-1]*U[k,k]))+UX[k,k]/U[k+1,k+1]*(-U[k,k+1]/U[k,k])+UX[k,k+1]/U[k+1,k+1]  # this is dot(view(UX,k,k-1:k+1), U[k-1:k+1,k-1:k+1] \ ek)
     end
-end
+end

test/test_choleskyQR.jl

@@ -205,6 +205,15 @@ import LazyArrays: AbstractCachedMatrix, resizedata!
             @test JqrQ[1:20,1:20] ≈ F[1:20,1:20]
             @test JqrQ[50:70,50:70] ≈ F[50:70,50:70]
         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]
+        end
     end
 
     @testset "ConvertedOP" begin