Skip to content

Commit 186a218

Browse files
dlfivefiftydlfivefifty
committed
Support caching
1 parent 0780e9a commit 186a218

File tree

5 files changed

+98
-66
lines changed

5 files changed

+98
-66
lines changed

‎Project.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,7 @@ FillArrays = "1.0"
2626
InfiniteArrays = "0.15"
2727
InfiniteRandomArrays = "0.2"
2828
Infinities = "0.1"
29-
LazyArrays = "2.3"
29+
LazyArrays = "2.5"
3030
LazyBandedMatrices = "0.11"
3131
LinearAlgebra = "1"
3232
MatrixFactorizations = "3.0"

‎src/InfiniteLinearAlgebra.jl‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,7 @@ import Base.Broadcast: BroadcastStyle, Broadcasted, broadcasted
1414

1515
import ArrayLayouts: AbstractBandedLayout, AbstractQLayout, AdjQRPackedQLayout, CNoPivot, DenseColumnMajor, FillLayout,
1616
MatLdivVec, MatLmulMat, MatLmulVec, MemoryLayout, QRPackedQLayout, RangeCumsum, TriangularLayout,
17-
TridiagonalLayout, __qr, _factorize, _qr, check_mul_axes, colsupport,
17+
TridiagonalLayout, _qr_layout, factorize_layout, qr_layout, check_mul_axes, colsupport,
1818
diagonaldata, ldiv!, lmul!, mul, mulreduce, reflector!, reflectorApply!,
1919
rowsupport, sub_materialize, subdiagonaldata, sublayout, supdiagonaldata, transposelayout,
2020
triangulardata, triangularlayout, zero!, materialize!
@@ -40,7 +40,7 @@ import LazyArrays: AbstractCachedMatrix, AbstractCachedVector, AbstractLazyLayou
4040
CachedArray, CachedLayout, CachedMatrix, CachedVector, LazyArrayStyle, LazyLayout,
4141
LazyLayouts, LazyMatrix, LazyVector, AbstractPaddedLayout, PaddedColumns, _broadcast_sub_arguments,
4242
applybroadcaststyle, applylayout, arguments, cacheddata, paddeddata, resizedata!, simplifiable,
43-
simplify, islazy, islazy_layout, cache_getindex
43+
simplify, islazy, islazy_layout, cache_getindex, cache_layout
4444

4545
import LazyBandedMatrices: AbstractLazyBandedBlockBandedLayout, AbstractLazyBandedLayout, ApplyBandedLayout, BlockVec,
4646
BroadcastBandedLayout, KronTravBandedBlockBandedLayout, LazyBandedLayout,

‎src/banded/tridiagonalconjugation.jl‎

Lines changed: 43 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,11 @@ function upper_mul_tri_triview!(UX::Tridiagonal, U::BandedMatrix, X::Tridiagonal
2828
end
2929

3030
# populate first row of UX with UX
31-
function initiate_upper_mul_tri_triview!(UX, U, X)
31+
32+
initiate_upper_mul_tri_triview!(UX, U::UpperTriangular, X) = initiate_upper_mul_tri_triview!(UX, parent(U), X)
33+
initiate_upper_mul_tri_triview!(UX, U::CachedMatrix, X) = initiate_upper_mul_tri_triview!(UX, U.data, X)
34+
35+
function initiate_upper_mul_tri_triview!(UX, U::BandedMatrix, X)
3236
Xdl, Xd, Xdu = X.dl, X.d, X.du
3337
UXdl, UXd, UXdu = UX.dl, UX.d, UX.du
3438
Udat = U.data
@@ -46,7 +50,12 @@ function initiate_upper_mul_tri_triview!(UX, U, X)
4650
end
4751

4852
# fills in the rows kr of UX
49-
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])
53+
function main_upper_mul_tri_triview!(UX, U::CachedMatrix, X, kr, kwds...)
54+
resizedata!(U, kr[end], kr[end]+2)
55+
main_upper_mul_tri_triview!(UX, U.data, X, kr, kwds...)
56+
end
57+
58+
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])
5059
Xdl, Xd, Xdu = X.dl, X.d, X.du
5160
UXdl, UXd, UXdu = UX.dl, UX.d, UX.du
5261
Udat = U.data
@@ -109,8 +118,14 @@ function tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::BandedMatr
109118
finalize_tri_mul_invupper_triview!(Y, X, R, n, Rₖₖ, Rₖₖ₊₁)
110119
end
111120

112-
# populate first row of X/R
113-
function initiate_tri_mul_invupper_triview!(Y, X, R)
121+
# partially-populate first row of X/R
122+
# Ydu[k] is updated below
123+
function initiate_tri_mul_invupper_triview!(Y, X, R::CachedMatrix)
124+
resizedata!(R, 1, 2)
125+
initiate_tri_mul_invupper_triview!(Y, X, R.data)
126+
end
127+
128+
function initiate_tri_mul_invupper_triview!(Y, X, R::BandedMatrix)
114129
Xdl, Xd, Xdu = X.dl, X.d, X.du
115130
Ydl, Yd, Ydu = Y.dl, Y.d, Y.du
116131
Rdat = R.data
@@ -120,15 +135,21 @@ function initiate_tri_mul_invupper_triview!(Y, X, R)
120135
k = 1
121136
aₖ,bₖ = Xd[k], Xdu[k]
122137
Rₖₖ,Rₖₖ₊₁ = Rdat[u+1,k], Rdat[u,k+1] # R[1,1], R[1,2]
138+
123139
Yd[k] = aₖ/Rₖₖ
124140
Ydu[k] = bₖ - aₖ * Rₖₖ₊₁/Rₖₖ
125141

126142
Y, Rₖₖ, Rₖₖ₊₁
127143
end
128144

129145

130-
# populate rows kr of X/R
131-
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])
146+
# populate rows kr of X/R. Ydu[k] is wrong until next run.
147+
function main_tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::CachedMatrix, kr, kwds...)
148+
resizedata!(R, kr[end], kr[end]+1)
149+
main_tri_mul_invupper_triview!(Y, X, R.data, kr, kwds...)
150+
end
151+
152+
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)])
132153
Xdl, Xd, Xdu = X.dl, X.d, X.du
133154
Ydl, Yd, Ydu = Y.dl, Y.d, Y.du
134155
Rdat = R.data
@@ -187,11 +208,18 @@ function TridiagonalConjugationData(U, X, V)
187208
n_init = 100
188209
UX = Tridiagonal(Vector{T}(undef, n_init-1), Vector{T}(undef, n_init), Vector{T}(undef, n_init-1))
189210
Y = Tridiagonal(Vector{T}(undef, n_init-1), Vector{T}(undef, n_init), Vector{T}(undef, n_init-1))
211+
resizedata!(U, n_init, n_init)
212+
resizedata!(V, n_init, n_init)
190213
initiate_upper_mul_tri_triview!(UX, U, X) # fill-in 1st row
191-
return TridiagonalConjugationData(U, X, V, UX, Y, 1)
214+
initiate_tri_mul_invupper_triview!(Y, UX, V)
215+
return TridiagonalConjugationData(U, X, V, UX, Y, 0)
192216
end
193217

194-
TridiagonalConjugationData(U, X) = TridiagonalConjugationData(U, X, U)
218+
219+
function TridiagonalConjugationData(U, X)
220+
C = cache(U)
221+
TridiagonalConjugationData(C, X, C)
222+
end
195223

196224
copy(data::TridiagonalConjugationData) = TridiagonalConjugationData(copy(data.U), copy(data.X), copy(data.V), copy(data.UX), copy(data.Y), data.datasize)
197225

@@ -209,8 +237,13 @@ function resizedata!(data::TridiagonalConjugationData, n)
209237
resize!(data.Y.du, 2n)
210238
end
211239

212-
main_upper_mul_tri_triview!(data.UX, data.U, data.X, data.datasize+1:n)
213240

214-
data.datasize = n
241+
if n > data.datasize
242+
main_upper_mul_tri_triview!(data.UX, data.U, data.X, data.datasize+2:n+1)
243+
main_tri_mul_invupper_triview!(data.Y, data.UX, data.U, data.datasize+2:n+1) # need one extra as it updates first row
244+
data.datasize = n
245+
end
246+
247+
data
215248
end
216249

‎src/infqr.jl‎

Lines changed: 8 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -149,11 +149,14 @@ function adaptiveqr(A)
149149
QR(AdaptiveQRFactors(data), AdaptiveQRTau(data))
150150
end
151151

152-
_qr(::AbstractBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A)
153-
_qr(::AbstractAlmostBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A)
154-
__qr(_, ::NTuple{2,InfiniteCardinal{0}}, A) = adaptiveqr(A)
155-
_qr(::AbstractBlockBandedLayout, ::NTuple{2,InfiniteCardinal{0}}, A) = adaptiveqr(A)
156-
_factorize(::AbstractBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = qr(A)
152+
qr_layout(::AbstractBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A)
153+
qr_layout(::AbstractAlmostBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = adaptiveqr(A)
154+
__qr_layout(_, ::NTuple{2,InfiniteCardinal{0}}, A) = adaptiveqr(A)
155+
qr_layout(::AbstractBlockBandedLayout, ::NTuple{2,InfiniteCardinal{0}}, A) = adaptiveqr(A)
156+
factorize_layout(::AbstractBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = qr(A)
157+
158+
159+
cache_layout(::TriangularLayout{UPLO, UNIT, <:AdaptiveLayout}, A::AbstractMatrix) where {UPLO, UNIT} = A # already cached
157160

158161
partialqr!(F::QR, n) = partialqr!(F.factors, n)
159162
partialqr!(F::AdaptiveQRFactors, n) = partialqr!(F.data, n)

‎test/test_bidiagonalconjugation.jl‎

Lines changed: 44 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@ using InfiniteLinearAlgebra: BidiagonalConjugation, OneToInf
33
using ArrayLayouts: supdiagonaldata, subdiagonaldata, diagonaldata
44
using LinearAlgebra
55
using LazyArrays: LazyLayout
6+
using BandedMatrices: _BandedMatrix
67

78
@testset "BidiagonalConjugation" begin
89
@test InfiniteLinearAlgebra._to_uplo('U') == 'U'
@@ -102,11 +103,25 @@ end
102103

103104

104105
@testset "TridiagonalConjugation" begin
105-
@testset "T -> U" begin
106-
R = BandedMatrices._BandedMatrix(Vcat(-Ones(1,∞)/2,
107-
Zeros(1,∞),
108-
Hcat(Ones(1,1),Ones(1,∞)/2)), ℵ₀, 0,2)
109-
X_T = LazyBandedMatrices.Tridiagonal(Vcat(1.0, Fill(1/2,∞)), Zeros(∞), Fill(1/2,∞))
106+
for (R,X_T) in (
107+
# T -> U
108+
(_BandedMatrix(Vcat(-Ones(1,∞)/2, Zeros(1,∞), Hcat(Ones(1,1),Ones(1,∞)/2)), ℵ₀, 0,2),
109+
LazyBandedMatrices.Tridiagonal(Vcat(1.0, Fill(1/2,∞)), Zeros(∞), Fill(1/2,∞))),
110+
# P -> C^(3/2)
111+
(_BandedMatrix(Vcat((-1 ./ (1:2:∞))',
112+
Zeros(1,∞),
113+
(1 ./ (1:2:∞))'), ℵ₀, 0,2),
114+
LazyBandedMatrices.Tridiagonal((1:∞) ./ (1:2:∞), Zeros(∞), (1:∞) ./ (3:2:∞))),
115+
# P^(1,0) -> P^(2,0)
116+
(_BandedMatrix(Vcat(Zeros(1,∞), # extra band since code assumes two bands
117+
(-(0:∞) ./ (2:2:∞))',
118+
((2:∞) ./ (2:2:∞))'), ℵ₀, 0,2),
119+
LazyBandedMatrices.Tridiagonal((2:∞) ./ (3:2:∞), -1 ./ ((1:2:∞) .* (3:2:∞)), (1:∞) ./ (3:2:∞))),
120+
# P -> C^(5/2)
121+
(_BandedMatrix(Vcat((-3 ./ (3:2:∞))', Zeros(1,∞), (3 ./ (3:2:∞))'), ℵ₀, 0,2) *
122+
_BandedMatrix(Vcat((-1 ./ (1:2:∞))', Zeros(1,∞), (1 ./ (1:2:∞))'), ℵ₀, 0,2),
123+
LazyBandedMatrices.Tridiagonal((1:∞) ./ (1:2:∞), Zeros(∞), (1:∞) ./ (3:2:∞)))
124+
)
110125
n = 1000
111126
@time U = V = R[1:n,1:n];
112127
@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
116131
@time Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U)
117132
@test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y
118133

119-
InfiniteLinearAlgebra.TridiagonalConjugationData(U, X, U)
120-
end
121-
@testset "P -> Ultraspherical(3/2)" begin
122-
R = BandedMatrices._BandedMatrix(Vcat((-1 ./ (1:2:∞))',
123-
Zeros(1,∞),
124-
(1 ./ (1:2:∞))'), ℵ₀, 0,2)
125-
X_P = LazyBandedMatrices.Tridiagonal((1:∞) ./ (1:2:∞), Zeros(∞), (1:∞) ./ (3:2:∞))
126-
n = 1000
127-
@time U = V = R[1:n,1:n]
128-
@time X = Tridiagonal(Vector(X_P.dl[1:n-1]), Vector(X_P.d[1:n]), Vector(X_P.du[1:n-1]))
129-
@time UX = InfiniteLinearAlgebra.upper_mul_tri_triview(U, X)
130-
@test Tridiagonal(U*X) ≈ UX
131-
# U*X*inv(U) only depends on Tridiagonal(U*X)
132-
@time Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U)
133-
@test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y
134-
end
134+
data = InfiniteLinearAlgebra.TridiagonalConjugationData(R, X_T)
135+
@test data.UX[1,:] ≈ UX[1,1:100]
136+
InfiniteLinearAlgebra.resizedata!(data, 1)
137+
@test data.UX[1:2,:] ≈ UX[1:2,1:100]
138+
@test data.Y[1,:] ≈ Y[1,1:100]
139+
InfiniteLinearAlgebra.resizedata!(data, 2)
140+
@test data.UX[1:3,:] ≈ UX[1:3,1:100]
141+
@test data.Y[1:2,:] ≈ Y[1:2,1:100]
142+
InfiniteLinearAlgebra.resizedata!(data, 3)
143+
@test data.UX[1:4,:] ≈ UX[1:4,1:100]
144+
@test data.Y[1:3,:] ≈ Y[1:3,1:100]
145+
InfiniteLinearAlgebra.resizedata!(data, 1000)
146+
@test data.UX[1:999,1:999] ≈ UX[1:999,1:999]
147+
@test data.Y[1:999,1:999] ≈ Y[1:999,1:999]
135148

136-
@testset "Jacobi(1,0) -> Jacobi(2,0)" begin
137-
R = BandedMatrices._BandedMatrix(Vcat(Zeros(1,∞), # extra band since code assumes two bands
138-
(-(0:∞) ./ (2:2:∞))',
139-
((2:∞) ./ (2:2:∞))'), ℵ₀, 0,2)
140-
X_P = LazyBandedMatrices.Tridiagonal((2:∞) ./ (3:2:∞), -1 ./ ((1:2:∞) .* (3:2:∞)), (1:∞) ./ (3:2:∞))
141-
n = 1000
142-
@time U = V = R[1:n,1:n]
143-
@time X = Tridiagonal(Vector(X_P.dl[1:n-1]), Vector(X_P.d[1:n]), Vector(X_P.du[1:n-1]))
144-
@time UX = InfiniteLinearAlgebra.upper_mul_tri_triview(U, X)
145-
@test Tridiagonal(U*X) ≈ UX
146-
# U*X*inv(U) only depends on Tridiagonal(U*X)
147-
@time Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U)
148-
@test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y
149+
data = InfiniteLinearAlgebra.TridiagonalConjugationData(R, X_T)
150+
InfiniteLinearAlgebra.resizedata!(data, 1000)
151+
@test data.UX[1:999,1:999] ≈ UX[1:999,1:999]
152+
@test data.Y[1:999,1:999] ≈ Y[1:999,1:999]
149153
end
150-
151-
@testset "Legendre() -> Jacobi(5/2)" begin
152-
R = BandedMatrices._BandedMatrix(Vcat((-3 ./ (3:2:∞))', Zeros(1,∞), (3 ./ (3:2:∞))'), ℵ₀, 0,2) *
153-
BandedMatrices._BandedMatrix(Vcat((-1 ./ (1:2:∞))', Zeros(1,∞), (1 ./ (1:2:∞))'), ℵ₀, 0,2)
154-
X_P = LazyBandedMatrices.Tridiagonal((1:∞) ./ (1:2:∞), Zeros(∞), (1:∞) ./ (3:2:∞))
155-
156-
n = 1000
157-
@time U = V = R[1:n,1:n]
158-
@time X = Tridiagonal(Vector(X_P.dl[1:n-1]), Vector(X_P.d[1:n]), Vector(X_P.du[1:n-1]))
159-
@time UX = InfiniteLinearAlgebra.upper_mul_tri_triview(U, X)
160-
@test Tridiagonal(U*X) ≈ UX
161-
# U*X*inv(U) only depends on Tridiagonal(U*X)
162-
@time Y = InfiniteLinearAlgebra.tri_mul_invupper_triview(UX, U)
163-
@test Tridiagonal(U*X / U) ≈ Tridiagonal(UX / U) ≈ Y
154+
155+
@testset "Cholesky" begin
156+
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))
157+
R = cholesky(M).U
158+
X_T = LazyBandedMatrices.Tridiagonal(Vcat(1/sqrt(2), Fill(1/2,∞)), Zeros(∞), Vcat(1/sqrt(2), Fill(1/2,∞)))
159+
data = InfiniteLinearAlgebra.TridiagonalConjugationData(R, X_T)
164160
end
165161
end

0 commit comments

Comments
 (0)