Skip to content

Commit ee98f40

Browse files
trial support of libfasttransforms
no binarybuilder yet
1 parent a80acce commit ee98f40

File tree

12 files changed

+305
-62
lines changed

12 files changed

+305
-62
lines changed

.gitignore

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,2 +1,3 @@
11
docs/build/
22
docs/site/
3+
FastTransforms

src/FastTransforms.jl

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -106,6 +106,8 @@ jac2jac(x...) = th_jac2jac(x...)
106106
plan_jac2jac(x...) = th_jac2jacplan(x...)
107107
plan_ultra2ultra(x...) = th_ultra2ultraplan(x...)
108108

109+
include("libfasttransforms.jl")
110+
109111
include("hierarchical.jl")
110112
include("SphericalHarmonics/SphericalHarmonics.jl")
111113
include("TriangularHarmonics/TriangularHarmonics.jl")

src/SphericalHarmonics/SphericalHarmonics.jl

Lines changed: 3 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,15 +1,13 @@
1-
abstract type SphericalHarmonicPlan{T} end
2-
3-
function *(P::SphericalHarmonicPlan, X::AbstractMatrix)
1+
function *(P::HarmonicPlan, X::AbstractMatrix)
42
mul!(zero(X), P, X)
53
end
64

75
if VERSION < v"0.7-"
8-
function \(P::SphericalHarmonicPlan, X::AbstractMatrix)
6+
function \(P::HarmonicPlan, X::AbstractMatrix)
97
At_mul_B!(zero(X), P, X)
108
end
119
else
12-
function \(P::SphericalHarmonicPlan, X::AbstractMatrix)
10+
function \(P::HarmonicPlan, X::AbstractMatrix)
1311
mul!(zero(X), transpose(P), X)
1412
end
1513
end

src/SphericalHarmonics/fastplan.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
struct FastSphericalHarmonicPlan{T} <: SphericalHarmonicPlan{T}
1+
struct FastSphericalHarmonicPlan{T} <: HarmonicPlan{T}
22
RP::RotationPlan{T}
33
BF::Vector{Butterfly{T}}
44
p1::NormalizedLegendreToChebyshevPlan{T}

src/SphericalHarmonics/slowplan.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -308,7 +308,7 @@ end
308308

309309

310310

311-
struct SlowSphericalHarmonicPlan{T} <: SphericalHarmonicPlan{T}
311+
struct SlowSphericalHarmonicPlan{T} <: HarmonicPlan{T}
312312
RP::RotationPlan{T}
313313
p1::NormalizedLegendreToChebyshevPlan{T}
314314
p2::NormalizedLegendre1ToChebyshev2Plan{T}

src/SphericalHarmonics/thinplan.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@ const LAYERSKELETON = 64
22

33
checklayer(J::Int) = J÷LAYERSKELETON == J/LAYERSKELETON
44

5-
struct ThinSphericalHarmonicPlan{T} <: SphericalHarmonicPlan{T}
5+
struct ThinSphericalHarmonicPlan{T} <: HarmonicPlan{T}
66
RP::RotationPlan{T}
77
BF::Vector{Butterfly{T}}
88
p1::NormalizedLegendreToChebyshevPlan{T}

src/TriangularHarmonics/TriangularHarmonics.jl

Lines changed: 0 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,3 @@
1-
abstract type TriangularHarmonicPlan{T} end
2-
3-
function *(P::TriangularHarmonicPlan, X::AbstractMatrix)
4-
mul!(zero(X), P, X)
5-
end
6-
7-
function \(P::TriangularHarmonicPlan, X::AbstractMatrix)
8-
At_mul_B!(zero(X), P, X)
9-
end
10-
111
include("trifunctions.jl")
122
include("slowplan.jl")
133

src/TriangularHarmonics/slowplan.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -97,7 +97,7 @@ else
9797
lmul!(transpose(parent(Pc)), A)
9898
end
9999

100-
struct SlowTriangularHarmonicPlan{T} <: TriangularHarmonicPlan{T}
100+
struct SlowTriangularHarmonicPlan{T} <: HarmonicPlan{T}
101101
RP::TriRotationPlan{T}
102102
p::NormalizedLegendreToChebyshevPlan{T}
103103
pinv::ChebyshevToNormalizedLegendrePlan{T}
@@ -169,6 +169,6 @@ else
169169
tri_zero_spurious_modes!(At_mul_B!(RP, Y))
170170
end
171171

172-
LinearAlgebra.mul!(Y::Matrix, SPc::Adjoint{T,<:SlowTriangularHarmonicPlan}, X::Matrix) where T =
172+
LinearAlgebra.mul!(Y::Matrix, SPc::Adjoint{T,<:SlowTriangularHarmonicPlan}, X::Matrix) where T =
173173
mul!(Y, transpose(parent(SPc)), X)
174174
end

src/chebyshevtransform.jl

Lines changed: 42 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -121,48 +121,48 @@ end
121121
# Matrix inputs
122122
#
123123
#
124-
# function chebyshevtransform!(X::AbstractMatrix{T}; kind::Integer=1) where T<:fftwNumber
125-
# if kind == 1
126-
# if size(X) == (1,1)
127-
# X
128-
# else
129-
# X=r2r!(X,REDFT10)
130-
# X[:,1]/=2;X[1,:]/=2;
131-
# lmul!(1/(size(X,1)*size(X,2)),X)
132-
# end
133-
# elseif kind == 2
134-
# if size(X) == (1,1)
135-
# X
136-
# else
137-
# X=r2r!(X,REDFT00)
138-
# lmul!(1/((size(X,1)-1)*(size(X,2)-1)),X)
139-
# X[:,1]/=2;X[:,end]/=2
140-
# X[1,:]/=2;X[end,:]/=2
141-
# X
142-
# end
143-
# end
144-
# end
124+
function chebyshevtransform!(X::AbstractMatrix{T}; kind::Integer=1) where T<:fftwNumber
125+
if kind == 1
126+
if size(X) == (1,1)
127+
X
128+
else
129+
X=FFTW.r2r!(X,REDFT10)
130+
X[:,1]/=2;X[1,:]/=2;
131+
lmul!(1/(size(X,1)*size(X,2)),X)
132+
end
133+
elseif kind == 2
134+
if size(X) == (1,1)
135+
X
136+
else
137+
X=FFTW.r2r!(X,REDFT00)
138+
lmul!(1/((size(X,1)-1)*(size(X,2)-1)),X)
139+
X[:,1]/=2;X[:,end]/=2
140+
X[1,:]/=2;X[end,:]/=2
141+
X
142+
end
143+
end
144+
end
145145
#
146-
# function ichebyshevtransform!(X::AbstractMatrix{T}; kind::Integer=1) where T<:fftwNumber
147-
# if kind == 1
148-
# if size(X) == (1,1)
149-
# X
150-
# else
151-
# X[1,:]*=2;X[:,1]*=2
152-
# X = r2r(X,REDFT01)
153-
# lmul!(0.25, X)
154-
# end
155-
# elseif kind == 2
156-
# if size(X) == (1,1)
157-
# X
158-
# else
159-
# X[1,:]*=2;X[end,:]*=2;X[:,1]*=2;X[:,end]*=2
160-
# X=chebyshevtransform!(X;kind=kind)
161-
# X[1,:]*=2;X[end,:]*=2;X[:,1]*=2;X[:,end]*=2
162-
# lmul!((size(X,1)-1)*(size(X,2)-1)/4,X)
163-
# end
164-
# end
165-
# end
146+
function ichebyshevtransform!(X::AbstractMatrix{T}; kind::Integer=1) where T<:fftwNumber
147+
if kind == 1
148+
if size(X) == (1,1)
149+
X
150+
else
151+
X[1,:]*=2;X[:,1]*=2
152+
X = FFTW.r2r(X,REDFT01)
153+
lmul!(0.25, X)
154+
end
155+
elseif kind == 2
156+
if size(X) == (1,1)
157+
X
158+
else
159+
X[1,:]*=2;X[end,:]*=2;X[:,1]*=2;X[:,end]*=2
160+
X=chebyshevtransform!(X;kind=kind)
161+
X[1,:]*=2;X[end,:]*=2;X[:,1]*=2;X[:,end]*=2
162+
lmul!((size(X,1)-1)*(size(X,2)-1)/4,X)
163+
end
164+
end
165+
end
166166
#
167167

168168

@@ -294,7 +294,7 @@ function chebyshevpoints(::Type{T}, n::Integer; kind::Int=1) where T<:Number
294294
if n == 1
295295
zeros(T,1)
296296
else
297-
T[sinpi((n-2k-one(T))/(2n-2)) for k=0:n-1]
297+
T[sinpi((n-2k-one(T))/(2n-2)) for k=0:n-1]
298298
end
299299
else
300300
throw(ArgumentError("kind $kind not a valid kind of Chebyshev points"))

src/libfasttransforms.jl

Lines changed: 160 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,160 @@
1+
# This file shows how to call `libfasttransforms` from Julia.
2+
3+
# Step 1: In this repository, `git clone -b v0.1.0 https://github.com/MikaelSlevinsky/FastTransforms.git`
4+
5+
# Step 2: use a version of gcc that supports OpenMP: on OS X, this means using a
6+
# version of `gcc` from Homebrew, `brew install gcc`; on linux, `gcc-4.6` and up should work.
7+
# `export CC=gcc-"the-right-version"`.
8+
9+
# Step 3: get the remaining dependencies: On OS X, either `brew install openblas`
10+
# or change the Make.inc to use `BLAS = APPLEBLAS` instead of `BLAS = OPENBLAS`.
11+
# Furthermore, `brew install fftw`. For linux, see the `Travis.yml` file.
12+
13+
# Step 4: run `make` and check the tests by running `./test_drivers 3 3 0`.
14+
# All the errors should be roughly on the order of machine precision.
15+
16+
if VERSION < v"0.7-"
17+
using Base.Libdl
18+
else
19+
using Libdl
20+
end
21+
22+
import Base: unsafe_convert
23+
24+
const libfasttransforms = joinpath(dirname(@__DIR__), "FastTransforms", "libfasttransforms")
25+
26+
export plan_sph2fourier, plan_sph_synthesis, plan_sph_analysis, plan_sphv_synthesis, plan_sphv_analysis,
27+
sph2fourier!, fourier2sph!, sphv2fourier!, fourier2sphv!,
28+
sph_synthesis!, sph_analysis!, sphv_synthesis!, sphv_analysis!,
29+
plan_tri2cheb, plan_tri_synthesis, plan_tri_analysis,
30+
tri2cheb!, cheb2tri!, tri_synthesis!, tri_analysis!,
31+
plan_disk2cxf, plan_disk_synthesis, plan_disk_analysis,
32+
disk2cxf!, cxf2disk!, disk_synthesis!, disk_analysis!
33+
34+
35+
struct ft_plan_struct end
36+
37+
abstract type HarmonicPlan{T} end
38+
39+
for P in (:SphericalHarmonicPlan, :TriangularHarmonicPlan, :DiskHarmonicPlan)
40+
@eval begin
41+
mutable struct $P{T} <: HarmonicPlan{T}
42+
plan::Ptr{ft_plan_struct}
43+
function $P{T}(plan::Ptr{ft_plan_struct}) where T
44+
p = new(plan)
45+
finalizer(destroy_harmonic_plan, p)
46+
p
47+
end
48+
end
49+
end
50+
end
51+
52+
abstract type FTFFTWPlan{T} end
53+
54+
mutable struct SphereFFTWPlan{T} <: FTFFTWPlan{T}
55+
plan::Ptr{ft_plan_struct}
56+
function SphereFFTWPlan{T}(plan::Ptr{ft_plan_struct}) where T
57+
p = new(plan)
58+
finalizer(destroy_sphere_fftw_plan, p)
59+
p
60+
end
61+
end
62+
63+
mutable struct TriangleFFTWPlan{T} <: FTFFTWPlan{T}
64+
plan::Ptr{ft_plan_struct}
65+
function TriangleFFTWPlan{T}(plan::Ptr{ft_plan_struct}) where T
66+
p = new(plan)
67+
finalizer(destroy_triangle_fftw_plan, p)
68+
p
69+
end
70+
end
71+
72+
mutable struct DiskFFTWPlan{T} <: FTFFTWPlan{T}
73+
plan::Ptr{ft_plan_struct}
74+
function DiskFFTWPlan{T}(plan::Ptr{ft_plan_struct}) where T
75+
p = new(plan)
76+
finalizer(destroy_disk_fftw_plan, p)
77+
p
78+
end
79+
end
80+
81+
unsafe_convert(::Type{Ptr{ft_plan_struct}}, P::HarmonicPlan{Float64}) = P.plan
82+
unsafe_convert(::Type{Ptr{ft_plan_struct}}, P::FTFFTWPlan{Float64}) = P.plan
83+
84+
if Libdl.find_library(libfasttransforms) libfasttransforms
85+
set_num_threads(n::Int) = ccall((:ft_set_num_threads, libfasttransforms), Nothing, (Int, ), n)
86+
87+
destroy_harmonic_plan(P::HarmonicPlan{Float64}) = ccall((:ft_destroy_harmonic_plan, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, ), P)
88+
89+
function plan_sph2fourier(n::Int)
90+
plan = ccall((:ft_plan_sph2fourier, libfasttransforms), Ptr{ft_plan_struct}, (Int, ), n)
91+
return SphericalHarmonicPlan{Float64}(plan)
92+
end
93+
sph2fourier!(P::SphericalHarmonicPlan{Float64}, A::Matrix{Float64}) = ccall((:ft_execute_sph2fourier, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
94+
fourier2sph!(P::SphericalHarmonicPlan{Float64}, A::Matrix{Float64}) = ccall((:ft_execute_fourier2sph, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
95+
sphv2fourier!(P::SphericalHarmonicPlan{Float64}, A::Matrix{Float64}) = ccall((:ft_execute_sphv2fourier, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
96+
fourier2sphv!(P::SphericalHarmonicPlan{Float64}, A::Matrix{Float64}) = ccall((:ft_execute_fourier2sphv, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
97+
98+
function plan_tri2cheb(n::Int, α::Float64, β::Float64, γ::Float64)
99+
plan = ccall((:ft_plan_tri2cheb, libfasttransforms), Ptr{ft_plan_struct}, (Int, Float64, Float64, Float64), n, α, β, γ)
100+
return TriangularHarmonicPlan{Float64}(plan)
101+
end
102+
tri2cheb!(P::TriangularHarmonicPlan{Float64}, A::Matrix{Float64}) = ccall((:ft_execute_tri2cheb, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
103+
cheb2tri!(P::TriangularHarmonicPlan{Float64}, A::Matrix{Float64}) = ccall((:ft_execute_cheb2tri, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
104+
105+
function plan_disk2cxf(n::Int)
106+
plan = ccall((:ft_plan_disk2cxf, libfasttransforms), Ptr{ft_plan_struct}, (Int, ), n)
107+
return DiskHarmonicPlan{Float64}(plan)
108+
end
109+
disk2cxf!(P::DiskHarmonicPlan{Float64}, A::Matrix{Float64}) = ccall((:ft_execute_disk2cxf, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
110+
cxf2disk!(P::DiskHarmonicPlan{Float64}, A::Matrix{Float64}) = ccall((:ft_execute_cxf2disk, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
111+
112+
113+
destroy_sphere_fftw_plan(P::SphereFFTWPlan{Float64}) = ccall((:ft_destroy_sphere_fftw_plan, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, ), P)
114+
destroy_triangle_fftw_plan(P::TriangleFFTWPlan{Float64}) = ccall((:ft_destroy_triangle_fftw_plan, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, ), P)
115+
destroy_disk_fftw_plan(P::DiskFFTWPlan{Float64}) = ccall((:ft_destroy_disk_fftw_plan, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, ), P)
116+
117+
function plan_sph_synthesis(n::Int, m::Int)
118+
plan = ccall((:ft_plan_sph_synthesis, libfasttransforms), Ptr{ft_plan_struct}, (Int, Int), n, m)
119+
SphereFFTWPlan{Float64}(plan)
120+
end
121+
function plan_sph_analysis(n::Int, m::Int)
122+
plan = ccall((:ft_plan_sph_analysis, libfasttransforms), Ptr{ft_plan_struct}, (Int, Int), n, m)
123+
SphereFFTWPlan{Float64}(plan)
124+
end
125+
function plan_sphv_synthesis(n::Int, m::Int)
126+
plan = ccall((:ft_plan_sphv_synthesis, libfasttransforms), Ptr{ft_plan_struct}, (Int, Int), n, m)
127+
SphereFFTWPlan{Float64}(plan)
128+
end
129+
function plan_sphv_analysis(n::Int, m::Int)
130+
plan = ccall((:ft_plan_sphv_analysis, libfasttransforms), Ptr{ft_plan_struct}, (Int, Int), n, m)
131+
SphereFFTWPlan{Float64}(plan)
132+
end
133+
134+
sph_synthesis!(P::SphereFFTWPlan{Float64}, A::Matrix{Float64}) = ccall((:ft_execute_sph_synthesis, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
135+
sph_analysis!(P::SphereFFTWPlan{Float64}, A::Matrix{Float64}) = ccall((:ft_execute_sph_analysis, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
136+
sphv_synthesis!(P::SphereFFTWPlan{Float64}, A::Matrix{Float64}) = ccall((:ft_execute_sphv_synthesis, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
137+
sphv_analysis!(P::SphereFFTWPlan{Float64}, A::Matrix{Float64}) = ccall((:ft_execute_sphv_analysis, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
138+
139+
function plan_tri_synthesis(n::Int, m::Int)
140+
plan = ccall((:ft_plan_tri_synthesis, libfasttransforms), Ptr{ft_plan_struct}, (Int, Int), n, m)
141+
TriangleFFTWPlan{Float64}(plan)
142+
end
143+
function plan_tri_analysis(n::Int, m::Int)
144+
plan = ccall((:ft_plan_tri_analysis, libfasttransforms), Ptr{ft_plan_struct}, (Int, Int), n, m)
145+
TriangleFFTWPlan{Float64}(plan)
146+
end
147+
tri_synthesis!(P::TriangleFFTWPlan{Float64}, A::Matrix{Float64}) = ccall((:ft_execute_tri_synthesis, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
148+
tri_analysis!(P::TriangleFFTWPlan{Float64}, A::Matrix{Float64}) = ccall((:ft_execute_tri_analysis, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
149+
150+
function plan_disk_synthesis(n::Int, m::Int)
151+
plan = ccall((:ft_plan_disk_synthesis, libfasttransforms), Ptr{ft_plan_struct}, (Int, Int), n, m)
152+
DiskFFTWPlan{Float64}(plan)
153+
end
154+
function plan_disk_analysis(n::Int, m::Int)
155+
plan = ccall((:ft_plan_disk_analysis, libfasttransforms), Ptr{ft_plan_struct}, (Int, Int), n, m)
156+
DiskFFTWPlan{Float64}(plan)
157+
end
158+
disk_synthesis!(P::DiskFFTWPlan{Float64}, A::Matrix{Float64}) = ccall((:ft_execute_disk_synthesis, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
159+
disk_analysis!(P::DiskFFTWPlan{Float64}, A::Matrix{Float64}) = ccall((:ft_execute_disk_analysis, libfasttransforms), Nothing, (Ptr{ft_plan_struct}, Ptr{Float64}, Int, Int), P, A, size(A, 1), size(A, 2))
160+
end

0 commit comments

Comments
 (0)