Skip to content

Commit acfefb6

Browse files
MikaelSlevinskyMikaelSlevinsky
committed
add associated jac2jac, update references
1 parent cfae853 commit acfefb6

5 files changed

+98
-30
lines changed

‎Project.toml

+2-2
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
name = "FastTransforms"
22
uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c"
3-
version = "0.11.0"
3+
version = "0.11.1"
44

55
[deps]
66
AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
@@ -25,7 +25,7 @@ BinaryProvider = "0.5"
2525
DSP = "0.6"
2626
FFTW = "1"
2727
FastGaussQuadrature = "0.4"
28-
FastTransforms_jll = "0.4.0"
28+
FastTransforms_jll = "0.4.1"
2929
FillArrays = "0.9, 0.10"
3030
Reexport = "0.2"
3131
SpecialFunctions = "0.10, 1"

‎README.md

+5-3
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,7 @@ julia> using FastTransforms, LinearAlgebra
1919

2020
## Fast orthogonal polynomial transforms
2121

22-
The 33 orthogonal polynomial transforms are listed in `FastTransforms.kind2string.(0:32)`. Univariate transforms may be planned with the standard normalization or with orthonormalization. For multivariate transforms, the standard normalization may be too severe for floating-point computations, so it is omitted. Here are two examples:
22+
The 34 orthogonal polynomial transforms are listed in `FastTransforms.kind2string.(0:33)`. Univariate transforms may be planned with the standard normalization or with orthonormalization. For multivariate transforms, the standard normalization may be too severe for floating-point computations, so it is omitted. Here are two examples:
2323

2424
### The Chebyshev--Legendre transform
2525

@@ -159,6 +159,8 @@ julia> @time norm(ipaduatransform(paduatransform(v)) - v)/norm(v)
159159

160160
[1] D. Ruiz—Antolín and A. Townsend. <a href="https://doi.org/10.1137/17M1134822">A nonuniform fast Fourier transform based on low rank approximation</a>, *SIAM J. Sci. Comput.*, **40**:A529–A547, 2018.
161161

162-
[2] R. M. Slevinsky. <a href="https://doi.org/10.1016/j.acha.2017.11.001">Fast and backward stable transforms between spherical harmonic expansions and bivariate Fourier series</a>, *Appl. Comput. Harmon. Anal.*, **47**:585—606, 2019.
162+
[2] S. Olver, R. M. Slevinsky, and A. Townsend. <a href="https://doi.org/10.1017/S0962492920000045">Fast algorithms using orthogonal polynomials</a>, *Acta Numerica*, **29**:573—699, 2020.
163163

164-
[3] R. M. Slevinsky, <a href="https://arxiv.org/abs/1711.07866">Conquering the pre-computation in two-dimensional harmonic polynomial transforms</a>, arXiv:1711.07866, 2017.
164+
[3] R. M. Slevinsky. <a href="https://doi.org/10.1016/j.acha.2017.11.001">Fast and backward stable transforms between spherical harmonic expansions and bivariate Fourier series</a>, *Appl. Comput. Harmon. Anal.*, **47**:585—606, 2019.
165+
166+
[4] R. M. Slevinsky, <a href="https://arxiv.org/abs/1711.07866">Conquering the pre-computation in two-dimensional harmonic polynomial transforms</a>, arXiv:1711.07866, 2017.

‎src/FastTransforms.jl

+2-2
Original file line numberDiff line numberDiff line change
@@ -34,13 +34,13 @@ import LinearAlgebra: mul!, lmul!, ldiv!
3434

3535
export leg2cheb, cheb2leg, ultra2ultra, jac2jac,
3636
lag2lag, jac2ultra, ultra2jac, jac2cheb,
37-
cheb2jac, ultra2cheb, cheb2ultra,
37+
cheb2jac, ultra2cheb, cheb2ultra, associatedjac2jac,
3838
sph2fourier, sphv2fourier, disk2cxf, rectdisk2cheb, tri2cheb, tet2cheb,
3939
fourier2sph, fourier2sphv, cxf2disk, cheb2rectdisk, cheb2tri, cheb2tet
4040

4141
export plan_leg2cheb, plan_cheb2leg, plan_ultra2ultra, plan_jac2jac,
4242
plan_lag2lag, plan_jac2ultra, plan_ultra2jac, plan_jac2cheb,
43-
plan_cheb2jac, plan_ultra2cheb, plan_cheb2ultra,
43+
plan_cheb2jac, plan_ultra2cheb, plan_cheb2ultra, plan_associatedjac2jac,
4444
plan_sph2fourier, plan_sph_synthesis, plan_sph_analysis,
4545
plan_sphv2fourier, plan_sphv_synthesis, plan_sphv_analysis,
4646
plan_disk2cxf, plan_disk_synthesis, plan_disk_analysis,

‎src/libfasttransforms.jl

+79-23
Original file line numberDiff line numberDiff line change
@@ -116,28 +116,29 @@ const JAC2CHEB = 7
116116
const CHEB2JAC = 8
117117
const ULTRA2CHEB = 9
118118
const CHEB2ULTRA = 10
119-
const SPHERE = 11
120-
const SPHEREV = 12
121-
const DISK = 13
122-
const RECTDISK = 14
123-
const TRIANGLE = 15
124-
const TETRAHEDRON = 16
125-
const SPINSPHERE = 17
126-
const SPHERESYNTHESIS = 18
127-
const SPHEREANALYSIS = 19
128-
const SPHEREVSYNTHESIS = 20
129-
const SPHEREVANALYSIS = 21
130-
const DISKSYNTHESIS = 22
131-
const DISKANALYSIS = 23
132-
const RECTDISKSYNTHESIS = 24
133-
const RECTDISKANALYSIS = 25
134-
const TRIANGLESYNTHESIS = 26
135-
const TRIANGLEANALYSIS = 27
136-
const TETRAHEDRONSYNTHESIS = 28
137-
const TETRAHEDRONANALYSIS = 29
138-
const SPINSPHERESYNTHESIS = 30
139-
const SPINSPHEREANALYSIS = 31
140-
const SPHERICALISOMETRY = 32
119+
const ASSOCIATEDJAC2JAC = 11
120+
const SPHERE = 12
121+
const SPHEREV = 13
122+
const DISK = 14
123+
const RECTDISK = 15
124+
const TRIANGLE = 16
125+
const TETRAHEDRON = 17
126+
const SPINSPHERE = 18
127+
const SPHERESYNTHESIS = 19
128+
const SPHEREANALYSIS = 20
129+
const SPHEREVSYNTHESIS = 21
130+
const SPHEREVANALYSIS = 22
131+
const DISKSYNTHESIS = 23
132+
const DISKANALYSIS = 24
133+
const RECTDISKSYNTHESIS = 25
134+
const RECTDISKANALYSIS = 26
135+
const TRIANGLESYNTHESIS = 27
136+
const TRIANGLEANALYSIS = 28
137+
const TETRAHEDRONSYNTHESIS = 29
138+
const TETRAHEDRONANALYSIS = 30
139+
const SPINSPHERESYNTHESIS = 31
140+
const SPINSPHEREANALYSIS = 32
141+
const SPHERICALISOMETRY = 33
141142

142143

143144
let k2s = Dict(LEG2CHEB => "Legendre--Chebyshev",
@@ -151,6 +152,7 @@ let k2s = Dict(LEG2CHEB => "Legendre--Chebyshev",
151152
CHEB2JAC => "Chebyshev--Jacobi",
152153
ULTRA2CHEB => "ultraspherical--Chebyshev",
153154
CHEB2ULTRA => "Chebyshev--ultraspherical",
155+
ASSOCIATEDJAC2JAC => "Associated Jacobi--Jacobi",
154156
SPHERE => "Spherical harmonic--Fourier",
155157
SPHEREV => "Spherical vector field--Fourier",
156158
DISK => "Zernike--Chebyshev×Fourier",
@@ -244,6 +246,8 @@ unsafe_convert(::Type{Ptr{mpfr_t}}, p::FTPlan) = unsafe_convert(Ptr{mpfr_t}, p.p
244246
destroy_plan(p::FTPlan{Float32, 1}) = ccall((:ft_destroy_tb_eigen_FMMf, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
245247
destroy_plan(p::FTPlan{Float64, 1}) = ccall((:ft_destroy_tb_eigen_FMM, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
246248
destroy_plan(p::FTPlan{BigFloat, 1}) = ccall((:ft_mpfr_destroy_plan, libfasttransforms), Cvoid, (Ptr{mpfr_t}, Cint), p, p.n)
249+
destroy_plan(p::FTPlan{Float32, 1, ASSOCIATEDJAC2JAC}) = ccall((:ft_destroy_btb_eigen_FMMf, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
250+
destroy_plan(p::FTPlan{Float64, 1, ASSOCIATEDJAC2JAC}) = ccall((:ft_destroy_btb_eigen_FMM, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
247251
destroy_plan(p::FTPlan{Float64, 2}) = ccall((:ft_destroy_harmonic_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
248252
destroy_plan(p::FTPlan{Float64, 3}) = ccall((:ft_destroy_tetrahedral_harmonic_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
249253
destroy_plan(p::FTPlan{Complex{Float64}, 2, SPINSPHERE}) = ccall((:ft_destroy_spin_harmonic_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
@@ -307,7 +311,7 @@ unsafe_convert(::Type{Ptr{mpfr_t}}, p::TransposeFTPlan{T, FTPlan{T, N, K}}) wher
307311

308312
for f in (:leg2cheb, :cheb2leg, :ultra2ultra, :jac2jac,
309313
:lag2lag, :jac2ultra, :ultra2jac, :jac2cheb,
310-
:cheb2jac, :ultra2cheb, :cheb2ultra,
314+
:cheb2jac, :ultra2cheb, :cheb2ultra, :associatedjac2jac,
311315
:sph2fourier, :sphv2fourier, :disk2cxf,
312316
:rectdisk2cheb, :tri2cheb, :tet2cheb)
313317
plan_f = Symbol("plan_", f)
@@ -385,6 +389,11 @@ function plan_cheb2ultra(::Type{Float32}, n::Integer, λ; normcheb::Bool=false,
385389
return FTPlan{Float32, 1, CHEB2ULTRA}(plan, n)
386390
end
387391

392+
function plan_associatedjac2jac(::Type{Float32}, n::Integer, c::Integer, α, β, γ, δ; norm1::Bool=false, norm2::Bool=false)
393+
plan = ccall((:ft_plan_associated_jacobi_to_jacobif, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Cint, Float32, Float32, Float32, Float32), norm1, norm2, n, c, α, β, γ, δ)
394+
return FTPlan{Float32, 1, ASSOCIATEDJAC2JAC}(plan, n)
395+
end
396+
388397

389398
function plan_leg2cheb(::Type{Float64}, n::Integer; normleg::Bool=false, normcheb::Bool=false)
390399
plan = ccall((:ft_plan_legendre_to_chebyshev, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint), normleg, normcheb, n)
@@ -441,6 +450,11 @@ function plan_cheb2ultra(::Type{Float64}, n::Integer, λ; normcheb::Bool=false,
441450
return FTPlan{Float64, 1, CHEB2ULTRA}(plan, n)
442451
end
443452

453+
function plan_associatedjac2jac(::Type{Float64}, n::Integer, c::Integer, α, β, γ, δ; norm1::Bool=false, norm2::Bool=false)
454+
plan = ccall((:ft_plan_associated_jacobi_to_jacobi, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Cint, Float64, Float64, Float64, Float64), norm1, norm2, n, c, α, β, γ, δ)
455+
return FTPlan{Float64, 1, ASSOCIATEDJAC2JAC}(plan, n)
456+
end
457+
444458

445459
function plan_leg2cheb(::Type{BigFloat}, n::Integer; normleg::Bool=false, normcheb::Bool=false)
446460
plan = ccall((:ft_mpfr_plan_legendre_to_chebyshev, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Clong, Int32), normleg, normcheb, n, precision(BigFloat), Base.MPFR.ROUNDING_MODE[])
@@ -664,6 +678,27 @@ for (fJ, fC, elty) in ((:lmul!, :ft_bfmvf, :Float32),
664678
end
665679
end
666680

681+
for (fJ, fC, elty) in ((:lmul!, :ft_bbbfmvf, :Float32),
682+
(:lmul!, :ft_bbbfmv , :Float64))
683+
@eval begin
684+
function $fJ(p::FTPlan{$elty, 1, ASSOCIATEDJAC2JAC}, x::Vector{$elty})
685+
checksize(p, x)
686+
ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Cint, Ptr{ft_plan_struct}, Ptr{$elty}), 'N', '2', '1', p, x)
687+
return x
688+
end
689+
function $fJ(p::AdjointFTPlan{$elty, FTPlan{$elty, 1, ASSOCIATEDJAC2JAC}}, x::Vector{$elty})
690+
checksize(p, x)
691+
ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Cint, Ptr{ft_plan_struct}, Ptr{$elty}), 'T', '1', '2', p, x)
692+
return x
693+
end
694+
function $fJ(p::TransposeFTPlan{$elty, FTPlan{$elty, 1, ASSOCIATEDJAC2JAC}}, x::Vector{$elty})
695+
checksize(p, x)
696+
ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Cint, Ptr{ft_plan_struct}, Ptr{$elty}), 'T', '1', '2', p, x)
697+
return x
698+
end
699+
end
700+
end
701+
667702
for (fJ, fC) in ((:lmul!, :ft_mpfr_trmv_ptr),
668703
(:ldiv!, :ft_mpfr_trsv_ptr))
669704
@eval begin
@@ -708,6 +743,27 @@ for (fJ, fC, elty) in ((:lmul!, :ft_bfmmf, :Float32),
708743
end
709744
end
710745

746+
for (fJ, fC, elty) in ((:lmul!, :ft_bbbfmmf, :Float32),
747+
(:lmul!, :ft_bbbfmm , :Float64))
748+
@eval begin
749+
function $fJ(p::FTPlan{$elty, 1, ASSOCIATEDJAC2JAC}, x::Matrix{$elty})
750+
checksize(p, x)
751+
ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Cint, Ptr{ft_plan_struct}, Ptr{$elty}, Cint, Cint), 'N', '2', '1', p, x, size(x, 1), size(x, 2))
752+
return x
753+
end
754+
function $fJ(p::AdjointFTPlan{$elty, FTPlan{$elty, 1, ASSOCIATEDJAC2JAC}}, x::Matrix{$elty})
755+
checksize(p, x)
756+
ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Cint, Ptr{ft_plan_struct}, Ptr{$elty}, Cint, Cint), 'T', '1', '2', p, x, size(x, 1), size(x, 2))
757+
return x
758+
end
759+
function $fJ(p::TransposeFTPlan{$elty, FTPlan{$elty, 1, ASSOCIATEDJAC2JAC}}, x::Matrix{$elty})
760+
checksize(p, x)
761+
ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Cint, Ptr{ft_plan_struct}, Ptr{$elty}, Cint, Cint), 'T', '1', '2', p, x, size(x, 1), size(x, 2))
762+
return x
763+
end
764+
end
765+
end
766+
711767
for (fJ, fC) in ((:lmul!, :ft_mpfr_trmm_ptr),
712768
(:ldiv!, :ft_mpfr_trsm_ptr))
713769
@eval begin

‎test/libfasttransformstests.jl

+10
Original file line numberDiff line numberDiff line change
@@ -75,6 +75,16 @@ FastTransforms.set_num_threads(ceil(Int, Base.Sys.CPU_THREADS/2))
7575
end
7676
end
7777

78+
for T in (Float32, Float64, Complex{Float32}, Complex{Float64})
79+
x = T(1)./(1:n)
80+
Id = Matrix{T}(I, n, n)
81+
p = plan_associatedjac2jac(Id, 1, α, β, γ, δ)
82+
V = p*I
83+
@test V ≈ p*Id
84+
y = p*x
85+
@test V\y ≈ x
86+
end
87+
7888
function test_nd_plans(p, ps, pa, A)
7989
B = copy(A)
8090
C = ps*(p*A)

0 commit comments

Comments
 (0)