|
| 1 | +""" |
| 2 | +cholesky_jacobimatrix(w, P) |
| 3 | +
|
| 4 | +returns the Jacobi matrix `X` associated to a quasi-matrix of polynomials |
| 5 | +orthogonal with respect to `w(x) w_p(x)` where `w_p(x)` is the weight of the polynomials in `P`. |
| 6 | +
|
| 7 | +The resulting polynomials are orthonormal on the same domain as `P`. The supplied `P` must be normalized. Accepted inputs are `w` as a function or `W` as an infinite matrix representing multiplication with the function `w` on the basis `P`. |
| 8 | +
|
| 9 | +An optional bool can be supplied, i.e. `cholesky_jacobimatrix(sqrtw, P, false)` to disable checks of symmetry for the weight multiplication matrix and orthonormality for the basis (use with caution). |
| 10 | +""" |
| 11 | +function cholesky_jacobimatrix(w::Function, P::OrthogonalPolynomial, checks::Bool = true) |
| 12 | + checks && !(P isa Normalized) && error("Polynomials must be orthonormal.") |
| 13 | + W = Symmetric(P \ (w.(axes(P,1)) .* P)) # Compute weight multiplication via Clenshaw |
| 14 | + return cholesky_jacobimatrix(W, P, false) # At this point checks already passed or were entered as false, no need to recheck |
| 15 | +end |
| 16 | +function cholesky_jacobimatrix(W::AbstractMatrix, P::OrthogonalPolynomial, checks::Bool = true) |
| 17 | + checks && !(P isa Normalized) && error("Polynomials must be orthonormal.") |
| 18 | + checks && !(W isa Symmetric) && error("Weight modification matrix must be symmetric.") |
| 19 | + return SymTridiagonal(CholeskyJacobiBands{:dv}(W,P),CholeskyJacobiBands{:ev}(W,P)) |
| 20 | +end |
| 21 | + |
| 22 | +# The generated Jacobi operators are symmetric tridiagonal, so we store their data in cached bands |
| 23 | +mutable struct CholeskyJacobiBands{dv,T} <: AbstractCachedVector{T} |
| 24 | + data::Vector{T} # store band entries, :dv for diagonal, :ev for off-diagonal |
| 25 | + U::UpperTriangular{T} # store upper triangular conversion matrix (needed to extend available entries) |
| 26 | + UX::ApplyArray{T} # store U*X, where X is the Jacobi matrix of the original P (needed to extend available entries) |
| 27 | + datasize::Int # size of so-far computed block |
| 28 | +end |
| 29 | + |
| 30 | +# Computes the initial data for the Jacobi operator bands |
| 31 | +function CholeskyJacobiBands{:dv}(W, P::OrthogonalPolynomial{T}) where T |
| 32 | + U = cholesky(W).U |
| 33 | + X = jacobimatrix(P) |
| 34 | + UX = ApplyArray(*,U,X) |
| 35 | + dv = zeros(T,2) # compute a length 2 vector on first go |
| 36 | + dv[1] = dot(view(UX,1,1), U[1,1] \ [one(T)]) |
| 37 | + dv[2] = dot(view(UX,2,1:2), U[1:2,1:2] \ [zero(T); one(T)]) |
| 38 | + return CholeskyJacobiBands{:dv,T}(dv, U, UX, 2) |
| 39 | +end |
| 40 | +function CholeskyJacobiBands{:ev}(W, P::OrthogonalPolynomial{T}) where T |
| 41 | + U = cholesky(W).U |
| 42 | + X = jacobimatrix(P) |
| 43 | + UX = ApplyArray(*,U,X) |
| 44 | + ev = zeros(T,2) # compute a length 2 vector on first go |
| 45 | + ev[1] = dot(view(UX,1,1:2), U[1:2,1:2] \ [zero(T); one(T)]) |
| 46 | + ev[2] = dot(view(UX,2,1:3), U[1:3,1:3] \ [zeros(T,2); one(T)]) |
| 47 | + return CholeskyJacobiBands{:ev,T}(ev, U, UX, 2) |
| 48 | +end |
| 49 | + |
| 50 | +size(::CholeskyJacobiBands) = (ℵ₀,) # Stored as an infinite cached vector |
| 51 | + |
| 52 | +# Resize and filling functions for cached implementation |
| 53 | +function resizedata!(K::CholeskyJacobiBands, nm::Integer) |
| 54 | + νμ = K.datasize |
| 55 | + if nm > νμ |
| 56 | + resize!(K.data,nm) |
| 57 | + cache_filldata!(K, νμ:nm) |
| 58 | + K.datasize = nm |
| 59 | + end |
| 60 | + K |
| 61 | +end |
| 62 | +function cache_filldata!(J::CholeskyJacobiBands{:dv,T}, inds::UnitRange{Int}) where T |
| 63 | + # pre-fill U and UX to prevent expensive step-by-step filling in of cached U and UX in the loop |
| 64 | + getindex(J.U,inds[end]+1,inds[end]+1) |
| 65 | + getindex(J.UX,inds[end]+1,inds[end]+1) |
| 66 | + |
| 67 | + ek = [zero(T); one(T)] |
| 68 | + @inbounds for k in inds |
| 69 | + J.data[k] = dot(view(J.UX,k,k-1:k), J.U[k-1:k,k-1:k] \ ek) |
| 70 | + end |
| 71 | +end |
| 72 | +function cache_filldata!(J::CholeskyJacobiBands{:ev, T}, inds::UnitRange{Int}) where T |
| 73 | + # pre-fill U and UX to prevent expensive step-by-step filling in of cached U and UX in the loop |
| 74 | + getindex(J.U,inds[end]+1,inds[end]+1) |
| 75 | + getindex(J.UX,inds[end]+1,inds[end]+1) |
| 76 | + |
| 77 | + ek = [zeros(T,2); one(T)] |
| 78 | + @inbounds for k in inds |
| 79 | + J.data[k] = dot(view(J.UX,k,k-1:k+1), J.U[k-1:k+1,k-1:k+1] \ ek) |
| 80 | + end |
| 81 | +end |
| 82 | + |
| 83 | + |
| 84 | +""" |
| 85 | +qr_jacobimatrix(sqrtw, P) |
| 86 | +
|
| 87 | +returns the Jacobi matrix `X` associated to a quasi-matrix of polynomials |
| 88 | +orthogonal with respect to `w(x) w_p(x)` where `w_p(x)` is the weight of the polynomials in `P`. |
| 89 | +
|
| 90 | +The resulting polynomials are orthonormal on the same domain as `P`. The supplied `P` must be normalized. Accepted inputs for `sqrtw` are the square root of the weight modification as a function or `sqrtW` as an infinite matrix representing multiplication with the function `sqrt(w)` on the basis `P`. |
| 91 | +
|
| 92 | +An optional bool can be supplied, i.e. `qr_jacobimatrix(sqrtw, P, false)` to disable checks of symmetry for the weight multiplication matrix and orthonormality for the basis (use with caution). |
| 93 | +""" |
| 94 | +function qr_jacobimatrix(sqrtw::Function, P::OrthogonalPolynomial, checks::Bool = true) |
| 95 | + checks && !(P isa Normalized) && error("Polynomials must be orthonormal.") |
| 96 | + sqrtW = (P \ (sqrtw.(axes(P,1)) .* P)) # Compute weight multiplication via Clenshaw |
| 97 | + return qr_jacobimatrix(sqrtW, P, false) # At this point checks already passed or were entered as false, no need to recheck |
| 98 | +end |
| 99 | +function qr_jacobimatrix(sqrtW::AbstractMatrix, P::OrthogonalPolynomial, checks::Bool = true) |
| 100 | + checks && !(P isa Normalized) && error("Polynomials must be orthonormal.") |
| 101 | + checks && !(sqrtW isa Symmetric) && error("Weight modification matrix must be symmetric.") |
| 102 | + K = SymTridiagonal(QRJacobiBands{:dv}(sqrtW,P),QRJacobiBands{:ev}(sqrtW,P)) |
| 103 | + return K |
| 104 | +end |
| 105 | + |
| 106 | +# The generated Jacobi operators are symmetric tridiagonal, so we store their data in cached bands |
| 107 | +mutable struct QRJacobiBands{dv,T} <: AbstractCachedVector{T} |
| 108 | + data::Vector{T} # store band entries, :dv for diagonal, :ev for off-diagonal |
| 109 | + U::ApplyArray{T} # store upper triangular conversion matrix (needed to extend available entries) |
| 110 | + UX::ApplyArray{T} # store U*X, where X is the Jacobi matrix of the original P (needed to extend available entries) |
| 111 | + datasize::Int # size of so-far computed block |
| 112 | +end |
| 113 | + |
| 114 | +# Computes the initial data for the Jacobi operator bands |
| 115 | +function QRJacobiBands{:dv}(sqrtW, P::OrthogonalPolynomial{T}) where T |
| 116 | + U = qr(sqrtW).R |
| 117 | + U = ApplyArray(*,Diagonal(sign.(view(U,band(0)))),U) |
| 118 | + X = jacobimatrix(P) |
| 119 | + UX = ApplyArray(*,U,X) |
| 120 | + dv = zeros(T,2) # compute a length 2 vector on first go |
| 121 | + dv[1] = dot(view(UX,1,1), U[1,1] \ [one(T)]) |
| 122 | + dv[2] = dot(view(UX,2,1:2), U[1:2,1:2] \ [zero(T); one(T)]) |
| 123 | + return QRJacobiBands{:dv,T}(dv, U, UX, 2) |
| 124 | +end |
| 125 | +function QRJacobiBands{:ev}(sqrtW, P::OrthogonalPolynomial{T}) where T |
| 126 | + U = qr(sqrtW).R |
| 127 | + U = ApplyArray(*,Diagonal(sign.(view(U,band(0)))),U) |
| 128 | + X = jacobimatrix(P) |
| 129 | + UX = ApplyArray(*,U,X) |
| 130 | + ev = zeros(T,2) # compute a length 2 vector on first go |
| 131 | + ev[1] = dot(view(UX,1,1:2), U[1:2,1:2] \ [zero(T); one(T)]) |
| 132 | + ev[2] = dot(view(UX,2,1:3), U[1:3,1:3] \ [zeros(T,2); one(T)]) |
| 133 | + return QRJacobiBands{:ev,T}(ev, U, UX, 2) |
| 134 | +end |
| 135 | + |
| 136 | +size(::QRJacobiBands) = (ℵ₀,) # Stored as an infinite cached vector |
| 137 | + |
| 138 | +# Resize and filling functions for cached implementation |
| 139 | +function resizedata!(K::QRJacobiBands, nm::Integer) |
| 140 | + νμ = K.datasize |
| 141 | + if nm > νμ |
| 142 | + resize!(K.data,nm) |
| 143 | + cache_filldata!(K, νμ:nm) |
| 144 | + K.datasize = nm |
| 145 | + end |
| 146 | + K |
| 147 | +end |
| 148 | +function cache_filldata!(J::QRJacobiBands{:dv,T}, inds::UnitRange{Int}) where T |
| 149 | + # pre-fill U and UX to prevent expensive step-by-step filling in of cached U and UX in the loop |
| 150 | + getindex(J.U,inds[end]+1,inds[end]+1) |
| 151 | + getindex(J.UX,inds[end]+1,inds[end]+1) |
| 152 | + |
| 153 | + ek = [zero(T); one(T)] |
| 154 | + @inbounds for k in inds |
| 155 | + J.data[k] = dot(view(J.UX,k,k-1:k), J.U[k-1:k,k-1:k] \ ek) |
| 156 | + end |
| 157 | +end |
| 158 | +function cache_filldata!(J::QRJacobiBands{:ev, T}, inds::UnitRange{Int}) where T |
| 159 | + # pre-fill U and UX to prevent expensive step-by-step filling in of cached U and UX in the loop |
| 160 | + getindex(J.U,inds[end]+1,inds[end]+1) |
| 161 | + getindex(J.UX,inds[end]+1,inds[end]+1) |
| 162 | + |
| 163 | + ek = [zeros(T,2); one(T)] |
| 164 | + @inbounds for k in inds |
| 165 | + J.data[k] = dot(view(J.UX,k,k-1:k+1), J.U[k-1:k+1,k-1:k+1] \ ek) |
| 166 | + end |
| 167 | +end |
0 commit comments