From 5c38eff3a194fbc7a7568a221802e9035951432d Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Wed, 29 Jan 2025 07:31:49 +0000 Subject: [PATCH 01/18] Add Chebyshev example --- test/test_bidiagonalconjugation.jl | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl index 510d180..e18eb32 100644 --- a/test/test_bidiagonalconjugation.jl +++ b/test/test_bidiagonalconjugation.jl @@ -4,7 +4,7 @@ using ArrayLayouts: supdiagonaldata, subdiagonaldata, diagonaldata using LinearAlgebra using LazyArrays: LazyLayout -@testset "BidiagonalConjugationData" begin +@testset "BidiagonalConjugation" begin @test InfiniteLinearAlgebra._to_uplo('U') == 'U' @test InfiniteLinearAlgebra._to_uplo('L') == 'L' @test_throws ArgumentError InfiniteLinearAlgebra._to_uplo('a') @@ -87,4 +87,12 @@ using LazyArrays: LazyLayout @test B'[3000:3005, 2993:3006] ≈ A[2993:3006, 3000:3005]' end end + + @testset "Chebyshev" begin + R0 = BandedMatrix(0 => [1; Fill(0.5,∞)], 2 => Fill(-0.5,∞)) + D0 = BandedMatrix(1 => 1:∞) + R1 = BandedMatrix(0 => 1 ./ (1:∞), 2 => -1 ./ (3:∞)) + + BidiagonalConjugation(R1, D0, R0, :U) + end end From 1a59fa866d80d3dbaa37925ba51af7cdfae48f9e Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Wed, 29 Jan 2025 12:03:08 +0000 Subject: [PATCH 02/18] Start TridiagonalConjugation --- test/test_bidiagonalconjugation.jl | 39 ++++++++++++++++++++++++++++-- 1 file changed, 37 insertions(+), 2 deletions(-) diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl index e18eb32..7209c0f 100644 --- a/test/test_bidiagonalconjugation.jl +++ b/test/test_bidiagonalconjugation.jl @@ -89,10 +89,45 @@ using LazyArrays: LazyLayout end @testset "Chebyshev" begin - R0 = BandedMatrix(0 => [1; Fill(0.5,∞)], 2 => Fill(-0.5,∞)) + 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:∞)) - BidiagonalConjugation(R1, D0, R0, :U) + B = BidiagonalConjugation(R0', D0', R1', :L)' end end + + +""" +tridiagonal(A, B) == Tridiagonal(A*B) +""" +function tridiagonalmul(A, B) + T = promote_type(eltype(A), eltype(B)) + UX = Tridiagonal(Vector{T}(undef, n-1), Vector{T}(undef, n), Vector{T}(undef, n-1)) + + for j = 1:n-1 + UX.d[j] = U.data[3,j]*X.d[j] + U.data[2,j]*X.dl[j] + end + UX.d[n] = U.data[3,n]*X.d[n] + + for j = 1:n-1 + UX.dl[j] = U.data[3,j]*X.d[j] + U.data[2,j]*X.dl[j] + 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 U*X + T = Float64 + +end \ No newline at end of file From e64a85451f0170f8598e3be5f59f5493689049c2 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Thu, 30 Jan 2025 17:04:33 +0000 Subject: [PATCH 03/18] Fast Tridiagonal(U*X) --- test/test_bidiagonalconjugation.jl | 90 ++++++++++++++++++++++++++---- 1 file changed, 78 insertions(+), 12 deletions(-) diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl index 7209c0f..646befe 100644 --- a/test/test_bidiagonalconjugation.jl +++ b/test/test_bidiagonalconjugation.jl @@ -102,32 +102,98 @@ end """ -tridiagonal(A, B) == Tridiagonal(A*B) +upper_mul_tri_triview(A, B) == Tridiagonal(A*B) where A is Upper triangular BandedMatrix and B is """ -function tridiagonalmul(A, B) - T = promote_type(eltype(A), eltype(B)) - UX = Tridiagonal(Vector{T}(undef, n-1), Vector{T}(undef, n), Vector{T}(undef, n-1)) +function upper_mul_tri_triview(U, X) + T = promote_type(eltype(U), eltype(X)) + n = size(U,1) + upper_mul_tri_triview!(Tridiagonal(Vector{T}(undef, n-1), Vector{T}(undef, n), Vector{T}(undef, n-1)), U, X) +end + + +# function upper_mul_tri_triview!(UX::Tridiagonal, U::BandedMatrix, X::Tridiagonal) +# n = size(U,1) +# @inbounds for j = 1:n-1 +# UX.d[j] = U.data[3,j]*X.d[j] + U.data[2,j]*X.dl[j] +# end +# UX.d[n] = U.data[3,n]*X.d[n] + +# @inbounds for j = 1:n-1 +# UX.dl[j] = U.data[3,j+1]*X.dl[j] +# end + +# @inbounds for j = 1:n-2 +# UX.du[j] = U.data[3,j]*X.du[j] + U.data[2,j+1]*X.d[j+1] + U.data[1,j+2]*X.dl[j+1] +# end + +# UX.du[n-1] = U.data[3,n-1]*X.du[n-1] + U.data[2,n]*X.d[n] + +# UX +# 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] + UX.d[1] = Uⱼⱼ*Xⱼⱼ + Uⱼⱼ₊₁*Xⱼ₊₁ⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j] + Xⱼⱼ₊₁, Xⱼⱼ, Xⱼ₊₁ⱼ, Xⱼⱼ₋₁ = X.du[1], X.d[2], X.dl[2], Xⱼ₊₁ⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j] + UX.du[1] = Uⱼⱼ*Xⱼⱼ₊₁ + Uⱼⱼ₊₁*Xⱼⱼ + Uⱼⱼ₊₂*Xⱼ₊₁ⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+1]*X[j+1,j] + + @inbounds for j = 2:n-2 + Uⱼⱼ, Uⱼⱼ₊₁, Uⱼⱼ₊₂ = U.data[3,j], U.data[2,j+1], U.data[1,j+2] # U[j,j], U[j,j+1], U[j,j+2] + UX.dl[j-1] = Uⱼⱼ*Xⱼⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1] + UX.d[j] = Uⱼⱼ*Xⱼⱼ + Uⱼⱼ₊₁*Xⱼ₊₁ⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j] + Xⱼⱼ₊₁, Xⱼⱼ, Xⱼ₊₁ⱼ, Xⱼⱼ₋₁ = X.du[j], X.d[j+1], X.dl[j+1], Xⱼ₊₁ⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j] + UX.du[j] = Uⱼⱼ*Xⱼⱼ₊₁ + Uⱼⱼ₊₁*Xⱼⱼ + Uⱼⱼ₊₂*Xⱼ₊₁ⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+2]*X[j+2,j+1] + end + + j = n-1 + Uⱼⱼ, Uⱼⱼ₊₁ = U.data[3,j], U.data[2,j+1] # U[j,j], U[j,j+1] + UX.dl[j-1] = Uⱼⱼ*Xⱼⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1] + UX.d[j] = Uⱼⱼ*Xⱼⱼ + Uⱼⱼ₊₁*Xⱼ₊₁ⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j] + Xⱼⱼ₊₁, Xⱼⱼ, Xⱼⱼ₋₁ = X.du[j], X.d[j+1], Xⱼ₊₁ⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j] + UX.du[j] = Uⱼⱼ*Xⱼⱼ₊₁ + Uⱼⱼ₊₁*Xⱼⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+2]*X[j+2,j+1] + + j = n + Uⱼⱼ = U.data[3,j] # U[j,j] + UX.dl[j-1] = Uⱼⱼ*Xⱼⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1] + UX.d[j] = Uⱼⱼ*Xⱼⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j] + + UX +end + +tri_mul_invupper_triview(UX::Tridiagonal, R::BandedMatrix) = tri_mul_invupper_triview!(similar(UX, promote_type(eltype(UX), eltype(R))), UX, R) - for j = 1:n-1 - UX.d[j] = U.data[3,j]*X.d[j] + U.data[2,j]*X.dl[j] +function tri_mul_invupper_triview!(X, UX, R) + @inbounds for j = 1:n-1 + UX.d[j] = UX.d[j]/ U.data[3,j]*X.d[j] + U.data[2,j]*X.dl[j] end UX.d[n] = U.data[3,n]*X.d[n] - for j = 1:n-1 - UX.dl[j] = U.data[3,j]*X.d[j] + U.data[2,j]*X.dl[j] + @inbounds for j = 1:n-1 + UX.dl[j] = U.data[3,j+1]*X.dl[j] + end + + @inbounds for j = 1:n-2 + UX.du[j] = U.data[3,j]*X.du[j] + U.data[2,j+1]*X.d[j+1] + U.data[1,j+2]*X.dl[j+1] end + + UX.du[n-1] = U.data[3,n-1]*X.du[n-1] + U.data[2,n]*X.d[n] + + UX 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])) + 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 U*X - T = Float64 + @test Tridiagonal(U*X) ≈ upper_mul_tri_triview(U, X) end \ No newline at end of file From 955d096d44a79e48f38dcc3ed51122e622cbe19f Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Thu, 30 Jan 2025 22:26:23 +0000 Subject: [PATCH 04/18] Fast finite tri-conj --- test/test_bidiagonalconjugation.jl | 80 +++++++++++++++--------------- 1 file changed, 39 insertions(+), 41 deletions(-) diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl index 646befe..1b29976 100644 --- a/test/test_bidiagonalconjugation.jl +++ b/test/test_bidiagonalconjugation.jl @@ -102,37 +102,17 @@ end """ -upper_mul_tri_triview(A, B) == Tridiagonal(A*B) where A is Upper triangular BandedMatrix and B is +upper_mul_tri_triview(U, X) == Tridiagonal(U*X) where U is Upper triangular BandedMatrix and X is Tridiagonal """ -function upper_mul_tri_triview(U, X) +function upper_mul_tri_triview(U::BandedMatrix, X::Tridiagonal) T = promote_type(eltype(U), eltype(X)) n = size(U,1) upper_mul_tri_triview!(Tridiagonal(Vector{T}(undef, n-1), Vector{T}(undef, n), Vector{T}(undef, n-1)), U, X) end - -# function upper_mul_tri_triview!(UX::Tridiagonal, U::BandedMatrix, X::Tridiagonal) -# n = size(U,1) -# @inbounds for j = 1:n-1 -# UX.d[j] = U.data[3,j]*X.d[j] + U.data[2,j]*X.dl[j] -# end -# UX.d[n] = U.data[3,n]*X.d[n] - -# @inbounds for j = 1:n-1 -# UX.dl[j] = U.data[3,j+1]*X.dl[j] -# end - -# @inbounds for j = 1:n-2 -# UX.du[j] = U.data[3,j]*X.du[j] + U.data[2,j+1]*X.d[j+1] + U.data[1,j+2]*X.dl[j+1] -# end - -# UX.du[n-1] = U.data[3,n-1]*X.du[n-1] + U.data[2,n]*X.d[n] - -# UX -# 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] @@ -163,25 +143,42 @@ function upper_mul_tri_triview!(UX::Tridiagonal, U::BandedMatrix, X::Tridiagonal UX end -tri_mul_invupper_triview(UX::Tridiagonal, R::BandedMatrix) = tri_mul_invupper_triview!(similar(UX, promote_type(eltype(UX), eltype(R))), UX, R) -function tri_mul_invupper_triview!(X, UX, R) - @inbounds for j = 1:n-1 - UX.d[j] = UX.d[j]/ U.data[3,j]*X.d[j] + U.data[2,j]*X.dl[j] +# X*R^{-1} = X*[1/R₁₁ -R₁₂/(R₁₁R₂₂) -R₁₃/R₂₂ … +# 0 1/R₂₂ -R₂₃/R₃₃ +# 1/R₃₃ + +tri_mul_invupper_triview(X::Tridiagonal, R::BandedMatrix) = tri_mul_invupper_triview!(similar(X, promote_type(eltype(X), eltype(R))), X, R) + +function tri_mul_invupper_triview!(Y, X, R) + n = size(X,1) + k = 1 + Xₖₖ,Xₖₖ₊₁ = X.d[k], X.du[k] + 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ₖₖ + + for k = 2:n-1 + Xₖₖ₋₁,Xₖₖ,Xₖₖ₊₁ = X.dl[k-1], X.d[k], X.du[k] + Y.dl[k-1] = Xₖₖ₋₁/Rₖₖ + Y.d[k] = Xₖₖ-Xₖₖ₋₁*Rₖₖ₊₁/Rₖₖ + Y.du[k] = Xₖₖ₋₁/Rₖₖ + Rₖₖ,Rₖₖ₊₁,Rₖ₋₁ₖ₊₁,Rₖ₋₁ₖ = R.data[3,k], R.data[2,k+1],R.data[1,k+1],Rₖₖ₊₁ # R[2,2], R[2,3], R[1,3] + Y.d[k] /= Rₖₖ + Y.du[k-1] /= Rₖₖ + Y.du[k] *= Rₖ₋₁ₖ*Rₖₖ₊₁/Rₖₖ - Rₖ₋₁ₖ₊₁ + Y.du[k] += Xₖₖ₊₁ - Xₖₖ * Rₖₖ₊₁ / Rₖₖ end - UX.d[n] = U.data[3,n]*X.d[n] - @inbounds for j = 1:n-1 - UX.dl[j] = U.data[3,j+1]*X.dl[j] - end + k = n + Xₖₖ₋₁,Xₖₖ = X.dl[k-1], X.d[k] + Y.dl[k-1] = Xₖₖ₋₁/Rₖₖ + Y.d[k] = Xₖₖ-Xₖₖ₋₁*Rₖₖ₊₁/Rₖₖ + Rₖₖ = R.data[3,k] # R[2,2], R[2,3], R[1,3] + Y.d[k] /= Rₖₖ + Y.du[k-1] /= Rₖₖ - @inbounds for j = 1:n-2 - UX.du[j] = U.data[3,j]*X.du[j] + U.data[2,j+1]*X.d[j+1] + U.data[1,j+2]*X.dl[j+1] - end - - UX.du[n-1] = U.data[3,n-1]*X.du[n-1] + U.data[2,n]*X.d[n] - - UX + Y end @@ -193,7 +190,8 @@ end 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])); - - @test Tridiagonal(U*X) ≈ upper_mul_tri_triview(U, X) - + @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) end \ No newline at end of file From 4389524dd13fd9388140997bf71248cb811b23fb Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Thu, 30 Jan 2025 22:29:56 +0000 Subject: [PATCH 05/18] Update test_bidiagonalconjugation.jl --- test/test_bidiagonalconjugation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl index 1b29976..eda098b 100644 --- a/test/test_bidiagonalconjugation.jl +++ b/test/test_bidiagonalconjugation.jl @@ -158,7 +158,7 @@ function tri_mul_invupper_triview!(Y, X, R) Y.d[k] = Xₖₖ/Rₖₖ Y.du[k] = Xₖₖ₊₁ - Xₖₖ * Rₖₖ₊₁/Rₖₖ - for k = 2:n-1 + @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ₖₖ Y.d[k] = Xₖₖ-Xₖₖ₋₁*Rₖₖ₊₁/Rₖₖ From 0325f286f8009147b1b3c5ed7114602d4a8b63ce Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Fri, 31 Jan 2025 17:02:16 +0000 Subject: [PATCH 06/18] Update test_bidiagonalconjugation.jl --- test/test_bidiagonalconjugation.jl | 92 ++++++++++++++++++++++++------ 1 file changed, 73 insertions(+), 19 deletions(-) diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl index eda098b..6d7cfd2 100644 --- a/test/test_bidiagonalconjugation.jl +++ b/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 \ No newline at end of file From 50d30ade82456bcf8d01d69bddaf50b94fbe91ff Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sat, 1 Feb 2025 09:22:18 +0000 Subject: [PATCH 07/18] Allow general bandwiddths in R --- src/InfiniteLinearAlgebra.jl | 1 + src/banded/tridiagonalconjugation.jl | 141 ++++++++++++++++++++++++ test/test_bidiagonalconjugation.jl | 156 ++++++--------------------- 3 files changed, 176 insertions(+), 122 deletions(-) create mode 100644 src/banded/tridiagonalconjugation.jl diff --git a/src/InfiniteLinearAlgebra.jl b/src/InfiniteLinearAlgebra.jl index d7ab620..d6f7336 100644 --- a/src/InfiniteLinearAlgebra.jl +++ b/src/InfiniteLinearAlgebra.jl @@ -125,5 +125,6 @@ include("infqr.jl") include("inful.jl") include("infcholesky.jl") include("banded/bidiagonalconjugation.jl") +include("banded/tridiagonalconjugation.jl") end # module diff --git a/src/banded/tridiagonalconjugation.jl b/src/banded/tridiagonalconjugation.jl new file mode 100644 index 0000000..8107d5d --- /dev/null +++ b/src/banded/tridiagonalconjugation.jl @@ -0,0 +1,141 @@ +""" +upper_mul_tri_triview(U, X) == Tridiagonal(U*X) where U is Upper triangular BandedMatrix and X is Tridiagonal +""" +function upper_mul_tri_triview(U::BandedMatrix, X::Tridiagonal) + T = promote_type(eltype(U), eltype(X)) + n = size(U,1) + upper_mul_tri_triview!(Tridiagonal(Vector{T}(undef, n-1), Vector{T}(undef, n), Vector{T}(undef, n-1)), U, X) +end + +function upper_mul_tri_triview!(UX::Tridiagonal, U::BandedMatrix, X::Tridiagonal) + n = size(UX,1) + + + Xdl, Xd, Xdu = X.dl, X.d, X.du + UXdl, UXd, UXdu = UX.dl, UX.d, UX.du + Udat = U.data + + l,u = bandwidths(U) + + @assert size(U) == (n,n) + @assert l == 0 && u ≥ 2 + # Tridiagonal bands can be resized + @assert length(Xdl)+1 == length(Xd) == length(Xdu)+1 == length(UXdl)+1 == length(UXd) == length(UXdu)+1 == n + + j = 1 + aⱼ, cⱼ = Xd[1], Xdl[1] + Uⱼⱼ, Uⱼⱼ₊₁, Uⱼⱼ₊₂ = Udat[u+1,1], Udat[u,2], Udat[u-1,3] # U[j,j], U[j,j+1], U[j,j+2] + UXd[1] = Uⱼⱼ*aⱼ + Uⱼⱼ₊₁*cⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j] + bⱼ, aⱼ, cⱼ, cⱼ₋₁ = Xdu[1], Xd[2], Xdl[2], cⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j] + UXdu[1] = Uⱼⱼ*bⱼ + Uⱼⱼ₊₁*aⱼ + Uⱼⱼ₊₂*cⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+1]*X[j+1,j] + + @inbounds for j = 2:n-2 + Uⱼⱼ, Uⱼⱼ₊₁, Uⱼⱼ₊₂ = Udat[u+1,j], Udat[u,j+1], Udat[u-1,j+2] # U[j,j], U[j,j+1], U[j,j+2] + UXdl[j-1] = Uⱼⱼ*cⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1] + UXd[j] = Uⱼⱼ*aⱼ + Uⱼⱼ₊₁*cⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j] + bⱼ, aⱼ, cⱼ, cⱼ₋₁ = Xdu[j], Xd[j+1], Xdl[j+1], cⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j] + UXdu[j] = Uⱼⱼ*bⱼ + Uⱼⱼ₊₁*aⱼ + Uⱼⱼ₊₂*cⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+2]*X[j+2,j+1] + end + + j = n-1 + Uⱼⱼ, Uⱼⱼ₊₁ = Udat[u+1,j], Udat[u,j+1] # U[j,j], U[j,j+1] + UXdl[j-1] = Uⱼⱼ*cⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1] + UXd[j] = Uⱼⱼ*aⱼ + Uⱼⱼ₊₁*cⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j] + bⱼ, aⱼ, cⱼ₋₁ = Xdu[j], Xd[j+1], cⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j] + UXdu[j] = Uⱼⱼ*bⱼ + Uⱼⱼ₊₁*aⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+2]*X[j+2,j+1] + + j = n + Uⱼⱼ = Udat[u+1,j] # U[j,j] + UXdl[j-1] = Uⱼⱼ*cⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1] + UXd[j] = Uⱼⱼ*aⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j] + + UX +end + + +# X*R^{-1} = X*[1/R₁₁ -R₁₂/(R₁₁R₂₂) -R₁₃/R₂₂ … +# 0 1/R₂₂ -R₂₃/R₃₃ +# 1/R₃₃ + +tri_mul_invupper_triview(X::Tridiagonal, R::BandedMatrix) = tri_mul_invupper_triview!(similar(X, promote_type(eltype(X), eltype(R))), X, R) + +function tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::BandedMatrix) + n = size(X,1) + Xdl, Xd, Xdu = X.dl, X.d, X.du + Ydl, Yd, Ydu = Y.dl, Y.d, Y.du + Rdat = R.data + + l,u = bandwidths(R) + + @assert size(R) == (n,n) + @assert l == 0 && u ≥ 2 + # Tridiagonal bands can be resized + @assert length(Xdl)+1 == length(Xd) == length(Xdu)+1 == length(Ydl)+1 == length(Yd) == length(Ydu)+1 == n + + + k = 1 + aₖ,bₖ = Xd[k], Xdu[k] + Rₖₖ,Rₖₖ₊₁ = Rdat[u+1,k], Rdat[u,k+1] # R[1,1], R[1,2] + Yd[k] = aₖ/Rₖₖ + Ydu[k] = bₖ - aₖ * Rₖₖ₊₁/Rₖₖ + + @inbounds for k = 2:n-1 + cₖ₋₁,aₖ,bₖ = Xdl[k-1], Xd[k], Xdu[k] + Ydl[k-1] = cₖ₋₁/Rₖₖ + Yd[k] = aₖ-cₖ₋₁*Rₖₖ₊₁/Rₖₖ + Ydu[k] = cₖ₋₁/Rₖₖ + Rₖₖ,Rₖₖ₊₁,Rₖ₋₁ₖ₊₁,Rₖ₋₁ₖ = Rdat[u+1,k], Rdat[u,k+1],Rdat[u-1,k+1],Rₖₖ₊₁ # R[2,2], R[2,3], R[1,3] + Yd[k] /= Rₖₖ + Ydu[k-1] /= Rₖₖ + Ydu[k] *= Rₖ₋₁ₖ*Rₖₖ₊₁/Rₖₖ - Rₖ₋₁ₖ₊₁ + Ydu[k] += bₖ - aₖ * Rₖₖ₊₁ / Rₖₖ + end + + k = n + cₖ₋₁,aₖ = Xdl[k-1], Xd[k] + Ydl[k-1] = cₖ₋₁/Rₖₖ + Yd[k] = aₖ-cₖ₋₁*Rₖₖ₊₁/Rₖₖ + Rₖₖ = Rdat[u+1,k] # R[2,2], R[2,3], R[1,3] + Yd[k] /= Rₖₖ + Ydu[k-1] /= 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 + diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl index 6d7cfd2..bfbf481 100644 --- a/test/test_bidiagonalconjugation.jl +++ b/test/test_bidiagonalconjugation.jl @@ -101,124 +101,6 @@ using LazyArrays: LazyLayout end -""" -upper_mul_tri_triview(U, X) == Tridiagonal(U*X) where U is Upper triangular BandedMatrix and X is Tridiagonal -""" -function upper_mul_tri_triview(U::BandedMatrix, X::Tridiagonal) - T = promote_type(eltype(U), eltype(X)) - n = size(U,1) - upper_mul_tri_triview!(Tridiagonal(Vector{T}(undef, n-1), Vector{T}(undef, n), Vector{T}(undef, n-1)), U, X) -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] - UX.d[1] = Uⱼⱼ*Xⱼⱼ + Uⱼⱼ₊₁*Xⱼ₊₁ⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j] - Xⱼⱼ₊₁, Xⱼⱼ, Xⱼ₊₁ⱼ, Xⱼⱼ₋₁ = X.du[1], X.d[2], X.dl[2], Xⱼ₊₁ⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j] - UX.du[1] = Uⱼⱼ*Xⱼⱼ₊₁ + Uⱼⱼ₊₁*Xⱼⱼ + Uⱼⱼ₊₂*Xⱼ₊₁ⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+1]*X[j+1,j] - - @inbounds for j = 2:n-2 - Uⱼⱼ, Uⱼⱼ₊₁, Uⱼⱼ₊₂ = U.data[3,j], U.data[2,j+1], U.data[1,j+2] # U[j,j], U[j,j+1], U[j,j+2] - UX.dl[j-1] = Uⱼⱼ*Xⱼⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1] - UX.d[j] = Uⱼⱼ*Xⱼⱼ + Uⱼⱼ₊₁*Xⱼ₊₁ⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j] - Xⱼⱼ₊₁, Xⱼⱼ, Xⱼ₊₁ⱼ, Xⱼⱼ₋₁ = X.du[j], X.d[j+1], X.dl[j+1], Xⱼ₊₁ⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j] - UX.du[j] = Uⱼⱼ*Xⱼⱼ₊₁ + Uⱼⱼ₊₁*Xⱼⱼ + Uⱼⱼ₊₂*Xⱼ₊₁ⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+2]*X[j+2,j+1] - end - - j = n-1 - Uⱼⱼ, Uⱼⱼ₊₁ = U.data[3,j], U.data[2,j+1] # U[j,j], U[j,j+1] - UX.dl[j-1] = Uⱼⱼ*Xⱼⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1] - UX.d[j] = Uⱼⱼ*Xⱼⱼ + Uⱼⱼ₊₁*Xⱼ₊₁ⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j] - Xⱼⱼ₊₁, Xⱼⱼ, Xⱼⱼ₋₁ = X.du[j], X.d[j+1], Xⱼ₊₁ⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j] - UX.du[j] = Uⱼⱼ*Xⱼⱼ₊₁ + Uⱼⱼ₊₁*Xⱼⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+2]*X[j+2,j+1] - - j = n - Uⱼⱼ = U.data[3,j] # U[j,j] - UX.dl[j-1] = Uⱼⱼ*Xⱼⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1] - UX.d[j] = Uⱼⱼ*Xⱼⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j] - - UX -end - - -# X*R^{-1} = X*[1/R₁₁ -R₁₂/(R₁₁R₂₂) -R₁₃/R₂₂ … -# 0 1/R₂₂ -R₂₃/R₃₃ -# 1/R₃₃ - -tri_mul_invupper_triview(X::Tridiagonal, R::BandedMatrix) = tri_mul_invupper_triview!(similar(X, promote_type(eltype(X), eltype(R))), X, R) - -function tri_mul_invupper_triview!(Y, X, R) - n = size(X,1) - k = 1 - Xₖₖ,Xₖₖ₊₁ = X.d[k], X.du[k] - 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ₖₖ - Y.d[k] = Xₖₖ-Xₖₖ₋₁*Rₖₖ₊₁/Rₖₖ - Y.du[k] = Xₖₖ₋₁/Rₖₖ - Rₖₖ,Rₖₖ₊₁,Rₖ₋₁ₖ₊₁,Rₖ₋₁ₖ = R.data[3,k], R.data[2,k+1],R.data[1,k+1],Rₖₖ₊₁ # R[2,2], R[2,3], R[1,3] - Y.d[k] /= Rₖₖ - Y.du[k-1] /= Rₖₖ - Y.du[k] *= Rₖ₋₁ₖ*Rₖₖ₊₁/Rₖₖ - Rₖ₋₁ₖ₊₁ - Y.du[k] += Xₖₖ₊₁ - Xₖₖ * Rₖₖ₊₁ / Rₖₖ - end - - k = n - Xₖₖ₋₁,Xₖₖ = X.dl[k-1], X.d[k] - Y.dl[k-1] = Xₖₖ₋₁/Rₖₖ - Y.d[k] = Xₖₖ-Xₖₖ₋₁*Rₖₖ₊₁/Rₖₖ - Rₖₖ = R.data[3,k] # R[2,2], R[2,3], R[1,3] - Y.d[k] /= Rₖₖ - Y.du[k-1] /= 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 @testset "T -> U" begin R = BandedMatrices._BandedMatrix(Vcat(-Ones(1,∞)/2, @@ -228,10 +110,10 @@ end 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) + @time UX = InfiniteLinearAlgebra.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) + @time Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U) @test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y end @testset "P -> Ultraspherical(3/2)" begin @@ -242,10 +124,40 @@ end 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) + @time UX = InfiniteLinearAlgebra.upper_mul_tri_triview(U, X) + @test Tridiagonal(U*X) ≈ UX + # U*X*inv(U) only depends on Tridiagonal(U*X) + @time Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U) + @test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y + end + + @testset "Jacobi(1,0) -> Jacobi(2,0)" begin + R = BandedMatrices._BandedMatrix(Vcat(Zeros(1,∞), # extra band since code assumes two bands + (-(0:∞) ./ (2:2:∞))', + ((2:∞) ./ (2:2:∞))'), ℵ₀, 0,2) + X_P = LazyBandedMatrices.Tridiagonal((2:∞) ./ (3:2:∞), -1 ./ ((1:2:∞) .* (3:2:∞)), (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 = InfiniteLinearAlgebra.upper_mul_tri_triview(U, X) + @test Tridiagonal(U*X) ≈ UX + # U*X*inv(U) only depends on Tridiagonal(U*X) + @time Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U) + @test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y + end + + @testset "Legendre() -> Jacobi(5/2)" begin + R = BandedMatrices._BandedMatrix(Vcat((-3 ./ (3:2:∞))', Zeros(1,∞), (3 ./ (3:2:∞))'), ℵ₀, 0,2) * + 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 = InfiniteLinearAlgebra.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) + @time Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U) @test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y end end \ No newline at end of file From b8c29cf35d29bfb4e1f1d0703576c52eabe0565f Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sat, 1 Feb 2025 09:52:01 +0000 Subject: [PATCH 08/18] split out pieces --- src/banded/tridiagonalconjugation.jl | 50 ++++++++++++++++++++++++---- 1 file changed, 44 insertions(+), 6 deletions(-) diff --git a/src/banded/tridiagonalconjugation.jl b/src/banded/tridiagonalconjugation.jl index 8107d5d..62246c1 100644 --- a/src/banded/tridiagonalconjugation.jl +++ b/src/banded/tridiagonalconjugation.jl @@ -4,7 +4,9 @@ upper_mul_tri_triview(U, X) == Tridiagonal(U*X) where U is Upper triangular Band function upper_mul_tri_triview(U::BandedMatrix, X::Tridiagonal) T = promote_type(eltype(U), eltype(X)) n = size(U,1) - upper_mul_tri_triview!(Tridiagonal(Vector{T}(undef, n-1), Vector{T}(undef, n), Vector{T}(undef, n-1)), U, X) + UX = Tridiagonal(Vector{T}(undef, n-1), Vector{T}(undef, n), Vector{T}(undef, n-1)) + + upper_mul_tri_triview!(UX, U, X) end function upper_mul_tri_triview!(UX::Tridiagonal, U::BandedMatrix, X::Tridiagonal) @@ -22,6 +24,19 @@ function upper_mul_tri_triview!(UX::Tridiagonal, U::BandedMatrix, X::Tridiagonal # Tridiagonal bands can be resized @assert length(Xdl)+1 == length(Xd) == length(Xdu)+1 == length(UXdl)+1 == length(UXd) == length(UXdu)+1 == n + UX, bⱼ, aⱼ, cⱼ, cⱼ₋₁ = initiate_upper_mul_tri_triview!(UX, U, X) + UX, bⱼ, aⱼ, cⱼ, cⱼ₋₁ = main_upper_mul_tri_triview!(UX, U, X, 2:n-2, bⱼ, aⱼ, cⱼ, cⱼ₋₁) + finalize_upper_mul_tri_triview!(UX, U, X, n-1, bⱼ, aⱼ, cⱼ, cⱼ₋₁) +end + + +function initiate_upper_mul_tri_triview!(UX, U, X) + Xdl, Xd, Xdu = X.dl, X.d, X.du + UXdl, UXd, UXdu = UX.dl, UX.d, UX.du + Udat = U.data + + l,u = bandwidths(U) + j = 1 aⱼ, cⱼ = Xd[1], Xdl[1] Uⱼⱼ, Uⱼⱼ₊₁, Uⱼⱼ₊₂ = Udat[u+1,1], Udat[u,2], Udat[u-1,3] # U[j,j], U[j,j+1], U[j,j+2] @@ -29,7 +44,17 @@ function upper_mul_tri_triview!(UX::Tridiagonal, U::BandedMatrix, X::Tridiagonal bⱼ, aⱼ, cⱼ, cⱼ₋₁ = Xdu[1], Xd[2], Xdl[2], cⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j] UXdu[1] = Uⱼⱼ*bⱼ + Uⱼⱼ₊₁*aⱼ + Uⱼⱼ₊₂*cⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+1]*X[j+1,j] - @inbounds for j = 2:n-2 + UX, bⱼ, aⱼ, cⱼ, cⱼ₋₁ +end + + +function main_upper_mul_tri_triview!(UX, U, X, jr, bⱼ, aⱼ, cⱼ, cⱼ₋₁) + Xdl, Xd, Xdu = X.dl, X.d, X.du + UXdl, UXd, UXdu = UX.dl, UX.d, UX.du + Udat = U.data + l,u = bandwidths(U) + + @inbounds for j = jr Uⱼⱼ, Uⱼⱼ₊₁, Uⱼⱼ₊₂ = Udat[u+1,j], Udat[u,j+1], Udat[u-1,j+2] # U[j,j], U[j,j+1], U[j,j+2] UXdl[j-1] = Uⱼⱼ*cⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1] UXd[j] = Uⱼⱼ*aⱼ + Uⱼⱼ₊₁*cⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j] @@ -37,14 +62,22 @@ function upper_mul_tri_triview!(UX::Tridiagonal, U::BandedMatrix, X::Tridiagonal UXdu[j] = Uⱼⱼ*bⱼ + Uⱼⱼ₊₁*aⱼ + Uⱼⱼ₊₂*cⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+2]*X[j+2,j+1] end - j = n-1 + UX, bⱼ, aⱼ, cⱼ, cⱼ₋₁ +end + +function finalize_upper_mul_tri_triview!(UX, U, X, j, bⱼ, aⱼ, cⱼ, cⱼ₋₁) + Xdl, Xd, Xdu = X.dl, X.d, X.du + UXdl, UXd, UXdu = UX.dl, UX.d, UX.du + Udat = U.data + l,u = bandwidths(U) + Uⱼⱼ, Uⱼⱼ₊₁ = Udat[u+1,j], Udat[u,j+1] # U[j,j], U[j,j+1] UXdl[j-1] = Uⱼⱼ*cⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1] UXd[j] = Uⱼⱼ*aⱼ + Uⱼⱼ₊₁*cⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j] bⱼ, aⱼ, cⱼ₋₁ = Xdu[j], Xd[j+1], cⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j] UXdu[j] = Uⱼⱼ*bⱼ + Uⱼⱼ₊₁*aⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+2]*X[j+2,j+1] - j = n + j += 1 Uⱼⱼ = Udat[u+1,j] # U[j,j] UXdl[j-1] = Uⱼⱼ*cⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1] UXd[j] = Uⱼⱼ*aⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j] @@ -132,8 +165,13 @@ function resizedata!(data::TridiagonalConjugationData, n) 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) + resize!(data.UX.dl, 2n) + resize!(data.UX.d, 2n + 1) + resize!(data.UX.du, 2n) + + resize!(data.Y.dl, 2n) + resize!(data.Y.d, 2n + 1) + resize!(data.Y.du, 2n) end From e815efe1d75007ad0a562a13131627fc4e84ccd5 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sat, 1 Feb 2025 22:30:20 +0000 Subject: [PATCH 09/18] adaptively populate UX --- src/banded/tridiagonalconjugation.jl | 79 +++++++++++++++------------- test/test_bidiagonalconjugation.jl | 2 + 2 files changed, 45 insertions(+), 36 deletions(-) diff --git a/src/banded/tridiagonalconjugation.jl b/src/banded/tridiagonalconjugation.jl index 62246c1..f1043bb 100644 --- a/src/banded/tridiagonalconjugation.jl +++ b/src/banded/tridiagonalconjugation.jl @@ -24,9 +24,9 @@ function upper_mul_tri_triview!(UX::Tridiagonal, U::BandedMatrix, X::Tridiagonal # Tridiagonal bands can be resized @assert length(Xdl)+1 == length(Xd) == length(Xdu)+1 == length(UXdl)+1 == length(UXd) == length(UXdu)+1 == n - UX, bⱼ, aⱼ, cⱼ, cⱼ₋₁ = initiate_upper_mul_tri_triview!(UX, U, X) - UX, bⱼ, aⱼ, cⱼ, cⱼ₋₁ = main_upper_mul_tri_triview!(UX, U, X, 2:n-2, bⱼ, aⱼ, cⱼ, cⱼ₋₁) - finalize_upper_mul_tri_triview!(UX, U, X, n-1, bⱼ, aⱼ, cⱼ, cⱼ₋₁) + UX, bₖ, aₖ, cₖ, cₖ₋₁ = initiate_upper_mul_tri_triview!(UX, U, X) + UX, bₖ, aₖ, cₖ, cₖ₋₁ = main_upper_mul_tri_triview!(UX, U, X, 2:n-2, bₖ, aₖ, cₖ, cₖ₋₁) + finalize_upper_mul_tri_triview!(UX, U, X, n-1, bₖ, aₖ, cₖ, cₖ₋₁) end @@ -37,50 +37,50 @@ function initiate_upper_mul_tri_triview!(UX, U, X) l,u = bandwidths(U) - j = 1 - aⱼ, cⱼ = Xd[1], Xdl[1] - Uⱼⱼ, Uⱼⱼ₊₁, Uⱼⱼ₊₂ = Udat[u+1,1], Udat[u,2], Udat[u-1,3] # U[j,j], U[j,j+1], U[j,j+2] - UXd[1] = Uⱼⱼ*aⱼ + Uⱼⱼ₊₁*cⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j] - bⱼ, aⱼ, cⱼ, cⱼ₋₁ = Xdu[1], Xd[2], Xdl[2], cⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j] - UXdu[1] = Uⱼⱼ*bⱼ + Uⱼⱼ₊₁*aⱼ + Uⱼⱼ₊₂*cⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+1]*X[j+1,j] + k = 1 + aₖ, cₖ = Xd[1], Xdl[1] + Uₖₖ, Uₖₖ₊₁, Uₖₖ₊₂ = Udat[u+1,1], Udat[u,2], Udat[u-1,3] # U[k,k], U[k,k+1], U[k,k+2] + UXd[1] = Uₖₖ*aₖ + Uₖₖ₊₁*cₖ # UX[k,k] = U[k,k]*X[k,k] + U[k,k+1]*X[k+1,k] + bₖ, aₖ, cₖ, cₖ₋₁ = Xdu[1], Xd[2], Xdl[2], cₖ # X[k,k+1], X[k+1,k+1], X[k+2,k+1], X[k+1,k] + UXdu[1] = Uₖₖ*bₖ + Uₖₖ₊₁*aₖ + Uₖₖ₊₂*cₖ # UX[k,k+1] = U[k,k]*X[k,k+1] + U[k,k+1]*X[k+1,k+1] + U[k,k+1]*X[k+1,k] - UX, bⱼ, aⱼ, cⱼ, cⱼ₋₁ + UX, bₖ, aₖ, cₖ, cₖ₋₁ end - -function main_upper_mul_tri_triview!(UX, U, X, jr, bⱼ, aⱼ, cⱼ, cⱼ₋₁) +# fills in the rows kr of UX +function main_upper_mul_tri_triview!(UX, U, X, kr, bₖ=X.du[kr[1]-1], aₖ=X.d[kr[1]], cₖ=X.dl[kr[1]], cₖ₋₁=X.du[kr[1]-1]) Xdl, Xd, Xdu = X.dl, X.d, X.du UXdl, UXd, UXdu = UX.dl, UX.d, UX.du Udat = U.data l,u = bandwidths(U) - @inbounds for j = jr - Uⱼⱼ, Uⱼⱼ₊₁, Uⱼⱼ₊₂ = Udat[u+1,j], Udat[u,j+1], Udat[u-1,j+2] # U[j,j], U[j,j+1], U[j,j+2] - UXdl[j-1] = Uⱼⱼ*cⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1] - UXd[j] = Uⱼⱼ*aⱼ + Uⱼⱼ₊₁*cⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j] - bⱼ, aⱼ, cⱼ, cⱼ₋₁ = Xdu[j], Xd[j+1], Xdl[j+1], cⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j] - UXdu[j] = Uⱼⱼ*bⱼ + Uⱼⱼ₊₁*aⱼ + Uⱼⱼ₊₂*cⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+2]*X[j+2,j+1] + for k = kr + Uₖₖ, Uₖₖ₊₁, Uₖₖ₊₂ = Udat[u+1,k], Udat[u,k+1], Udat[u-1,k+2] # U[k,k], U[k,k+1], U[k,k+2] + UXdl[k-1] = Uₖₖ*cₖ₋₁ # UX[k,k-1] = U[k,k]*X[k,k-1] + UXd[k] = Uₖₖ*aₖ + Uₖₖ₊₁*cₖ # UX[k,k] = U[k,k]*X[k,k] + U[k,k+1]*X[k+1,k] + bₖ, aₖ, cₖ, cₖ₋₁ = Xdu[k], Xd[k+1], Xdl[k+1], cₖ # X[k,k+1], X[k+1,k+1], X[k+2,k+1], X[k+1,k] + UXdu[k] = Uₖₖ*bₖ + Uₖₖ₊₁*aₖ + Uₖₖ₊₂*cₖ # UX[k,k+1] = U[k,k]*X[k,k+1] + U[k,k+1]*X[k+1,k+1] + U[k,k+2]*X[k+2,k+1] end - UX, bⱼ, aⱼ, cⱼ, cⱼ₋₁ + UX, bₖ, aₖ, cₖ, cₖ₋₁ end -function finalize_upper_mul_tri_triview!(UX, U, X, j, bⱼ, aⱼ, cⱼ, cⱼ₋₁) +function finalize_upper_mul_tri_triview!(UX, U, X, k, bₖ, aₖ, cₖ, cₖ₋₁) Xdl, Xd, Xdu = X.dl, X.d, X.du UXdl, UXd, UXdu = UX.dl, UX.d, UX.du Udat = U.data l,u = bandwidths(U) - Uⱼⱼ, Uⱼⱼ₊₁ = Udat[u+1,j], Udat[u,j+1] # U[j,j], U[j,j+1] - UXdl[j-1] = Uⱼⱼ*cⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1] - UXd[j] = Uⱼⱼ*aⱼ + Uⱼⱼ₊₁*cⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j] - bⱼ, aⱼ, cⱼ₋₁ = Xdu[j], Xd[j+1], cⱼ # X[j,j+1], X[j+1,j+1], X[j+2,j+1], X[j+1,j] - UXdu[j] = Uⱼⱼ*bⱼ + Uⱼⱼ₊₁*aⱼ # UX[j,j+1] = U[j,j]*X[j,j+1] + U[j,j+1]*X[j+1,j+1] + U[j,j+2]*X[j+2,j+1] + Uₖₖ, Uₖₖ₊₁ = Udat[u+1,k], Udat[u,k+1] # U[k,k], U[k,k+1] + UXdl[k-1] = Uₖₖ*cₖ₋₁ # UX[k,k-1] = U[k,k]*X[k,k-1] + UXd[k] = Uₖₖ*aₖ + Uₖₖ₊₁*cₖ # UX[k,k] = U[k,k]*X[k,k] + U[k,k+1]*X[k+1,k] + bₖ, aₖ, cₖ₋₁ = Xdu[k], Xd[k+1], cₖ # X[k,k+1], X[k+1,k+1], X[k+2,k+1], X[k+1,k] + UXdu[k] = Uₖₖ*bₖ + Uₖₖ₊₁*aₖ # UX[k,k+1] = U[k,k]*X[k,k+1] + U[k,k+1]*X[k+1,k+1] + U[k,k+2]*X[k+2,k+1] - j += 1 - Uⱼⱼ = Udat[u+1,j] # U[j,j] - UXdl[j-1] = Uⱼⱼ*cⱼ₋₁ # UX[j,j-1] = U[j,j]*X[j,j-1] - UXd[j] = Uⱼⱼ*aⱼ # UX[j,j] = U[j,j]*X[j,j] + U[j,j+1]*X[j+1,j] + k += 1 + Uₖₖ = Udat[u+1,k] # U[k,k] + UXdl[k-1] = Uₖₖ*cₖ₋₁ # UX[k,k-1] = U[k,k]*X[k,k-1] + UXd[k] = Uₖₖ*aₖ # UX[k,k] = U[k,k]*X[k,k] + U[k,k+1]*X[k+1,k] UX end @@ -152,19 +152,24 @@ mutable struct TridiagonalConjugationData{T} 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) +function TridiagonalConjugationData(U, X, V) + T = promote_type(typeof(inv(V[1, 1])), eltype(U), eltype(X)) # include inv so that we can't get Ints + n_init = 100 + UX = Tridiagonal(Vector{T}(undef, n_init-1), Vector{T}(undef, n_init), Vector{T}(undef, n_init-1)) + Y = Tridiagonal(Vector{T}(undef, n_init-1), Vector{T}(undef, n_init), Vector{T}(undef, n_init-1)) + initiate_upper_mul_tri_triview!(UX, U, X) # fill-in 1st row + return TridiagonalConjugationData(U, X, V, UX, Y, 1) end +TridiagonalConjugationData(U, X) = TridiagonalConjugationData(U, X, U) + 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) + n ≤ data.datasize && return data + + if n > length(data.UX.d) # 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) @@ -174,6 +179,8 @@ function resizedata!(data::TridiagonalConjugationData, n) resize!(data.Y.du, 2n) end + main_upper_mul_tri_triview!(data.UX, data.U, data.X, data.datasize+1:n) + data.datasize = n end diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl index bfbf481..e706295 100644 --- a/test/test_bidiagonalconjugation.jl +++ b/test/test_bidiagonalconjugation.jl @@ -115,6 +115,8 @@ end # U*X*inv(U) only depends on Tridiagonal(U*X) @time Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U) @test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y + + InfiniteLinearAlgebra.TridiagonalConjugationData(U, X, U) end @testset "P -> Ultraspherical(3/2)" begin R = BandedMatrices._BandedMatrix(Vcat((-1 ./ (1:2:∞))', From 0780e9a1afe365645725329dc877a8c2aa0ed68c Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sun, 2 Feb 2025 14:29:02 +0000 Subject: [PATCH 10/18] split tri_mul_invupper in 3 --- src/banded/tridiagonalconjugation.jl | 66 ++++++++++++++++++++-------- 1 file changed, 48 insertions(+), 18 deletions(-) diff --git a/src/banded/tridiagonalconjugation.jl b/src/banded/tridiagonalconjugation.jl index f1043bb..96a384c 100644 --- a/src/banded/tridiagonalconjugation.jl +++ b/src/banded/tridiagonalconjugation.jl @@ -1,22 +1,20 @@ -""" -upper_mul_tri_triview(U, X) == Tridiagonal(U*X) where U is Upper triangular BandedMatrix and X is Tridiagonal -""" + +# upper_mul_tri_triview(U, X) == Tridiagonal(U*X) where U is Upper triangular BandedMatrix and X is Tridiagonal function upper_mul_tri_triview(U::BandedMatrix, X::Tridiagonal) T = promote_type(eltype(U), eltype(X)) n = size(U,1) UX = Tridiagonal(Vector{T}(undef, n-1), Vector{T}(undef, n), Vector{T}(undef, n-1)) - + upper_mul_tri_triview!(UX, U, X) end function upper_mul_tri_triview!(UX::Tridiagonal, U::BandedMatrix, X::Tridiagonal) n = size(UX,1) - + Xdl, Xd, Xdu = X.dl, X.d, X.du UXdl, UXd, UXdu = UX.dl, UX.d, UX.du - Udat = U.data - + l,u = bandwidths(U) @assert size(U) == (n,n) @@ -29,7 +27,7 @@ function upper_mul_tri_triview!(UX::Tridiagonal, U::BandedMatrix, X::Tridiagonal finalize_upper_mul_tri_triview!(UX, U, X, n-1, bₖ, aₖ, cₖ, cₖ₋₁) end - +# populate first row of UX with UX function initiate_upper_mul_tri_triview!(UX, U, X) Xdl, Xd, Xdu = X.dl, X.d, X.du UXdl, UXd, UXdu = UX.dl, UX.d, UX.du @@ -65,6 +63,7 @@ function main_upper_mul_tri_triview!(UX, U, X, kr, bₖ=X.du[kr[1]-1], aₖ=X.d[ UX, bₖ, aₖ, cₖ, cₖ₋₁ end +# populate rows k and k+1 of UX, assuming we are at the bottom-right function finalize_upper_mul_tri_triview!(UX, U, X, k, bₖ, aₖ, cₖ, cₖ₋₁) Xdl, Xd, Xdu = X.dl, X.d, X.du UXdl, UXd, UXdu = UX.dl, UX.d, UX.du @@ -92,43 +91,74 @@ end tri_mul_invupper_triview(X::Tridiagonal, R::BandedMatrix) = tri_mul_invupper_triview!(similar(X, promote_type(eltype(X), eltype(R))), X, R) + function tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::BandedMatrix) n = size(X,1) Xdl, Xd, Xdu = X.dl, X.d, X.du Ydl, Yd, Ydu = Y.dl, Y.d, Y.du - Rdat = R.data - + l,u = bandwidths(R) - + @assert size(R) == (n,n) @assert l == 0 && u ≥ 2 # Tridiagonal bands can be resized @assert length(Xdl)+1 == length(Xd) == length(Xdu)+1 == length(Ydl)+1 == length(Yd) == length(Ydu)+1 == n - - + + UX, Rₖₖ, Rₖₖ₊₁ = initiate_tri_mul_invupper_triview!(Y, X, R) + UX, Rₖₖ, Rₖₖ₊₁ = main_tri_mul_invupper_triview!(Y, X, R, 2:n-1, Rₖₖ, Rₖₖ₊₁) + finalize_tri_mul_invupper_triview!(Y, X, R, n, Rₖₖ, Rₖₖ₊₁) +end + +# populate first row of X/R +function initiate_tri_mul_invupper_triview!(Y, X, R) + Xdl, Xd, Xdu = X.dl, X.d, X.du + Ydl, Yd, Ydu = Y.dl, Y.d, Y.du + Rdat = R.data + + l,u = bandwidths(R) + k = 1 aₖ,bₖ = Xd[k], Xdu[k] Rₖₖ,Rₖₖ₊₁ = Rdat[u+1,k], Rdat[u,k+1] # R[1,1], R[1,2] Yd[k] = aₖ/Rₖₖ Ydu[k] = bₖ - aₖ * Rₖₖ₊₁/Rₖₖ - @inbounds for k = 2:n-1 + Y, Rₖₖ, Rₖₖ₊₁ +end + + +# populate rows kr of X/R +function main_tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::BandedMatrix, kr, Rₖₖ=R[first(kr),first(kr)], Rₖₖ₊₁=R[first(kr),first(kr)+1]) + Xdl, Xd, Xdu = X.dl, X.d, X.du + Ydl, Yd, Ydu = Y.dl, Y.d, Y.du + Rdat = R.data + l,u = bandwidths(R) + + @inbounds for k = kr cₖ₋₁,aₖ,bₖ = Xdl[k-1], Xd[k], Xdu[k] Ydl[k-1] = cₖ₋₁/Rₖₖ Yd[k] = aₖ-cₖ₋₁*Rₖₖ₊₁/Rₖₖ Ydu[k] = cₖ₋₁/Rₖₖ - Rₖₖ,Rₖₖ₊₁,Rₖ₋₁ₖ₊₁,Rₖ₋₁ₖ = Rdat[u+1,k], Rdat[u,k+1],Rdat[u-1,k+1],Rₖₖ₊₁ # R[2,2], R[2,3], R[1,3] + Rₖₖ,Rₖₖ₊₁,Rₖ₋₁ₖ₊₁,Rₖ₋₁ₖ = Rdat[u+1,k], Rdat[u,k+1],Rdat[u-1,k+1],Rₖₖ₊₁ # R[k,k], R[k,k+1], R[k-1,k] Yd[k] /= Rₖₖ Ydu[k-1] /= Rₖₖ Ydu[k] *= Rₖ₋₁ₖ*Rₖₖ₊₁/Rₖₖ - Rₖ₋₁ₖ₊₁ Ydu[k] += bₖ - aₖ * Rₖₖ₊₁ / Rₖₖ end + Y, Rₖₖ, Rₖₖ₊₁ +end + - k = n +# populate row k of X/R, assuming we are at the bottom-right +function finalize_tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::BandedMatrix, k, Rₖₖ=R[k-1,k-1], Rₖₖ₊₁=R[k-1,k]) + Xdl, Xd, Xdu = X.dl, X.d, X.du + Ydl, Yd, Ydu = Y.dl, Y.d, Y.du + Rdat = R.data + l,u = bandwidths(R) cₖ₋₁,aₖ = Xdl[k-1], Xd[k] Ydl[k-1] = cₖ₋₁/Rₖₖ Yd[k] = aₖ-cₖ₋₁*Rₖₖ₊₁/Rₖₖ - Rₖₖ = Rdat[u+1,k] # R[2,2], R[2,3], R[1,3] + Rₖₖ = Rdat[u+1,k] # R[k,k] Yd[k] /= Rₖₖ Ydu[k-1] /= Rₖₖ @@ -173,7 +203,7 @@ function resizedata!(data::TridiagonalConjugationData, n) resize!(data.UX.dl, 2n) resize!(data.UX.d, 2n + 1) resize!(data.UX.du, 2n) - + resize!(data.Y.dl, 2n) resize!(data.Y.d, 2n + 1) resize!(data.Y.du, 2n) From 186a218f2f5b6ac4f5feea462e239abaed2c28b7 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 3 Feb 2025 22:13:33 +0000 Subject: [PATCH 11/18] Support caching --- Project.toml | 2 +- src/InfiniteLinearAlgebra.jl | 4 +- src/banded/tridiagonalconjugation.jl | 53 +++++++++++++--- src/infqr.jl | 13 ++-- test/test_bidiagonalconjugation.jl | 92 +++++++++++++--------------- 5 files changed, 98 insertions(+), 66 deletions(-) diff --git a/Project.toml b/Project.toml index af9c2b2..6b5c318 100644 --- a/Project.toml +++ b/Project.toml @@ -26,7 +26,7 @@ FillArrays = "1.0" InfiniteArrays = "0.15" InfiniteRandomArrays = "0.2" Infinities = "0.1" -LazyArrays = "2.3" +LazyArrays = "2.5" LazyBandedMatrices = "0.11" LinearAlgebra = "1" MatrixFactorizations = "3.0" diff --git a/src/InfiniteLinearAlgebra.jl b/src/InfiniteLinearAlgebra.jl index d6f7336..70d8425 100644 --- a/src/InfiniteLinearAlgebra.jl +++ b/src/InfiniteLinearAlgebra.jl @@ -14,7 +14,7 @@ import Base.Broadcast: BroadcastStyle, Broadcasted, broadcasted import ArrayLayouts: AbstractBandedLayout, AbstractQLayout, AdjQRPackedQLayout, CNoPivot, DenseColumnMajor, FillLayout, MatLdivVec, MatLmulMat, MatLmulVec, MemoryLayout, QRPackedQLayout, RangeCumsum, TriangularLayout, - TridiagonalLayout, __qr, _factorize, _qr, check_mul_axes, colsupport, + TridiagonalLayout, _qr_layout, factorize_layout, qr_layout, check_mul_axes, colsupport, diagonaldata, ldiv!, lmul!, mul, mulreduce, reflector!, reflectorApply!, rowsupport, sub_materialize, subdiagonaldata, sublayout, supdiagonaldata, transposelayout, triangulardata, triangularlayout, zero!, materialize! @@ -40,7 +40,7 @@ import LazyArrays: AbstractCachedMatrix, AbstractCachedVector, AbstractLazyLayou CachedArray, CachedLayout, CachedMatrix, CachedVector, LazyArrayStyle, LazyLayout, LazyLayouts, LazyMatrix, LazyVector, AbstractPaddedLayout, PaddedColumns, _broadcast_sub_arguments, applybroadcaststyle, applylayout, arguments, cacheddata, paddeddata, resizedata!, simplifiable, - simplify, islazy, islazy_layout, cache_getindex + simplify, islazy, islazy_layout, cache_getindex, cache_layout import LazyBandedMatrices: AbstractLazyBandedBlockBandedLayout, AbstractLazyBandedLayout, ApplyBandedLayout, BlockVec, BroadcastBandedLayout, KronTravBandedBlockBandedLayout, LazyBandedLayout, diff --git a/src/banded/tridiagonalconjugation.jl b/src/banded/tridiagonalconjugation.jl index 96a384c..1e082d9 100644 --- a/src/banded/tridiagonalconjugation.jl +++ b/src/banded/tridiagonalconjugation.jl @@ -28,7 +28,11 @@ function upper_mul_tri_triview!(UX::Tridiagonal, U::BandedMatrix, X::Tridiagonal end # populate first row of UX with UX -function initiate_upper_mul_tri_triview!(UX, U, X) + +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) + +function initiate_upper_mul_tri_triview!(UX, U::BandedMatrix, X) Xdl, Xd, Xdu = X.dl, X.d, X.du UXdl, UXd, UXdu = UX.dl, UX.d, UX.du Udat = U.data @@ -46,7 +50,12 @@ function initiate_upper_mul_tri_triview!(UX, U, X) end # fills in the rows kr of UX -function main_upper_mul_tri_triview!(UX, U, X, kr, bₖ=X.du[kr[1]-1], aₖ=X.d[kr[1]], cₖ=X.dl[kr[1]], cₖ₋₁=X.du[kr[1]-1]) +function main_upper_mul_tri_triview!(UX, U::CachedMatrix, X, kr, kwds...) + resizedata!(U, kr[end], kr[end]+2) + main_upper_mul_tri_triview!(UX, U.data, X, kr, kwds...) +end + +function main_upper_mul_tri_triview!(UX, U::BandedMatrix, X, kr, bₖ=X.du[kr[1]-1], aₖ=X.d[kr[1]], cₖ=X.dl[kr[1]], cₖ₋₁=X.dl[kr[1]-1]) Xdl, Xd, Xdu = X.dl, X.d, X.du UXdl, UXd, UXdu = UX.dl, UX.d, UX.du Udat = U.data @@ -109,8 +118,14 @@ function tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::BandedMatr finalize_tri_mul_invupper_triview!(Y, X, R, n, Rₖₖ, Rₖₖ₊₁) end -# populate first row of X/R -function initiate_tri_mul_invupper_triview!(Y, X, R) +# partially-populate first row of X/R +# Ydu[k] is updated below +function initiate_tri_mul_invupper_triview!(Y, X, R::CachedMatrix) + resizedata!(R, 1, 2) + initiate_tri_mul_invupper_triview!(Y, X, R.data) +end + +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 Rdat = R.data @@ -120,6 +135,7 @@ function initiate_tri_mul_invupper_triview!(Y, X, R) k = 1 aₖ,bₖ = Xd[k], Xdu[k] Rₖₖ,Rₖₖ₊₁ = Rdat[u+1,k], Rdat[u,k+1] # R[1,1], R[1,2] + Yd[k] = aₖ/Rₖₖ Ydu[k] = bₖ - aₖ * Rₖₖ₊₁/Rₖₖ @@ -127,8 +143,13 @@ function initiate_tri_mul_invupper_triview!(Y, X, R) end -# populate rows kr of X/R -function main_tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::BandedMatrix, kr, Rₖₖ=R[first(kr),first(kr)], Rₖₖ₊₁=R[first(kr),first(kr)+1]) +# 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...) + resizedata!(R, kr[end], kr[end]+1) + main_tri_mul_invupper_triview!(Y, X, R.data, kr, kwds...) +end + +function main_tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::BandedMatrix, kr, Rₖₖ=R[first(kr)-1,first(kr)-1], Rₖₖ₊₁=R[first(kr)-1,first(kr)]) Xdl, Xd, Xdu = X.dl, X.d, X.du Ydl, Yd, Ydu = Y.dl, Y.d, Y.du Rdat = R.data @@ -187,11 +208,18 @@ function TridiagonalConjugationData(U, X, V) n_init = 100 UX = Tridiagonal(Vector{T}(undef, n_init-1), Vector{T}(undef, n_init), Vector{T}(undef, n_init-1)) Y = Tridiagonal(Vector{T}(undef, n_init-1), Vector{T}(undef, n_init), Vector{T}(undef, n_init-1)) + resizedata!(U, n_init, n_init) + resizedata!(V, n_init, n_init) initiate_upper_mul_tri_triview!(UX, U, X) # fill-in 1st row - return TridiagonalConjugationData(U, X, V, UX, Y, 1) + initiate_tri_mul_invupper_triview!(Y, UX, V) + return TridiagonalConjugationData(U, X, V, UX, Y, 0) end -TridiagonalConjugationData(U, X) = TridiagonalConjugationData(U, X, U) + +function TridiagonalConjugationData(U, X) + C = cache(U) + TridiagonalConjugationData(C, X, C) +end copy(data::TridiagonalConjugationData) = TridiagonalConjugationData(copy(data.U), copy(data.X), copy(data.V), copy(data.UX), copy(data.Y), data.datasize) @@ -209,8 +237,13 @@ function resizedata!(data::TridiagonalConjugationData, n) resize!(data.Y.du, 2n) end - main_upper_mul_tri_triview!(data.UX, data.U, data.X, data.datasize+1:n) - data.datasize = n + if n > data.datasize + main_upper_mul_tri_triview!(data.UX, data.U, data.X, data.datasize+2:n+1) + main_tri_mul_invupper_triview!(data.Y, data.UX, data.U, data.datasize+2:n+1) # need one extra as it updates first row + data.datasize = n + end + + data end diff --git a/src/infqr.jl b/src/infqr.jl index 5d86f28..76de794 100644 --- a/src/infqr.jl +++ b/src/infqr.jl @@ -149,11 +149,14 @@ function adaptiveqr(A) QR(AdaptiveQRFactors(data), AdaptiveQRTau(data)) end -_qr(::AbstractBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A) -_qr(::AbstractAlmostBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A) -__qr(_, ::NTuple{2,InfiniteCardinal{0}}, A) = adaptiveqr(A) -_qr(::AbstractBlockBandedLayout, ::NTuple{2,InfiniteCardinal{0}}, A) = adaptiveqr(A) -_factorize(::AbstractBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = qr(A) +qr_layout(::AbstractBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A) +qr_layout(::AbstractAlmostBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A) +__qr_layout(_, ::NTuple{2,InfiniteCardinal{0}}, A) = adaptiveqr(A) +qr_layout(::AbstractBlockBandedLayout, ::NTuple{2,InfiniteCardinal{0}}, A) = adaptiveqr(A) +factorize_layout(::AbstractBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = qr(A) + + +cache_layout(::TriangularLayout{UPLO, UNIT, <:AdaptiveLayout}, A::AbstractMatrix) where {UPLO, UNIT} = A # already cached partialqr!(F::QR, n) = partialqr!(F.factors, n) partialqr!(F::AdaptiveQRFactors, n) = partialqr!(F.data, n) diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl index e706295..a9c2307 100644 --- a/test/test_bidiagonalconjugation.jl +++ b/test/test_bidiagonalconjugation.jl @@ -3,6 +3,7 @@ using InfiniteLinearAlgebra: BidiagonalConjugation, OneToInf using ArrayLayouts: supdiagonaldata, subdiagonaldata, diagonaldata using LinearAlgebra using LazyArrays: LazyLayout +using BandedMatrices: _BandedMatrix @testset "BidiagonalConjugation" begin @test InfiniteLinearAlgebra._to_uplo('U') == 'U' @@ -102,11 +103,25 @@ end @testset "TridiagonalConjugation" begin - @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,∞)) + for (R,X_T) in ( + # T -> U + (_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,∞))), + # P -> C^(3/2) + (_BandedMatrix(Vcat((-1 ./ (1:2:∞))', + Zeros(1,∞), + (1 ./ (1:2:∞))'), ℵ₀, 0,2), + LazyBandedMatrices.Tridiagonal((1:∞) ./ (1:2:∞), Zeros(∞), (1:∞) ./ (3:2:∞))), + # P^(1,0) -> P^(2,0) + (_BandedMatrix(Vcat(Zeros(1,∞), # extra band since code assumes two bands + (-(0:∞) ./ (2:2:∞))', + ((2:∞) ./ (2:2:∞))'), ℵ₀, 0,2), + LazyBandedMatrices.Tridiagonal((2:∞) ./ (3:2:∞), -1 ./ ((1:2:∞) .* (3:2:∞)), (1:∞) ./ (3:2:∞))), + # P -> C^(5/2) + (_BandedMatrix(Vcat((-3 ./ (3:2:∞))', Zeros(1,∞), (3 ./ (3:2:∞))'), ℵ₀, 0,2) * + _BandedMatrix(Vcat((-1 ./ (1:2:∞))', Zeros(1,∞), (1 ./ (1:2:∞))'), ℵ₀, 0,2), + LazyBandedMatrices.Tridiagonal((1:∞) ./ (1:2:∞), Zeros(∞), (1:∞) ./ (3: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])); @@ -116,50 +131,31 @@ end @time Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U) @test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y - InfiniteLinearAlgebra.TridiagonalConjugationData(U, X, U) - 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 = InfiniteLinearAlgebra.upper_mul_tri_triview(U, X) - @test Tridiagonal(U*X) ≈ UX - # U*X*inv(U) only depends on Tridiagonal(U*X) - @time Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U) - @test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y - end + data = InfiniteLinearAlgebra.TridiagonalConjugationData(R, X_T) + @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] - @testset "Jacobi(1,0) -> Jacobi(2,0)" begin - R = BandedMatrices._BandedMatrix(Vcat(Zeros(1,∞), # extra band since code assumes two bands - (-(0:∞) ./ (2:2:∞))', - ((2:∞) ./ (2:2:∞))'), ℵ₀, 0,2) - X_P = LazyBandedMatrices.Tridiagonal((2:∞) ./ (3:2:∞), -1 ./ ((1:2:∞) .* (3:2:∞)), (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 = InfiniteLinearAlgebra.upper_mul_tri_triview(U, X) - @test Tridiagonal(U*X) ≈ UX - # U*X*inv(U) only depends on Tridiagonal(U*X) - @time Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U) - @test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y + data = InfiniteLinearAlgebra.TridiagonalConjugationData(R, X_T) + 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 - - @testset "Legendre() -> Jacobi(5/2)" begin - R = BandedMatrices._BandedMatrix(Vcat((-3 ./ (3:2:∞))', Zeros(1,∞), (3 ./ (3:2:∞))'), ℵ₀, 0,2) * - 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 = InfiniteLinearAlgebra.upper_mul_tri_triview(U, X) - @test Tridiagonal(U*X) ≈ UX - # U*X*inv(U) only depends on Tridiagonal(U*X) - @time Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U) - @test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y + + @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) end end \ No newline at end of file From 04b6ffe181d69eae77550c5cae7ee68afb4a4f40 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 4 Feb 2025 22:20:09 +0000 Subject: [PATCH 12/18] Cholesky R --- src/banded/tridiagonalconjugation.jl | 15 +++++++++++++-- src/infcholesky.jl | 4 ++++ test/test_bidiagonalconjugation.jl | 21 ++++++++++++++++++++- 3 files changed, 37 insertions(+), 3 deletions(-) diff --git a/src/banded/tridiagonalconjugation.jl b/src/banded/tridiagonalconjugation.jl index 1e082d9..27ba3ce 100644 --- a/src/banded/tridiagonalconjugation.jl +++ b/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 diff --git a/src/infcholesky.jl b/src/infcholesky.jl index 6292579..7324b4b 100644 --- a/src/infcholesky.jl +++ b/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] diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl index a9c2307..463c5b2 100644 --- a/test/test_bidiagonalconjugation.jl +++ b/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 \ No newline at end of file From 9cd942735421c376abbbd986f1ca987f138c73dc Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Wed, 5 Feb 2025 19:57:40 +0000 Subject: [PATCH 13/18] Add Sym/TridiagonalConjugation --- src/banded/tridiagonalconjugation.jl | 34 +++++++++++++++++++++- test/test_bidiagonalconjugation.jl | 42 ++++++++++++++++++---------- 2 files changed, 60 insertions(+), 16 deletions(-) diff --git a/src/banded/tridiagonalconjugation.jl b/src/banded/tridiagonalconjugation.jl index 27ba3ce..8e00fe4 100644 --- a/src/banded/tridiagonalconjugation.jl +++ b/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 \ No newline at end of file diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl index 463c5b2..8ff176b 100644 --- a/test/test_bidiagonalconjugation.jl +++ b/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 \ No newline at end of file From dba640ae88bdcd89b5fc2d2991bf966735b5c812 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Wed, 5 Feb 2025 22:31:16 +0000 Subject: [PATCH 14/18] allow SymTridiagonal --- src/InfiniteLinearAlgebra.jl | 5 +++++ src/banded/tridiagonalconjugation.jl | 12 ++++++------ src/infqr.jl | 4 ++-- 3 files changed, 13 insertions(+), 8 deletions(-) diff --git a/src/InfiniteLinearAlgebra.jl b/src/InfiniteLinearAlgebra.jl index 70d8425..2e60767 100644 --- a/src/InfiniteLinearAlgebra.jl +++ b/src/InfiniteLinearAlgebra.jl @@ -46,6 +46,11 @@ import LazyBandedMatrices: AbstractLazyBandedBlockBandedLayout, AbstractLazyBand BroadcastBandedLayout, KronTravBandedBlockBandedLayout, LazyBandedLayout, _block_interlace_axes, _krontrav_axes, krontravargs +const StructuredLayoutTypes{Lay} = Union{SymmetricLayout{Lay}, HermitianLayout{Lay}, TriangularLayout{'L','N',Lay}, TriangularLayout{'U','N',Lay}, TriangularLayout{'L','U',Lay}, TriangularLayout{'U','U',Lay}} + +const BandedLayouts = Union{AbstractBandedLayout, StructuredLayoutTypes{<:AbstractBandedLayout}} + + import LinearAlgebra: AbstractQ, AdjointQ, AdjOrTrans, factorize, matprod, qr import MatrixFactorizations: AdjQLPackedQLayout, LayoutQ, QL, QLPackedQ, QLPackedQLayout, QR, QRPackedQ, diff --git a/src/banded/tridiagonalconjugation.jl b/src/banded/tridiagonalconjugation.jl index 8e00fe4..4cfe5d5 100644 --- a/src/banded/tridiagonalconjugation.jl +++ b/src/banded/tridiagonalconjugation.jl @@ -31,10 +31,10 @@ 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) +initiate_upper_mul_tri_triview!(UX, U::Union{AdaptiveCholeskyFactors,AdaptiveQRFactors}, 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 + Xdl, Xd, Xdu = subdiagonaldata(X), diagonaldata(X), supdiagonaldata(X) UXdl, UXd, UXdu = UX.dl, UX.d, UX.du Udat = U.data @@ -53,13 +53,13 @@ end # fills in the rows kr of UX 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...) +function main_upper_mul_tri_triview!(UX, U::Union{CachedMatrix,AdaptiveCholeskyFactors,AdaptiveQRFactors}, X, kr, kwds...) resizedata!(U, kr[end], kr[end]+2) main_upper_mul_tri_triview!(UX, U.data, X, kr, kwds...) end -function main_upper_mul_tri_triview!(UX, U::BandedMatrix, X, kr, bₖ=X.du[kr[1]-1], aₖ=X.d[kr[1]], cₖ=X.dl[kr[1]], cₖ₋₁=X.dl[kr[1]-1]) - Xdl, Xd, Xdu = X.dl, X.d, X.du +function main_upper_mul_tri_triview!(UX, U::BandedMatrix, X, kr, bₖ=X[kr[1]-1,kr[1]], aₖ=X[kr[1],kr[1]], cₖ=X[kr[1]+1,kr[1]], cₖ₋₁=X[kr[1],kr[1]-1]) + Xdl, Xd, Xdu = subdiagonaldata(X), diagonaldata(X), supdiagonaldata(X) UXdl, UXd, UXdu = UX.dl, UX.d, UX.du Udat = U.data l,u = bandwidths(U) @@ -128,7 +128,7 @@ 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) +function initiate_tri_mul_invupper_triview!(Y, X, R::Union{AdaptiveCholeskyFactors,AdaptiveQRFactors}) resizedata!(R, 1, 2) initiate_tri_mul_invupper_triview!(Y, X, R.data.data) end diff --git a/src/infqr.jl b/src/infqr.jl index 76de794..4f908f9 100644 --- a/src/infqr.jl +++ b/src/infqr.jl @@ -149,11 +149,11 @@ function adaptiveqr(A) QR(AdaptiveQRFactors(data), AdaptiveQRTau(data)) end -qr_layout(::AbstractBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A) +qr_layout(::BandedLayouts, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A) qr_layout(::AbstractAlmostBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A) __qr_layout(_, ::NTuple{2,InfiniteCardinal{0}}, A) = adaptiveqr(A) qr_layout(::AbstractBlockBandedLayout, ::NTuple{2,InfiniteCardinal{0}}, A) = adaptiveqr(A) -factorize_layout(::AbstractBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = qr(A) +factorize_layout(::BandedLayouts, ::NTuple{2,OneToInf{Int}}, A) = qr(A) cache_layout(::TriangularLayout{UPLO, UNIT, <:AdaptiveLayout}, A::AbstractMatrix) where {UPLO, UNIT} = A # already cached From b98359ce36d67ba1c398b76266813c826dd4e882 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Thu, 6 Feb 2025 07:30:50 +0000 Subject: [PATCH 15/18] support QR --- src/banded/tridiagonalconjugation.jl | 17 ++++++++++++++--- src/infcholesky.jl | 2 +- src/infqr.jl | 7 +++++++ 3 files changed, 22 insertions(+), 4 deletions(-) diff --git a/src/banded/tridiagonalconjugation.jl b/src/banded/tridiagonalconjugation.jl index 4cfe5d5..ce3a80d 100644 --- a/src/banded/tridiagonalconjugation.jl +++ b/src/banded/tridiagonalconjugation.jl @@ -53,11 +53,17 @@ end # fills in the rows kr of UX 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,AdaptiveQRFactors}, 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 +function main_upper_mul_tri_triview!(UX, U::AdaptiveQRFactors, X, kr, kwds...) + resizedata!(U, kr[end], kr[end]+2) + main_upper_mul_tri_triview!(UX, U.data.data, X, kr, kwds...) +end + + function main_upper_mul_tri_triview!(UX, U::BandedMatrix, X, kr, bₖ=X[kr[1]-1,kr[1]], aₖ=X[kr[1],kr[1]], cₖ=X[kr[1]+1,kr[1]], cₖ₋₁=X[kr[1],kr[1]-1]) Xdl, Xd, Xdu = subdiagonaldata(X), diagonaldata(X), supdiagonaldata(X) UXdl, UXd, UXdu = UX.dl, UX.d, UX.du @@ -160,6 +166,11 @@ function main_tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::Union main_tri_mul_invupper_triview!(Y, X, R.data, kr, kwds...) end +function main_tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::AdaptiveQRFactors, kr, kwds...) + resizedata!(R, kr[end], kr[end]+1) + main_tri_mul_invupper_triview!(Y, X, R.data.data, kr, kwds...) +end + function main_tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::BandedMatrix, kr, Rₖₖ=R[first(kr)-1,first(kr)-1], Rₖₖ₊₁=R[first(kr)-1,first(kr)]) Xdl, Xd, Xdu = X.dl, X.d, X.du Ydl, Yd, Ydu = Y.dl, Y.d, Y.du @@ -217,8 +228,8 @@ end function TridiagonalConjugationData(U, X, V) T = promote_type(typeof(inv(V[1, 1])), eltype(U), eltype(X)) # include inv so that we can't get Ints n_init = 100 - UX = Tridiagonal(Vector{T}(undef, n_init-1), Vector{T}(undef, n_init), Vector{T}(undef, n_init-1)) - Y = Tridiagonal(Vector{T}(undef, n_init-1), Vector{T}(undef, n_init), Vector{T}(undef, n_init-1)) + UX = Tridiagonal(zeros(T, n_init-1), zeros(T, n_init), zeros(T, n_init-1)) # zeros for BigFLoat + Y = Tridiagonal(zeros(T, n_init-1), zeros(T, n_init), zeros(T, n_init-1)) resizedata!(U, n_init, n_init) resizedata!(V, n_init, n_init) initiate_upper_mul_tri_triview!(UX, U, X) # fill-in 1st row diff --git a/src/infcholesky.jl b/src/infcholesky.jl index 7324b4b..f382aaf 100644 --- a/src/infcholesky.jl +++ b/src/infcholesky.jl @@ -12,7 +12,7 @@ const SymmetricBandedLayouts = Union{SymTridiagonalLayout,SymmetricLayout{<:Abst function AdaptiveCholeskyFactors(::SymmetricBandedLayouts, S::AbstractMatrix{T}) where T A = parent(S) l,u = bandwidths(A) - data = BandedMatrix{T}(undef,(0,0),(l,u)) # pad super + data = BandedMatrix(Zeros{T}(0,0),(l,u)) # pad super AdaptiveCholeskyFactors(CachedArray(data,A), 0) end AdaptiveCholeskyFactors(A::AbstractMatrix{T}) where T = AdaptiveCholeskyFactors(MemoryLayout(A), A) diff --git a/src/infqr.jl b/src/infqr.jl index 4f908f9..2fb9f30 100644 --- a/src/infqr.jl +++ b/src/infqr.jl @@ -87,6 +87,7 @@ end partialqr!(F::AdaptiveQRData{<:Any,<:BlockSkylineMatrix}, n::Int) = partialqr!(F, findblock(axes(F.data,2), n)) + struct AdaptiveQRFactors{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}} <: LayoutMatrix{T} data::AdaptiveQRData{T,DM,M} @@ -124,6 +125,11 @@ function getindex(F::AdaptiveQRFactors, k::Int, j::Int) F.data.data[k,j] end +function resizedata!(F::AdaptiveQRFactors, k::Int, j::Int) + partialqr!(F.data, j) + F +end + colsupport(F::QRPackedQ{<:Any,<:AdaptiveQRFactors}, j) = 1:last(colsupport(F.factors, j)) rowsupport(F::QRPackedQ{<:Any,<:AdaptiveQRFactors}, j) = first(rowsupport(F.factors, j)):size(F,2) @@ -150,6 +156,7 @@ function adaptiveqr(A) end qr_layout(::BandedLayouts, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A) +qr_layout(::AbstractBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A) qr_layout(::AbstractAlmostBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A) __qr_layout(_, ::NTuple{2,InfiniteCardinal{0}}, A) = adaptiveqr(A) qr_layout(::AbstractBlockBandedLayout, ::NTuple{2,InfiniteCardinal{0}}, A) = adaptiveqr(A) From 9b4e57d485d71bceaeedc35e92e1379e6c393df0 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Thu, 6 Feb 2025 07:49:33 +0000 Subject: [PATCH 16/18] Update tridiagonalconjugation.jl --- src/banded/tridiagonalconjugation.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/banded/tridiagonalconjugation.jl b/src/banded/tridiagonalconjugation.jl index ce3a80d..0661041 100644 --- a/src/banded/tridiagonalconjugation.jl +++ b/src/banded/tridiagonalconjugation.jl @@ -228,8 +228,8 @@ end function TridiagonalConjugationData(U, X, V) T = promote_type(typeof(inv(V[1, 1])), eltype(U), eltype(X)) # include inv so that we can't get Ints n_init = 100 - UX = Tridiagonal(zeros(T, n_init-1), zeros(T, n_init), zeros(T, n_init-1)) # zeros for BigFLoat - Y = Tridiagonal(zeros(T, n_init-1), zeros(T, n_init), zeros(T, n_init-1)) + UX = Tridiagonal(Vector{T}(undef, n_init-1), Vector{T}(undef, n_init), Vector{T}(undef, n_init-1)) + Y = Tridiagonal(Vector{T}(undef, n_init-1), Vector{T}(undef, n_init), Vector{T}(undef, n_init-1)) resizedata!(U, n_init, n_init) resizedata!(V, n_init, n_init) initiate_upper_mul_tri_triview!(UX, U, X) # fill-in 1st row From 21811a78315bc284d1830f1d250c666872211764 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sat, 8 Feb 2025 08:02:22 +0000 Subject: [PATCH 17/18] Support bidiagonal --- src/banded/tridiagonalconjugation.jl | 12 ++++++------ test/test_bidiagonalconjugation.jl | 8 ++++++-- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/src/banded/tridiagonalconjugation.jl b/src/banded/tridiagonalconjugation.jl index 0661041..9e46d3f 100644 --- a/src/banded/tridiagonalconjugation.jl +++ b/src/banded/tridiagonalconjugation.jl @@ -18,7 +18,7 @@ function upper_mul_tri_triview!(UX::Tridiagonal, U::BandedMatrix, X::Tridiagonal l,u = bandwidths(U) @assert size(U) == (n,n) - @assert l == 0 && u ≥ 2 + @assert l ≥ 0 # Tridiagonal bands can be resized @assert length(Xdl)+1 == length(Xd) == length(Xdu)+1 == length(UXdl)+1 == length(UXd) == length(UXdu)+1 == n @@ -42,7 +42,7 @@ function initiate_upper_mul_tri_triview!(UX, U::BandedMatrix, X) k = 1 aₖ, cₖ = Xd[1], Xdl[1] - Uₖₖ, Uₖₖ₊₁, Uₖₖ₊₂ = Udat[u+1,1], Udat[u,2], Udat[u-1,3] # U[k,k], U[k,k+1], U[k,k+2] + Uₖₖ, Uₖₖ₊₁, Uₖₖ₊₂ = Udat[u+1,1], Udat[u,2], (u > 1 ? Udat[u-1,3] : zero(eltype(Udat))) # U[k,k], U[k,k+1], U[k,k+2] UXd[1] = Uₖₖ*aₖ + Uₖₖ₊₁*cₖ # UX[k,k] = U[k,k]*X[k,k] + U[k,k+1]*X[k+1,k] bₖ, aₖ, cₖ, cₖ₋₁ = Xdu[1], Xd[2], Xdl[2], cₖ # X[k,k+1], X[k+1,k+1], X[k+2,k+1], X[k+1,k] UXdu[1] = Uₖₖ*bₖ + Uₖₖ₊₁*aₖ + Uₖₖ₊₂*cₖ # UX[k,k+1] = U[k,k]*X[k,k+1] + U[k,k+1]*X[k+1,k+1] + U[k,k+1]*X[k+1,k] @@ -71,7 +71,7 @@ function main_upper_mul_tri_triview!(UX, U::BandedMatrix, X, kr, bₖ=X[kr[1]-1, l,u = bandwidths(U) for k = kr - Uₖₖ, Uₖₖ₊₁, Uₖₖ₊₂ = Udat[u+1,k], Udat[u,k+1], Udat[u-1,k+2] # U[k,k], U[k,k+1], U[k,k+2] + Uₖₖ, Uₖₖ₊₁, Uₖₖ₊₂ = Udat[u+1,k], Udat[u,k+1], (u > 1 ? Udat[u-1,k+2] : zero(eltype(Udat))) # U[k,k], U[k,k+1], U[k,k+2] UXdl[k-1] = Uₖₖ*cₖ₋₁ # UX[k,k-1] = U[k,k]*X[k,k-1] UXd[k] = Uₖₖ*aₖ + Uₖₖ₊₁*cₖ # UX[k,k] = U[k,k]*X[k,k] + U[k,k+1]*X[k+1,k] bₖ, aₖ, cₖ, cₖ₋₁ = Xdu[k], Xd[k+1], Xdl[k+1], cₖ # X[k,k+1], X[k+1,k+1], X[k+2,k+1], X[k+1,k] @@ -118,7 +118,7 @@ function tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::BandedMatr l,u = bandwidths(R) @assert size(R) == (n,n) - @assert l == 0 && u ≥ 2 + @assert l ≥ 0 && u ≥ 1 # Tridiagonal bands can be resized @assert length(Xdl)+1 == length(Xd) == length(Xdu)+1 == length(Ydl)+1 == length(Yd) == length(Ydu)+1 == n @@ -177,12 +177,12 @@ function main_tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::Bande Rdat = R.data l,u = bandwidths(R) - @inbounds for k = kr + for k = kr cₖ₋₁,aₖ,bₖ = Xdl[k-1], Xd[k], Xdu[k] Ydl[k-1] = cₖ₋₁/Rₖₖ Yd[k] = aₖ-cₖ₋₁*Rₖₖ₊₁/Rₖₖ Ydu[k] = cₖ₋₁/Rₖₖ - Rₖₖ,Rₖₖ₊₁,Rₖ₋₁ₖ₊₁,Rₖ₋₁ₖ = Rdat[u+1,k], Rdat[u,k+1],Rdat[u-1,k+1],Rₖₖ₊₁ # R[k,k], R[k,k+1], R[k-1,k] + Rₖₖ,Rₖₖ₊₁,Rₖ₋₁ₖ₊₁,Rₖ₋₁ₖ = Rdat[u+1,k], Rdat[u,k+1],(u > 1 ? Rdat[u-1,k+1] : zero(eltype(Rdat))),Rₖₖ₊₁ # R[k,k], R[k,k+1], R[k-1,k+1] Yd[k] /= Rₖₖ Ydu[k-1] /= Rₖₖ Ydu[k] *= Rₖ₋₁ₖ*Rₖₖ₊₁/Rₖₖ - Rₖ₋₁ₖ₊₁ diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl index 8ff176b..0cd14ba 100644 --- a/test/test_bidiagonalconjugation.jl +++ b/test/test_bidiagonalconjugation.jl @@ -117,14 +117,18 @@ end (-(0:∞) ./ (2:2:∞))', ((2:∞) ./ (2:2:∞))'), ℵ₀, 0,2), LazyBandedMatrices.Tridiagonal((2:∞) ./ (3:2:∞), -1 ./ ((1:2:∞) .* (3:2:∞)), (1:∞) ./ (3:2:∞))), + (_BandedMatrix(Vcat( + (-(0:∞) ./ (2:2:∞))', + ((2:∞) ./ (2:2:∞))'), ℵ₀, 0,1), + LazyBandedMatrices.Tridiagonal((2:∞) ./ (3:2:∞), -1 ./ ((1:2:∞) .* (3:2:∞)), (1:∞) ./ (3:2:∞))) # P -> C^(5/2) (_BandedMatrix(Vcat((-3 ./ (3:2:∞))', Zeros(1,∞), (3 ./ (3:2:∞))'), ℵ₀, 0,2) * _BandedMatrix(Vcat((-1 ./ (1:2:∞))', Zeros(1,∞), (1 ./ (1:2:∞))'), ℵ₀, 0,2), LazyBandedMatrices.Tridiagonal((1:∞) ./ (1:2:∞), Zeros(∞), (1:∞) ./ (3: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 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 = InfiniteLinearAlgebra.upper_mul_tri_triview(U, X) @test Tridiagonal(U*X) ≈ UX # U*X*inv(U) only depends on Tridiagonal(U*X) From 332dde84ac5abb8febf483547b36699cf1557d96 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sat, 8 Feb 2025 10:18:37 +0000 Subject: [PATCH 18/18] fix tests --- src/infqr.jl | 3 ++- test/test_bidiagonalconjugation.jl | 14 +++++++------- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/infqr.jl b/src/infqr.jl index 2fb9f30..d7b42ac 100644 --- a/src/infqr.jl +++ b/src/infqr.jl @@ -158,9 +158,10 @@ end qr_layout(::BandedLayouts, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A) qr_layout(::AbstractBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A) qr_layout(::AbstractAlmostBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A) -__qr_layout(_, ::NTuple{2,InfiniteCardinal{0}}, A) = adaptiveqr(A) +_qr_layout(_, ::NTuple{2,InfiniteCardinal{0}}, A) = adaptiveqr(A) qr_layout(::AbstractBlockBandedLayout, ::NTuple{2,InfiniteCardinal{0}}, A) = adaptiveqr(A) factorize_layout(::BandedLayouts, ::NTuple{2,OneToInf{Int}}, A) = qr(A) +factorize_layout(::AbstractBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = qr(A) cache_layout(::TriangularLayout{UPLO, UNIT, <:AdaptiveLayout}, A::AbstractMatrix) where {UPLO, UNIT} = A # already cached diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl index 0cd14ba..52f528e 100644 --- a/test/test_bidiagonalconjugation.jl +++ b/test/test_bidiagonalconjugation.jl @@ -120,19 +120,19 @@ end (_BandedMatrix(Vcat( (-(0:∞) ./ (2:2:∞))', ((2:∞) ./ (2:2:∞))'), ℵ₀, 0,1), - LazyBandedMatrices.Tridiagonal((2:∞) ./ (3:2:∞), -1 ./ ((1:2:∞) .* (3:2:∞)), (1:∞) ./ (3:2:∞))) + LazyBandedMatrices.Tridiagonal((2:∞) ./ (3:2:∞), -1 ./ ((1:2:∞) .* (3:2:∞)), (1:∞) ./ (3:2:∞))), # P -> C^(5/2) (_BandedMatrix(Vcat((-3 ./ (3:2:∞))', Zeros(1,∞), (3 ./ (3:2:∞))'), ℵ₀, 0,2) * _BandedMatrix(Vcat((-1 ./ (1:2:∞))', Zeros(1,∞), (1 ./ (1:2:∞))'), ℵ₀, 0,2), LazyBandedMatrices.Tridiagonal((1:∞) ./ (1:2:∞), Zeros(∞), (1:∞) ./ (3: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 = InfiniteLinearAlgebra.upper_mul_tri_triview(U, X) + U = V = R[1:n,1:n] + X = Tridiagonal(Vector(X_T.dl[1:n-1]), Vector(X_T.d[1:n]), Vector(X_T.du[1:n-1])) + UX = InfiniteLinearAlgebra.upper_mul_tri_triview(U, X) @test Tridiagonal(U*X) ≈ UX # U*X*inv(U) only depends on Tridiagonal(U*X) - @time Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U) + Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U) @test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y data = TridiagonalConjugationData(R, X_T) @@ -162,8 +162,8 @@ end X_T = LazyBandedMatrices.Tridiagonal(Vcat(1/sqrt(2), Fill(1/2,∞)), Zeros(∞), Vcat(1/sqrt(2), Fill(1/2,∞))) 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])); + U = V = R[1:n,1:n]; + 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]