1+ """
2+ Store a diagonally-scaled Toeplitz∘Hankel matrix:
3+ DL(T∘H)DR
4+ where the Hankel matrix `H` is non-negative definite. This allows a Cholesky decomposition in 𝒪(K²N) operations and 𝒪(KN) storage, K = log N log ɛ⁻¹.
5+ """
6+ struct ToeplitzHankelPlan{S, N, M, N1, TP<: ToeplitzPlan{S,N1} } <: Plan{S}
7+ T:: NTuple{M,TP}
8+ L:: NTuple{M,Matrix{S}}
9+ R:: NTuple{M,Matrix{S}}
10+ tmp:: Array{S,N1}
11+ dims:: NTuple{M,Int}
12+ function ToeplitzHankelPlan {S,N,M,N1,TP} (T:: NTuple{M,TP} , L, R, dims) where {S,TP,N,N1,M}
13+ tmp = Array {S} (undef, max .(size .(T)... )... )
14+ new {S,N,M,N1,TP} (T, L, R, tmp, dims)
15+ end
16+ ToeplitzHankelPlan {S,N,M,N1,TP} (T:: NTuple{M,TP} , L, R, dims:: Int ) where {S,TP,N,N1,M} =
17+ ToeplitzHankelPlan {S,N,M,N1,TP} (T, L, R, (dims,))
18+ end
19+
20+ ToeplitzHankelPlan (T:: ToeplitzPlan{S,2} , L:: Matrix , R:: Matrix , dims= 1 ) where S =
21+ ToeplitzHankelPlan {S, 1, 1, 2, typeof(T)} ((T,), (L,), (R,), dims)
22+
23+ ToeplitzHankelPlan (T:: ToeplitzPlan{S,3} , L:: Matrix , R:: Matrix , dims) where S =
24+ ToeplitzHankelPlan {S, 2, 1,3, typeof(T)} ((T,), (L,), (R,), dims)
25+
26+ ToeplitzHankelPlan (T:: NTuple{2,TP} , L:: Tuple , R:: Tuple , dims) where {S,TP<: ToeplitzPlan{S,3} } =
27+ ToeplitzHankelPlan {S, 2,2,3, TP} (T, L, R, dims)
28+
29+
30+ function * (P:: ToeplitzHankelPlan{<:Any,1} , v:: AbstractVector )
31+ (R,),(L,),(T,),tmp = P. R,P. L,P. T,P. tmp
32+ tmp .= R .* v
33+ T * tmp
34+ tmp .= L .* tmp
35+ sum! (v, tmp)
36+ end
37+
38+ function _th_applymul1! (v, T, L, R, tmp)
39+ N = size (R,2 )
40+ m,n = size (v)
41+ tmp[1 : m,1 : n,1 : N] .= reshape (R,size (R,1 ),1 ,N) .* v
42+ T * view (tmp,1 : m,1 : n,1 : N)
43+ view (tmp,1 : m,1 : n,1 : N) .*= reshape (L,size (L,1 ),1 ,N)
44+ sum! (v, view (tmp,1 : m,1 : n,1 : N))
45+ end
46+
47+ function _th_applymul2! (v, T, L, R, tmp)
48+ N = size (R,2 )
49+ m,n = size (v)
50+ tmp[1 : m,1 : n,1 : N] .= reshape (R,1 ,size (R,1 ),N) .* v
51+ T * view (tmp,1 : m,1 : n,1 : N)
52+ view (tmp,1 : m,1 : n,1 : N) .*= reshape (L,1 ,size (L,1 ),N)
53+ sum! (v, view (tmp,1 : m,1 : n,1 : N))
54+ end
55+
56+
57+ function * (P:: ToeplitzHankelPlan{<:Any,2,1} , v:: AbstractMatrix )
58+ (R,),(L,),(T,),tmp = P. R,P. L,P. T,P. tmp
59+ if P. dims == (1 ,)
60+ _th_applymul1! (v, T, L, R, tmp)
61+ else
62+ _th_applymul2! (v, T, L, R, tmp)
63+ end
64+ v
65+ end
66+
67+ function * (P:: ToeplitzHankelPlan{<:Any,2,2} , v:: AbstractMatrix )
68+ (R1,R2),(L1,L2),(T1,T2),tmp = P. R,P. L,P. T,P. tmp
69+
70+ _th_applymul1! (v, T1, L1, R1, tmp)
71+ _th_applymul2! (v, T2, L2, R2, tmp)
72+
73+ v
74+ end
75+
76+ # partial cholesky for a Hankel matrix
77+
78+ function hankel_partialchol (v:: Vector{T} ) where T
79+ # Assumes positive definite
80+ σ = T[]
81+ n = isempty (v) ? 0 : (length (v)+ 2 ) ÷ 2
82+ C = Matrix {T} (undef, n, n)
83+ d = v[1 : 2 : end ] # diag of H
84+ @assert length (v) ≥ 2 n- 1
85+ reltol = maximum (abs,d)* eps (T)* log (n)
86+ r = 0
87+ for k = 1 : n
88+ mx,idx = findmax (d)
89+ if mx ≤ reltol break end
90+ push! (σ, inv (mx))
91+ C[:,k] .= view (v,idx: n+ idx- 1 )
92+ for j = 1 : k- 1
93+ nCjidxσj = - C[idx,j]* σ[j]
94+ LinearAlgebra. axpy! (nCjidxσj, view (C,:,j), view (C,:,k))
95+ end
96+ @inbounds for p= 1 : n
97+ d[p] -= C[p,k]^ 2 / mx
98+ end
99+ r += 1
100+ end
101+ for k= 1 : length (σ) rmul! (view (C,:,k), sqrt (σ[k])) end
102+ C[:,1 : r]
103+ end
104+
105+
106+ # Diagonally-scaled Toeplitz∘Hankel polynomial transforms
107+
108+
109+
110+ struct ChebyshevToLegendrePlanTH{TH}
111+ toeplitzhankel:: TH
112+ end
113+
114+ function * (P:: ChebyshevToLegendrePlanTH , v:: AbstractVector{S} ) where S
115+ n = length (v)
116+ ret = zero (S)
117+ @inbounds for k = 1 : 2 : n
118+ ret += - v[k]/ (k* (k- 2 ))
119+ end
120+ v[1 ] = ret
121+ P. toeplitzhankel* view (v,2 : n)
122+ v
123+ end
124+
125+ function _cheb2leg_rescale1! (V:: AbstractMatrix{S} ) where S
126+ m,n = size (V)
127+ for j = 1 : n
128+ ret = zero (S)
129+ @inbounds for k = 1 : 2 : m
130+ ret += - V[k,j]/ (k* (k- 2 ))
131+ end
132+ V[1 ,j] = ret
133+ end
134+ V
135+ end
136+
137+
138+ function * (P:: ChebyshevToLegendrePlanTH , V:: AbstractMatrix )
139+ m,n = size (V)
140+ dims = P. toeplitzhankel. dims
141+ if dims == (1 ,)
142+ _cheb2leg_rescale1! (V)
143+ P. toeplitzhankel* view (V,2 : m,:)
144+ elseif dims == (2 ,)
145+ _cheb2leg_rescale1! (transpose (V))
146+ P. toeplitzhankel* view (V,:,2 : n)
147+ else
148+ @assert dims == (1 ,2 )
149+ (R1,R2),(L1,L2),(T1,T2),tmp = P. toeplitzhankel. R,P. toeplitzhankel. L,P. toeplitzhankel. T,P. toeplitzhankel. tmp
150+
151+ _cheb2leg_rescale1! (V)
152+ _th_applymul1! (view (V,2 : m,:), T1, L1, R1, tmp)
153+ _cheb2leg_rescale1! (transpose (V))
154+ _th_applymul2! (view (V,:,2 : n), T2, L2, R2, tmp)
155+ end
156+ V
157+ end
158+
159+
160+
161+ function _leg2chebTH_TLC (:: Type{S} , mn, d) where S
162+ n = mn[d]
163+ λ = Λ .(0 : half (real (S)): n- 1 )
164+ t = zeros (S,n)
165+ t[1 : 2 : end ] .= 2 .* view (λ, 1 : 2 : n) ./ π
166+ C = hankel_partialchol (λ)
167+ T = plan_uppertoeplitz! (t, (mn... , size (C,2 )), d)
168+ L = copy (C)
169+ L[1 ,:] ./= 2
170+ T,L,C
171+ end
172+
173+ function _leg2chebuTH_TLC (:: Type{S} , mn, d) where {S}
174+ n = mn[d]
175+ S̃ = real (S)
176+ λ = Λ .(0 : half (S̃): n- 1 )
177+ t = zeros (S,n)
178+ t[1 : 2 : end ] = λ[1 : 2 : n]. / (((1 : 2 : n). - 2 ))
179+ h = λ./ ((1 : 2 n- 1 ). + 1 )
180+ C = hankel_partialchol (h)
181+ T = plan_uppertoeplitz! (- 2 t/ π, (length (t), size (C,2 )), 1 )
182+ (T, (1 : n) .* C, C)
183+ end
184+
185+
186+ for f in (:leg2cheb , :leg2chebu )
187+ plan = Symbol (" plan_th_" " !"
188+ TLC = Symbol (" _" " TH_TLC"
189+ @eval begin
190+ $ plan (:: Type{S} , mn:: Tuple , dims:: Int ) where {S} = ToeplitzHankelPlan ($ TLC (S, mn, dims)... , dims)
191+
192+ function $plan (:: Type{S} , mn:: NTuple{2,Int} , dims:: NTuple{2,Int} ) where {S}
193+ @assert dims == (1 ,2 )
194+ T1,L1,C1 = $ TLC (S, mn, 1 )
195+ T2,L2,C2 = $ TLC (S, mn, 2 )
196+ ToeplitzHankelPlan ((T1,T2), (L1,L2), (C1,C2), dims)
197+ end
198+ end
199+ end
200+
201+ _sub_dim_by_one (d) = ()
202+ _sub_dim_by_one (d, m, n... ) = (isone (d) ? m- 1 : m, _sub_dim_by_one (d- 1 , n... )... )
203+
204+ function _cheb2legTH_TLC (:: Type{S} , mn, d) where S
205+ n = mn[d]
206+ t = zeros (S,n- 1 )
207+ S̃ = real (S)
208+ if n > 1
209+ t[1 : 2 : end ] = Λ .(0 : one (S̃): div (n- 2 ,2 ), - half (S̃), one (S̃))
210+ end
211+ h = Λ .(1 : half (S̃): n- 1 , zero (S̃), 3 half (S̃))
212+ DL = (3 half (S̃): n- half (S̃))
213+ DR = - (one (S̃): n- one (S̃)). / 4
214+ C = hankel_partialchol (h)
215+ T = plan_uppertoeplitz! (t, (_sub_dim_by_one (d, mn... )... , size (C,2 )), d)
216+ T, DL .* C, DR .* C
217+ end
218+
219+ plan_th_cheb2leg! (:: Type{S} , mn:: Tuple , dims:: Int ) where {S} = ChebyshevToLegendrePlanTH (ToeplitzHankelPlan (_cheb2legTH_TLC (S, mn, dims)... , dims))
220+
221+ function plan_th_cheb2leg! (:: Type{S} , mn:: NTuple{2,Int} , dims:: NTuple{2,Int} ) where {S}
222+ @assert dims == (1 ,2 )
223+ T1,L1,C1 = _cheb2legTH_TLC (S, mn, 1 )
224+ T2,L2,C2 = _cheb2legTH_TLC (S, mn, 2 )
225+ ChebyshevToLegendrePlanTH (ToeplitzHankelPlan ((T1,T2), (L1,L2), (C1,C2), dims))
226+ end
227+
228+ function plan_th_ultra2ultra! (:: Type{S} , (n,):: Tuple{Int} , λ₁, λ₂) where {S}
229+ @assert abs (λ₁- λ₂) < 1
230+ S̃ = real (S)
231+ DL = (zero (S̃): n- one (S̃)) .+ λ₂
232+ jk = 0 : half (S̃): n- 1
233+ t = zeros (S,n)
234+ t[1 : 2 : n] = Λ .(jk,λ₁- λ₂,one (S̃))[1 : 2 : n]
235+ h = Λ .(jk,λ₁,λ₂+ one (S̃))
236+ lmul! (gamma (λ₂)/ gamma (λ₁),h)
237+ C = hankel_partialchol (h)
238+ T = plan_uppertoeplitz! (lmul! (inv (gamma (λ₁- λ₂)),t), (length (t), size (C,2 )), 1 )
239+ ToeplitzHankelPlan (T, DL .* C, C)
240+ end
241+
242+ function alternatesign! (v)
243+ @inbounds for k = 2 : 2 : length (v)
244+ v[k] = - v[k]
245+ end
246+ v
247+ end
248+
249+ function plan_th_jac2jac! (:: Type{S} , (n,), α, β, γ, δ) where {S}
250+ if β == δ
251+ @assert abs (α- γ) < 1
252+ @assert α+ β > - 1
253+ jk = 0 : n- 1
254+ DL = (2 jk .+ γ .+ β .+ 1 ). * Λ .(jk,γ+ β+ 1 ,β+ 1 )
255+ t = convert (AbstractVector{S}, Λ .(jk, α- γ,1 ))
256+ h = Λ .(0 : 2 n- 2 ,α+ β+ 1 ,γ+ β+ 2 )
257+ DR = Λ .(jk,β+ 1 ,α+ β+ 1 )./ gamma (α- γ)
258+ C = hankel_partialchol (h)
259+ T = plan_uppertoeplitz! (t, (length (t), size (C,2 )), 1 )
260+ elseif α == γ
261+ jk = 0 : n- 1
262+ DL = (2 jk .+ δ .+ α .+ 1 ). * Λ .(jk,δ+ α+ 1 ,α+ 1 )
263+ h = Λ .(0 : 2 n- 2 ,α+ β+ 1 ,δ+ α+ 2 )
264+ DR = Λ .(jk,α+ 1 ,α+ β+ 1 )./ gamma (β- δ)
265+ t = alternatesign! (convert (AbstractVector{S}, Λ .(jk,β- δ,1 )))
266+ C = hankel_partialchol (h)
267+ T = plan_uppertoeplitz! (t, (length (t), size (C,2 )), 1 )
268+ else
269+ throw (ArgumentError (" Cannot create Toeplitz dot Hankel, use a sequence of plans."
270+ end
271+
272+ ToeplitzHankelPlan (T, DL .* C, DR .* C)
273+ end
274+
275+ for f in (:th_leg2cheb , :th_cheb2leg , :th_leg2chebu )
276+ plan = Symbol (" plan_" " !"
277+ @eval begin
278+ $ plan (:: Type{S} , mn:: NTuple{N,Int} , dims:: UnitRange ) where {N,S} = $ plan (S, mn, tuple (dims... ))
279+ $ plan (:: Type{S} , mn:: Tuple{Int} , dims:: Tuple{Int} = (1 ,)) where {S} = $ plan (S, mn, dims... )
280+ $ plan (:: Type{S} , (m,n):: NTuple{2,Int} ) where {S} = $ plan (S, (m,n), (1 ,2 ))
281+ $ plan (arr:: AbstractArray{T} , dims... ) where T = $ plan (T, size (arr), dims... )
282+ $ f (v, dims... ) = $ plan (eltype (v), size (v), dims... )* copy (v)
283+ end
284+ end
285+
286+ th_ultra2ultra (v, λ₁, λ₂, dims... ) = plan_th_ultra2ultra! (eltype (v),size (v),λ₁,λ₂, dims... )* copy (v)
287+ th_jac2jac (v, α, β, γ, δ, dims... ) = plan_th_jac2jac! (eltype (v),size (v),α,β,γ,δ, dims... )* copy (v)
0 commit comments