Cholesky R

Author: dlfivefifty

src/banded/tridiagonalconjugation.jl

@@ -31,6 +31,7 @@ end
 
 initiate_upper_mul_tri_triview!(UX, U::UpperTriangular, X) = initiate_upper_mul_tri_triview!(UX, parent(U), X)
 initiate_upper_mul_tri_triview!(UX, U::CachedMatrix, X) = initiate_upper_mul_tri_triview!(UX, U.data, X)
+initiate_upper_mul_tri_triview!(UX, U::AdaptiveCholeskyFactors, X) = initiate_upper_mul_tri_triview!(UX, U.data.data, X)
 
 function initiate_upper_mul_tri_triview!(UX, U::BandedMatrix, X)
     Xdl, Xd, Xdu = X.dl, X.d, X.du
@@ -50,7 +51,9 @@ function initiate_upper_mul_tri_triview!(UX, U::BandedMatrix, X)
 end
 
 # fills in the rows kr of UX
-function main_upper_mul_tri_triview!(UX, U::CachedMatrix, X, kr, kwds...)
+main_upper_mul_tri_triview!(UX, U::UpperTriangular, X, kr, kwds...) = main_upper_mul_tri_triview!(UX, parent(U), X, kr, kwds...)
+
+function main_upper_mul_tri_triview!(UX, U::Union{CachedMatrix,AdaptiveCholeskyFactors}, X, kr, kwds...)
     resizedata!(U, kr[end], kr[end]+2)
     main_upper_mul_tri_triview!(UX, U.data, X, kr, kwds...)
 end
@@ -125,6 +128,13 @@ function initiate_tri_mul_invupper_triview!(Y, X, R::CachedMatrix)
     initiate_tri_mul_invupper_triview!(Y, X, R.data)
 end
 
+function initiate_tri_mul_invupper_triview!(Y, X, R::AdaptiveCholeskyFactors)
+    resizedata!(R, 1, 2)
+    initiate_tri_mul_invupper_triview!(Y, X, R.data.data)
+end
+
+initiate_tri_mul_invupper_triview!(Y, X, R::UpperTriangular) = initiate_tri_mul_invupper_triview!(Y, X, parent(R))
+
 function initiate_tri_mul_invupper_triview!(Y, X, R::BandedMatrix)
     Xdl, Xd, Xdu = X.dl, X.d, X.du
     Ydl, Yd, Ydu = Y.dl, Y.d, Y.du
@@ -144,7 +154,8 @@ end
 
 
 # populate rows kr of X/R. Ydu[k] is wrong until next run.
-function main_tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::CachedMatrix, kr, kwds...)
+main_tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::UpperTriangular, kr, kwds...) = main_tri_mul_invupper_triview!(Y, X, parent(R), kr, kwds...)
+function main_tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::Union{AdaptiveCholeskyFactors,CachedMatrix}, kr, kwds...)
     resizedata!(R, kr[end], kr[end]+1)
     main_tri_mul_invupper_triview!(Y, X, R.data, kr, kwds...)
 end

src/infcholesky.jl

@@ -41,6 +41,10 @@ function partialcholesky!(F::AdaptiveCholeskyFactors{T,<:BandedMatrix}, n::Int)
     F
 end
 
+resizedata!(F::AdaptiveCholeskyFactors, m, n) = partialcholesky!(F, n) # support cache interface
+
+resizedata!(R::UpperTriangular{<:Any,<:AdaptiveCholeskyFactors}, m...) = resizedata!(parent(R), m...)
+
 function getindex(F::AdaptiveCholeskyFactors, k::Int, j::Int)
     partialcholesky!(F, max(k,j))
     F.data.data[k,j]

test/test_bidiagonalconjugation.jl

@@ -156,6 +156,25 @@ end
         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 = InfiniteLinearAlgebra.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)
+        @test data.UX[1:2,:] ≈ UX[1:2,1:100]
+        @test data.Y[1,:] ≈ Y[1,1:100]
+        InfiniteLinearAlgebra.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)
+        @test data.UX[1:4,:] ≈ UX[1:4,1:100]
+        @test data.Y[1:3,:] ≈ Y[1:3,1:100]
+        InfiniteLinearAlgebra.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
 end