1
1
# This currently takes the weight multiplication operator as input.
2
2
# I will probably change this to take the weight function instead.
3
- function cholesky_jacobimatrix (W)
3
+ function cholesky_jacobimatrix (W:: Symmetric )
4
4
bands = CholeskyJacobiBands (W) # the cached array only needs to store two bands bc of symmetry
5
5
return SymTridiagonal (bands[1 ,:],bands[2 ,:])
6
6
end
7
7
8
8
# The generated Jacobi operators are symmetric tridiagonal, so we store their data as two bands
9
9
mutable struct CholeskyJacobiBands{T} <: AbstractCachedMatrix{T}
10
- data
11
- U
12
- X
10
+ data:: Matrix{T}
11
+ U:: UpperTriangular
12
+ X:: Symmetric{T}
13
13
datasize:: Int
14
14
array
15
15
end
18
18
symmjacobim (J:: SymTridiagonal ) = Symmetric (BandedMatrix (0 => J. dv, 1 => J. ev))
19
19
20
20
# Computes the initial data for the Jacobi operator bands
21
- function CholeskyJacobiBands (W)
22
- T = eltype (W)
21
+ function CholeskyJacobiBands (W:: Symmetric{T} ) where T
23
22
U = cholesky (W). U
24
23
X = symmjacobim (jacobimatrix (Normalized (Legendre ()[affine (zero (T).. one (T),Inclusion (- one (T).. one (T))),:])))
25
24
dat = zeros (T,2 ,10 )
@@ -37,15 +36,15 @@ function resizedata!(K::CholeskyJacobiBands, nm::Integer)
37
36
νμ = K. datasize
38
37
if nm > νμ
39
38
olddata = copy (K. data)
40
- K. data = similar (K. data, nm , nm)
39
+ K. data = similar (K. data, 2 , nm)
41
40
K. data[axes (olddata)... ] = olddata
42
41
inds = νμ:nm
43
42
cache_filldata! (K, inds)
44
43
K. datasize = nm
45
44
end
46
45
K
47
46
end
48
- function cache_filldata! (J:: CholeskyJacobiBands , inds)
47
+ function cache_filldata! (J:: CholeskyJacobiBands , inds:: UnitRange{Int} )
49
48
for k in inds
50
49
J. data[1 ,k] = (J. U * (J. X * (J. U \ [zeros (k- 1 ); 1 ; zeros (∞)])))[k]
51
50
J. data[2 ,k] = (J. U * (J. X * (J. U \ [zeros (k); 1 ; zeros (∞)])))[k]
@@ -56,7 +55,7 @@ function getindex(K::CholeskyJacobiBands, k::Integer, j::Integer)
56
55
resizedata! (K, max (k,j))
57
56
K. data[k, j]
58
57
end
59
- function getindex (K:: CholeskyJacobiBands , kr:: Integer , jr:: UnitRange{Integer } )
58
+ function getindex (K:: CholeskyJacobiBands , kr:: Integer , jr:: UnitRange{Int } )
60
59
resizedata! (K, maximum (jr))
61
60
K. data[kr, jr]
62
61
end
0 commit comments