Skip to content

Commit cc10528

Browse files
MikaelSlevinskyMikaelSlevinsky
authored
add modified classical op ccalls (#173)
1 parent ee31e98 commit cc10528

File tree

5 files changed

+146
-33
lines changed

5 files changed

+146
-33
lines changed

‎Project.toml‎

Lines changed: 2 additions & 2 deletions
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.13.9"
3+
version = "0.14.0"
44

55
[deps]
66
AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
@@ -23,7 +23,7 @@ ArrayLayouts = "0.4, 0.5, 0.6, 0.7, 0.8"
2323
DSP = "0.6, 0.7"
2424
FFTW = "1"
2525
FastGaussQuadrature = "0.4"
26-
FastTransforms_jll = "0.5.4"
26+
FastTransforms_jll = "0.6.0"
2727
FillArrays = "0.9, 0.10, 0.11, 0.12, 0.13"
2828
Reexport = "0.2, 1.0"
2929
SpecialFunctions = "0.10, 1, 2"

‎docs/src/dev.md‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -40,34 +40,34 @@ This lets the developer experiment with new features through `ccall`ing into ble
4040
To get from a C library release to a Julia package release, the developer needs to update Yggdrasil's [build_tarballs.jl](https://github.com/JuliaPackaging/Yggdrasil/blob/master/F/FastTransforms/build_tarballs.jl) script for the new version and its 256-bit SHA. On macOS, the SHA can be found by:
4141
4242
```julia
43-
shell> curl https://codeload.github.com/MikaelSlevinsky/FastTransforms/tar.gz/v0.5.0 --output FastTransforms-0.5.0.tar.gz
43+
shell> curl https://codeload.github.com/MikaelSlevinsky/FastTransforms/tar.gz/v0.6.0 --output FastTransforms.tar.gz
4444
% Total % Received % Xferd Average Speed Time Time Time Current
4545
Dload Upload Total Spent Left Speed
46-
100 156k 100 156k 0 0 349k 0 --:--:-- --:--:-- --:--:-- 348k
46+
100 162k 0 162k 0 0 252k 0 --:--:-- --:--:-- --:--:-- 252k
4747

48-
shell> shasum -a 256 FastTransforms-0.5.0.tar.gz
49-
9556d0037bd5348a33f15ad6100e32053b6e22cab16a97c504f30d6c52fd0efd FastTransforms-0.5.0.tar.gz
48+
shell> shasum -a 256 FastTransforms.tar.gz
49+
ae2db2fa808ca17c5dc5ac25b079eba2dbe598d061b9b4e14c948680870abc3c FastTransforms.tar.gz
5050

51-
shell> rm -f FastTransforms-0.5.0.tar.gz
51+
shell> rm -f FastTransforms.tar.gz
5252

5353
```
5454
5555
Using [SHA.jl](https://github.com/JuliaCrypto/SHA.jl), the SHA can also be found by:
5656
5757
```julia
58-
shell> curl https://codeload.github.com/MikaelSlevinsky/FastTransforms/tar.gz/v0.5.0 --output FastTransforms-0.5.0.tar.gz
58+
shell> curl https://codeload.github.com/MikaelSlevinsky/FastTransforms/tar.gz/v0.6.0 --output FastTransforms.tar.gz
5959
% Total % Received % Xferd Average Speed Time Time Time Current
6060
Dload Upload Total Spent Left Speed
6161
100 156k 0 156k 0 0 443k 0 --:--:-- --:--:-- --:--:-- 443k
6262

6363
julia> using SHA
6464

65-
julia> open("FastTransforms-0.5.0.tar.gz") do f
65+
julia> open("FastTransforms.tar.gz") do f
6666
bytes2hex(sha256(f))
6767
end
68-
"9556d0037bd5348a33f15ad6100e32053b6e22cab16a97c504f30d6c52fd0efd"
68+
"ae2db2fa808ca17c5dc5ac25b079eba2dbe598d061b9b4e14c948680870abc3c"
6969

70-
shell> rm -f FastTransforms-0.5.0.tar.gz
70+
shell> rm -f FastTransforms.tar.gz
7171

7272
```
7373

‎src/FastTransforms.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -34,12 +34,14 @@ import LinearAlgebra: mul!, lmul!, ldiv!
3434
export leg2cheb, cheb2leg, ultra2ultra, jac2jac,
3535
lag2lag, jac2ultra, ultra2jac, jac2cheb,
3636
cheb2jac, ultra2cheb, cheb2ultra, associatedjac2jac,
37+
modifiedjac2jac, modifiedlag2lag, modifiedherm2herm,
3738
sph2fourier, sphv2fourier, disk2cxf, rectdisk2cheb, tri2cheb, tet2cheb,
3839
fourier2sph, fourier2sphv, cxf2disk, cheb2rectdisk, cheb2tri, cheb2tet
3940

4041
export plan_leg2cheb, plan_cheb2leg, plan_ultra2ultra, plan_jac2jac,
4142
plan_lag2lag, plan_jac2ultra, plan_ultra2jac, plan_jac2cheb,
4243
plan_cheb2jac, plan_ultra2cheb, plan_cheb2ultra, plan_associatedjac2jac,
44+
plan_modifiedjac2jac, plan_modifiedlag2lag, plan_modifiedherm2herm,
4345
plan_sph2fourier, plan_sph_synthesis, plan_sph_analysis,
4446
plan_sphv2fourier, plan_sphv_synthesis, plan_sphv_analysis,
4547
plan_disk2cxf, plan_disk_synthesis, plan_disk_analysis,

‎src/libfasttransforms.jl‎

Lines changed: 106 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -115,28 +115,31 @@ const CHEB2JAC = 8
115115
const ULTRA2CHEB = 9
116116
const CHEB2ULTRA = 10
117117
const ASSOCIATEDJAC2JAC = 11
118-
const SPHERE = 12
119-
const SPHEREV = 13
120-
const DISK = 14
121-
const RECTDISK = 15
122-
const TRIANGLE = 16
123-
const TETRAHEDRON = 17
124-
const SPINSPHERE = 18
125-
const SPHERESYNTHESIS = 19
126-
const SPHEREANALYSIS = 20
127-
const SPHEREVSYNTHESIS = 21
128-
const SPHEREVANALYSIS = 22
129-
const DISKSYNTHESIS = 23
130-
const DISKANALYSIS = 24
131-
const RECTDISKSYNTHESIS = 25
132-
const RECTDISKANALYSIS = 26
133-
const TRIANGLESYNTHESIS = 27
134-
const TRIANGLEANALYSIS = 28
135-
const TETRAHEDRONSYNTHESIS = 29
136-
const TETRAHEDRONANALYSIS = 30
137-
const SPINSPHERESYNTHESIS = 31
138-
const SPINSPHEREANALYSIS = 32
139-
const SPHERICALISOMETRY = 33
118+
const MODIFIEDJAC2JAC = 12
119+
const MODIFIEDLAG2LAG = 13
120+
const MODIFIEDHERM2HERM = 14
121+
const SPHERE = 15
122+
const SPHEREV = 16
123+
const DISK = 17
124+
const RECTDISK = 18
125+
const TRIANGLE = 19
126+
const TETRAHEDRON = 20
127+
const SPINSPHERE = 21
128+
const SPHERESYNTHESIS = 22
129+
const SPHEREANALYSIS = 23
130+
const SPHEREVSYNTHESIS = 24
131+
const SPHEREVANALYSIS = 25
132+
const DISKSYNTHESIS = 26
133+
const DISKANALYSIS = 27
134+
const RECTDISKSYNTHESIS = 28
135+
const RECTDISKANALYSIS = 29
136+
const TRIANGLESYNTHESIS = 30
137+
const TRIANGLEANALYSIS = 31
138+
const TETRAHEDRONSYNTHESIS = 32
139+
const TETRAHEDRONANALYSIS = 33
140+
const SPINSPHERESYNTHESIS = 34
141+
const SPINSPHEREANALYSIS = 35
142+
const SPHERICALISOMETRY = 36
140143

141144

142145
let k2s = Dict(LEG2CHEB => "Legendre--Chebyshev",
@@ -151,6 +154,9 @@ let k2s = Dict(LEG2CHEB => "Legendre--Chebyshev",
151154
ULTRA2CHEB => "ultraspherical--Chebyshev",
152155
CHEB2ULTRA => "Chebyshev--ultraspherical",
153156
ASSOCIATEDJAC2JAC => "Associated Jacobi--Jacobi",
157+
MODIFIEDJAC2JAC => "Modified Jacobi--Jacobi",
158+
MODIFIEDLAG2LAG => "Modified Laguerre--Laguerre",
159+
MODIFIEDHERM2HERM => "Modified Hermite--Hermite",
154160
SPHERE => "Spherical harmonic--Fourier",
155161
SPHEREV => "Spherical vector field--Fourier",
156162
DISK => "Zernike--Chebyshev×Fourier",
@@ -266,6 +272,9 @@ destroy_plan(p::FTPlan{Float64, 1}) = ccall((:ft_destroy_tb_eigen_FMM, libfasttr
266272
destroy_plan(p::FTPlan{BigFloat, 1}) = ccall((:ft_mpfr_destroy_plan, libfasttransforms), Cvoid, (Ptr{mpfr_t}, Cint), p, p.n)
267273
destroy_plan(p::FTPlan{Float32, 1, ASSOCIATEDJAC2JAC}) = ccall((:ft_destroy_btb_eigen_FMMf, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
268274
destroy_plan(p::FTPlan{Float64, 1, ASSOCIATEDJAC2JAC}) = ccall((:ft_destroy_btb_eigen_FMM, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
275+
destroy_plan(p::FTPlan{Float64, 1, MODIFIEDJAC2JAC}) = ccall((:ft_destroy_modified_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
276+
destroy_plan(p::FTPlan{Float64, 1, MODIFIEDLAG2LAG}) = ccall((:ft_destroy_modified_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
277+
destroy_plan(p::FTPlan{Float64, 1, MODIFIEDHERM2HERM}) = ccall((:ft_destroy_modified_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
269278
destroy_plan(p::FTPlan{Float64}) = ccall((:ft_destroy_harmonic_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
270279
destroy_plan(p::FTPlan{Complex{Float64}, 2, SPINSPHERE}) = ccall((:ft_destroy_spin_harmonic_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
271280
destroy_plan(p::FTPlan{Float64, 2, SPHERESYNTHESIS}) = ccall((:ft_destroy_sphere_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
@@ -383,6 +392,7 @@ end
383392
for f in (:leg2cheb, :cheb2leg, :ultra2ultra, :jac2jac,
384393
:lag2lag, :jac2ultra, :ultra2jac, :jac2cheb,
385394
:cheb2jac, :ultra2cheb, :cheb2ultra, :associatedjac2jac,
395+
:modifiedjac2jac, :modifiedlag2lag, :modifiedherm2herm,
386396
:sph2fourier, :sphv2fourier, :disk2cxf,
387397
:rectdisk2cheb, :tri2cheb, :tet2cheb)
388398
plan_f = Symbol("plan_", f)
@@ -526,6 +536,36 @@ function plan_associatedjac2jac(::Type{Float64}, n::Integer, c::Integer, α, β,
526536
return FTPlan{Float64, 1, ASSOCIATEDJAC2JAC}(plan, n)
527537
end
528538

539+
function plan_modifiedjac2jac(::Type{Float64}, n::Integer, α, β, w::Vector{Float64}; verbose::Bool=false)
540+
#plan_modifiedjac2jac(Float64, n, α, β, w, Vector{Float64}(undef, 0); verbose=verbose)
541+
plan = ccall((:ft_plan_modified_jacobi_to_jacobi, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Float64, Float64, Cint, Ptr{Float64}, Cint, Ptr{Float64}, Cint), n, α, β, length(w), w, 0, C_NULL, verbose)
542+
return FTPlan{Float64, 1, MODIFIEDJAC2JAC}(plan, n)
543+
end
544+
545+
function plan_modifiedjac2jac(::Type{Float64}, n::Integer, α, β, u::Vector{Float64}, v::Vector{Float64}; verbose::Bool=false)
546+
plan = ccall((:ft_plan_modified_jacobi_to_jacobi, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Float64, Float64, Cint, Ptr{Float64}, Cint, Ptr{Float64}, Cint), n, α, β, length(u), u, length(v), v, verbose)
547+
return FTPlan{Float64, 1, MODIFIEDJAC2JAC}(plan, n)
548+
end
549+
550+
function plan_modifiedlag2lag(::Type{Float64}, n::Integer, α, w::Vector{Float64}; verbose::Bool=false)
551+
plan = ccall((:ft_plan_modified_laguerre_to_laguerre, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Float64, Cint, Ptr{Float64}, Cint, Ptr{Float64}, Cint), n, α, length(w), w, 0, C_NULL, verbose)
552+
return FTPlan{Float64, 1, MODIFIEDLAG2LAG}(plan, n)
553+
end
554+
555+
function plan_modifiedlag2lag(::Type{Float64}, n::Integer, α, u::Vector{Float64}, v::Vector{Float64}; verbose::Bool=false)
556+
plan = ccall((:ft_plan_modified_laguerre_to_laguerre, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Float64, Cint, Ptr{Float64}, Cint, Ptr{Float64}, Cint), n, α, length(u), u, length(v), v, verbose)
557+
return FTPlan{Float64, 1, MODIFIEDLAG2LAG}(plan, n)
558+
end
559+
560+
function plan_modifiedherm2herm(::Type{Float64}, n::Integer, w::Vector{Float64}; verbose::Bool=false)
561+
plan = ccall((:ft_plan_modified_hermite_to_hermite, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Ptr{Float64}, Cint, Ptr{Float64}, Cint), n, length(w), w, 0, C_NULL, verbose)
562+
return FTPlan{Float64, 1, MODIFIEDHERM2HERM}(plan, n)
563+
end
564+
565+
function plan_modifiedherm2herm(::Type{Float64}, n::Integer, u::Vector{Float64}, v::Vector{Float64}; verbose::Bool=false)
566+
plan = ccall((:ft_plan_modified_hermite_to_hermite, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Ptr{Float64}, Cint, Ptr{Float64}, Cint), n, length(u), u, length(v), v, verbose)
567+
return FTPlan{Float64, 1, MODIFIEDHERM2HERM}(plan, n)
568+
end
529569

530570
function plan_leg2cheb(::Type{BigFloat}, n::Integer; normleg::Bool=false, normcheb::Bool=false)
531571
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[])
@@ -789,6 +829,28 @@ for (fJ, fC, elty) in ((:lmul!, :ft_bbbfmvf, :Float32),
789829
end
790830
end
791831

832+
for (fJ, fC, elty) in ((:lmul!, :ft_mpmv, :Float64),
833+
(:ldiv!, :ft_mpsv, :Float64))
834+
@eval begin
835+
ModifiedFTPlan = Union{FTPlan{$elty, 1, MODIFIEDJAC2JAC}, FTPlan{$elty, 1, MODIFIEDLAG2LAG}, FTPlan{$elty, 1, MODIFIEDHERM2HERM}}
836+
function $fJ(p::ModifiedFTPlan, x::Vector{$elty})
837+
checksize(p, x)
838+
ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Ptr{ft_plan_struct}, Ptr{$elty}), 'N', p, x)
839+
return x
840+
end
841+
function $fJ(p::AdjointFTPlan{$elty, ModifiedFTPlan}, x::Vector{$elty})
842+
checksize(p, x)
843+
ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Ptr{ft_plan_struct}, Ptr{$elty}), 'T', p, x)
844+
return x
845+
end
846+
function $fJ(p::TransposeFTPlan{$elty, ModifiedFTPlan}, x::Vector{$elty})
847+
checksize(p, x)
848+
ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Ptr{ft_plan_struct}, Ptr{$elty}), 'T', p, x)
849+
return x
850+
end
851+
end
852+
end
853+
792854
for (fJ, fC) in ((:lmul!, :ft_mpfr_trmv_ptr),
793855
(:ldiv!, :ft_mpfr_trsv_ptr))
794856
@eval begin
@@ -854,6 +916,28 @@ for (fJ, fC, elty) in ((:lmul!, :ft_bbbfmmf, :Float32),
854916
end
855917
end
856918

919+
for (fJ, fC, elty) in ((:lmul!, :ft_mpmm, :Float64),
920+
(:ldiv!, :ft_mpsm, :Float64))
921+
@eval begin
922+
ModifiedFTPlan = Union{FTPlan{$elty, 1, MODIFIEDJAC2JAC}, FTPlan{$elty, 1, MODIFIEDLAG2LAG}, FTPlan{$elty, 1, MODIFIEDHERM2HERM}}
923+
function $fJ(p::ModifiedFTPlan, x::Matrix{$elty})
924+
checksize(p, x)
925+
ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Ptr{ft_plan_struct}, Ptr{$elty}, Cint, Cint), 'N', p, x, size(x, 1), size(x, 2))
926+
return x
927+
end
928+
function $fJ(p::AdjointFTPlan{$elty, ModifiedFTPlan}, x::Matrix{$elty})
929+
checksize(p, x)
930+
ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Ptr{ft_plan_struct}, Ptr{$elty}, Cint, Cint), 'T', p, x, size(x, 1), size(x, 2))
931+
return x
932+
end
933+
function $fJ(p::TransposeFTPlan{$elty, ModifiedFTPlan}, x::Matrix{$elty})
934+
checksize(p, x)
935+
ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Ptr{ft_plan_struct}, Ptr{$elty}, Cint, Cint), 'T', p, x, size(x, 1), size(x, 2))
936+
return x
937+
end
938+
end
939+
end
940+
857941
for (fJ, fC) in ((:lmul!, :ft_mpfr_trmm_ptr),
858942
(:ldiv!, :ft_mpfr_trsm_ptr))
859943
@eval begin

‎test/libfasttransformstests.jl‎

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -85,6 +85,33 @@ FastTransforms.ft_set_num_threads(ceil(Int, Base.Sys.CPU_THREADS/2))
8585
@test V\y ≈ x
8686
end
8787

88+
@testset "Modified classical orthonormal polynomial transforms" begin
89+
(n, α, β) = (16, 0, 0)
90+
P1 = plan_modifiedjac2jac(Float64, n, α, β, [0.9428090415820636, -0.32659863237109055, -0.42163702135578396, 0.2138089935299396]) # u1(x) = (1-x)^2*(1+x)
91+
P2 = plan_modifiedjac2jac(Float64, n, α, β, [0.9428090415820636, -0.32659863237109055, -0.42163702135578396, 0.2138089935299396], [1.4142135623730951]) # u2(x) = (1-x)^2*(1+x)
92+
P3 = plan_modifiedjac2jac(Float64, n, α, β, [-0.9428090415820636, 0.32659863237109055, 0.42163702135578396, -0.2138089935299396], [-5.185449728701348, 0.0, 0.42163702135578374]) # u3(x) = -(1-x)^2*(1+x), v3(x) = -(2-x)*(2+x)
93+
P4 = plan_modifiedjac2jac(Float64, n, α+2, β+1, [1.1547005383792517], [4.387862045841156, 0.1319657758147716, -0.20865621238292037]) # v4(x) = (2-x)*(2+x)
94+
95+
@test P1*I ≈ P2*I
96+
@test P1\I ≈ P2\I
97+
@test P3*I ≈ P2*(P4*I)
98+
@test P3\I ≈ P4\(P2\I)
99+
100+
P5 = plan_modifiedlag2lag(Float64, n, α, [2.0, -4.0, 2.0]) # u5(x) = x^2
101+
P6 = plan_modifiedlag2lag(Float64, n, α, [2.0, -4.0, 2.0], [1.0]) # u6(x) = x^2
102+
P7 = plan_modifiedlag2lag(Float64, n, α, [2.0, -4.0, 2.0], [7.0, -7.0, 2.0]) # u7(x) = x^2, v7(x) = (1+x)*(2+x)
103+
P8 = plan_modifiedlag2lag(Float64, n, α+2, [sqrt(2.0)], [sqrt(1058.0), -sqrt(726.0), sqrt(48.0)]) # v8(x) = (1+x)*(2+x)
104+
105+
@test P5*I ≈ P6*I
106+
@test P5\I ≈ P6\I
107+
@test P7*I ≈ P6*(P8*I)
108+
@test P7\I ≈ P8\(P6\I)
109+
110+
P9 = plan_modifiedherm2herm(Float64, n, [2.995504568550877, 0.0, 3.7655850551068593, 0.0, 1.6305461589167827], [2.995504568550877, 0.0, 3.7655850551068593, 0.0, 1.6305461589167827]) # u9(x) = 1+x^2+x^4, v9(x) = 1+x^2+x^4
111+
112+
@test P9*I ≈ P9\I
113+
end
114+
88115
function test_nd_plans(p, ps, pa, A)
89116
B = copy(A)
90117
C = ps*(p*A)

0 commit comments

Comments
 (0)