Skip to content

Commit 84517d2

Browse files
dlfivefiftydlfivefifty
authored
TridiagonalConjugation (#194)
* Add Chebyshev example * Start TridiagonalConjugation * Fast Tridiagonal(U*X) * Fast finite tri-conj * Update test_bidiagonalconjugation.jl * Update test_bidiagonalconjugation.jl * Allow general bandwiddths in R * split out pieces * adaptively populate UX * split tri_mul_invupper in 3 * Support caching * Cholesky R * Add Sym/TridiagonalConjugation * allow SymTridiagonal * support QR * Update tridiagonalconjugation.jl * Support bidiagonal * fix tests
1 parent 0155e48 commit 84517d2

File tree

6 files changed

+446
-16
lines changed

6 files changed

+446
-16
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: 8 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,12 +40,17 @@ 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,
4747
_block_interlace_axes, _krontrav_axes, krontravargs
4848

49+
const StructuredLayoutTypes{Lay} = Union{SymmetricLayout{Lay}, HermitianLayout{Lay}, TriangularLayout{'L','N',Lay}, TriangularLayout{'U','N',Lay}, TriangularLayout{'L','U',Lay}, TriangularLayout{'U','U',Lay}}
50+
51+
const BandedLayouts = Union{AbstractBandedLayout, StructuredLayoutTypes{<:AbstractBandedLayout}}
52+
53+
4954
import LinearAlgebra: AbstractQ, AdjointQ, AdjOrTrans, factorize, matprod, qr
5055

5156
import MatrixFactorizations: AdjQLPackedQLayout, LayoutQ, QL, QLPackedQ, QLPackedQLayout, QR, QRPackedQ,
@@ -125,5 +130,6 @@ include("infqr.jl")
125130
include("inful.jl")
126131
include("infcholesky.jl")
127132
include("banded/bidiagonalconjugation.jl")
133+
include("banded/tridiagonalconjugation.jl")
128134

129135
end # module

‎src/banded/tridiagonalconjugation.jl‎

Lines changed: 303 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,303 @@
1+
2+
# upper_mul_tri_triview(U, X) == Tridiagonal(U*X) where U is Upper triangular BandedMatrix and X is Tridiagonal
3+
function upper_mul_tri_triview(U::BandedMatrix, X::Tridiagonal)
4+
T = promote_type(eltype(U), eltype(X))
5+
n = size(U,1)
6+
UX = Tridiagonal(Vector{T}(undef, n-1), Vector{T}(undef, n), Vector{T}(undef, n-1))
7+
8+
upper_mul_tri_triview!(UX, U, X)
9+
end
10+
11+
function upper_mul_tri_triview!(UX::Tridiagonal, U::BandedMatrix, X::Tridiagonal)
12+
n = size(UX,1)
13+
14+
15+
Xdl, Xd, Xdu = X.dl, X.d, X.du
16+
UXdl, UXd, UXdu = UX.dl, UX.d, UX.du
17+
18+
l,u = bandwidths(U)
19+
20+
@assert size(U) == (n,n)
21+
@assert l ≥ 0
22+
# Tridiagonal bands can be resized
23+
@assert length(Xdl)+1 == length(Xd) == length(Xdu)+1 == length(UXdl)+1 == length(UXd) == length(UXdu)+1 == n
24+
25+
UX, bₖ, aₖ, cₖ, cₖ₋₁ = initiate_upper_mul_tri_triview!(UX, U, X)
26+
UX, bₖ, aₖ, cₖ, cₖ₋₁ = main_upper_mul_tri_triview!(UX, U, X, 2:n-2, bₖ, aₖ, cₖ, cₖ₋₁)
27+
finalize_upper_mul_tri_triview!(UX, U, X, n-1, bₖ, aₖ, cₖ, cₖ₋₁)
28+
end
29+
30+
# populate first row of UX with UX
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+
initiate_upper_mul_tri_triview!(UX, U::Union{AdaptiveCholeskyFactors,AdaptiveQRFactors}, X) = initiate_upper_mul_tri_triview!(UX, U.data.data, X)
35+
36+
function initiate_upper_mul_tri_triview!(UX, U::BandedMatrix, X)
37+
Xdl, Xd, Xdu = subdiagonaldata(X), diagonaldata(X), supdiagonaldata(X)
38+
UXdl, UXd, UXdu = UX.dl, UX.d, UX.du
39+
Udat = U.data
40+
41+
l,u = bandwidths(U)
42+
43+
k = 1
44+
aₖ, cₖ = Xd[1], Xdl[1]
45+
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]
46+
UXd[1] = Uₖₖ*aₖ + Uₖₖ₊₁*cₖ # UX[k,k] = U[k,k]*X[k,k] + U[k,k+1]*X[k+1,k]
47+
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]
48+
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]
49+
50+
UX, bₖ, aₖ, cₖ, cₖ₋₁
51+
end
52+
53+
# fills in the rows kr of UX
54+
main_upper_mul_tri_triview!(UX, U::UpperTriangular, X, kr, kwds...) = main_upper_mul_tri_triview!(UX, parent(U), X, kr, kwds...)
55+
56+
function main_upper_mul_tri_triview!(UX, U::Union{CachedMatrix,AdaptiveCholeskyFactors}, X, kr, kwds...)
57+
resizedata!(U, kr[end], kr[end]+2)
58+
main_upper_mul_tri_triview!(UX, U.data, X, kr, kwds...)
59+
end
60+
61+
function main_upper_mul_tri_triview!(UX, U::AdaptiveQRFactors, X, kr, kwds...)
62+
resizedata!(U, kr[end], kr[end]+2)
63+
main_upper_mul_tri_triview!(UX, U.data.data, X, kr, kwds...)
64+
end
65+
66+
67+
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])
68+
Xdl, Xd, Xdu = subdiagonaldata(X), diagonaldata(X), supdiagonaldata(X)
69+
UXdl, UXd, UXdu = UX.dl, UX.d, UX.du
70+
Udat = U.data
71+
l,u = bandwidths(U)
72+
73+
for k = kr
74+
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]
75+
UXdl[k-1] = Uₖₖ*cₖ₋₁ # UX[k,k-1] = U[k,k]*X[k,k-1]
76+
UXd[k] = Uₖₖ*aₖ + Uₖₖ₊₁*cₖ # UX[k,k] = U[k,k]*X[k,k] + U[k,k+1]*X[k+1,k]
77+
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]
78+
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]
79+
end
80+
81+
UX, bₖ, aₖ, cₖ, cₖ₋₁
82+
end
83+
84+
# populate rows k and k+1 of UX, assuming we are at the bottom-right
85+
function finalize_upper_mul_tri_triview!(UX, U, X, k, bₖ, aₖ, cₖ, cₖ₋₁)
86+
Xdl, Xd, Xdu = X.dl, X.d, X.du
87+
UXdl, UXd, UXdu = UX.dl, UX.d, UX.du
88+
Udat = U.data
89+
l,u = bandwidths(U)
90+
91+
Uₖₖ, Uₖₖ₊₁ = Udat[u+1,k], Udat[u,k+1] # U[k,k], U[k,k+1]
92+
UXdl[k-1] = Uₖₖ*cₖ₋₁ # UX[k,k-1] = U[k,k]*X[k,k-1]
93+
UXd[k] = Uₖₖ*aₖ + Uₖₖ₊₁*cₖ # UX[k,k] = U[k,k]*X[k,k] + U[k,k+1]*X[k+1,k]
94+
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]
95+
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]
96+
97+
k += 1
98+
Uₖₖ = Udat[u+1,k] # U[k,k]
99+
UXdl[k-1] = Uₖₖ*cₖ₋₁ # UX[k,k-1] = U[k,k]*X[k,k-1]
100+
UXd[k] = Uₖₖ*aₖ # UX[k,k] = U[k,k]*X[k,k] + U[k,k+1]*X[k+1,k]
101+
102+
UX
103+
end
104+
105+
106+
# X*R^{-1} = X*[1/R₁₁ -R₁₂/(R₁₁R₂₂) -R₁₃/R₂₂ …
107+
# 0 1/R₂₂ -R₂₃/R₃₃
108+
# 1/R₃₃
109+
110+
tri_mul_invupper_triview(X::Tridiagonal, R::BandedMatrix) = tri_mul_invupper_triview!(similar(X, promote_type(eltype(X), eltype(R))), X, R)
111+
112+
113+
function tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::BandedMatrix)
114+
n = size(X,1)
115+
Xdl, Xd, Xdu = X.dl, X.d, X.du
116+
Ydl, Yd, Ydu = Y.dl, Y.d, Y.du
117+
118+
l,u = bandwidths(R)
119+
120+
@assert size(R) == (n,n)
121+
@assert l ≥ 0 && u ≥ 1
122+
# Tridiagonal bands can be resized
123+
@assert length(Xdl)+1 == length(Xd) == length(Xdu)+1 == length(Ydl)+1 == length(Yd) == length(Ydu)+1 == n
124+
125+
UX, Rₖₖ, Rₖₖ₊₁ = initiate_tri_mul_invupper_triview!(Y, X, R)
126+
UX, Rₖₖ, Rₖₖ₊₁ = main_tri_mul_invupper_triview!(Y, X, R, 2:n-1, Rₖₖ, Rₖₖ₊₁)
127+
finalize_tri_mul_invupper_triview!(Y, X, R, n, Rₖₖ, Rₖₖ₊₁)
128+
end
129+
130+
# partially-populate first row of X/R
131+
# Ydu[k] is updated below
132+
function initiate_tri_mul_invupper_triview!(Y, X, R::CachedMatrix)
133+
resizedata!(R, 1, 2)
134+
initiate_tri_mul_invupper_triview!(Y, X, R.data)
135+
end
136+
137+
function initiate_tri_mul_invupper_triview!(Y, X, R::Union{AdaptiveCholeskyFactors,AdaptiveQRFactors})
138+
resizedata!(R, 1, 2)
139+
initiate_tri_mul_invupper_triview!(Y, X, R.data.data)
140+
end
141+
142+
initiate_tri_mul_invupper_triview!(Y, X, R::UpperTriangular) = initiate_tri_mul_invupper_triview!(Y, X, parent(R))
143+
144+
function initiate_tri_mul_invupper_triview!(Y, X, R::BandedMatrix)
145+
Xdl, Xd, Xdu = X.dl, X.d, X.du
146+
Ydl, Yd, Ydu = Y.dl, Y.d, Y.du
147+
Rdat = R.data
148+
149+
l,u = bandwidths(R)
150+
151+
k = 1
152+
aₖ,bₖ = Xd[k], Xdu[k]
153+
Rₖₖ,Rₖₖ₊₁ = Rdat[u+1,k], Rdat[u,k+1] # R[1,1], R[1,2]
154+
155+
Yd[k] = aₖ/Rₖₖ
156+
Ydu[k] = bₖ - aₖ * Rₖₖ₊₁/Rₖₖ
157+
158+
Y, Rₖₖ, Rₖₖ₊₁
159+
end
160+
161+
162+
# populate rows kr of X/R. Ydu[k] is wrong until next run.
163+
main_tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::UpperTriangular, kr, kwds...) = main_tri_mul_invupper_triview!(Y, X, parent(R), kr, kwds...)
164+
function main_tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::Union{AdaptiveCholeskyFactors,CachedMatrix}, kr, kwds...)
165+
resizedata!(R, kr[end], kr[end]+1)
166+
main_tri_mul_invupper_triview!(Y, X, R.data, kr, kwds...)
167+
end
168+
169+
function main_tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::AdaptiveQRFactors, kr, kwds...)
170+
resizedata!(R, kr[end], kr[end]+1)
171+
main_tri_mul_invupper_triview!(Y, X, R.data.data, kr, kwds...)
172+
end
173+
174+
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)])
175+
Xdl, Xd, Xdu = X.dl, X.d, X.du
176+
Ydl, Yd, Ydu = Y.dl, Y.d, Y.du
177+
Rdat = R.data
178+
l,u = bandwidths(R)
179+
180+
for k = kr
181+
cₖ₋₁,aₖ,bₖ = Xdl[k-1], Xd[k], Xdu[k]
182+
Ydl[k-1] = cₖ₋₁/Rₖₖ
183+
Yd[k] = aₖ-cₖ₋₁*Rₖₖ₊₁/Rₖₖ
184+
Ydu[k] = cₖ₋₁/Rₖₖ
185+
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]
186+
Yd[k] /= Rₖₖ
187+
Ydu[k-1] /= Rₖₖ
188+
Ydu[k] *= Rₖ₋₁ₖ*Rₖₖ₊₁/Rₖₖ - Rₖ₋₁ₖ₊₁
189+
Ydu[k] += bₖ - aₖ * Rₖₖ₊₁ / Rₖₖ
190+
end
191+
Y, Rₖₖ, Rₖₖ₊₁
192+
end
193+
194+
195+
# populate row k of X/R, assuming we are at the bottom-right
196+
function finalize_tri_mul_invupper_triview!(Y::Tridiagonal, X::Tridiagonal, R::BandedMatrix, k, Rₖₖ=R[k-1,k-1], Rₖₖ₊₁=R[k-1,k])
197+
Xdl, Xd, Xdu = X.dl, X.d, X.du
198+
Ydl, Yd, Ydu = Y.dl, Y.d, Y.du
199+
Rdat = R.data
200+
l,u = bandwidths(R)
201+
cₖ₋₁,aₖ = Xdl[k-1], Xd[k]
202+
Ydl[k-1] = cₖ₋₁/Rₖₖ
203+
Yd[k] = aₖ-cₖ₋₁*Rₖₖ₊₁/Rₖₖ
204+
Rₖₖ = Rdat[u+1,k] # R[k,k]
205+
Yd[k] /= Rₖₖ
206+
Ydu[k-1] /= Rₖₖ
207+
208+
Y
209+
end
210+
"""
211+
TridiagonalConjugationData(U, X, V, Y)
212+
213+
caches the infinite dimensional Tridiagonal(U*X/V)
214+
in the tridiagonal matrix `Y`
215+
"""
216+
217+
mutable struct TridiagonalConjugationData{T}
218+
const U::AbstractMatrix{T}
219+
const X::AbstractMatrix{T}
220+
const V::AbstractMatrix{T}
221+
222+
const UX::Tridiagonal{T,Vector{T}} # cache Tridiagonal(U*X)
223+
const Y::Tridiagonal{T,Vector{T}} # cache Tridiagonal(U*X/V)
224+
225+
datasize::Int
226+
end
227+
228+
function TridiagonalConjugationData(U, X, V)
229+
T = promote_type(typeof(inv(V[1, 1])), eltype(U), eltype(X)) # include inv so that we can't get Ints
230+
n_init = 100
231+
UX = Tridiagonal(Vector{T}(undef, n_init-1), Vector{T}(undef, n_init), Vector{T}(undef, n_init-1))
232+
Y = Tridiagonal(Vector{T}(undef, n_init-1), Vector{T}(undef, n_init), Vector{T}(undef, n_init-1))
233+
resizedata!(U, n_init, n_init)
234+
resizedata!(V, n_init, n_init)
235+
initiate_upper_mul_tri_triview!(UX, U, X) # fill-in 1st row
236+
initiate_tri_mul_invupper_triview!(Y, UX, V)
237+
return TridiagonalConjugationData(U, X, V, UX, Y, 0)
238+
end
239+
240+
241+
function TridiagonalConjugationData(U, X)
242+
C = cache(U)
243+
TridiagonalConjugationData(C, X, C)
244+
end
245+
246+
copy(data::TridiagonalConjugationData) = TridiagonalConjugationData(copy(data.U), copy(data.X), copy(data.V), copy(data.UX), copy(data.Y), data.datasize)
247+
248+
249+
function resizedata!(data::TridiagonalConjugationData, n)
250+
n ≤ data.datasize && return data
251+
252+
if n ≥ length(data.UX.dl) # Avoid O(n²) growing. Note min(length(dv), length(ev)) == length(ev)
253+
resize!(data.UX.dl, 2n)
254+
resize!(data.UX.d, 2n + 1)
255+
resize!(data.UX.du, 2n)
256+
257+
resize!(data.Y.dl, 2n)
258+
resize!(data.Y.d, 2n + 1)
259+
resize!(data.Y.du, 2n)
260+
end
261+
262+
263+
if n > data.datasize
264+
main_upper_mul_tri_triview!(data.UX, data.U, data.X, data.datasize+2:n+1)
265+
main_tri_mul_invupper_triview!(data.Y, data.UX, data.U, data.datasize+2:n+1) # need one extra as it updates first row
266+
data.datasize = n
267+
end
268+
269+
data
270+
end
271+
272+
struct TridiagonalConjugationBand{T} <: LazyVector{T}
273+
data::TridiagonalConjugationData{T}
274+
diag::Symbol
275+
end
276+
277+
size(P::TridiagonalConjugationBand) = (ℵ₀,)
278+
resizedata!(A::TridiagonalConjugationBand, n) = resizedata!(A.data, n)
279+
280+
function _triconj_getindex(C::TridiagonalConjugationBand, I)
281+
resizedata!(C, maximum(I)+1)
282+
getfield(C.data.Y, C.diag)[I]
283+
end
284+
285+
getindex(A::TridiagonalConjugationBand, I::Integer) = _triconj_getindex(A, I)
286+
getindex(A::TridiagonalConjugationBand, I::AbstractVector) = _triconj_getindex(A, I)
287+
getindex(K::TridiagonalConjugationBand, k::AbstractInfUnitRange{<:Integer}) = view(K, k)
288+
getindex(K::SubArray{<:Any,1,<:TridiagonalConjugationBand}, k::AbstractInfUnitRange{<:Integer}) = view(K, k)
289+
290+
copy(A::TridiagonalConjugationBand) = A # immutable
291+
292+
293+
const TridiagonalConjugation{T} = Tridiagonal{T, TridiagonalConjugationBand{T}}
294+
const SymTridiagonalConjugation{T} = SymTridiagonal{T, TridiagonalConjugationBand{T}}
295+
function TridiagonalConjugation(R, X, Y...)
296+
data = TridiagonalConjugationData(R, X, Y...)
297+
Tridiagonal(TridiagonalConjugationBand(data, :dl), TridiagonalConjugationBand(data, :d), TridiagonalConjugationBand(data, :du))
298+
end
299+
300+
function SymTridiagonalConjugation(R, X, Y...)
301+
data = TridiagonalConjugationData(R, X, Y...)
302+
SymTridiagonal(TridiagonalConjugationBand(data, :d), TridiagonalConjugationBand(data, :du))
303+
end

‎src/infcholesky.jl‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ const SymmetricBandedLayouts = Union{SymTridiagonalLayout,SymmetricLayout{<:Abst
1212
function AdaptiveCholeskyFactors(::SymmetricBandedLayouts, S::AbstractMatrix{T}) where T
1313
A = parent(S)
1414
l,u = bandwidths(A)
15-
data = BandedMatrix{T}(undef,(0,0),(l,u)) # pad super
15+
data = BandedMatrix(Zeros{T}(0,0),(l,u)) # pad super
1616
AdaptiveCholeskyFactors(CachedArray(data,A), 0)
1717
end
1818
AdaptiveCholeskyFactors(A::AbstractMatrix{T}) where T = AdaptiveCholeskyFactors(MemoryLayout(A), A)
@@ -41,6 +41,10 @@ function partialcholesky!(F::AdaptiveCholeskyFactors{T,<:BandedMatrix}, n::Int)
4141
F
4242
end
4343

44+
resizedata!(F::AdaptiveCholeskyFactors, m, n) = partialcholesky!(F, n) # support cache interface
45+
46+
resizedata!(R::UpperTriangular{<:Any,<:AdaptiveCholeskyFactors}, m...) = resizedata!(parent(R), m...)
47+
4448
function getindex(F::AdaptiveCholeskyFactors, k::Int, j::Int)
4549
partialcholesky!(F, max(k,j))
4650
F.data.data[k,j]

0 commit comments

Comments
 (0)