Add Sym/TridiagonalConjugation

Author: dlfivefifty

src/banded/tridiagonalconjugation.jl

@@ -238,7 +238,7 @@ copy(data::TridiagonalConjugationData) = TridiagonalConjugationData(copy(data.U)
 function resizedata!(data::TridiagonalConjugationData, n)
     n ≤ data.datasize && return data
 
-    if n > length(data.UX.d) # Avoid O(n²) growing. Note min(length(dv), length(ev)) == length(ev)
+    if n ≥ length(data.UX.dl) # Avoid O(n²) growing. Note min(length(dv), length(ev)) == length(ev)
         resize!(data.UX.dl, 2n)
         resize!(data.UX.d, 2n + 1)
         resize!(data.UX.du, 2n)
@@ -258,3 +258,35 @@ function resizedata!(data::TridiagonalConjugationData, n)
     data
 end
 
+struct TridiagonalConjugationBand{T} <: LazyVector{T}
+    data::TridiagonalConjugationData{T}
+    diag::Symbol
+end
+
+size(P::TridiagonalConjugationBand) = (ℵ₀,)
+resizedata!(A::TridiagonalConjugationBand, n) = resizedata!(A.data, n)
+
+function _triconj_getindex(C::TridiagonalConjugationBand, I)
+    resizedata!(C, maximum(I)+1)
+    getfield(C.data.Y, C.diag)[I]
+end
+
+getindex(A::TridiagonalConjugationBand, I::Integer) = _triconj_getindex(A, I)
+getindex(A::TridiagonalConjugationBand, I::AbstractVector) = _triconj_getindex(A, I)
+getindex(K::TridiagonalConjugationBand, k::AbstractInfUnitRange{<:Integer}) = view(K, k)
+getindex(K::SubArray{<:Any,1,<:TridiagonalConjugationBand}, k::AbstractInfUnitRange{<:Integer}) = view(K, k)
+
+copy(A::TridiagonalConjugationBand) = A # immutable
+
+
+const TridiagonalConjugation{T} = Tridiagonal{T, TridiagonalConjugationBand{T}}
+const SymTridiagonalConjugation{T} = SymTridiagonal{T, TridiagonalConjugationBand{T}}
+function TridiagonalConjugation(R, X, Y...)
+    data = TridiagonalConjugationData(R, X, Y...)
+    Tridiagonal(TridiagonalConjugationBand(data, :dl), TridiagonalConjugationBand(data, :d), TridiagonalConjugationBand(data, :du))
+end
+
+function SymTridiagonalConjugation(R, X, Y...)
+    data = TridiagonalConjugationData(R, X, Y...)
+    SymTridiagonal(TridiagonalConjugationBand(data, :d), TridiagonalConjugationBand(data, :du))
+end

test/test_bidiagonalconjugation.jl

@@ -1,5 +1,5 @@
 using InfiniteLinearAlgebra, InfiniteRandomArrays, BandedMatrices, LazyArrays, LazyBandedMatrices, InfiniteArrays, ArrayLayouts, Test
-using InfiniteLinearAlgebra: BidiagonalConjugation, OneToInf
+using InfiniteLinearAlgebra: BidiagonalConjugation, SymTridiagonalConjugation, TridiagonalConjugation, OneToInf, TridiagonalConjugationData, resizedata!, TridiagonalConjugationBand
 using ArrayLayouts: supdiagonaldata, subdiagonaldata, diagonaldata
 using LinearAlgebra
 using LazyArrays: LazyLayout
@@ -131,50 +131,62 @@ end
         @time Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U)
         @test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y
 
-        data = InfiniteLinearAlgebra.TridiagonalConjugationData(R, X_T)
+        data = TridiagonalConjugationData(R, X_T)
         @test data.UX[1,:] ≈ UX[1,1:100]
-        InfiniteLinearAlgebra.resizedata!(data, 1)
+        resizedata!(data, 1)
         @test data.UX[1:2,:] ≈ UX[1:2,1:100]
         @test data.Y[1,:] ≈ Y[1,1:100]
-        InfiniteLinearAlgebra.resizedata!(data, 2)
+        resizedata!(data, 2)
         @test data.UX[1:3,:] ≈ UX[1:3,1:100]
         @test data.Y[1:2,:] ≈ Y[1:2,1:100]
-        InfiniteLinearAlgebra.resizedata!(data, 3)
+        resizedata!(data, 3)
         @test data.UX[1:4,:] ≈ UX[1:4,1:100]
         @test data.Y[1:3,:] ≈ Y[1:3,1:100]
-        InfiniteLinearAlgebra.resizedata!(data, 1000)
+        resizedata!(data, 1000)
         @test data.UX[1:999,1:999] ≈ UX[1:999,1:999]
         @test data.Y[1:999,1:999] ≈ Y[1:999,1:999]
 
-        data = InfiniteLinearAlgebra.TridiagonalConjugationData(R, X_T)
-        InfiniteLinearAlgebra.resizedata!(data, 1000)
+        data = TridiagonalConjugationData(R, X_T)
+        resizedata!(data, 1000)
         @test data.UX[1:999,1:999] ≈ UX[1:999,1:999]
         @test data.Y[1:999,1:999] ≈ Y[1:999,1:999]
     end
-    
+
     @testset "Cholesky" begin
         M = Symmetric(_BandedMatrix(Vcat(Hcat(Fill(-1/(2sqrt(2)),1,3), Fill(-1/4,1,∞)), Zeros(1,∞), Hcat([0.5 0.25], Fill(0.5,1,∞))), ∞, 0, 2))
         R = cholesky(M).U
         X_T = LazyBandedMatrices.Tridiagonal(Vcat(1/sqrt(2), Fill(1/2,∞)), Zeros(∞), Vcat(1/sqrt(2), Fill(1/2,∞)))
-        data = InfiniteLinearAlgebra.TridiagonalConjugationData(R, X_T);
+        data = TridiagonalConjugationData(R, X_T);
         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]));
         UX = Tridiagonal(U*X)
         Y = Tridiagonal(UX / U)
         @test data.UX[1,:] ≈ UX[1,1:100]
-        InfiniteLinearAlgebra.resizedata!(data, 1)
+        resizedata!(data, 1)
         @test data.UX[1:2,:] ≈ UX[1:2,1:100]
         @test data.Y[1,:] ≈ Y[1,1:100]
-        InfiniteLinearAlgebra.resizedata!(data, 2)
+        resizedata!(data, 2)
         @test data.UX[1:3,:] ≈ UX[1:3,1:100]
         @test data.Y[1:2,:] ≈ Y[1:2,1:100]
-        InfiniteLinearAlgebra.resizedata!(data, 3)
+        resizedata!(data, 3)
         @test data.UX[1:4,:] ≈ UX[1:4,1:100]
         @test data.Y[1:3,:] ≈ Y[1:3,1:100]
-        InfiniteLinearAlgebra.resizedata!(data, 1000)
+        resizedata!(data, 1000)
         @test data.UX[1:999,1:999] ≈ UX[1:999,1:999]
         @test data.Y[1:999,1:999] ≈ Y[1:999,1:999]
-        
+    end
+
+    @testset "TridiagonalConjugationBand" begin
+        R,X = (_BandedMatrix(Vcat(-Ones(1,∞)/2, Zeros(1,∞), Hcat(Ones(1,1),Ones(1,∞)/2)), ℵ₀, 0,2),
+                LazyBandedMatrices.Tridiagonal(Vcat(1.0, Fill(1/2,∞)), Zeros(∞), Fill(1/2,∞)))
+
+        Y = TridiagonalConjugation(R, X)
+        n = 100_000
+        @test Y[n,n+1] ≈ 1/2
+
+        Y = SymTridiagonalConjugation(R, X)
+        n = 100_000
+        @test Y[n,n+1] ≈ 1/2
     end
 end