Update test_bidiagonalconjugation.jl

Author: dlfivefifty

test/test_bidiagonalconjugation.jl

@@ -22,7 +22,7 @@ using LazyArrays: LazyLayout
         V2 = brand(∞, 0, 1)
         A2 = LazyBandedMatrices.Bidiagonal(Fill(0.2, ∞), 2.0 ./ (1.0 .+ (1:∞)), 'L') # LinearAlgebra.Bidiagonal not playing nice for this case
         X2 = InfRandBidiagonal('L')
-        U2 = X2 * V2 * ApplyArray(inv, A2) 
+        U2 = X2 * V2 * ApplyArray(inv, A2)
         B2 = BidiagonalConjugation(U2, X2, V2, :L);
 
         for (A, B, uplo) in ((A1, B1, 'U'), (A2, B2, 'L'))
@@ -54,12 +54,12 @@ using LazyArrays: LazyLayout
             @test LazyBandedMatrices.Bidiagonal(B) === LazyBandedMatrices.Bidiagonal(B.dv, B.ev, Symbol(uplo))
             @test B[1:10, 1:10] ≈ A[1:10, 1:10]
             @test B[230, 230] ≈ A[230, 230]
-            @test B[102, 102] ≈ A[102, 102] # make sure we compute intermediate columns correctly when skipping 
+            @test B[102, 102] ≈ A[102, 102] # make sure we compute intermediate columns correctly when skipping
             @test B[band(0)][1:100] == B.dv[1:100]
             if uplo == 'U'
                 @test B[band(1)][1:100] == B.ev[1:100]
-                # @test B[band(-1)][1:100] == zeros(100) # This test requires that we define a 
-                # convert(::Type{BidiagonalConjugationBand{T}}, ::Zeros{V, 1, Tuple{OneToInf{Int}}}) where {T, V} method, 
+                # @test B[band(-1)][1:100] == zeros(100) # This test requires that we define a
+                # convert(::Type{BidiagonalConjugationBand{T}}, ::Zeros{V, 1, Tuple{OneToInf{Int}}}) where {T, V} method,
                 # which we probably don't need beyond this test
             else
                 @test B[band(-1)][1:100] == B.ev[1:100]
@@ -74,7 +74,7 @@ using LazyArrays: LazyLayout
             @test (B*I)[1:100, 1:100] ≈ B[1:100, 1:100]
             # @test (B*Diagonal(1:∞))[1:100, 1:100] ≈ B[1:100, 1:100] * Diagonal(1:100) # Uncomment once https://github.com/JuliaLinearAlgebra/ArrayLayouts.jl/pull/241 is registered
 
-            # Pointwise tests 
+            # Pointwise tests
             for i in 1:10
                 for j in 1:10
                     @test B[i, j] ≈ A[i, j]
@@ -92,7 +92,7 @@ using LazyArrays: LazyLayout
         R0 = BandedMatrices._BandedMatrix(Vcat(-Ones(1,∞)/2,
                                     Zeros(1,∞),
                                     Hcat(Ones(1,1),Ones(1,∞)/2)), ℵ₀, 0,2)
-        
+
         D0 = BandedMatrix(1 => 1:∞)
         R1 = BandedMatrix(0 => 1 ./ (1:∞), 2 => -1 ./ (3:∞))
 
@@ -112,7 +112,7 @@ end
 
 function upper_mul_tri_triview!(UX::Tridiagonal, U::BandedMatrix, X::Tridiagonal)
     n = size(U,1)
-    
+
     j = 1
     Xⱼⱼ, Xⱼ₊₁ⱼ = X.d[1], X.dl[1]
     Uⱼⱼ, Uⱼⱼ₊₁, Uⱼⱼ₊₂ =  U.data[3,1], U.data[2,2],  U.data[1,3] # U[j,j], U[j,j+1], U[j,j+2]
@@ -157,7 +157,7 @@ function tri_mul_invupper_triview!(Y, X, R)
     Rₖₖ,Rₖₖ₊₁ = R.data[3,k], R.data[2,k+1] # R[1,1], R[1,2]
     Y.d[k] = Xₖₖ/Rₖₖ
     Y.du[k] = Xₖₖ₊₁ - Xₖₖ * Rₖₖ₊₁/Rₖₖ
-    
+
     @inbounds for k = 2:n-1
         Xₖₖ₋₁,Xₖₖ,Xₖₖ₊₁ = X.dl[k-1], X.d[k], X.du[k]
         Y.dl[k-1] = Xₖₖ₋₁/Rₖₖ
@@ -180,18 +180,72 @@ function tri_mul_invupper_triview!(Y, X, R)
 
     Y
 end
+"""
+    TridiagonalConjugationData(U, X, V, Y)
+
+caches the infinite dimensional Tridiagonal(U*X/V)
+in the tridiagonal matrix `Y`
+"""
+
+mutable struct TridiagonalConjugationData{T}
+    const U::AbstractMatrix{T}
+    const X::AbstractMatrix{T}
+    const V::AbstractMatrix{T}
+
+    const UX::Tridiagonal{T,Vector{T}} # cache Tridiagonal(U*X)
+    const Y::Tridiagonal{T,Vector{T}} # cache Tridiagonal(U*X/V)
+
+    datasize::Int
+end
+
+function TridiagonalConjugationData(U, X, V, uplo::Char)
+    T = promote_type(typeof(inv(V[1, 1])), eltype(U), eltype(C)) # include inv so that we can't get Ints
+    return BidiagonalConjugationData(U, X, V, Tridiagonal(T[], T[], T[]), Tridiagonal(T[], T[], T[]), 0)
+end
+
+copy(data::TridiagonalConjugationData) = TridiagonalConjugationData(copy(data.U), copy(data.X), copy(data.V), copy(data.UX), copy(data.Y), data.datasize)
+
+
+function resizedata!(data::TridiagonalConjugationData, n)
+    n ≤ 0 && return data
+    n = max(v, n)
+    dv, ev = data.dv, data.ev
+    if n > length(ev) # Avoid O(n²) growing. Note min(length(dv), length(ev)) == length(ev)
+        resize!(dv, 2n + 1)
+        resize!(ev, 2n)
+    end
+
+
+end
 
 
 @testset "TridiagonalConjugation" begin
-    R0 = BandedMatrices._BandedMatrix(Vcat(-Ones(1,∞)/2,
-                                    Zeros(1,∞),
-                                    Hcat(Ones(1,1),Ones(1,∞)/2)), ℵ₀, 0,2)
-    X_T = LazyBandedMatrices.Tridiagonal(Vcat(1.0, Fill(1/2,∞)), Zeros(∞), Fill(1/2,∞))
-
-    n = 1000; @time U = V = R0[1:n,1:n];
-    @time X = Tridiagonal(Vector(X_T.dl[1:n-1]), Vector(X_T.d[1:n]), Vector(X_T.du[1:n-1]));
-    @time UX = upper_mul_tri_triview(U, X)
-    @test Tridiagonal(U*X) ≈  UX
-    # U*X*inv(U) only depends on Tridiagonal(U*X)
-    @test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ tri_mul_invupper_triview(UX, U)
+    @testset "T -> U" begin
+        R = BandedMatrices._BandedMatrix(Vcat(-Ones(1,∞)/2,
+                                        Zeros(1,∞),
+                                        Hcat(Ones(1,1),Ones(1,∞)/2)), ℵ₀, 0,2)
+        X_T = LazyBandedMatrices.Tridiagonal(Vcat(1.0, Fill(1/2,∞)), Zeros(∞), Fill(1/2,∞))
+        n = 1000
+        @time U = V = R[1:n,1:n];
+        @time X = Tridiagonal(Vector(X_T.dl[1:n-1]), Vector(X_T.d[1:n]), Vector(X_T.du[1:n-1]));
+        @time UX = upper_mul_tri_triview(U, X)
+        @test Tridiagonal(U*X) ≈  UX
+        # U*X*inv(U) only depends on Tridiagonal(U*X)
+        @time Y = tri_mul_invupper_triview(UX, U)
+        @test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y
+    end
+    @testset "P -> Ultraspherical(3/2)" begin
+        R = BandedMatrices._BandedMatrix(Vcat((-1 ./ (1:2:∞))',
+                                         Zeros(1,∞),
+                                         (1 ./ (1:2:∞))'), ℵ₀, 0,2)
+        X_P = LazyBandedMatrices.Tridiagonal((1:∞) ./ (1:2:∞), Zeros(∞), (1:∞) ./ (3:2:∞))
+        n = 1000
+        @time U = V = R[1:n,1:n]
+        @time X = Tridiagonal(Vector(X_P.dl[1:n-1]), Vector(X_P.d[1:n]), Vector(X_P.du[1:n-1]))
+        @time UX = upper_mul_tri_triview(U, X)
+        @test Tridiagonal(U*X) ≈  UX
+        # U*X*inv(U) only depends on Tridiagonal(U*X)
+        @time Y = tri_mul_invupper_triview(UX, U)
+        @test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y
+    end
 end