From 77e7910a422abcaf86d2c593df92ae41f2c127c7 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 17 Oct 2023 23:11:50 +0100 Subject: [PATCH] jac2jac matrix --- src/toeplitzhankel.jl | 47 +++++++++++++++++++++++++++---------- test/toeplitzhankeltests.jl | 11 ++++++++- 2 files changed, 44 insertions(+), 14 deletions(-) diff --git a/src/toeplitzhankel.jl b/src/toeplitzhankel.jl index ad3ff062..75e4ef88 100644 --- a/src/toeplitzhankel.jl +++ b/src/toeplitzhankel.jl @@ -279,14 +279,6 @@ end # th_ultra2ultra ### -plan_th_ultra2ultra!(::Type{S}, mn, λ₁, λ₂, dims::Int) where {S} = ToeplitzHankelPlan(_ultra2ultraTH_TLC(S, mn, λ₁, λ₂, dims)..., dims) - -function plan_th_ultra2ultra!(::Type{S}, mn::NTuple{2,Int}, λ₁, λ₂, dims::NTuple{2,Int}) where {S} - @assert dims == (1,2) - T1,L1,C1 = _ultra2ultraTH_TLC(S, mn, λ₁, λ₂, 1) - T2,L2,C2 = _ultra2ultraTH_TLC(S, mn, λ₁, λ₂, 2) - ToeplitzHankelPlan((T1,T2), (L1,L2), (C1,C2), dims) -end function _ultra2ultraTH_TLC(::Type{S}, mn, λ₁, λ₂, d) where {S} n = mn[d] @@ -303,6 +295,16 @@ function _ultra2ultraTH_TLC(::Type{S}, mn, λ₁, λ₂, d) where {S} T, DL .* C, C end +plan_th_ultra2ultra!(::Type{S}, mn, λ₁, λ₂, dims::Int) where {S} = ToeplitzHankelPlan(_ultra2ultraTH_TLC(S, mn, λ₁, λ₂, dims)..., dims) + +function plan_th_ultra2ultra!(::Type{S}, mn::NTuple{2,Int}, λ₁, λ₂, dims::NTuple{2,Int}) where {S} + @assert dims == (1,2) + T1,L1,C1 = _ultra2ultraTH_TLC(S, mn, λ₁, λ₂, 1) + T2,L2,C2 = _ultra2ultraTH_TLC(S, mn, λ₁, λ₂, 2) + ToeplitzHankelPlan((T1,T2), (L1,L2), (C1,C2), dims) +end + + ### # th_jac2jac ### @@ -314,7 +316,8 @@ function alternatesign!(v) v end -function plan_th_jac2jac!(::Type{S}, (n,), α, β, γ, δ) where {S} +function _jac2jacTH_TLC(::Type{S}, mn, α, β, γ, δ, d) where {S} + n = mn[d] if β == δ @assert abs(α-γ) < 1 @assert α+β > -1 @@ -324,7 +327,7 @@ function plan_th_jac2jac!(::Type{S}, (n,), α, β, γ, δ) where {S} h = Λ.(0:2n-2,α+β+1,γ+β+2) DR = Λ.(jk,β+1,α+β+1)./gamma(α-γ) C = hankel_partialchol(h) - T = plan_uppertoeplitz!(t, (length(t), size(C,2)), 1) + T = plan_uppertoeplitz!(t, (mn..., size(C,2)), d) elseif α == γ jk = 0:n-1 DL = (2jk .+ δ .+ α .+ 1).*Λ.(jk,δ+α+1,α+1) @@ -332,14 +335,28 @@ function plan_th_jac2jac!(::Type{S}, (n,), α, β, γ, δ) where {S} DR = Λ.(jk,α+1,α+β+1)./gamma(β-δ) t = alternatesign!(convert(AbstractVector{S}, Λ.(jk,β-δ,1))) C = hankel_partialchol(h) - T = plan_uppertoeplitz!(t, (length(t), size(C,2)), 1) + T = plan_uppertoeplitz!(t, (mn..., size(C,2)), d) else throw(ArgumentError("Cannot create Toeplitz dot Hankel, use a sequence of plans.")) end - ToeplitzHankelPlan(T, DL .* C, DR .* C) + (T, DL .* C, DR .* C) +end + +plan_th_jac2jac!(::Type{S}, mn, α, β, γ, δ, dims::Int) where {S} = ToeplitzHankelPlan(_jac2jacTH_TLC(S, mn, α, β, γ, δ, dims)..., dims) + +function plan_th_jac2jac!(::Type{S}, mn::NTuple{2,Int}, α, β, γ, δ, dims::NTuple{2,Int}) where {S} + @assert dims == (1,2) + T1,L1,C1 = _jac2jacTH_TLC(S, mn, α, β, γ, δ, 1) + T2,L2,C2 = _jac2jacTH_TLC(S, mn, α, β, γ, δ, 2) + ToeplitzHankelPlan((T1,T2), (L1,L2), (C1,C2), dims) end + +### +# other routines +### + for f in (:th_leg2cheb, :th_cheb2leg, :th_leg2chebu) plan = Symbol("plan_", f, "!") @eval begin @@ -357,4 +374,8 @@ plan_th_ultra2ultra!(::Type{S}, (m,n)::NTuple{2,Int}, λ₁, λ₂) where {S} = plan_th_ultra2ultra!(arr::AbstractArray{T}, λ₁, λ₂, dims...) where T = plan_th_ultra2ultra!(T, size(arr), λ₁, λ₂, dims...) th_ultra2ultra(v, λ₁, λ₂, dims...) = plan_th_ultra2ultra!(eltype(v), size(v), λ₁, λ₂, dims...)*copy(v) -th_jac2jac(v, α, β, γ, δ, dims...) = plan_th_jac2jac!(eltype(v),size(v),α,β,γ,δ, dims...)*copy(v) \ No newline at end of file +plan_th_jac2jac!(::Type{S}, mn::NTuple{N,Int}, α, β, γ, δ, dims::UnitRange) where {N,S} = plan_th_jac2jac!(S, mn, α, β, γ, δ, tuple(dims...)) +plan_th_jac2jac!(::Type{S}, mn::Tuple{Int}, α, β, γ, δ, dims::Tuple{Int}=(1,)) where {S} = plan_th_jac2jac!(S, mn, α, β, γ, δ, dims...) +plan_th_jac2jac!(::Type{S}, (m,n)::NTuple{2,Int}, α, β, γ, δ) where {S} = plan_th_jac2jac!(S, (m,n), α, β, γ, δ, (1,2)) +plan_th_jac2jac!(arr::AbstractArray{T}, α, β, γ, δ, dims...) where T = plan_th_jac2jac!(T, size(arr), α, β, γ, δ, dims...) +th_jac2jac(v, α, β, γ, δ, dims...) = plan_th_jac2jac!(eltype(v), size(v), α, β, γ, δ, dims...)*copy(v) \ No newline at end of file diff --git a/test/toeplitzhankeltests.jl b/test/toeplitzhankeltests.jl index 55d08a58..5d194c69 100644 --- a/test/toeplitzhankeltests.jl +++ b/test/toeplitzhankeltests.jl @@ -1,7 +1,7 @@ using FastTransforms, Test import FastTransforms: th_leg2cheb, th_cheb2leg, th_leg2chebu, th_ultra2ultra,th_jac2jac, th_leg2chebu, lib_leg2cheb, lib_cheb2leg, lib_ultra2ultra, lib_jac2jac, - plan_th_cheb2leg!, plan_th_leg2chebu!, plan_th_leg2cheb!, plan_th_ultra2ultra! + plan_th_cheb2leg!, plan_th_leg2chebu!, plan_th_leg2cheb!, plan_th_ultra2ultra!, plan_th_jac2jac! @testset "ToeplitzHankel" begin for x in ([1.0], [1.0,2,3,4,5], [1.0+im,2-3im,3+4im,4-5im,5+10im], collect(1.0:1000)) @@ -48,6 +48,15 @@ import FastTransforms: th_leg2cheb, th_cheb2leg, th_leg2chebu, th_ultra2ultra,th @test th_ultra2ultra(X, 0.1, 0.6) == plan_th_ultra2ultra!(X, 0.1, 0.6, 1:2)*copy(X) @test th_ultra2ultra(th_ultra2ultra(X, 0.1, 0.6), 0.6, 0.1) ≈ X + + @test th_jac2jac(X, 0.1, 0.6, 0.1, 0.8, 1) ≈ hcat([jac2jac(X[:,j], 0.1, 0.6, 0.1, 0.8) for j=1:size(X,2)]...) + @test th_jac2jac(X, 0.1, 0.6, 0.1, 0.8, 2) ≈ vcat([permutedims(jac2jac(X[k,:], 0.1, 0.6, 0.1, 0.8)) for k=1:size(X,1)]...) + @test th_jac2jac(X, 0.1, 0.6, 0.1, 0.8) ≈ th_jac2jac(th_jac2jac(X, 0.1, 0.6, 0.1, 0.8, 1), 0.1, 0.6, 0.1, 0.8, 2) + + @test th_jac2jac(X, 0.1, 0.6, 0.1, 0.8) == plan_th_jac2jac!(X, 0.1, 0.6, 0.1, 0.8, 1:2)*copy(X) + @test th_jac2jac(X, 0.1, 0.6, 0.1, 0.8) == plan_th_jac2jac!(X, 0.1, 0.6, 0.1, 0.8, 1:2)*copy(X) + + @test th_jac2jac(th_jac2jac(X, 0.1, 0.6, 0.1, 0.8), 0.1, 0.8, 0.1, 0.6) ≈ X end @testset "BigFloat" begin